M.Hiroi's Home Page

Common Lisp Programming

お気楽 Common Lisp プログラミング入門 : 自作ライブラリ編

[ Common Lisp | library ]

ntheory

ntheory は素数など「整数論 (Number Theroy)」に関連した関数を集めたライブラリです。

●インストール

アーカイブファイル minlib.tar.gz をインストールする、または プログラムリスト にある以下の 4 つのファイルを、~/common-lisp/ 以下の適当なサブディレクトリ (たとえば ntheory など) に配置してください。

●仕様

●説明

●簡単なテスト

* (asdf:test-system :ntheory)
; compiling file ... 略 ...

----- test start -----

(PRIME-NTH 0)
=> 2 OK

(PRIME-NTH 24)
=> 97 OK

(PRIME-NTH 40)
=> 179 OK

(LNTH 0 *TWIN-PRIMES*)
=> (3 5) OK

(LNTH 8 *TWIN-PRIMES*)
=> (101 103) OK

(LNTH 21 *TWIN-PRIMES*)
=> (419 421) OK

(PRIME-COUNT 100)
=> 25 OK

(PRIME-COUNT 1000)
=> 168 OK

(PRIME-COUNT 10000)
=> 1229 OK

(PRIME-NEXT 97)
=> 101 OK

(PRIME-NEXT 1000)
=> 1009 OK

(PRIME-NEXT 10000)
=> 10007 OK

(PRIME-PREV 2)
=> NIL OK

(PRIME-PREV 3)
=> 2 OK

(PRIME-PREV 101)
=> 97 OK

(PRIME-PREV 1000)
=> 997 OK

(PRIME-PREV 10000)
=> 9973 OK

(PRIME-RANGE 2 100)
=> (2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97) OK

(PRIME-RANGE 1000 1100)
=> (1009 1013 1019 1021 1031 1033 1039 1049 1051 1061 1063 1069 1087 1091 1093
    1097) OK

(SIEVE 100)
=> #(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97) OK

(LET ((XS (SIEVE 541)))
  (AREF XS (1- (LENGTH XS))))
=> 541 OK

(SEGMENTED-SIEVE 2 100)
=> #(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97) OK

(SEGMENTED-SIEVE 2147483648 2147483904)
=> #(2147483659 2147483693 2147483713 2147483743 2147483777 2147483783
     2147483813 2147483857 2147483867 2147483869 2147483887 2147483893) OK

(MODPOW 4 13 497)
=> 445 OK

(MODPOW 5 4321 10)
=> 5 OK

(MODPOW 7 654321 10)
=> 7 OK

(FERMAT-TEST 4759123129)
=> T OK

(FERMAT-TEST 4759123141)
=> NIL OK

(FERMAT-TEST 4759123151)
=> T OK

(MILLER-RABIN-TEST 4759123129)
=> T OK

(MILLER-RABIN-TEST 4759123141)
=> NIL OK

(MILLER-RABIN-TEST 4759123151)
=> T OK

(MAPCAR #'FERMAT-TEST
        (MAPCAR (LAMBDA (X) (1- (EXPT 2 X)))
                '(1 2 3 4 5 6 7 8 9 10 13 31 32 61)))
=> (NIL T T NIL T NIL T NIL NIL NIL T T NIL T) OK

(MAPCAR #'MILLER-RABIN-TEST
        (MAPCAR (LAMBDA (X) (1- (EXPT 2 X)))
                '(1 2 3 4 5 6 7 8 9 10 13 31 32 61)))
=> (NIL T T NIL T NIL T NIL NIL NIL T T NIL T) OK

(MAPCAR #'LUCAS-LEHMER-TEST '(1 2 3 4 5 6 7 8 9 10 13 31 32 61))
=> (NIL T T NIL T NIL T NIL NIL NIL T T NIL T) OK

(FACTORIZATION 6)
=> ((2 . 1) (3 . 1)) OK

(FACTORIZATION 12345678)
=> ((2 . 1) (3 . 2) (47 . 1) (14593 . 1)) OK

(FACTORIZATION 123456789)
=> ((3 . 2) (3607 . 1) (3803 . 1)) OK

(FACTORIZATION 1234567890)
=> ((2 . 1) (3 . 2) (5 . 1) (3607 . 1) (3803 . 1)) OK

(FACTORIZATION 1111111111)
=> ((11 . 1) (41 . 1) (271 . 1) (9091 . 1)) OK

(FACTORS-TO-NUMBER '((2 . 1) (3 . 2) (47 . 1) (14593 . 1)))
=> 12345678 OK

(FACTORS-TO-NUMBER '((3 . 2) (3607 . 1) (3803 . 1)))
=> 123456789 OK

(FACTORS-TO-NUMBER '((2 . 1) (3 . 2) (5 . 1) (3607 . 1) (3803 . 1)))
=> 1234567890 OK

(FACTORS-TO-NUMBER (FACTORIZATION (1- (EXPT 2 64))))
=> 18446744073709551615 OK

(FACTORS-TO-NUMBER (FACTORIZATION (1+ (EXPT 2 64))))
=> 18446744073709551617 OK

(DIVISOR-COUNT 6)
=> 4 OK

(DIVISOR-COUNT 12345678)
=> 24 OK

(DIVISOR-COUNT 123456789)
=> 12 OK

(DIVISOR-COUNT 1234567890)
=> 48 OK

(DIVISOR-COUNT 1111111111)
=> 16 OK

(DIVISOR-TOTAL 6)
=> 12 OK

(DIVISOR-TOTAL 12345678)
=> 27319968 OK

(DIVISOR-TOTAL 123456789)
=> 178422816 OK

(DIVISOR-TOTAL 1234567890)
=> 3211610688 OK

(DIVISOR-TOTAL 1111111111)
=> 1246404096 OK

(DIVISOR-SIGMA 6 0)
=> 4 OK

(DIVISOR-SIGMA 6 1)
=> 12 OK

(DIVISOR-SIGMA 6 2)
=> 50 OK

(DIVISOR-SIGMA 1111111111 0)
=> 16 OK

(DIVISOR-SIGMA 1111111111 1)
=> 1246404096 OK

(DIVISOR-SIGMA 1111111111 2)
=> 1245528410223519376 OK

(DIVISORS 6)
=> (1 2 3 6) OK

(DIVISORS 12345678)
=> (1 2 3 6 9 18 47 94 141 282 423 846 14593 29186 43779 87558 131337 262674
    685871 1371742 2057613 4115226 6172839 12345678) OK

(DIVISORS 123456789)
=> (1 3 9 3607 3803 10821 11409 32463 34227 13717421 41152263 123456789) OK

(DIVISORS 1234567890)
=> (1 2 3 5 6 9 10 15 18 30 45 90 3607 3803 7214 7606 10821 11409 18035 19015
    21642 22818 32463 34227 36070 38030 54105 57045 64926 68454 108210 114090
    162315 171135 324630 342270 13717421 27434842 41152263 68587105 82304526
    123456789 137174210 205761315 246913578 411522630 617283945 1234567890) OK

(DIVISORS 1111111111)
=> (1 11 41 271 451 2981 9091 11111 100001 122221 372731 2463661 4100041
    27100271 101010101 1111111111) OK

(MAPCAR (LAMBDA (X) (TOTIENT X)) '(1 2 3 4 5 10 20 30 40 50 97))
=> (1 1 2 2 4 4 8 8 16 20 96) OK

(MAPCAR (LAMBDA (X) (MOBIUS X)) '(1 2 3 4 5 6 7 8 9 10))
=> (1 -1 -1 0 -1 1 -1 0 0 1) OK

(MAPCAR (LAMBDA (X) (ORDER-N X 11)) '(2 3 4 5 6 7 8 9 10))
=> (10 5 5 5 10 10 10 5 2) OK

(MAPCAR (LAMBDA (X) (ORDER-N X 13)) '(2 3 4 5 6 7 8 9 10 11 12))
=> (12 3 6 4 12 12 4 3 6 12 2) OK

(MAPCAR (LAMBDA (X) (PRIMITIVE-ROOT-P X 11)) '(2 3 4 5 6 7 8 9 10))
=> (T NIL NIL NIL T T T NIL NIL) OK

(MAPCAR (LAMBDA (X) (PRIMITIVE-ROOT-P X 13)) '(2 3 4 5 6 7 8 9 10 11 12))
=> (T NIL NIL NIL T T NIL NIL NIL T NIL) OK

(LET ((A NIL))
  (PRIMITIVE-ROOT (LAMBDA (X) (PUSH X A)) 11)
  A)
=> (8 7 6 2) OK

(LET ((A NIL))
  (PRIMITIVE-ROOT (LAMBDA (X) (PUSH X A)) 13)
  A)
=> (11 7 6 2) OK

(LET ((A NIL))
  (PRIMITIVE-ROOT (LAMBDA (X) (PUSH X A)) 41)
  A)
=> (35 34 30 29 28 26 24 22 19 17 15 13 12 11 7 6) OK

(DIGITS 1234567890)
=> (1 2 3 4 5 6 7 8 9 0) OK

(DIGITS 255 2)
=> (1 1 1 1 1 1 1 1) OK

(DIGITS 256 2)
=> (1 0 0 0 0 0 0 0 0) OK

(DIGITS 255 8)
=> (3 7 7) OK

(DIGITS 256 8)
=> (4 0 0) OK

(DIGITS 256 16)
=> (1 0 0) OK

(DIGITS 65535 16)
=> (15 15 15 15) OK

(MAPCAR #'PERFECT-NUMBER-P '(6 10 28 100 496 1000 8128 10000))
=> (T NIL T NIL T NIL T NIL) OK

(MAPCAR #'ABUNDANT-NUMBER-P '(6 10 28 100 496 1000 8128 10000))
=> (NIL NIL NIL T NIL T NIL T) OK

(MAPCAR #'DEFICIENT-NUMBER-P '(6 10 28 100 496 999 8128 9999))
=> (NIL T NIL NIL NIL T NIL T) OK

(MAPCAR #'AMICABLE-NUMBERS-P '(220 222 6232 6234) '(284 286 6368 6370))
=> (T NIL T NIL) OK

(EXTENDED-EUCLID 11 3)
=> (1 -1 4) OK

(EXTENDED-EUCLID 4321 1234)
=> (1 309 -1082) OK

(EXTENDED-EUCLID 12357 100102)
=> (1 30127 -3719) OK

(MAPCAR (LAMBDA (X) (POLYGONAL-NUMBER 3 X)) '(1 2 3 4 5 6 7 8 9 10))
=> (1 3 6 10 15 21 28 36 45 55) OK

(MAPCAR (LAMBDA (X) (POLYGONAL-NUMBER 4 X)) '(1 2 3 4 5 6 7 8 9 10))
=> (1 4 9 16 25 36 49 64 81 100) OK

(MAPCAR (LAMBDA (X) (POLYGONAL-NUMBER 5 X)) '(1 2 3 4 5 6 7 8 9 10))
=> (1 5 12 22 35 51 70 92 117 145) OK

(POLYGONAL-NUMBER-P 3 (POLYGONAL-NUMBER 3 100))
=> T OK

(POLYGONAL-NUMBER-P 3 (1+ (POLYGONAL-NUMBER 3 100)))
=> NIL OK

(POLYGONAL-NUMBER-P 4 (POLYGONAL-NUMBER 4 100))
=> T OK

(POLYGONAL-NUMBER-P 4 (1+ (POLYGONAL-NUMBER 4 100)))
=> NIL OK

(POLYGONAL-NUMBER-P 5 (POLYGONAL-NUMBER 5 100))
=> T OK

(POLYGONAL-NUMBER-P 5 (1+ (POLYGONAL-NUMBER 5 100)))
=> NIL OK

(MAPCAR #'FERMAT-NUMBER '(0 1 2 3 4 5))
=> (3 5 17 257 65537 4294967297) OK

(BINOMIAL-COEFFICIENTS 2)
=> (1 2 1) OK

(BINOMIAL-COEFFICIENTS 10)
=> (1 10 45 120 210 252 210 120 45 10 1) OK

(CRT '((2 3) (4 5)))
=> (14 15) OK

(CRT '((20 30) (40 50)))
=> (140 150) OK

(CRT '((2 3) (3 5) (2 7)))
=> (23 105) OK

(CRT '((20 30) (30 50) (20 70)))
=> (230 1050) OK

(CRT '((2 3) (3 5) (2 7) (4 11)))
=> (653 1155) OK

(CRT '((20 30) (30 50) (20 70) (40 110)))
=> (6530 11550) OK

(MAPCAR (LAMBDA (A) (LEGENDRE-SYMBOL A 5)) '(0 1 2 3 4))
=> (0 1 -1 -1 1) OK

(MAPCAR (LAMBDA (A) (LEGENDRE-SYMBOL A 7)) '(0 1 2 3 4 5 6))
=> (0 1 1 -1 1 -1 -1) OK

(MAPCAR (LAMBDA (A) (LEGENDRE-SYMBOL A 11)) '(0 1 2 3 4 5 6 7 8 9 10))
=> (0 1 -1 1 1 1 -1 -1 -1 1 -1) OK

(MAPCAR (LAMBDA (A) (JACOBI-SYMBOL A 9)) '(0 1 2 3 4 5 6 7 8))
=> (0 1 1 0 1 1 0 1 1) OK

(MAPCAR (LAMBDA (A) (JACOBI-SYMBOL A 11)) '(0 1 2 3 4 5 6 7 8 9 10))
=> (0 1 -1 1 1 1 -1 -1 -1 1 -1) OK

(MAPCAR (LAMBDA (A) (JACOBI-SYMBOL A 15)) '(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14))
=> (0 1 1 0 1 0 0 -1 1 0 0 -1 0 -1 -1) OK

(MAPCAR (LAMBDA (A) (MODSQRT A 11)) '(1 2 3 4 5 6 7 8 9 10))
=> (1 NIL 5 9 4 NIL NIL NIL 3 NIL) OK

(MAPCAR (LAMBDA (A) (MODSQRT A 13)) '(1 2 3 4 5 6 7 8 9 10 11 12))
=> (1 NIL 9 11 NIL NIL NIL NIL 3 7 NIL 8) OK

(MAPCAR (LAMBDA (A) (MODSQRT A 17)) '(1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16))
=> (1 6 NIL 2 NIL NIL NIL 12 14 NIL NIL NIL 8 NIL 7 4) OK

(CONT-FRAC 355/113)
=> (3 7 16) OK

(CONT-FRAC 81201/56660)
=> (1 2 3 4 5 6 7 8) OK

(CONT-FRAC-FLOAT (SQRT 2.0d0))
=> (1 2 2 2 2 2 2 2) OK

(CONT-FRAC-FLOAT (SQRT 3.0d0) K 16)
=> (1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1) OK

(CONT-FRAC-REDUCE '(3 7 16))
=> 355/113 OK

(CONT-FRAC-REDUCE '(1 2 3 4 5 6 7 8))
=> 81201/56660 OK

(CONT-FRAC-REDUCE '(1 2 2 2 2 2 2 2 2))
=> 1393/985 OK

(CONT-FRAC-REDUCE '(1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2))
=> 51409/29681 OK

(MULTIPLE-VALUE-LIST (PELLS-EQUATION 7))
=> (8 3) OK

(MULTIPLE-VALUE-LIST (PELLS-EQUATION 61))
=> (1766319049 226153980) OK

(MULTIPLE-VALUE-LIST (PELLS-EQUATION 109))
=> (158070671986249 15140424455100) OK

(REPEAT-DECIMAL 1/7)
=> ((0) (1 4 2 8 5 7)) OK

(REPEAT-DECIMAL 355/113)
=> ((3)
    (1 4 1 5 9 2 9 2 0 3 5 3 9 8 2 3 0 0 8 8 4 9 5 5 7 5 2 2 1 2 3 8 9 3 8 0 5
     3 0 9 7 3 4 5 1 3 2 7 4 3 3 6 2 8 3 1 8 5 8 4 0 7 0 7 9 6 4 6 0 1 7 6 9 9
     1 1 5 0 4 4 2 4 7 7 8 7 6 1 0 6 1 9 4 6 9 0 2 6 5 4 8 6 7 2 5 6 6 3 7 1 6
     8)) OK

(APPLY #'FROM-REPEAT-DECIMAL '((0) (1 4 2 8 5 7)))
=> 1/7 OK

(APPLY #'FROM-REPEAT-DECIMAL
       '((3)
         (1 4 1 5 9 2 9 2 0 3 5 3 9 8 2 3 0 0 8 8 4 9 5 5 7 5 2 2 1 2 3 8 9 3 8
          0 5 3 0 9 7 3 4 5 1 3 2 7 4 3 3 6 2 8 3 1 8 5 8 4 0 7 0 7 9 6 4 6 0 1
          7 6 9 9 1 1 5 0 4 4 2 4 7 7 8 7 6 1 0 6 1 9 4 6 9 0 2 6 5 4 8 6 7 2 5
          6 6 3 7 1 6 8)))
=> 355/113 OK

(CONT-FRAC-ROOT 2)
=> (1 2) OK

(CONT-FRAC-ROOT 7)
=> (2 1 1 1 4) OK

(CONT-FRAC-ROOT 97)
=> (9 1 5 1 1 1 1 1 1 5 1 18) OK

(CONT-FRAC-ROOT 111)
=> (10 1 1 6 1 1 20) OK

----- test end -----
TEST: 135
OK: 135
NG: 0
ERR: 0
T

●簡単なサンプルプログラム

●プログラムリスト

;;;
;;; ntheory.lisp : Number Theory (整数論)
;;;
;;; Copyright (c) 2023-2024 Makoto Hiroi
;;;
;;; Released under the MIT license
;;; https://opensource.org/license/mit/
;;;
(provide :ntheory)
(defpackage :ntheory (:use :cl :combination :lazy))
(in-package :ntheory)
(export '(*primes*
          *twin-primes*
          prime-nth
          prime-count
          prime-next
          prime-prev
          prime-range
          sieve
          segmented-sieve
          modpow
          fermat-test
          miller-rabin-test
          lucas-lehmer-test
          factorization
          factors-to-number
          divisor-count
          divisor-total
          divisor-sigma
          divisors
          totient
          mobius
          order-n
          primitive-root-p
          primitive-root
          digits
          perfect-number-p
          abundant-number-p
          deficient-number-p
          amicable-numbers-p
          extended-euclid
          polygonal-number
          polygonal-number-p
          fermat-number
          binomial-coefficients
          legendre-symbol
          jacobi-symbol
          modsqrt
          cont-frac-reduce
          cont-frac
          cont-frac-float
          cont-frac-root
          pells-equation
          repeat-decimal
          from-repeat-decimal
          ))

;;; 素数 (遅延ストリーム版)
(declaim (ftype (function (integer list) t) primes-from))
(defvar *primes* (lcons 2 (lcons 3 (lcons 5 (primes-from 7 '#1=(4 2 4 2 4 6 2 6 . #1#))))))

(defun primep (n &optional (ps (lcdr *primes*)))
  (let ((p (lcar ps)))
    (cond
     ((> (* p p) n) t)
     ((zerop (mod n p)) nil)
     (t (primep n (lcdr ps))))))

(defun primes-from (n ks)
  (if (primep n)
      (lcons n (primes-from (+ n (car ks)) (cdr ks)))
    (primes-from (+ n (car ks)) (cdr ks))))

;;; 双子素数
(defvar *twin-primes*
  (lfilter (lambda (xs) (= (apply #'- xs) -2)) (lzip *primes* (lcdr *primes*))))

;;; n 番目の素数を返す
(defun prime-nth (n) (lnth n *primes*))

;;; n までの素数を個数を求める
(defun prime-count (n &optional (ps *primes*) (c 0))
  (if (< n (lcar ps))
      c
    (prime-count n (lcdr ps) (1+ c))))

;;; n の次の素数を求める
(defun prime-next (n &optional (ps *primes*))
  (lcar (ldrop-while (lambda (p) (<= p n)) ps)))

;;; n の前の素数を求める
(defun prime-prev-sub (n &optional (ps *primes*))
  (if (and (< (lcar ps) n)
           (<= n (lcar (lcdr ps))))
      (lcar ps)
    (prime-prev-sub n (lcdr ps))))

(defun prime-prev (n)
  (when (< 2 n) (prime-prev-sub n)))

;;; 指定した範囲の素数を求める
(defun prime-range (low high)
  (tolist (ltake-while (lambda (p) (<= p high))
                       (ldrop-while (lambda (p) (< p low)) *primes*))))

;;; エラトステネスの篩
(defun sieve (n)
  ;; 先頭から 3, 5, 7, ...
  (do ((ns (make-array (floor n 2) :element-type '(unsigned-byte 1) :initial-element 0))
       (ps (make-array 1 :fill-pointer t :adjustable t :initial-contents '(2)))
       (x 3 (+ x 2)))
      ((> x n) ps)
      (let ((y (floor (- x 3) 2)))
        (when (zerop (sbit ns y))
          (vector-push-extend x ps)
          (when (<= (* x x) n)
            (do ((y (+ y x) (+ y x)))
                ((>= y (length ns)))
                (setf (sbit ns y) 1)))))))

;;; 区間ふるい low 以上 high 以下
(defun segmented-sieve (low high)
  (when (< low 2) (setf low 2))
  (let* ((size (1+ (- high low)))
         (ns (make-array size :element-type '(unsigned-byte 1) :initial-element 0)))
    (map nil
         (lambda (p)
           (do* ((s0 (- (* (ceiling low p) p) low))
                 (s1 (if (< p low) s0 (+ s0 p)) (+ s1 p)))
                ((>= s1 size))
                (setf (sbit ns s1) 1)))
         (sieve (floor (sqrt high))))
    (do ((rs (make-array 1 :fill-pointer 0 :adjustable t))
         (i 0 (1+ i)))
        ((>= i size) rs)
        (when (zerop (sbit ns i))
          (vector-push-extend (+ low i) rs)))))

;;; べき剰余
;;; (mod (expt b e) m) を高速に求める
(defun modpow (b e m)
  (do ((r 1)
       (e e (ash e -1))
       (b b (mod (* b b) m)))
      ((zerop e) r)
      (when (plusp (logand e 1))
        (setf r (mod (* r b) m)))))

;;; フェルマーテスト
(defun fermat-test (p &optional (k 20))
  (cond
   ((< p 2) nil)
   ((= p 2) t)
   ((zerop (logand p 1)) nil)
   (t
    (do ((k k (1- k)))
        ((zerop k) t)
        (let ((a (+ (random (- p 2)) 2)))
          (when (or (/= (gcd p a) 1)
                    (/= (modpow a (1- p) p) 1))
            (return nil)))))))

;;; ミラーラビンテスト

;;; n を (2 ** c) * d に分解する
(defun factor2 (n &optional (c 0))
  (if (= (logand n 1) 1)
      (values c n)
    (factor2 (ash n -1) (1+ c))))

(defun miller-rabin-sub (p s d ks)
  (dolist (a ks t)
    (when (>= a p) (return t))
    (unless (= (modpow a d p) 1)
      (dotimes (r s (return-from miller-rabin-sub nil))
        (when (= (modpow a (* (expt 2 r) d) p) (1- p))
          (return t))))))

(defun miller-rabin-test (p &optional (k 20))
  (cond
   ((< p 2) nil)
   ((= p 2) t)
   ((zerop (logand p 1)) nil)
   (t
    (multiple-value-bind
     (s d)
     (factor2 (1- p))
     (cond
      ((< p 4759123141)
       (miller-rabin-sub p s d '(2 7 61)))
      ((<= p #x10000000000000000)
       (miller-rabin-sub p s d '(2 325 9375 28178 450775 9780504 1795265022)))
      (t
       (do ((ks nil))
           ((>= (length ks) k) (miller-rabin-sub p s d ks))
           (pushnew (1+ (random (- p 2))) ks))))))))

;;; リュカ-レーマーテスト
(defun lucas-lehmer-test (p)
  (cond
   ((< p 2) nil)
   ((= p 2) t)
   (t
    (do ((m (1- (expt 2 p)))
         (x 4 (mod (- (* x x) 2) m))
         (i (- p 2) (1- i)))
        ((zerop i) (zerop (mod x m)))))))

;;;
;;; 素因数分解
;;;
(defun factor-sub (n m &optional (c 0))
  (if (zerop (mod n m))
      (factor-sub (/ n m) m (1+ c))
    (values c n)))

;;; end まで試し割り (wheel factorization)
(defun trial (n &key (end 10000))
  (let ((x 2) (a nil))
    (dolist (k '(1 2 2 . #1=(4 2 4 2 4 6 2 6 . #1#)))
      (cond
       ((= n 1)
        (return (values 1 a)))
       ((< n (* x x))
        (return (values 1 (cons (cons n 1) a))))
       ((> x end)
        (return (values n a)))
       (t
        (multiple-value-bind
         (c m)
         (factor-sub n x)
         (when (plusp c)
           (push (cons x c) a))
         (setf n m)
         (incf x k)))))))

;;;
;;; ポラード・ロー素因数分解法 (Pollard's rho algorithm)
;;;
(defun gen (x n c) (mod (+ (* x x) c) n))

;;; ブレントの変形版
(defun find-factor-brent (n c a &optional (limit 1000000))
  (do* ((xi (gen a n c))
        (xj (gen xi n c))
        (m 2048) (k 1)
        (inc 3) (cnt 0))
       ((>= cnt limit))
       (do ((l 0) (w 1) (xjs xj))
           ((>= l k))
           (setf w 1)
           (dotimes (i (min (- k l) m))
             (setf xj (gen xj n c)
                   w  (mod (* (abs (- xj xi)) w) n))
             (incf l))
           (let ((d (gcd w n)))
             (cond
              ((< 1 d n)
               (return-from find-factor-brent d))
              ((= d n)
               (loop
                (setf xjs (gen xjs n c)
                      d (gcd (abs (- xjs xj)) n))
                (when (< 1 d)
                  (return-from find-factor-brent (if (< d n) d))))))))
       (setf xi xj)
       (dotimes (i inc) (setf xj (gen xj n c)))
       (incf cnt k)
       (setf k (* k 2))
       (incf inc k)))

;;;
;;; ポラードの p-1 法
;;;
(defun find-factor-pm1 (n a &optional (limit 1000000))
  (do ((m 4096)
       (k 2))
      ((>= k limit))
      (let ((a1 a) (k1 k) (w 1))
        (setf w 1)
        (dotimes (i m)
          (setf a (modpow a k n))
          (when (zerop a)
            (return-from find-factor-pm1 nil))
          (setf w (mod (* w (1- a)) n))
          (incf k))
        (let ((d (gcd w n)))
          (cond
           ((< 1 d n)
            (return-from find-factor-pm1 d))
           ((= d n)
            (loop
             (setf a1 (modpow a1 k1 n))
             (incf k1)
             (setf d (gcd (1- a1) n))
             (when (< 1 d)
               (return-from find-factor-pm1 (if (< d n) d))))))))))

;;; ロー法 + p-1 法
(defun find-factor (n &optional (limit 1000000))
  (let ((p (find-factor-brent n 1 2 limit)))
    (unless p
      (setf p (find-factor-pm1 n 3 limit)))
    (when p
      (when (miller-rabin-test p)
        (return-from find-factor p))
      (let ((q (/ n p)))
        (when (miller-rabin-test q)
          (return-from find-factor q)))
      (find-factor p limit))))

;;; 試し割りのあとロー法と p-1 法で素因数分解する
(defun factorization (n &optional (limit 1000000))
  (multiple-value-bind
   (n ps)
   (trial n)
   (loop
    (cond
     ((= n 1)
      (return (values (nreverse ps) 1)))
     ((miller-rabin-test n)
      (return (values (nreverse (cons (cons n 1) ps)) 1)))
     (t
      (let ((p (find-factor n limit)))
        (unless p
          (return (values (nreverse ps) n)))
        (multiple-value-bind
         (c m)
         (factor-sub n p)
         (push (cons p c) ps)
         (setf n m))))))))

;;; 素因数分解したリストから数を求める
(defun factors-to-number (xs)
  (reduce (lambda (a x) (* a (expt (car x) (cdr x)))) xs :initial-value 1))

;;;
;;; 約数
;;;

;;; 個数
(defun divisor-count (n)
  (multiple-value-bind
   (ps m)
   (factorization n)
   (when (= m 1)
     (reduce (lambda (a x) (* a (+ 1 (cdr x))))
             ps
             :initial-value 1))))

;;; 合計値
(defun divisor-total (n)
  (multiple-value-bind
   (ps m)
   (factorization n)
   (when (= m 1)
     (reduce (lambda (a x)
               ;; 等比数列の和 (p ** (q + 1) - 1) // (p - 1)
               (* a (/ (1- (expt (car x) (1+ (cdr x)))) (1- (car x)))))
             ps
             :initial-value 1))))

;;; 約数
(defun divisor-sub (p n &optional (a nil))
  (if (zerop n)
      (cons 1 a)
    (divisor-sub p (- n 1) (cons (expt p n) a))))

(defun divisors (n)
  (let ((x (factorization n)))
    (sort (reduce (lambda (a y)
                    (mapcar (lambda (xs) (apply #'* xs))
                            (product-set (divisor-sub (car y) (cdr y)) a)))
                  (cdr x)
                  :initial-value (divisor-sub (caar x) (cdar x)))
          #'<)))

;;; 約数関数 σ(n), n の約数の k 乗の総和
(defun divisor-sigma (n k)
  (if (zerop k)
      (divisor-count n)
    (multiple-value-bind
     (ps m)
     (factorization n)
     (when (= m 1)
       (reduce
        (lambda (a xs)
          (* a (/ (1- (expt (car xs) (* k (1+ (cdr xs)))))
                  (1- (expt (car xs) k)))))
        ps
        :initial-value 1)))))

;;; オイラーのトーシェント関数
(defun totient (n)
  (multiple-value-bind
   (xs m)
   (factorization n)
   (when (= m 1)
     (reduce (lambda (a ps) (/ (* a (1- (car ps))) (car ps)))
             xs
             :initial-value n))))

;;; メビウス関数
(defun mobius (n)
  (if (= n 1)
      1
    (multiple-value-bind
     (xs m)
     (factorization n)
     (when (= m 1)
       (if (some (lambda (ps) (< 1 (cdr ps))) xs)
           0
         (expt -1 (length xs)))))))

;;; 位数を求める
(defun order-n (a n)
  (when (= (gcd a n) 1)
    (let ((d (1- n)))
      (multiple-value-bind
       (xs m)
       (factorization d)
       (when (= m 1)
         (loop
          (dolist (ps xs (return-from order-n d))
            (let ((p (car ps)))
              (when (and (zerop (mod d p))
                         (= (modpow a (/ d p) n) 1))
                (setf d (/ d p))
                (return))))))))))

;;; 原始根の判定
(defun primitive-root-p (a n)
  (when (= (gcd a n) 1)
    (multiple-value-bind
     (xs m)
     (factorization (1- n))
     (when (= m 1)
       (notany (lambda (ps) (= (modpow a (/ (1- n) (car ps)) n) 1)) xs)))))

;;; 原始根を求める
(defun primitive-root (fn n)
  (multiple-value-bind
   (xs m)
   (factorization (1- n))
   (when (= m 1)
     (do ((a 2 (1+ a)))
         ((>= a n))
         (when (every (lambda (ps) (/= (modpow a (/ (1- n) (car ps)) n) 1)) xs)
           (funcall fn a))))))

;;; 正整数 n を base 進数の各桁ごとに分解
(defun digits (n &optional (base 10) (a nil))
  (if (zerop n)
      a
    (multiple-value-bind
     (p q)
     (truncate n base)
     (digits p base (cons q a)))))

;;; 完全数
(defun perfect-number-p (n)
  (= (- (divisor-total n) n) n))

;;; 過剰数
(defun abundant-number-p (n)
  (> (- (divisor-total n) n) n))

;;; 不足数
(defun deficient-number-p (n)
  (< (- (divisor-total n) n) n))

;;; 友愛数
(defun amicable-numbers-p (m n)
  (and (= (- (divisor-total m) m) n)
       (= (- (divisor-total n) n) m)))

;;; 拡張ユークリッドの互除法
(defun extended-euclid (a b)
  (do ((xs (list a 1 0))
       (ys (list b 0 1)))
      ((zerop (first ys)) xs)
      (multiple-value-bind
       (q z)
       (floor (first xs) (first ys))
       (psetf xs ys
              ys (list z
                       (- (second xs) (* q (second ys)))
                       (- (third xs) (* q (third ys))))))))

;;; 多角数
(defun polygonal-number (p n)
  (/ (- (* n n (- p 2)) (* (- p 4) n)) 2))

;;; m は p 角数か？
(defun polygonal-number-p (p m)
  (let* ((a (- p 2))
         (b (- p 4))
         (c (+ (* b b) (* 8 a m)))
         (d (isqrt c)))
    (and (= (* d d) c)
         (zerop (mod (+ b d) (* 2 a))))))

;;; フェルマー数
(defun fermat-number (n)
  (1+ (expt 2 (expt 2 n))))

;;; 二項係数 (パスカルの三角形の n 段目を返す)
(defun binomial-coefficients (n &optional (a '(1)))
  (if (zerop n)
      a
    (binomial-coefficients
     (1- n)
     (maplist (lambda (xs) (if (null (cdr xs)) 1 (+ (car xs) (cadr xs))))
              (cons 0 a)))))

;;; 中国剰余定理
(defun crt-sub (a m b n)
  (let ((zs (extended-euclid m n)))
    (when (zerop (mod (- b a) (car zs)))
      (let ((s (/ (- b a) (car zs)))
            (l (lcm m n)))
        (values (mod (+ (* s m (cadr zs)) a) l) l)))))

(defun crt (xs)
  (reduce (lambda (zs ys)
            (multiple-value-bind
             (z l)
             (crt-sub (car zs) (cadr zs) (car ys) (cadr ys))
             (when (null z) (return-from crt nil))
             (list z l)))
          xs))

;;; ルジャンドル記号
(defun legendre-symbol (a n)
  (when (or (= n 2)
            (not (miller-rabin-test n)))
    (error "n should be an odd prime"))
  (cond
   ((zerop (mod a n)) 0)
   ((= (modpow a (/ (1- n) 2) n) 1) 1)
   (t -1)))

;;; ヤコビ記号
(defun jacobi-symbol (a n)
  (when (or (minusp n) (evenp n))
    (error "n should be an odd positive integer"))
  (setf a (mod a n))
  (do ((k 1))
      ((<= a 0) (if (= n 1) k 0))
      ;; 第２法則
      (loop
       (unless (zerop (mod a 2)) (return))
       (setf a (/ a 2))
       (let ((r (mod n 8)))
         (when (or (= r 3) (= r 5))
           (setf k (- k)))))
      ;; 相互法則
      (psetf a n n a)
      (when (and (= (mod a 4) 3) (= (mod n 4) 3))
        (setf k (- k)))
      (setf a (mod a n))))

;;; 平方非剰余を一つ求める
(defun modsqrt-sub (n)
  (do ((i 2 (1+ i)))
      ((>= i n))
      (when (= (legendre-symbol i n) -1)
        (return i))))

;;; 合同式 x**2 ≡ a (mod n) を解く
(defun modsqrt (a n)
  (cond
   ((= (legendre-symbol a n) -1) nil)
   ((= (mod n 4) 3)
    (modpow a (/ (1+ n) 4) n))
   ((= (mod n 8) 5)
    (let ((x (modpow a (/ (+ n 3)  8) n)))
      (if (= (modpow x 2 n) a)
          x
        (mod (* x (modpow 2 (/ (1- n) 4) n)) n))))
   (t
    ;; Tonelli–Shanks のアルゴリズム
    (multiple-value-bind
     (s k)
     (factor2 (1- n))           ; n - 1 = 2^s * k
     (let* ((d (modsqrt-sub n)) ; 平方非剰余を選ぶ
            (a1 (modpow a k n))
            (d1 (modpow d k n))
            (m 0))
       (dotimes (i s)
         (when (= (modpow (* a1 (expt d1 m)) (expt 2 (- s i 1)) n) (1- n))
           (incf m (expt 2 i))))
       (mod (* (modpow a (/ (1+ k) 2) n) (modpow d1 (/ m 2) n)) n))))))

;;;
;;; 連分数
;;;

;;; 分数を正則連分数に変換
(defun cont-frac (x)
  (unless (rationalp x)
    (error "must be rational number"))
  (do ((a (numerator x))
       (b (denominator x))
       (xs nil))
      ((<= b 0) (nreverse xs))
      (multiple-value-bind
       (p q)
       (floor a b)
       (push p xs)
       (setf a b b q))))

;;; 正則連分数を有理数に戻す
(defun cont-frac-reduce (xs)
  (if (null (cdr xs))
      (car xs)
    (+ (car xs)
       (/ (cont-frac-reduce (cdr xs))))))

;;; 実数を正則連分数に展開する
(defun cont-frac-float (x &key (k 8) (e 1e-16))
  (do ((k k (1- k))
       (xs nil))
      ((zerop k) (nreverse xs))
      (multiple-value-bind
       (p q)
       (floor x)
       (push p xs)
       (when (< (abs q) e)
         (return (nreverse xs)))
       (setf x (/ q)))))

;;; √m の連分数展開
(defun cont-frac-root (m)
  (let* ((c0 (isqrt m))
         (cf (list c0)))
    (if (= (* c0 c0) m)
        cf
      (let ((a 1) (b 0) (c c0))
        (loop
         (setf b (- (* c a) b)
               a (floor (- m (* b b)) a)
               c (floor (+ c0 b) a))
         (push c cf)
         (when (= a 1)
           (return (nreverse cf))))))))

;;; ペル方程式 p^2 - m*q^2 = 1
(defun pells-equation (m)
  (let ((cf (butlast (cont-frac-root m))))
    (unless cf
      (error "pells-equation: domain error"))
    (let ((p (cont-frac-sub cf 1 (car cf)))
          (q (cont-frac-sub cf 0 1)))
      (if (= (- (* p p) (* m q q)) 1)
          (values p q)
        (values (+ (* p p) (* m q q)) (* 2 p q))))))

;;;
;;; 循環小数
;;;

;;; 分数を循環小数に変換する
(defun repeat-decimal (x)
  (unless (rationalp x)
    (error "must be rational number"))
  (do ((m (numerator x))
       (n (denominator x))
       (xs nil) (ys nil))
      (nil)
      (multiple-value-bind
       (p q)
       (floor m n)
       (cond
        ((zerop q)
         (push p ys)
         (return (list (nreverse ys) '(0))))
        (t
         (let ((i (position q xs)))
           (push p ys)
           (when (and i (>= i 0))
             (return (list (nreverse (subseq ys (1+ i)))
                           (nreverse (subseq ys 0 (1+ i))))))
           (push q xs)
           (setq m (* q 10))))))))


;;; 循環小数を分数に直す
(defun from-repeat-decimal (xs ys)
  (let ((p0 (reduce (lambda (a x) (+ (* a 10) x)) xs :initial-value 0))
        (q0 (expt 10 (1- (length xs))))
        (p1 (reduce (lambda (a x) (+ (* a 10) x)) ys :initial-value 0))
        (q1 (1- (expt 10 (length ys)))))
    (/ (+ (* q1 p0) p1) (* q0 q1))))

リスト : ntheory.asd

(defsystem :ntheory
  :description "Number Theory (整数論)"
  :version "0.2.10"
  :author "Makoto Hiroi"
  :license "MIT"
  :depends-on (:combination :lazy)
  :in-order-to ((test-op (test-op :ntheory_tst)))
  :components ((:file "ntheory")))

;;;
;;; ntheory_tst.lisp : ntheory のテスト
;;;
;;; Copyright (C) 2023-2024 Makoto Hiroi
;;;
;;; Released under the MIT license
;;; https://opensource.org/license/mit/
;;;
(provide :ntheory_tst)
(defpackage :ntheory_tst (:use :cl :mintst :ntheory :lazy))
(in-package :ntheory_tst)
(export '(test))

;;; 実行
(defun test ()
  (initial)
  (run (prime-nth 0) 2)
  (run (prime-nth 24) 97)
  (run (prime-nth 40) 179)
  (run (lnth 0 *twin-primes*) '(3 5))
  (run (lnth 8 *twin-primes*) '(101 103))
  (run (lnth 21 *twin-primes*) '(419 421))
  (run (prime-count 100) 25)
  (run (prime-count 1000) 168)
  (run (prime-count 10000) 1229)
  (run (prime-next 97) 101)
  (run (prime-next 1000) 1009)
  (run (prime-next 10000) 10007)
  (run (prime-prev 2) nil)
  (run (prime-prev 3) 2)
  (run (prime-prev 101) 97)
  (run (prime-prev 1000) 997)
  (run (prime-prev 10000) 9973)
  (run (prime-range 2 100)
       '(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97))
  (run (prime-range 1000 1100)
       '(1009 1013 1019 1021 1031 1033 1039 1049 1051 1061 1063 1069 1087 1091 1093 1097))
  (run (sieve 100)
       #(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97) :test #'equalp)
  (run (let ((xs (sieve 541))) (aref xs (1- (length xs)))) 541)
  (run (segmented-sieve 2 100)
       #(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97) :test #'equalp)
  (run (segmented-sieve #x80000000 #x80000100)
       #(2147483659 2147483693 2147483713 2147483743 2147483777 2147483783
         2147483813 2147483857 2147483867 2147483869 2147483887 2147483893) :test #'equalp)
  (run (modpow 4 13 497) 445)
  (run (modpow 5 4321 10) 5)
  (run (modpow 7 654321 10) 7)
  (run (fermat-test 4759123129) t)
  (run (fermat-test 4759123141) nil)
  (run (fermat-test 4759123151) t)
  (run (miller-rabin-test 4759123129) t)
  (run (miller-rabin-test 4759123141) nil)
  (run (miller-rabin-test 4759123151) t)
  (run (mapcar #'fermat-test (mapcar (lambda (x) (1- (expt 2 x))) '(1 2 3 4 5 6 7 8 9 10 13 31 32 61)))
       '(nil t t nil t nil t nil nil nil t t nil t))
  (run (mapcar #'miller-rabin-test (mapcar (lambda (x) (1- (expt 2 x))) '(1 2 3 4 5 6 7 8 9 10 13 31 32 61)))
       '(nil t t nil t nil t nil nil nil t t nil t))
  (run (mapcar #'lucas-lehmer-test '(1 2 3 4 5 6 7 8 9 10 13 31 32 61))
       '(nil t t nil t nil t nil nil nil t t nil t))
  (run (factorization 6) '((2 . 1) (3 . 1)))
  (run (factorization 12345678) '((2 . 1) (3 . 2) (47 . 1) (14593 . 1)))
  (run (factorization 123456789) '((3 . 2) (3607 . 1) (3803 . 1)))
  (run (factorization 1234567890) '((2 . 1) (3 . 2) (5 . 1) (3607 . 1) (3803 . 1)))
  (run (factorization 1111111111) '((11 . 1) (41 . 1) (271 . 1) (9091 . 1)))
  (run (factors-to-number '((2 . 1) (3 . 2) (47 . 1) (14593 . 1))) 12345678)
  (run (factors-to-number '((3 . 2) (3607 . 1) (3803 . 1))) 123456789)
  (run (factors-to-number '((2 . 1) (3 . 2) (5 . 1) (3607 . 1) (3803 . 1))) 1234567890)
  (run (factors-to-number (factorization (1- (expt 2 64)))) (1- (expt 2 64)))
  (run (factors-to-number (factorization (1+ (expt 2 64)))) (1+ (expt 2 64)))
  (run (divisor-count 6) 4)
  (run (divisor-count 12345678) 24)
  (run (divisor-count 123456789) 12)
  (run (divisor-count 1234567890) 48)
  (run (divisor-count 1111111111) 16)
  (run (divisor-total 6) 12)
  (run (divisor-total 12345678) 27319968)
  (run (divisor-total 123456789) 178422816)
  (run (divisor-total 1234567890) 3211610688)
  (run (divisor-total 1111111111) 1246404096)
  (run (divisor-sigma 6 0) 4)
  (run (divisor-sigma 6 1) 12)
  (run (divisor-sigma 6 2) 50)
  (run (divisor-sigma 1111111111 0) 16)
  (run (divisor-sigma 1111111111 1) 1246404096)
  (run (divisor-sigma 1111111111 2) 1245528410223519376)
  (run (divisors 6) '(1 2 3 6))
  (run (divisors 12345678)
       '(1 2 3 6 9 18 47 94 141 282 423 846 14593 29186 43779 87558 131337 262674 685871
         1371742 2057613 4115226 6172839 12345678))
  (run (divisors 123456789)
       '(1 3 9 3607 3803 10821 11409 32463 34227 13717421 41152263 123456789))
  (run (divisors 1234567890)
       '(1 2 3 5 6 9 10 15 18 30 45 90 3607 3803 7214 7606 10821 11409 18035 19015 21642
         22818 32463 34227 36070 38030 54105 57045 64926 68454 108210 114090 162315 171135
         324630 342270 13717421 27434842 41152263 68587105 82304526 123456789 137174210
         205761315 246913578 411522630 617283945 1234567890))
  (run (divisors 1111111111)
       '(1 11 41 271 451 2981 9091 11111 100001 122221 372731 2463661 4100041 27100271
           101010101 1111111111))
  (run (mapcar (lambda (x) (totient x)) '(1 2 3 4 5 10 20 30 40 50 97))
       '(1 1 2 2 4 4 8 8 16 20 96))
  (run (mapcar (lambda (x) (mobius x)) '(1 2 3 4 5 6 7 8 9 10))
       '(1 -1 -1 0 -1 1 -1 0 0 1))
  (run (mapcar (lambda (x) (order-n x 11)) '(2 3 4 5 6 7 8 9 10))
       '(10 5 5 5 10 10 10 5 2))
  (run (mapcar (lambda (x) (order-n x 13)) '(2 3 4 5 6 7 8 9 10 11 12))
       '(12 3 6 4 12 12 4 3 6 12 2))
  (run (mapcar (lambda (x) (primitive-root-p x 11)) '(2 3 4 5 6 7 8 9 10))
       '(T NIL NIL NIL T T T NIL NIL))
  (run (mapcar (lambda (x) (primitive-root-p x 13)) '(2 3 4 5 6 7 8 9 10 11 12))
       '(T NIL NIL NIL T T NIL NIL NIL T NIL))
  (run (let ((a nil)) (primitive-root (lambda (x) (push x a)) 11) a)
       '(8 7 6 2))
  (run (let ((a nil)) (primitive-root (lambda (x) (push x a)) 13) a)
       '(11 7 6 2))
  (run (let ((a nil)) (primitive-root (lambda (x) (push x a)) 41) a)
       '(35 34 30 29 28 26 24 22 19 17 15 13 12 11 7 6))
  (run (digits 1234567890) '(1 2 3 4 5 6 7 8 9 0))
  (run (digits 255 2) '(1 1 1 1 1 1 1 1))
  (run (digits 256 2) '(1 0 0 0 0 0 0 0 0))
  (run (digits 255 8) '(3 7 7))
  (run (digits 256 8) '(4 0 0))
  (run (digits 256 16) '(1 0 0))
  (run (digits 65535 16) '(15 15 15 15))
  (run (mapcar #'perfect-number-p '(6 10 28 100 496 1000 8128 10000))
       '(t nil t nil t nil t nil))
  (run (mapcar #'abundant-number-p '(6 10 28 100 496 1000 8128 10000))
       '(nil nil nil t nil t nil t))
  (run (mapcar #'deficient-number-p '(6 10 28 100 496 999 8128 9999))
       '(nil t nil nil nil t nil t))
  (run (mapcar #'amicable-numbers-p '(220 222 6232 6234) '(284 286 6368 6370))
       '(t nil t nil))
  (run (extended-euclid 11 3) '(1 -1 4))
  (run (extended-euclid 4321 1234) '(1 309 -1082))
  (run (extended-euclid 12357 100102) '(1 30127 -3719))
  (run (mapcar (lambda (x) (polygonal-number 3 x)) '(1 2 3 4 5 6 7 8 9 10))
       '(1 3 6 10 15 21 28 36 45 55))
  (run (mapcar (lambda (x) (polygonal-number 4 x)) '(1 2 3 4 5 6 7 8 9 10))
       '(1 4 9 16 25 36 49 64 81 100))
  (run (mapcar (lambda (x) (polygonal-number 5 x)) '(1 2 3 4 5 6 7 8 9 10))
       '(1 5 12 22 35 51 70 92 117 145))
  (run (polygonal-number-p 3 (polygonal-number 3 100)) t)
  (run (polygonal-number-p 3 (1+ (polygonal-number 3 100))) nil)
  (run (polygonal-number-p 4 (polygonal-number 4 100)) t)
  (run (polygonal-number-p 4 (1+ (polygonal-number 4 100))) nil)
  (run (polygonal-number-p 5 (polygonal-number 5 100)) t)
  (run (polygonal-number-p 5 (1+ (polygonal-number 5 100))) nil)
  (run (mapcar #'fermat-number '(0 1 2 3 4 5)) '(3 5 17 257 65537 4294967297))
  (run (binomial-coefficients 2) '(1 2 1))
  (run (binomial-coefficients 10) '(1 10 45 120 210 252 210 120 45 10 1))
  (run (crt '((2 3) (4 5))) '(14 15))
  (run (crt '((20 30) (40 50))) '(140 150))
  (run (crt '((2 3) (3 5) (2 7))) '(23 105))
  (run (crt '((20 30) (30 50) (20 70))) '(230 1050))
  (run (crt '((2 3) (3 5) (2 7) (4 11))) '(653 1155))
  (run (crt '((20 30) (30 50) (20 70) (40 110))) '(6530 11550))
  (run (mapcar (lambda (a) (legendre-symbol a 5)) '(0 1 2 3 4))
       '(0 1 -1 -1 1))
  (run (mapcar (lambda (a) (legendre-symbol a 7)) '(0 1 2 3 4 5 6))
       '(0 1 1 -1 1 -1 -1))
  (run (mapcar (lambda (a) (legendre-symbol a 11)) '(0 1 2 3 4 5 6 7 8 9 10))
       '(0 1 -1 1 1 1 -1 -1 -1 1 -1))
  (run (mapcar (lambda (a) (jacobi-symbol a 9)) '(0 1 2 3 4 5 6 7 8))
       '(0 1 1 0 1 1 0 1 1))
  (run (mapcar (lambda (a) (jacobi-symbol a 11)) '(0 1 2 3 4 5 6 7 8 9 10))
       '(0 1 -1 1 1 1 -1 -1 -1 1 -1))
  (run (mapcar (lambda (a) (jacobi-symbol a 15)) '(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14))
       '(0 1 1 0 1 0 0 -1 1 0 0 -1 0 -1 -1))
  (run (mapcar (lambda (a) (modsqrt a 11)) '(1 2 3 4 5 6 7 8 9 10))
       '(1 NIL 5 9 4 NIL NIL NIL 3 NIL))
  (run (mapcar (lambda (a) (modsqrt a 13)) '(1 2 3 4 5 6 7 8 9 10 11 12))
       '(1 NIL 9 11 NIL NIL NIL NIL 3 7 NIL 8))
  (run (mapcar (lambda (a) (modsqrt a 17)) '(1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16))
       '(1 6 NIL 2 NIL NIL NIL 12 14 NIL NIL NIL 8 NIL 7 4))
       '(1 6 NIL 2 NIL NIL NIL 12 14 NIL NIL NIL 8 NIL 7 4))
  (run (cont-frac 355/113) '(3 7 16))
  (run (cont-frac 81201/56660) '(1 2 3 4 5 6 7 8))
  (run (cont-frac-float (sqrt 2d0)) '(1 2 2 2 2 2 2 2))
  (run (cont-frac-float (sqrt 3d0) :k 16) '(1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1))
  (run (cont-frac-reduce '(3 7 16)) 355/113)
  (run (cont-frac-reduce '(1 2 3 4 5 6 7 8)) 81201/56660)
  (run (cont-frac-reduce '(1 2 2 2 2 2 2 2 2)) 1393/985)
  (run (cont-frac-reduce '(1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2)) 51409/29681)
  (run (multiple-value-list (pells-equation 7)) '(8 3))
  (run (multiple-value-list (pells-equation 61)) '(1766319049 226153980))
  (run (multiple-value-list (pells-equation 109)) '(158070671986249 15140424455100))
  (run (repeat-decimal 1/7) '((0) (1 4 2 8 5 7)))
  (run (repeat-decimal 355/113)
       '((3)
         (1 4 1 5 9 2 9 2 0 3 5 3 9 8 2 3 0 0 8 8 4 9 5 5 7 5 2 2 1 2 3 8 9 3 8 0 5 3 0
            9 7 3 4 5 1 3 2 7 4 3 3 6 2 8 3 1 8 5 8 4 0 7 0 7 9 6 4 6 0 1 7 6 9 9 1 1 5 0
            4 4 2 4 7 7 8 7 6 1 0 6 1 9 4 6 9 0 2 6 5 4 8 6 7 2 5 6 6 3 7 1 6 8)))
  (run (apply #'from-repeat-decimal '((0) (1 4 2 8 5 7))) 1/7)
  (run (apply #'from-repeat-decimal
       '((3)
         (1 4 1 5 9 2 9 2 0 3 5 3 9 8 2 3 0 0 8 8 4 9 5 5 7 5 2 2 1 2 3 8 9 3 8 0 5 3 0
            9 7 3 4 5 1 3 2 7 4 3 3 6 2 8 3 1 8 5 8 4 0 7 0 7 9 6 4 6 0 1 7 6 9 9 1 1 5 0
            4 4 2 4 7 7 8 7 6 1 0 6 1 9 4 6 9 0 2 6 5 4 8 6 7 2 5 6 6 3 7 1 6 8)))
       355/113)
  (run (cont-frac-root 2) '(1 2))
  (run (cont-frac-root 7) '(2 1 1 1 4))
  (run (cont-frac-root 97) '(9 1 5 1 1 1 1 1 1 5 1 18))
  (run (cont-frac-root 111) '(10 1 1 6 1 1 20))
  (final))

リスト : ntheory_tst.asd

(defsystem :ntheory_tst
  :description "test for ntheory"
  :version "0.2.10"
  :author "Makoto Hiroi"
  :license "MIT"
  :depends-on (:mintst :ntheory :lazy)
  :components ((:file "ntheory_tst"))
  :perform (test-op (o s) (symbol-call :ntheory_tst :test)))
Copyright (C) 2023-2024 Makoto Hiroi
All rights reserved.
[ Common Lisp | library ]