share/pslq/pslq_integer_relation.lisp

   1 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   2 ;;
   3 ;; Integer relations package
   4 ;;
   5 ;; Copyright: 2006, Andrej Vodopivec
   6 ;; Licence:   GPL
   7 ;;
   8 ;;   Integer relation package implements pslq algorithm to find integer
   9 ;;   relations between bigfloat numbers.
  10 ;;
  11 ;; Description:
  12 ;;
  13 ;;   For a given vector x of floating point numbers we want to find a
  14 ;;   vector of integers m such that m.x=0 (in given precision).
  15 ;;
  16 ;;   This implementation of PSLQ, to judge by the initialization of s
  17 ;;   (before it was changed, anyway) and the calculation of the bound M,
  18 ;;   may be derived from:
  19 ;;
  20 ;;     D.H. Bailey and S. Plouffe. "Recognizing Numerical Constants."
  21 ;;     http://www.cecm.sfu.ca/organics/papers/bailey/paper/html/paper.html
  22 ;;     See the section titled "The PSLQ Integer Relation Algorithm".
  23 ;;     The version at the above URL appears to have been published as
  24 ;;     part of a workshop proceedings on Dec. 11, 1995.
  25 ;;
  26 ;;     A PDF version, possibly identical, dated Dec. 15, 1995:
  27 ;;     https://www.davidhbailey.com/dhbpapers/recog.pdf
  28 ;;     which mentions:
  29 ;;     Canadian Mathematical Society, vol. 20 (1997), pg. 73-88.
  30 ;;
  31 ;;   Other sources, which might or might not have contributed to this implementation.
  32 ;;
  33 ;;     D.H.Bailey. "Integer relation detection,"
  34 ;;     Computing in Science and Engineering, Jan-Feb, 2000, pg. 24-28.
  35 ;;     Preprint published as "Integer Relation Detection and Lattice Reduction":
  36 ;;     https://www.davidhbailey.com/dhbpapers/pslq-cse.pdf
  37 ;;
  38 ;;     H.R.P. Ferguson, D.H. Bailey, S. Arno. "Analysis of PSLQ,
  39 ;;     an integer relation finding algorithm."
  40 ;;     Mathematics of Computation, vol. 68, no. 225 (Jan 1999), pg. 351-369.
  41 ;;     https://www.davidhbailey.com/dhbpapers/cpslq.pdf
  42 ;;
  43 ;;     Armin Straub. "A gentle introduction to PSLQ."
  44 ;;     https://arminstraub.com/downloads/math/pslq.pdf
  45 ;;
  46 ;;     Jingwei Chen, Damien Stehle, Gilles Villard.
  47 ;;     "A New View on HJLS and PSLQ: Sums and Projections of Lattices."
  48 ;;     Proceedings of ISSAC '13.
  49 ;;     https://arcnl.org/jchen/download/[CSV13].pdf
  50 ;;
  51 ;;     Yong Feng, Jingwei Chen, Wenyuan Wu.
  52 ;;     "The PSLQ algorithm for empirical data."
  53 ;;     Mathematics of Computation
  54 ;;     DOI: 10.1090/mcom/3356
  55 ;;     https://www.ams.org/journals/mcom/2019-88-317/S0025-5718-2018-03356-7/mcom3356_AM.pdf
  56 ;;
  57 ;;;;;;;;;;;;;;;;;;
  58
  59 (in-package :maxima)
  60 (macsyma-module pslq_integer_relation)
  61
  62 ($put '$pslq_integer_relation 1.0 '$version)
  63
  64 (defmfun $pslq_integer_relation (l)
  65   (if (not ($listp l))
  66       (merror (intl:gettext "pslq_integer_relation: argument must be a list; found: ~M") l))
  67   (let* ((n ($length l)) (l (cdr l)) (a (make-array n)) ($float2bf t))
  68     (loop for i from 0 to (1- n) do
  69          (let ((elm (car l)))
  70            (if ($floatnump elm)
  71                (setf (aref a i) elm)
  72                (progn
  73                  (setq elm ($bfloat elm))
  74                  (if (not ($freeof '$%i elm))
  75                      (return-from $pslq_integer_relation nil))
  76                  (if (not ($bfloatp elm))
  77                      (merror (intl:gettext "pslq_integer_relation: element can't be converted to bigfloat: ~M") elm))
  78                  (setf (aref a i) elm))))
  79          (setq l (cdr l)))
  80     (let ((ans (pslq-integer-relations a n)))
  81       (if ans
  82           `((mlist simp) ,@ans)
  83           nil))))
  84
  85 (defmvar $pslq_precision nil)
  86 (defmvar $pslq_threshold nil)
  87 (defmvar $pslq_depth nil)
  88 (defmvar $pslq_status 0)
  89
  90 (defmvar $pslq_fail_norm nil)
  91
  92 (defvar *pslq-debugging* nil)
  93
  94 (defun pslq-mabs (x)
  95   (if (mlsp x 0) (m- x) x))
  96
  97 (defun pslq-nearest-integer (x)
  98   (let ((nx ($entier x)))
  99     (if (mlsp 0.5 (m- x nx))
 100         (1+ nx)
 101         nx)))
 102
 103 (defun my-write-lisp-array (a s &rest args)
 104   (declare (ignore args))
 105   (let ((d (array-dimensions a)))
 106     (if (= (length d) 1)
 107       (progn (loop for i from 0 to (1- (first d))
 108                    do (format s " ~18a" (aref a i)))
 109              (write-char #\newline s))
 110       (progn (loop for i from 0 to (1- (first d))
 111                    do (loop for j from 0 to (1- (second d))
 112                             do (format s " ~18a" (aref a i j)))
 113                       (write-char #\newline s))
 114              (write-char #\newline s)))))
 115
 116 (defun pslq-integer-relations (x n)
 117   (let ((A (make-array `(,n ,n) :initial-element 0))
 118         (B (make-array `(,n ,n) :initial-element 0))
 119         (H (make-array `(,n ,(1- n)) :initial-element 0))
 120         (gamma (sqrt (/ 4 3)))
 121         (s (make-array `(,n) :initial-element 0))
 122         (y (make-array `(,n) :initial-element 0))
 123         ($pslq_precision (if $pslq_precision $pslq_precision (power 10 (- $fpprec 2))))
 124         ($pslq_threshold (if $pslq_threshold $pslq_threshold (power 10 (- 2 $fpprec))))
 125         ($pslq_depth (if $pslq_depth $pslq_depth (* 20 n)))
 126         (tt))
 127
 128     (when *pslq-debugging*
 129       (format t "PSLQ-INTEGER-RELATIONS: n = ~a, gamma = ~a, pslq_precision = ~a, pslq_threshold =  ~a, pslq_depth = ~a~%" n gamma $pslq_precision $pslq_threshold $pslq_depth))
 130
 131     ;; Initialize A, B and s
 132     (loop for k from 0 to (1- n) do
 133          (setf (aref A k k) 1)
 134          (setf (aref B k k) 1)
 135          (loop for j from k to (1- n) do
 136               (setf (aref s k) (m+ (aref s k) (power (aref x j) 2))))
 137          (setf (aref s k) ($sqrt (aref s k))))
 138     (setf tt (aref s 0))
 139
 140     ;; Initialize y, normalize s
 141     (loop for k from 0 to (1- n) do
 142          (setf (aref y k) (m// (aref x k) tt))
 143          (setf (aref s k) (m// (aref s k) tt)))
 144
 145     ;; Initialize H
 146     (loop for i from 0 to (1- n) do
 147          (if (< i (- n 1))
 148              (setf (aref H i i) (m// (aref s (1+ i)) (aref s i))))
 149          (loop for j from 0 to (1- i) do
 150               (setf (aref H i j) (m- (m// (m* (aref y i) (aref y j))
 151                                           (m* (aref s j) (aref s (1+ j))))))))
 152
 153 (when *pslq-debugging*
 154   (print "PSLQ-INTEGER-RELATIONS: just before initial reduce h:") (write-char #\newline)
 155   (print "A =") (write-char #\newline) (my-write-lisp-array A *terminal-io* '$comma 'text)
 156   (print "B =") (write-char #\newline) (my-write-lisp-array B *terminal-io* '$comma 'text)
 157   (print "s =") (write-char #\newline) (my-write-lisp-array s *terminal-io* '$comma 'text)
 158   (print "y =") (write-char #\newline) (my-write-lisp-array y *terminal-io* '$comma 'text)
 159   (print "H =") (write-char #\newline) (my-write-lisp-array H *terminal-io* '$comma 'text))
 160
 161     ;; Perform reduction on H, update A, B, y
 162     (loop for i from 1 to (- n 1) do
 163          (loop for j from (1- i) downto 0 do
 164               (setq tt (pslq-nearest-integer (m// (aref H i j) (aref H j j))))
 165               (setf (aref y j) (m+ (aref y j) (m* tt (aref y i))))
 166               (loop for k from 0 to j do
 167                     (let* ((bar (aref H i k)) (baz tt) (quux (aref H j k))
 168                            (foo (m- bar (m* baz quux))))
 169                       (when *pslq-debugging*
 170                         (format t "PSLQ-INTEGER-RELATIONS: assign ~a = ~a - ~a * ~a to H[~d, ~d]~%" foo bar baz quux i k))
 171                    (setf (aref H i k) foo)))
 172               (loop for k from 0 to (1- n) do
 173                    (setf (aref A i k) (m- (aref A i k) (m* tt (aref A j k))))
 174                    (setf (aref B k j) (m+ (aref B k j) (m* tt (aref B k i)))))))
 175
 176     (do ((r 1 (1+ r))) ((= r $pslq_depth))
 177       (let ((m 0) (mm 0) (s 1))
 178
 179 (when *pslq-debugging*
 180   (print "PSLQ-INTEGER-RELATIONS: just before bound check:") (write-char #\newline)
 181   (print "A =") (write-char #\newline) (my-write-lisp-array A *terminal-io* '$comma 'text)
 182   (print "B =") (write-char #\newline) (my-write-lisp-array B *terminal-io* '$comma 'text)
 183   (print "y =") (write-char #\newline) (my-write-lisp-array y *terminal-io* '$comma 'text)
 184   (print "H =") (write-char #\newline) (my-write-lisp-array H *terminal-io* '$comma 'text))
 185
 186         ;; Find the bound M
 187         (let ((maxNorm 0))
 188           (loop for j from 0 to (1- n) do
 189                (let ((absHj 0))
 190                  (loop for i from 0 to (- n 2) do
 191                       (if (mlsp absHj (pslq-mabs (aref H j i)))
 192                           (setq absHj (aref H j i))))
 193                  (if (mlsp maxNorm absHj)
 194                      (setq maxNorm absHj))))
 195           (setq $pslq_fail_norm (m// 1 maxNorm))
 196
 197           ;; Check to see if we have a relation
 198           (loop for j from 0 to (1- n) do
 199                (if (mlsp (pslq-mabs (aref y j)) $pslq_threshold)
 200                    (progn
 201                      (let ((ans ()))
 202                        (loop for i from 0 to (1- n) do
 203                             (setq ans (append ans `(,(aref B i j)))))
 204                        (setq $pslq_status 1)
 205                        (return-from pslq-integer-relations ans)))))
 206
 207           ;; Check to see if we exhausted the precision
 208           (loop for i from 0 to (1- n) do
 209                (loop for j from 0 to (1- n) do
 210                     (if (mlsp $pslq_precision (pslq-mabs (aref A i j)))
 211                         (progn
 212                           (setq $pslq_status 2)
 213                           (return-from pslq-integer-relations nil)))))
 214           )
 215
 216         ;; Find maximal value in H
 217         (loop for i from 0 to (m- n 2) do
 218              (setq s (* gamma s))
 219              (if (mlsp mm (m* s (pslq-mabs (aref H i i))))
 220                  (progn
 221                    (setf m i)
 222                    (setf mm (m* (expt gamma i) (pslq-mabs (aref H i i)))))))
 223
 224         ;; Swap entries in y
 225         (rotatef (aref y m) (aref y (1+ m)))
 226         ;; Swap rows in A and columns in B
 227         (loop for i from 0 to (- n 1) do
 228              (rotatef (aref A m i) (aref A (1+ m) i))
 229              (rotatef (aref B i m) (aref B i (1+ m))))
 230         ;; Swap rows in H
 231         (loop for i from 0 to (- n 2) do
 232              (rotatef (aref H m i) (aref H (1+ m) i)))
 233
 234         ;; Update H
 235         (if (< m (- n 2))
 236             (let* ((t0 ($sqrt (m+ (power (aref H m m) 2) (power (aref H m (1+ m)) 2))))
 237                    (t1 (m// (aref H m m) t0))
 238                    (t2 (m// (aref H m (1+ m)) t0))
 239                    (t3) (t4))
 240               (loop for i from m to (1- n) do
 241                    (setq t3 (aref H i m))
 242                    (setq t4 (aref H i (1+ m)))
 243                    (setf (aref H i m) (m+ (m* t1 t3) (m* t2 t4)))
 244                    (setf (aref H i (1+ m)) (m- (m* t1 t4) (m* t2 t3))))))
 245
 246         ;; Block reduction on H
 247         (loop for i from (1+ m) to (1- n) do
 248              (loop for j from (min (1- i) (1+ m)) downto 0 do
 249                   (setq tt (pslq-nearest-integer (m// (aref H i j) (aref H j j))))
 250                   (setf (aref y j) (m+ (aref y j) (m* tt (aref y i))))
 251                   (loop for k from 0 to j do
 252                        (setf (aref H i k) (m- (aref H i k) (m* tt (aref H j k)))))
 253                   (loop for k from 0 to (1- n) do
 254                        (setf (aref A i k) (m- (aref A i k) (m* tt (aref A j k))))
 255                        (setf (aref B k j) (m+ (aref B k j) (m* tt (aref B k i)))))))
 256     ))
 257     (setq $pslq_status 3)
 258     (return-from pslq-integer-relations nil) ))