ISLispの数学ライブラリ

はじめに

ISLispで書かれた初等数学のライブラリです。線形代数と初等整数論の理解のためにSchemeで書いたものをISLispに移植しました。EISLのコンパイラのテストも兼ねています。線形代数関連は添え字を数学と同様に１からスタートするものとしています。

コンパイラのテストが主目的で移植は不完全かもしれません。間違いがありましたらお知らせください。

コンパイル

インタプリタでも一応は動きますが、素因数分解になるとインタプリタでは無理があります。EISLでコンパイルしてお使いください。

コード

(defmacro when (test :rest body)
  `(if ,test (progn ,@body)))

(defmacro unless (test :rest body)
  `(if (not ,test) (progn ,@body)))

(defglobal *e* 2.718281828459045)
(defglobal *gamma* 0.57721566490153286060)

(defun oddp (n)
  (and (integerp n)
       (= (mod n 2) 1)))

(defun evenp (n)
  (and (integerp n)
       (= (mod n 2) 0)))

(defun zerop (n)
  (= n 0))

(defun square (n)
  (* n n))

(defun set-aref1 (val mat i j)
  (set-aref val mat (- i 1) (- j 1)))

(defun aref1 (mat i j)
  (aref mat (- i 1) (- j 1)))

(defun mult-scalar-mat (s mat)
  (let ((m (elt (array-dimensions mat) 0))
        (n (elt (array-dimensions mat) 1)))
    (for ((i 1 (+ i 1)))
         ((> i m) mat)
         (for ((j 1 (+ j 1)))
              ((> j n))
              (set-aref1 (* s (aref1 mat i j)) mat i j)))))

(defun matrix-ident (n)
  (let ((m (create-array (list n n) 0)))
    (for ((i 1 (+ i 1)))
         ((> i n) m)
         (set-aref1 1 m i i))))

(defun square-matrix-p (x)
  (let ((dim (array-dimensions x)))
    (and (= (length dim) 2)
         (= (elt dim 0)(elt dim 1)))))

(defun tr (x)
  (unless (square-matrix-p x)
    (error "tr require square matrix" x))
  (let ((l (elt (array-dimensions x) 0)))
    (for ((i 1 (+ i 1))
          (y 0))
         ((> i l) y)
          (setq y (+ (aref1 x i i) y)))))

;;位置(r,s)に対しての小行列
(defun sub-matrix (x r s)
  (let* ((m (elt (array-dimensions x) 0))
         (n (elt (array-dimensions x) 1))
         (y (create-array (- m 1)(- n 1) 0)))
    (for ((i 0 (+ i 1)))
         ((>= i m) y)
         (for ((j 0 (+ j 1)))
              ((>= j n))
              (cond ((and (< i r)(< j s)) (set-aref1 (aref1 x i j) y i j))
                    ((and (> i r)(< j s)) (set-aref1 (aref1 x i j) y (- i 1) j))
                    ((and (< i r)(> j s)) (set-aref1 (aref1 x i j) y i (- j 1)))
                    ((and (> i r)(> j s)) (set-aref1 (aref1 x i j) y (- i 1) (- j 1)))
                    ((and (= i r)(= j s)) nil))))))

(defun det (x)
  (unless (square-matrix-p x) 
    (error "det require square matrix" x))
  (let ((m (elt (array-dimensions x) 0)))
    (det1 x m)))



(defun det1 (x m)
  (if (= m 2)
      (- (* (aref1 x 1 1)(aref1 x 2 2))
         (* (aref1 x 1 2)(aref1 x 2 1)))
      (for ((i 1 (+ i 1))
            (y 0))
           ((> i m) y)
           (setq y (+ (* (sign (+ i 1))
                         (aref1 x i 1)
                         (det1 (sub-matrix x i 1) (- m 1)))
                      y)))))

(defun sign (x)
  (expt -1 x))

(defun transpose (x)
  (let* ((m (elt (array-dimensions x) 0))
         (n (elt (array-dimensions x) 1))
         (y (create-array (list n m) 0)))
    (for ((i 1 (+ i 1)))
         ((> i m) y)
         (for ((j 1 (+ j 1)))
              ((> j n))
              (set-aref1 (aref1 x i j) y j i)))))

;;逆行列　クラメルの公式。効率はよくないが5次元程度なら問題なし。
(defun inv (x)
  (unless (square-matrix-p x) 
    (error "inv require square matrix" x))
  (let ((m (elt (array-dimensions x) 0))
        (n (elt (array-dimensions x) 1)))
    (if (> m 2)
        (inv1 x m)
        (inv0 x m))))

(defun inv0 (x m)
  (let ((mat (create-array (list m m) 0))
        (d (det x)))
    (when (= d 0) (error "inv determinant is zero" x))
    (cond ((= m 1)
           (set-aref1 (quotient 1 d) mat 1 1) mat)
          (t (set-aref1 (aref1 x 2 2) mat 1 1)
             (set-aref1 (aref1 x 1 1) mat 2 2)
             (set-aref1 (- (aref1 x 1 2)) mat 1 2)
             (set-aref1 (- (aref1 x 2 1)) mat 2 1)
             (mult-scalar-mat (quotient 1 d) mat)))))

(defun inv1 (x m)
  (let ((d (det x)))
    (when (= d 0) (error "inv determinant is zero" x))
    (* (quotient 1 d) (transpose (inv2 x m)))))

(defun inv2 (x m)
  (let ((y (create-array (list m m) 0)))
    (for ((i 1 (+ i 1)))
        ((> i m) y)
        (for ((j 1 (+ j 1)))
            ((> j m))
            (set-aref1 y i j (* (sign (+ i j)) (det (sub-matrix x i j))))))))


;;リストlsに関数fを適用した値の総和を求める。
(defun sum (f ls)
  (if (null ls)
      0
      (+ (funcall f (car ls)) (sum f (cdr ls)))))

;;lsの各要素について関数fを適用してその積を求める。
(defun product (f ls)
  (if (null ls)
      1
      (* (funcall f (car ls)) (product f (cdr ls)))))


;;リストlsの要素すべてについて関数ｆが成り立つか？
(defun for-all (f ls)
  (cond ((null ls) t)
        ((not (funcall f (car ls))) nil)
        (t (for-all f (cdr ls)))))

;;リストlsの少なくとも1つに関数ｆが成り立つか？
(defun at-least (f ls)
  (cond ((null ls) nil)
        ((funcall f (car ls)) t)
        (t (at-least f (cdr ls)))))

;;ガウスの素数定理によりＸ以下の素数の概数を返す。
(defun gauss-primes (x)
  (quotient x (log x)))



;;; ;;ｍとｎが互いに素であればt そうでなければnil
(defun coprimep (m n)
  (= (gcd  m n) 1))

;;ｍがｎで割り切れるかどうか。割り切れればt そうでなければnil
;; n|m 相当
(defun divisiblep (m n)
  (and (integerp m)
       (integerp n)
       (= (mod m n) 0)))

;;ｍとｎが法ａで合同かどうか。合同ならt そうでなければnil
(defun eqmodp (m n a)
  (= (mod m a) (mod n a)))

;;素数判定 10^18より小さければ決定的判定アルゴリズム、大きければラビンミラー法で判定
(defun primep (n)
  (if (< n 1000000000000000000)
      (deterministic-prime-p n)
      (rabin-miller-p n)))

;;決定的方法でnが素数であるかどうかを返す。
(defun deterministic-prime-p (n)
  (labels ((iter (x y)
                 (cond ((> x y) t)
                       ((divisiblep n x) nil)
                       ((= x 2) (iter 3 y))
                       (t (iter (+ x 2) y)))))
    (if (< n 2)
        nil
        (iter 2 (sqrt n)))))


;;π(n) n以下の素数の数を返す。
(defun primepi (n)
  (labels ((iter (x y)
                 (cond ((> x n) y)
                       ((primep x) (iter (+ x 1) (+ y 1)))
                       (t (iter (+ x 1) y)))))
    (iter 2 0)))

;;約数の個数を返す。素因数分解をして計算している。τ(n)
(defun tau (n)
  (labels ((iter (ls m)
                 (if (null ls)
                     m
                     (iter (cdr ls) (* (+ (elt (elt ls 0) 1) 1) m)))))
    (if (= n 1)
        1
        (iter (factorize n) 1))))

;;リュービルのλ関数の下請け
(defun expt-1 (n)
  (if (oddp n) 
      -1
      1))

;;リュービルのλ関数本体
(defun liouville-lambda (n)
  (expt-1 (omega n)))

;;リュービルのλ関数の下請け
(defun omega (n)
  (if (= n 1)
      0
      (sum (lambda (x) (elt x 1)) (factorize n))))



;;リュービルのλ関数の性質を調べるためのもの
;;約数を与えるとその総和は平方数なら１、その他の場合は０。
(defun g (n)
  (sum #'liouville-lambda (divisors n)))


;;素因数分解約数の和を求める。
(defun sigma2 (ls)
  (let ((p (elt ls 0))
        (k (elt ls 1)))
    (quotient (- (expt p (+ k 1)) 1) (- p 1))))

;;約数の和を求める関数
;;σ(n) 相当
(defun sigma (n)
  (cond ((< n 1) nil)
        ((= n 1) 1)
        (t (product #'sigma2 (factorize n)))))


;;完全数かどうかをσ関数を利用して計算。
;;完全数ならt そうでなければnil
(defun perfectp (n)
  (= (sigma n) (* 2 n)))

;;メルセンヌ素数を計算して返す。
;;2^p -1 
(defun mersenne (p)
  (- (expt 2 p) 1))

;;ダブル完全数 
(defun double-perfect-number-p (n)
  (= (sigma n) (* 3 n)))

;;n以下のダブル完全数をリストにして返す。
(defun find-double-perfect (n)
  (labels ((iter (m ls)
                 (cond ((> m n) ls)
                       ((double-perfect-number-p m) (iter (+ m 1) (cons m ls)))
                       (t (iter (+ m 1) ls)))))
    (iter 1 '())))

;;フェルマの小定理を利用して確率的に素数判定をする。
;;素数ならt そうでなければnil
;;底を2～n-1まで乱数生成して10回試行する。
(defun fermatp (n)
  (labels ((iter (m)
           (cond ((< m 1) t)
                 ((not (= 1 (gaussmod (+ (random (- n 2)) 1) (- n 1) n))) nil)
                 (t (iter (- m 1))))))
    (iter 10)))

;;メルセンヌ数(2^p-1)に対するルカス-テスト
;;素数ならば#t そうでなければ#f
(defun lucasp (p)
  (labels ((iter (n i m)
                 (cond ((and (= i (- p 1)) (zerop (mod n m))) t)
                       ((and (= i (- p 1)) (not (zerop (mod n m)))) nil)
                       (t (iter (mod (- (expt n 2) 2) m) (+ i 1) m)))))
    (cond ((< p 2) nil)
          ((= p 2) t)
          (t (iter 4 1 (mersenne p))))))

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

;;ラビンミラーテスト
;;テスト対象のｎを n^k * q に分解する。
(defun rm1 (n)
  (labels ((iter (k q)
                 (if (oddp q)
                     (list k q)
                     (iter (+ k 1) (div q 2)))))
    (iter 0 (- n 1))))

;;ラビンミラーテスト条件１
;;合成数ならt 素数ならnil
(defun rm2 (a q n)
  (not (= (gaussmod a q n) 1)))

;;ラビンミラーテスト条件２
;;合成数ならt 素数ならnil
(defun rm3 (a k q n)
  (labels ((iter (i)
                 (cond ((>= i k) t)
                       ((= (gaussmod a (* (expt 2 i) q) n) -1) nil)
                       (t (iter (+ i 1))))))
    (iter 0)))


;;ラビンミラーテスト
;;nについて底aで条件１，２をテスト
;;合成数ならt 素数ならnil
(defun rm4 (n a)
  (let* ((ls (rm1 n))
         (k (elt ls 0))
         (q (elt ls 1)))
    (and (rm2 a q n)
         (rm3 a k q n))))

;;ラビンミラーテスト
;;関数本体
;;底を2～n-1(ただし32767以下)まで乱数で発生させ10回試行する。
;;素数ならt 合成数ならnil
;;擬素数である確率は 0.25^10 およそ0.000095% 
(defun rabin-miller-p (n)
  (labels ((iter (m)
                 (cond ((< m 1) nil)
                       ((rm4 n (+ (random (min (- n 2) 32767)) 1)) t)
                       (t (iter (- m 1))))))
    (if (= n 2)
        t
        (not (iter 10)))))



;;繰り返し二乗法によるmod計算で結果を法ｍとした場合 -m/2～m/2 で表す。
(defun gaussmod (a k m)
  (let ((k1 (expmod a k m)))
    (cond ((and (> k1 0) (> k1 (quotient m 2)) (< k1 m)) (- k1 m))
          ((and (< k1 0) (< k1 (- (quotient m 2))) (> k1 (- m))) (+ k1 m))
          (t k1))))

;;双子素数をn～n+mまで探す。
(defun twin-primes (n m)
  (labels ((iter (i j ls)
                 (cond ((> i j) (reverse ls))
                       ((and (primep i) (primep (+ i 2)))
                        (iter (+ i 2) j (cons (list i (+ i 2)) ls)))
                       (t (iter (+ i 2) j ls)))))
    (cond ((<= n 2) (iter 3 (+ n m) '()))
          ((evenp n) (iter (+ n 1) (+ n m) '()))
          (t (iter n (+ n m) '())))))



;;約数を求めてリストにして返す。
(defun divisors (n)
  (labels ((iter (m o ls)
                 (cond ((> m o) ls)
                       ((divisiblep n m) (iter (+ m 1) o (cons m ls)))
                       (t (iter (+ m 1) o ls)))))
    (cond ((not (integerp n)) (error "divisors require natural number" n))
          ((< n 1) (error "divisors require natural number" n))
          ((= n 1) '(1))
          (t (cons n (iter 1 (ceiling (quotient n 2))'()))))))


;;nを素因数分解する。指数形式ではなく単純に素数を並べたリストで返す。
;;n<0の場合には#f、n=0,n=1の場合には'(0),'(1)を返す。
(defun prime-factors (n)
  (labels ((iter (p x ls z)
                 (cond ((> p z) (cons x ls))
                       ((= (mod x p) 0) 
                        (let ((n1 (div x p)))
                          (iter 2 n1 (cons p ls) (isqrt n1))))
                       ((= p 2) (iter 3 x ls z))
                       (t (iter (+ p 2) x ls z)))))
    (cond ((< n 0) nil)
          ((< n 2) (list n))
          (t (iter 2 n '() (isqrt n))))))

;;nを素因数分解して標準形式にして返す。p^a + q^b + r^c ((p a)(q b)(r c))
(defun factorize (n)
  (labels ((iter (ls p n mult)
                 (cond ((null ls) (cons (list p n) mult))
                       ((= (car ls) p) (iter (cdr ls) p (+ n 1) mult))
                       (t (iter (cdr ls) (car ls) 1 (cons (list p n) mult))))))
    (let ((ls (prime-factors n)))
      (iter (cdr ls) (car ls) 1 '()))))

;;オイラーのφ関数
;;ｎ以下の数でｎと互いに素であるものの個数を返す。
;;素因数分解により計算している。 φ(n=p^a q^b r^c) = n(1-1/p)(1-1/q)(1-1/r)
(defun phi (n)
  (if (= n 1)
      1
      (convert 
        (* n (product 
               (lambda (ls) (- 1 (quotient 1 (elt ls 0))))
               (factorize n)))
        <integer>)))

;;原始根の判定
;;ｎが素数ｐを法として原始根であるなら#t 
;;素数なら必ず存在するが条件がそろっていなければ nilが返る。
(defun primitive-root-p (n p)
  (labels ((iter (i)
                 (cond ((>= i (- p 1)) t)
                       ((= (expmod n i p) 1) nil)
                       (t (iter (+ i 1))))))
    (and (iter 1)
         (= (expmod n (- p 1) p) 1))))

;;sicp
;;繰り返し二乗法によるmod計算。
;; a^n (mod m)を計算する。SICPより借用。
(defun expmod (a n m)
  (cond ((= 0 n) 1)
        ((evenp n)
         (mod (square (expmod a (div n 2) m)) m))
        (t
          (mod (* a (expmod a (- n 1) m)) m))))

;;素数ｐの最小の原始根を返す。
;;ｐの任意の原始根に成り立つ定理を試すのに一番小さな原始根を使うこととした。
;;計算が楽なので。
(defun primitive-root (p)
  (labels ((iter (n)
                 (cond ((> n p) nil)
                       ((primitive-root-p n p) n)
                       (t (iter (+ n 1))))))
    (iter 2)))

;;指数の計算
;;原始根ｒを底として素数ｐを法としたaに対する指数を求める
;;指数は必ず存在するが与えられた値が条件に合わなければnilが返る。
(defun ind (r a p)
  (labels ((iter (i)
                 (cond ((> i p) nil)
                       ((= (expmod r i p) a) i)
                       (t (iter (+ i 1))))))
    (iter 0)))



;;高度合成数
(defun highly-composite-number-p (n)
  (cond ((<= n 0) nil)
        ((= n 1) t)
        (t (> (tau n) (max-tau (- n 1) 0)))))

(defun max-tau (n m)
  (let ((x (tau n)))
    (cond ((= n 1) m)
          ((> x m) (max-tau (- n 1) x))
          (t (max-tau (- n 1) m)))))