src/ifactor.lisp

   1 ;;;;
   2 ;;;;                      ~*~  IFACTOR  ~*~
   3 ;;;;
   4 ;;;;  Maxima integer factorization package.
   5 ;;;;
   6 ;;;;  Copyright: 2005-2015 Andrej Vodopivec, Volker van Nek
   7 ;;;;  Licence  : GPL
   8 ;;;;
   9 ;;;;   - ifactors     : factorization of integers
  10 ;;;;   - primep       : probabilistic primality test
  11 ;;;;   - next_prime   : smallest prime > n
  12 ;;;;   - prev_prime   : greatest prime < n
  13 ;;;;   - primes       : list of primes
  14 ;;;;   - power_mod    : fast modular powers
  15 ;;;;   - inv_mod      : modular inverse
  16
  17
  18 (in-package :maxima)
  19 (macsyma-module ifactor)
  20
  21 (defmvar $save_primes nil "Save primes found." boolean)
  22
  23 (defmvar $primep_number_of_tests 25 "Number of Miller-Rabin tests." fixnum)
  24
  25 (defmvar $pollard_rho_limit 16000 "Limit for pollard-rho factorization depth." fixnum)
  26 (defmvar $pollard_pm1_limit 25000 "Limit for pollard-p1 factorization depth." fixnum)
  27
  28 (defmvar $pollard_rho_tests 5 "Number of pollard-rho rounds." fixnum)
  29 (defmvar $pollard_pm1_tests 3 "Number of pollard-p-1 rounds." fixnum)
  30
  31 (defmvar $pollard_rho_limit_step 1000 "Step for pollard-rho factorization limit." fixnum)
  32 (defmvar $pollard_pm1_limit_step 5000 "Step for pollard-p-1 factorization limit." fixnum)
  33
  34 (defmvar $ecm_number_of_curves 50 "Number of curves tried in one round of ecm." fixnum)
  35 (defmvar $ecm_limit       200  "Starting smootheness limit for ecm method." fixnum)
  36 (defmvar $ecm_max_limit 51199  "Maximum smootheness for ecm method." fixnum)
  37 (defmvar $ecm_limit_delta 200  "Increase smoothness limit for ecm method after each round." fixnum)
  38
  39 (defmvar $ifactor_verbose nil "Display factorization steps." boolean)
  40 (defmvar $factors_only nil  "Return a list of factors only." boolean)
  41
  42 (defun number-of-digits (n)
  43   (length (format nil "~d" n)))
  44
  45 ;;; List of primes up to *largest-small-prime*
  46
  47 (defvar *largest-small-prime*
  48   (car (last *small-primes*)))
  49
  50 ;;; List of numbers which have already been tested and are
  51 ;;; primes > *largest-small-prime* (only used if $save_primes is true!).
  52
  53 (defvar *large-primes* ())
  54
  55 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  56 ;;;                                                                       ;;;
  57 ;;;        ~*~  IMPLEMENTATION OF FACTORIZATION METHODS   ~*~             ;;;
  58 ;;;                                                                       ;;;
  59 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  60
  61 (defun get-factor-list (n)
  62   (when $ifactor_verbose
  63     (format t "~%Starting factorization of n = ~d~%" n))
  64   (init-prime-diffs 640000)
  65   (multiple-value-bind (large-part factor-list) (get-small-factors n)
  66     (if (> large-part 1)
  67         (append (convert-list (get-large-factors large-part)) factor-list)
  68         factor-list)))
  69
  70 (defmfun $ifactors (n)
  71   (unless (and (or (integerp n)
  72                    (and ($integerp n)
  73                         (setq n ($fix n)) ))
  74                (plusp n) )
  75     (merror (intl:gettext "ifactors: argument must be a positive integer; found: ~M") n))
  76   (let* (($intfaclim)
  77      (factor-list (get-factor-list n))
  78      (factor-list (if $factors_only
  79               (mapcar #'car factor-list)
  80               (mapcar #'(lambda (u) `((mlist simp) ,(car u) ,(cadr u))) factor-list))))
  81     ($sort `((mlist simp) ,@factor-list))))
  82
  83 ;; cfactor is the function used by maxima to factor integers
  84 ;; (used form outside ifactor package)
  85
  86 (defun cfactor (x)
  87   (cond ((null $factorflag) (return-from cfactor (list x 1)))
  88     ((floatp x) (rat-error "`factor' given floating arg"))
  89     ((pzerop x) (return-from cfactor (list (pzero) 1)))
  90     ((equal x -1) (return-from cfactor (list -1 1)))
  91     ((minusp x) (return-from cfactor (cons -1 (cons 1 (cfactor (- x))))))
  92     ((< x 2) (return-from cfactor (list x 1)))
  93     (t (let* ((factor-list (get-factor-list x))
  94           (factor-list (sort factor-list #'< :key #'car))
  95           (ans ()))
  96          (dolist (fac factor-list ans)
  97            (setq ans (cons (car fac) (cons (cadr fac) ans))))))))
  98
  99 ;;; we need to keep a list of differences between consecutive primes
 100 ;;; for trial division, for stage 2 of ecm and the big prime variation of Pollard p-1
 101
 102 (defvar *prime-diffs* (make-array 100000 :element-type 'fixnum :adjustable t :initial-element 0)
 103   "array of differences between consecutive primes")
 104
 105 (defvar *prime-diffs-limit* 1
 106   "biggest prime in *prime-diffs")
 107
 108 (defvar *prime-diffs-maxindex* 1
 109   "index of biggest valid prime difference")
 110
 111 (defvar *prime-diffs-maxdiff* 2
 112   "maximum difference between consecutive primes in *prime-diffs*")
 113
 114 ;;; factor out primes < *prime-diffs-limit* by trial division
 115 ;;; the array gets filled by a call to init-prime-diffs in get-factor-list
 116
 117 (defun get-small-factors (n)
 118   (when (= n 1)
 119     (return-from get-small-factors (values 1 nil)))
 120   (when (< n 4)                 ;n = 2 or 3
 121     (return-from get-small-factors (values 1 `((,n 1)))))
 122   (let (factors)
 123     ;; first divide off the even part
 124     (loop with deg = 0
 125        while (and (> n 1) (evenp n)) do
 126      (setq n (ash n -1))                ; divide n by 2
 127      (incf deg)                 ; and increment the exponent
 128        finally
 129      (when (plusp deg)
 130        (push `(2 ,deg) factors)
 131        (when $ifactor_verbose (format t "Factoring out 2: 2 (degree:~A)~%" deg))))
 132     (when (= n 1)
 133       (return-from get-small-factors (values 1 factors))) ; n was a power of 2
 134     ;; now use the *prime-diffs* array for trial-factoring
 135     (loop for i from 0 to *prime-diffs-maxindex*
 136        and d = 3 then (+ d (aref *prime-diffs* i))
 137        do
 138      (when (> (* d d) n)
 139        ;;(push `(,n 1) factors) replaced by workaround next line,
 140       (push (list n 1) factors) ;; see bug report 3510983 (van_nek, 2012-03-27)
 141        (when $ifactor_verbose  (format t "small prime cofactor: ~A~%" n))
 142        (return-from get-small-factors (values 1 factors)))
 143      (loop with deg = 0
 144         while (and (> n 1) (zerop (mod n d))) do
 145           (setq n (truncate n d))
 146           (incf deg)
 147         finally
 148           (when (plusp deg)
 149         (push `(,d ,deg) factors)
 150         (when $ifactor_verbose (format t "Factoring out small prime: ~A (degree:~A)~%" d deg))))
 151      (when (= n 1)
 152        (return-from get-small-factors (values 1 factors))))
 153     (return-from get-small-factors (values n factors))))
 154
 155 ;;; get-large-factors returns the list of factors of integer n (n has
 156 ;;; no small factor at this tage)
 157
 158 (defun get-large-factors (n)
 159   (if (primep n)
 160       (progn
 161     (when $ifactor_verbose (format t "========> Prime factor: ~d~%~%" n))
 162     (list n))
 163       (get-large-factors-1 n)))
 164
 165 (defun get-large-factors-1 (n)
 166   (let ((f (get-one-factor n)))
 167     (if (= f n)
 168     (progn
 169       (when $ifactor_verbose (format t "Warning: could not find factors of composite:~%~A~%" n))
 170       (list n))
 171     (append (get-large-factors f) (get-large-factors (/ n f))))))
 172
 173 (defun get-one-factor (n)
 174   (when $ifactor_verbose
 175     (format t "Factoring n = ~d~%" n))
 176   (let ((f nil)
 177     (lim_pollard $pollard_rho_limit)
 178     (lim_p-1 $pollard_pm1_limit)
 179     ($ecm_number_of_curves $ecm_number_of_curves))
 180
 181     ;; If $intfaclim is not false then we don't want to spend too much
 182     ;; time factoring integers so we return n and leave it
 183     ;; unfactored. The default value for $intfaclim is true, but most
 184     ;; functions which use integer factorization set it to false.
 185     (when $intfaclim
 186       (return-from get-one-factor n))
 187
 188     ;; first try known large primes
 189     (dolist (p *large-primes*)
 190       (when (zerop (mod n p))
 191         (return-from get-one-factor p)))
 192
 193     ;; try factoring smaller factors with pollard-rho
 194     (dotimes (i $pollard_rho_tests)
 195       (when $ifactor_verbose
 196     (format t "Pollard rho: round #~d of ~d (lim=~d)~%" (1+ i) $pollard_rho_tests lim_pollard))
 197       (setq f (get-one-factor-pollard n lim_pollard))
 198       (when (< 1 f n)
 199     (when $ifactor_verbose
 200       (format t "Pollard rho: found factor ~A (~d digits)~%" f (number-of-digits f)))
 201     (return-from get-one-factor f))
 202       (if (> lim_pollard 0)
 203       (incf lim_pollard $pollard_rho_limit_step)))
 204
 205     ;; now try factoring with pollards p-1 method
 206     (dotimes (i $pollard_pm1_tests)
 207       (when $ifactor_verbose
 208     (format t "Pollard p-1: round #~d of ~d (lim=~d)~%" (1+ i) $pollard_pm1_tests lim_p-1))
 209       (setq f (get-one-factor-p-1 n lim_p-1))
 210       (when (< 1 f n)
 211     (when $ifactor_verbose
 212       (format t "Pollard p-1: found factor ~A (~d digits)~%" f (number-of-digits f)))
 213     (return-from get-one-factor f))
 214       (when (plusp lim_pollard)
 215     (incf lim_p-1 $pollard_pm1_limit_step)))
 216
 217     ;; continue with ecm
 218     (do ()
 219     (nil)
 220       (setq f (get-one-factor-ecm n))
 221       (unless (null f)
 222     (return-from get-one-factor f))
 223       (setq $ecm_number_of_curves (+ $ecm_number_of_curves 50)))))
 224
 225 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 226 ;;;                                                                       ;;;
 227 ;;;   ~*~  IMPLEMENTATION OF POLLARDS-RHO FACTORIZATION METHOD   ~*~      ;;;
 228 ;;;                                                                       ;;;
 229 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 230
 231 ;; This is the Brent's variant of Pollard's rho method.
 232
 233 (defun get-one-factor-pollard (n lim)
 234   (let* ((x (+ (random (- n 3)) 2))
 235      (a (+ (random (- n 2)) 1))
 236      (b (+ (random (- n 5)) 1))
 237      (y x) (d 1) (r 2) (j 1) (k)
 238           (terms (integer-length n)) )
 239     (setq b (/ b (gcd a b)))
 240     (loop while (= d 1) do
 241      (setq y x)
 242      (incf j r)
 243      (dotimes (i r)
 244        (setq x (mod (+ (* a (mod (* x x) n)) b) n)))
 245      (setq k 0)
 246      (loop while (and (< k r) (equal d 1)) do
 247           (dotimes (i (min terms (- r k)))
 248         (setq x (mod (+ (* a (mod (* x x) n)) b) n))
 249         (setq d (mod (* d (- x y)) n)))
 250           (setq d (gcd d n))
 251           (incf k terms))
 252      (setq r (* 2 r))
 253      (when (< 0 lim j)
 254          (return-from get-one-factor-pollard d)))
 255     d))
 256
 257 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 258 ;;;                                                                       ;;;
 259 ;;;   ~*~  IMPLEMENTATION OF POLLARDS P-1 FACTORIZATION METHOD   ~*~      ;;;
 260 ;;;                                                                       ;;;
 261 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 262
 263 ;;; product of all primes low < p < high and all integers
 264 ;;; isqrt(low) < n < isqrt(high)
 265
 266 (defun ppexpo (low high)
 267   (declare (integer low high))
 268   (let ((x 1))
 269     (loop for i from (max 2 (1+ (isqrt low))) to (isqrt high) do
 270      (setq x (* i x)))
 271     (when (oddp low) (incf low))
 272     (loop for i from (1+ low) to high
 273        when (primep i) do (setq x (* i x)))
 274     x))
 275
 276 (defun big-prime-var (y n bound)
 277   (let ((x (make-array (1+ (ash *prime-diffs-maxdiff* -1)) :element-type 'integer :initial-element 1))
 278     (y2 (mod (* y y) n))
 279     (q 3)
 280     (count 1)
 281     (d 0)
 282     (z 0)
 283     (k 0))
 284     (loop for i from 1 to (ash *prime-diffs-maxdiff* -1) do
 285      (setf (aref x i) (mod (* y2 (aref x (1- i))) n)))
 286     (setq y (power-mod y q n))
 287     (setq z (1- y))
 288     (setq bound (min bound *prime-diffs-limit*))
 289     (loop for i from 0
 290      while (< q bound) do
 291      (setq k (aref *prime-diffs* i))
 292      (incf q k)
 293      (setq y (mod (* y (aref x (ash k -1))) n))
 294      (setq z (mod (* z (1- y)) n))
 295      (when (> (incf count) 1000)
 296        (setq count 0)
 297        (setq d (gcd z n))
 298        (when (> d 1)
 299          (return-from big-prime-var d))))
 300     d))
 301
 302 ;;; Pollard's p-1 factoring algorithm
 303 ;;; in general a prime factor p of x is found, if p-1 is
 304 ;;; a product of prime powers q^k <= lim
 305
 306 (defun get-one-factor-p-1 (n &optional (lim 16000))
 307   (declare (integer n lim))
 308   (let* ((base (+ 2 (random (- n 2))))
 309      (anz 256)
 310      (d (gcd base n)))
 311     (declare (fixnum anz)
 312          (integer base d))
 313     (when (< 1 d n) (return-from get-one-factor-p-1 d))
 314     (loop for n0 from 0 to (1- lim) by anz
 315        and ex = (ppexpo n0 (min lim (+ n0 anz))) do
 316        (setq base (power-mod base ex n))
 317        (when (<= base 1) (return-from get-one-factor-p-1 1))
 318        (setq d (gcd (1- base) n))
 319        (when (> d 1)
 320      (return-from get-one-factor-p-1 d)))
 321     (when (= d 1)
 322       (return-from get-one-factor-p-1 (big-prime-var base n *prime-diffs-limit*))))
 323   1)
 324
 325 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 326 ;;;                                                                       ;;;
 327 ;;;  ~*~  IMPLEMENTATION OF ELLIPTIC CURVE FACTORIZATION METHOD   ~*~     ;;;
 328 ;;;                                                                       ;;;
 329 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 330
 331 ;;; Elliptic curve factorization method as described in the paper
 332 ;;; R. Brent: "Factorization of the tenth and eleventh Fermat Number".
 333 ;;; The paper is in file rpb161tr.dvi.gz from
 334 ;;; ftp://ftp.comlab.ox.ac.uk/pub/Documents/techpapers/Richard.Brent/
 335 ;;;
 336 ;;; Based in part on the implementation from GAP4 FacInt package.
 337
 338 (defun init-prime-diffs (n)
 339   (when (> n *prime-diffs-limit*)
 340     (setq n (* 2 n))
 341     (when $ifactor_verbose (format t "Initializing prime diffs up to n=~d~%" n))
 342     (let ((sieve (make-array (1+ n) :element-type 'bit :initial-element 1)))
 343       (do ((p 3 ($next_prime p)))
 344       ((> p (isqrt n)))
 345     (do ((d (* 2 p) (+ d p)))
 346         ((> d n))
 347       (setf (sbit sieve d) 0)))
 348       (do ((q1 3)
 349        (i 0)
 350        (q2 5 (+ q2 2)))
 351       ((> q2 n) (setq *prime-diffs-maxindex* (1- i)))
 352     (when (= 1 (sbit sieve q2))
 353       (when (>= i (length *prime-diffs*))
 354         (setq *prime-diffs* (adjust-array *prime-diffs* (* 2 (length *prime-diffs*)) :element-type 'fixnum :initial-element 0))
 355         (when $ifactor_verbose
 356           (format t "init-prime-diffs: adjusting *prime-diffs* to size ~d~%" (* 2 (length *prime-diffs*)))))
 357       (setq *prime-diffs-limit* q2)
 358
 359       (let ((diff (- q2 q1)))
 360         (setf (aref *prime-diffs* i) diff)
 361         (when (> diff *prime-diffs-maxdiff*) (setq *prime-diffs-maxdiff* diff)))
 362       (setq q1 q2)
 363       (incf i))))))
 364
 365 ;;; modular inverse of a (modulus m)
 366 ;;;
 367 ;;; this implementation of the modular inverse is different to `invmod' in src/rat3a.lisp
 368 ;;; inv-mod returns a positive modulo or `nil' in case of a zero divisor
 369
 370 (defun inv-mod (a m)
 371   (let ((u1 0)(u2 m)(v1 1)(v2 (mod a m)) q r)
 372     (do ()((zerop v2)
 373             (if (= 1 u2) (mod u1 m) nil) )
 374       (multiple-value-setq (q r) (truncate u2 v2))
 375       (psetq u1 v1 v1 (- u1 (* q v1)))
 376       (setq u2 v2 v2 r) )))
 377
 378 (defmfun $inv_mod (a m)
 379   (unless (and (integerp a) (integerp m))
 380       (merror (intl:gettext "inv_mod: arguments must be integers; found: ~M, ~M") a m))
 381   (unless (= 0 a) (inv-mod a m)) )
 382
 383 ;;; computations on the elliptic curve:
 384 ;;; we use the elliptic curve in projective coordinates (x,y,z), but only
 385 ;;; use x and z
 386
 387 (defun ecm-product (q p1 p2 n)
 388   (let ((x1 (car p1)) (x2 (car p2))
 389     (z1 (cadr p1)) (z2 (cadr p2))
 390     (pr1) (pr2) (sq1) (sq2) (x3) (z3))
 391     (setq pr1 (mod (* (- x1 z1) (+ x2 z2)) n))
 392     (setq pr2 (mod (* (+ x1 z1) (- x2 z2)) n))
 393     (setq sq1 (mod (expt (+ pr1 pr2) 2) n))
 394     (setq sq2 (mod (expt (- pr1 pr2) 2) n))
 395     (setq x3 (mod (* (cadr q) sq1) n))
 396     (setq z3 (mod (* (car q) sq2) n))
 397     `(,x3 ,z3)))
 398
 399 (defun ecm-square (p n a)
 400   (let ((x1 (car p)) (z1 (cadr p))
 401     (x2) (z2) (sq1) (sq2) (f1) (f2))
 402     (setq sq1 (mod (* (+ x1 z1) (+ x1 z1)) n))
 403     (setq sq2 (mod (* (- x1 z1) (- x1 z1)) n))
 404     (setq f1 (- sq1 sq2))
 405     (setq f2 (mod (* a f1) n))
 406     (setq x2 (mod (* sq1 sq2) n))
 407     (setq z2 (mod (* f1 (+ sq2 f2)) n))
 408     `(,x2 ,z2)))
 409
 410 (defun ecm-power (base e n a)
 411   (let ((p base) (ptb (ecm-square base n a)) (l (integer-length e)))
 412     (do ((i (- l 2) (1- i)))
 413     ((< i 0))
 414       (if (logbitp i e)
 415       (progn
 416         (setq p (ecm-product base p ptb n))
 417         (setq ptb (ecm-square ptb n a)))
 418       (progn
 419         (setq ptb (ecm-product base p ptb n))
 420         (setq p (ecm-square p n a)))))
 421     p))
 422
 423 (defun ecm-factor-with-curve (n x z a lim1)
 424   (let ((g (gcd (- (* a a) 4) n)))
 425     (unless (= g 1) (return-from ecm-factor-with-curve g)))
 426   (setq a (mod (floor (+ a 2) 4) n))
 427   ;;
 428   ;; stage 1: compute p^M where M=p1^e1*...*pk^ek where
 429   ;;          p1,...,pk are primes < lim1 and ei=log[pi](n)
 430   ;;
 431   (let ((q 1)
 432     (last_q ($prev_prime lim1))
 433     (p `(,x ,z))
 434     (ex) (next_gcd)     (gcd_interval))
 435     (setq gcd_interval (floor lim1 4))
 436     (setq next_gcd gcd_interval)
 437     (do ()
 438     ((> q lim1))
 439       (setq q ($next_prime q))
 440       (setq ex (floor (log lim1) (log q)))
 441       (cond ((= q 2) (incf ex 2))
 442         ((= q 3) (incf ex)))
 443       (setq p (ecm-power p (expt q ex) n a))
 444       (when (>= q next_gcd)
 445     (let ((g (gcd (cadr p) n)))
 446       (when (< 1 g n)
 447         (when $ifactor_verbose
 448           (format t "ECM: found factor in stage 1: ~d (~d digits)~%" g (number-of-digits g)))
 449         (return-from ecm-factor-with-curve g))
 450       (setq next_gcd (min (+ next_gcd gcd_interval) last_q)))))
 451     ;;
 452     ;; stage 2: compute (p^M)^pi for each prime lim1<pi<lim2 (and some
 453     ;;          other exponents)
 454     ;;          Uses "Improved standard cotinuation".
 455     ;;
 456     (let* ((lim2 (* lim1 100))
 457        (power-after-1 p)
 458        (step-size (min (/ lim1 2) (isqrt (/ lim2 2))))
 459        (d-step-size (* 2 step-size))
 460        (power-table (make-array (+ 2 step-size)))
 461        (d-step-size-power (ecm-power power-after-1 d-step-size n a))
 462        (step-power (ecm-power power-after-1 (1+ d-step-size) n a))
 463        (last-step-power power-after-1)
 464        (step-pos 1)
 465        (q1 3)
 466        (prime-diffs-pos 0)
 467        (step-power-buff))
 468       (init-prime-diffs lim2)
 469       (setf (aref power-table 1) (ecm-square power-after-1 n a))
 470       (setf (aref power-table 2) (ecm-square (aref power-table 1) n a))
 471       (do ((i 3 (1+ i)))
 472       ((> i step-size))
 473     (setf (aref power-table i)
 474           (ecm-product (aref power-table (- i 2)) (aref power-table 1) (aref power-table (1- i)) n)))
 475       (do ()
 476       ((> step-pos (- lim2 d-step-size)))
 477     (let ((buff-prod 1)
 478           (q-limit (+ step-pos d-step-size))
 479           (power-table-pos (/ (- q1 step-pos) 2)))
 480       (when (zerop power-table-pos) ($error q1 step-pos))
 481       (do ()
 482           ((> q1 q-limit))
 483         (let* ((sp1 (car step-power))
 484            (sp2 (cadr step-power))
 485            (pp1 (car (aref power-table power-table-pos)))
 486            (pp2 (cadr (aref power-table power-table-pos)))
 487            (coord-diffs (mod (- (* sp1 pp2) (* sp2 pp1)) n)))
 488           (setq buff-prod (mod (* coord-diffs buff-prod) n)))
 489         (incf q1 (aref *prime-diffs* prime-diffs-pos))
 490         (incf power-table-pos (/ (aref *prime-diffs* prime-diffs-pos) 2))
 491         (incf prime-diffs-pos))
 492
 493       (let ((g (gcd n buff-prod)))
 494         (when (> g 1)
 495           (when $ifactor_verbose
 496         (format t "ECM: found factor in stage 2: ~d (~d digits)~%" g (number-of-digits g)))
 497           (return-from ecm-factor-with-curve g)))
 498
 499       (setq step-power-buff step-power)
 500       (setq step-power (ecm-product last-step-power     d-step-size-power step-power n))
 501       (setq last-step-power step-power-buff)
 502       (incf step-pos d-step-size))))
 503     nil))
 504
 505 (defun get-one-factor-ecm (n)
 506   (when (primep n) (return-from get-one-factor-ecm n))
 507   (let ((sigma (+ 6 (random (ash 1 20))))
 508     (x 0) (z 0) (u 0) (v 0) (a 0) (a1 0) (a2 0)
 509     (fact) (lim1 $ecm_limit) (a2_inv 0))
 510     (dotimes (i $ecm_number_of_curves)
 511       (setq u (mod (- (* sigma sigma) 5) n))
 512       (setq v (mod (* 4 sigma) n))
 513       (setq x (mod (expt u 3) n))
 514       (setq z (mod (expt v 3) n))
 515       (setq a1 (mod (* (expt (- v u) 3) (+ (* 3 u) v)) n))
 516       (setq a2 (mod (* 4 x v) n))
 517       (setq a2_inv (inv-mod a2 n))
 518       (when (null a2_inv)
 519     (return-from get-one-factor-ecm (gcd a2 n)))
 520       (setq a (mod (* a1 a2_inv) n))
 521       (setq sigma (max 6 (mod (+ (* sigma sigma) 1) n)))
 522       (when $ifactor_verbose
 523     (format t "ECM: trying with curve #~d of ~d (lim=~d)~%" (1+ i) $ecm_number_of_curves lim1))
 524       (setq fact (ecm-factor-with-curve n x z a lim1))
 525       (when (and fact (< fact n))
 526     (return-from get-one-factor-ecm fact))
 527       (setq lim1 (min (+ lim1 $ecm_limit_delta) $ecm_max_limit)))
 528     nil))
 529
 530
 531 ;;; convert (3 5 3 5 3 7) to ((3 3) (5 2) (7 1))
 532
 533 (defun convert-list (l)
 534   (labels ((convert-list-sub (e n l1 l2)
 535          (cond ((null l1)
 536             (cons (list e n) l2))
 537            ((= e (car l1))
 538             (convert-list-sub e (1+ n) (cdr l1) l2))
 539            (t (convert-list-sub (car l1) 1 (cdr l1) (cons `(,e ,n) l2))))))
 540     (let ((l1 (sort l #'>)))
 541       (convert-list-sub (car l1) 1 (rest l1) nil))))
 542
 543 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 544 ;;;                                                                       ;;;
 545 ;;;           ~*~  IMPLEMENTATION OF PRIMALITY TESTS  ~*~                 ;;;
 546 ;;;                                                                       ;;;
 547 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 548
 549 (defmfun $primep (n)
 550   (if (integerp n)
 551       (primep (abs n))
 552       (merror (intl:gettext "primep: argument must be an integer; found: ~M") n)))
 553
 554 (defun primep (n)
 555   (cond
 556     ((= n 1) nil)
 557     ((evenp n) (= n 2))
 558     ((<= n *largest-small-prime*) (when (member n *small-primes*) t))
 559     ((< n 9080191) (primep-small n '(31 73)))
 560     ((< n 4759123141) (primep-small n '(2 7 61)))
 561     ((< n 2152302898747) (primep-small n '(2 3 5 7 11)))
 562     ((< n 3474749660383) (primep-small n '(2 3 5 7 11 13)))
 563     ((< n 341550071728321) (primep-small n '(2 3 5 7 11 13 17)))
 564     ((< n 3825123056546413051) (primep-small n '(2 3 5 7 11 13 17 19 23)))
 565     ((< n 318665857834031151167461) (primep-small n '(2 3 5 7 11 13 17 19 23 29 31 37)))
 566     ((< n 3317044064679887385961981) (primep-small n '(2 3 5 7 11 13 17 19 23 29 31 37 41)))
 567     ((member n *large-primes*) t)
 568     (t (primep-prob n)) ))
 569
 570 ;;; A Miller-Rabin test is deterministic for small n if we test for small bases.
 571 ;;; Reference:
 572 ;;;  [1] G. Jaeschke, On Strong Pseudoprimes to Several Bases,
 573 ;;;         Math. Comp., 61 (1993), 915-926.
 574 ;;;  [2] http://primes.utm.edu/prove/prove2_3.html
 575 ;;;  [3] Jiang, Deng - Strong pseudoprimes to the first eight prime bases (2014)
 576 ;;;         Mathematics of Computation, Vol 83, Nr 290, Pages 2915--2924
 577 ;;;  [3] Sorenson, Webster - Strong Pseudoprimes to Twelve Prime Bases (2015)
 578 ;;;         arXiv:1509.00864v1 [math.NT]
 579
 580 (defun primep-small (n bases)
 581   (multiple-value-bind (q k) (miller-rabin-decomposition n)
 582     (dolist (x bases t)
 583       (unless (miller-rabin-kernel n q k x)
 584         (return-from primep-small nil) ))))
 585
 586 ;;; strong primality test:
 587 ;;;  - run $primep_number_of_tests times a Miller-Rabin test
 588 ;;;  - run one Lucas test
 589
 590 (defun primep-prob (n)
 591   ;; Miller-Rabin tests:
 592   (multiple-value-bind (q k) (miller-rabin-decomposition n)
 593     (dotimes (i $primep_number_of_tests)
 594       (unless (miller-rabin-kernel n q k)
 595         (return-from primep-prob nil) )))
 596   ;; Lucas test:
 597   (primep-lucas n) )
 598
 599
 600 ;;; Miller-Rabin (algorithm P from D. Knuth, TAOCP, 4.5.4)
 601 ;;;
 602 ;;;   - write n-1 = q*2^k (n,q odd, n > 2)
 603 ;;;   - x is a random number 1 < x < n
 604 ;;;   - n passes the test if
 605 ;;;         x^q = 1 (mod n)
 606 ;;;         or x^(q*2^j) = -1 (mod n) for some j = 0..k-1
 607 ;;;
 608 ;;; A prime number must pass this test.
 609 ;;; The probability of passing one test and not being a prime is less than 1/4.
 610
 611 ;; return values q,k with n-1 = q*2^k
 612 (defun miller-rabin-decomposition (n) ;; assume n > 2 (n-1 is even)
 613   (do ((k 1 (1+ k))
 614        (q (ash n -1) (ash q -1)) )
 615       ((logbitp 0 q) (values q k)) ))
 616 ;;
 617 ;; now assume n-1 = q*2^k, k >= 1
 618 (defun miller-rabin-kernel (n q k &optional x)
 619   (unless x
 620     (setq x (+ (random (- n 2)) 2)) )
 621   (let ((y (power-mod x q n)) ;; j = 0
 622         (minus1 (1- n)) )
 623     (if (or (= y 1) (= y minus1))
 624       t
 625       (do ((j 1 (1+ j)))
 626           ((= j k))
 627         (setq y (power-mod y 2 n))
 628         (when (= y minus1) (return t))
 629         (when (= y 1) (return)) )))) ;; n prime => last y must have been 1 or -1
 630
 631
 632 (defmfun $power_mod (b e m)
 633   (unless (and (integerp b) (integerp e) (integerp m))
 634     (merror (intl:gettext "power_mod: arguments must be integers; found: ~M, ~M, ~M") b e m) )
 635   (if (>= e 0)
 636     (power-mod b e m)
 637     (let ((inv (inv-mod b m)))
 638       (when inv
 639         (power-mod inv (- e) m) ))))
 640 ;;
 641 (defun power-mod (b e m)
 642   (declare (optimize (speed 3) (safety 0)))
 643   (cond
 644     ((zerop e)
 645       (mod 1 m) )
 646     ((typep e 'fixnum)
 647       (do ((res 1)) (())
 648         (when (logbitp 0 e)
 649           (setq res (mod (* res b) m))
 650           (when (= 1 e) (return res)) )
 651         (setq e (ash e -1)
 652               b (mod (* b b) m)) ))
 653     (t ;; sliding window variant:
 654       (let* ((l (integer-length e))
 655              (k (cond ((< l  65) 3)
 656                       ((< l 161) 4)
 657                       ((< l 385) 5)
 658                       ((< l 897) 6)
 659                       (t         7) ))
 660              (tab (power-mod-tab b k m))
 661              (res 1) s u tmp )
 662         (do ((i (1- l)))
 663             ((< i 0) res)
 664           (cond
 665             ((logbitp i e)
 666               (setq s (max (1+ (- i k)) 0))
 667               (do () ((logbitp s e)) (incf s))
 668               (setq tmp (1+ (- i s)))
 669               (dotimes (h tmp) (setq res (mod (* res res) m)))
 670               (setq u (ldb (byte tmp s) e))
 671               (unless (= u 0) (setq res (mod (* res (svref tab (ash u -1))) m)))
 672               (setq i (1- s)) )
 673             (t
 674               (setq res (mod (* res res) m))
 675               (decf i) )))))))
 676 ;;
 677 (defun power-mod-tab (b k m)
 678   (declare (optimize (speed 3) (safety 0)))
 679   (let* ((l (ash 1 (1- k)))
 680          (tab (make-array l :element-type 'integer :initial-element 1))
 681          (bi b)
 682          (bb (mod (* b b) m)) )
 683     (setf (svref tab 0) b)
 684     (do ((i 1 (1+ i)))
 685         ((= i l) tab)
 686       (setq bi (mod (* bi bb) m))
 687       (setf (svref tab i) bi) )))
 688
 689 ;;; primep-lucas:
 690 ;;;
 691 ;;;  Define: x^2-a*x+b, D=a^2-4*b; x1, x2 roots of x^2+a*x+b;
 692 ;;;  U[k]=(x1^k-x2^k)/(x1-x2), V[k]=x1^k+x2^k.
 693 ;;;
 694 ;;;  Lucas theorem: If p is an odd prime, gcd(p,b)=1 and jacobi(D,p)=-1,
 695 ;;;                 then p divides U[p+1].
 696 ;;;
 697 ;;;  We calculate U[p+1] for x^2-b*x+1 where jacobi(b^2-4,n)=-1
 698 ;;;  and test if p divides U[p+1].
 699
 700 (defun primep-lucas (n)
 701   (let (prmp (b 3))
 702     (loop while (not (= ($jacobi (- (* b b) 4) n) -1)) do
 703      (incf b))
 704     (setq prmp (zerop (lucas-sequence (1+ n) b n)))
 705     (when (and prmp $save_primes)
 706       (push n *large-primes*))
 707     prmp))
 708
 709 ;;; Get element U[p+1] of Lucas sequence for x^2-p*x+1.
 710 ;;;
 711 ;;; Uses algorithm from M. Joye and J.-J. Quisquater,
 712 ;;;                     Efficient computation of full Lucas sequences, 1996
 713
 714 (defun lucas-sequence (k p n)
 715   (let ((uh 1) (vl 2) (vh p) (s 0) l)
 716     (do ()
 717     ((logbitp 0 k))
 718       (setq k (ash k -1))
 719       (setq s (1+ s)))
 720
 721     (setq l (integer-length k))
 722
 723     (do ((j (1- l) (1- j)))
 724     ((= 0 j))
 725       (if (logbitp j k)
 726       (progn
 727         (setq uh (mod (* uh vh) n))
 728         (setq vl (mod (- (* vh vl) p) n))
 729         (setq vh (mod (- (* vh vh) 2) n)))
 730       (progn
 731         (setq uh (mod (1- (* uh vl)) n))
 732         (setq vh (mod (- (* vh vl) p) n))
 733         (setq vl (mod (- (* vl vl) 2) n)))))
 734
 735     (setq uh (mod (1- (* uh vl)) n))
 736     (setq vl (mod (- (* vh vl) p) n))
 737
 738     (dotimes (j s)
 739       (setq uh (mod (* uh vl) n))
 740       (setq vl (mod (- (* vl vl) 2) n)))
 741     uh))
 742
 743 ;;; first values of next_prime
 744 (defvar *next_prime_ar* #(0 2 3 5 5 7 7))
 745
 746 ;;; first values of prev_prime
 747 (defvar *prev_prime_ar* #(0 0 0 2 3 3 5 5 7 7 7 7))
 748
 749 ;;; gaps between numbers that are not multiples of 2,3,5,7
 750 (defvar deltaprimes_next
 751   '(1 10 9 8 7 6 5 4 3 2 1 2 1 4 3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 2
 752     1 6 5 4 3 2 1 4 3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 6 5 4 3 2 1 2 1
 753     6 5 4 3 2 1 4 3 2 1 2 1 6 5 4 3 2 1 4 3 2 1 6 5 4 3 2 1 8 7 6
 754     5 4 3 2 1 4 3 2 1 2 1 4 3 2 1 2 1 4 3 2 1 8 7 6 5 4 3 2 1 6 5
 755     4 3 2 1 4 3 2 1 6 5 4 3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 2 1 6 5 4
 756     3 2 1 6 5 4 3 2 1 4 3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 2 1 6 5 4 3
 757     2 1 4 3 2 1 2 1 4 3 2 1 2 1 10 9 8 7 6 5 4 3 2 1 2))
 758
 759 (defvar deltaprimes_prev
 760   '(-1 -2 -1 -2 -3 -4 -5 -6 -7 -8 -9 -10 -1 -2 -1 -2 -3 -4 -1 -2
 761     -1 -2 -3 -4 -1 -2 -3 -4 -5 -6 -1 -2 -1 -2 -3 -4 -5 -6 -1 -2 -3
 762     -4 -1 -2 -1 -2 -3 -4 -1 -2 -3 -4 -5 -6 -1 -2 -3 -4 -5 -6 -1 -2
 763     -1 -2 -3 -4 -5 -6 -1 -2 -3 -4 -1 -2 -1 -2 -3 -4 -5 -6 -1 -2 -3
 764     -4 -1 -2 -3 -4 -5 -6 -1 -2 -3 -4 -5 -6 -7 -8 -1 -2 -3 -4 -1 -2
 765     -1 -2 -3 -4 -1 -2 -1 -2 -3 -4 -1 -2 -3 -4 -5 -6 -7 -8 -1 -2 -3
 766     -4 -5 -6 -1 -2 -3 -4 -1 -2 -3 -4 -5 -6 -1 -2 -1 -2 -3 -4 -1 -2
 767     -3 -4 -5 -6 -1 -2 -1 -2 -3 -4 -5 -6 -1 -2 -3 -4 -5 -6 -1 -2 -3
 768     -4 -1 -2 -1 -2 -3 -4 -1 -2 -3 -4 -5 -6 -1 -2 -1 -2 -3 -4 -5 -6
 769     -1 -2 -3 -4 -1 -2 -1 -2 -3 -4 -1 -2 -1 -2 -3 -4 -5 -6 -7 -8 -9
 770     -10))
 771
 772 ;;; product of primes in [59..2897]
 773 (defvar bigprimemultiple 6805598092615180737440235028147472981586738014295015027644884201753964648883910180850814465749532893719128055374719237806417537893593625321589379773764981786235326314555704406245399180879758341371676681401881451390195684863765326592983982964414393796690715805513465774520452671995927595391575142047776807977863591126244782181086547150369260177339043045082132788709080989495477932949788444703905327686499493503904132269141007955089790798876488207574072278769735865653223865994494346936718462923487228576140267887355548289736131557613540186975875834980017431190021254898173201223012171417763388931502928376549397638685218312217808199405294916194758171476025904777185780125034583816795375331627264462778001498062163759312245245590800878057927864359433868165604228946307536835897173733369926842890411102870160854438921809703357774373318146115616129588245083207631664167515206143659538759733110973189757163548882116479710800109577584318611988710048552969742803870964125788279451564113232340649434743105271873797620278073136369295820926294656549976175331880139356684249842712956493849288710258349886914201056170180503844749859595207139766052196982574437241716274871254310342540993006427120762049161745282399431514257565489)
 774
 775 (defmfun $next_prime (n)
 776   (unless (and (integerp n))
 777     (merror (intl:gettext "next_prime: argument must be an integer; found: ~M") n))
 778   (cond ((< n 2) 2)
 779     ((<= n 6) (aref *next_prime_ar* n))
 780     ((< n 100000) (return-from $next_prime (next-prime-det n deltaprimes_next)))
 781     (t (next-prime-prob n deltaprimes_next))))
 782
 783 (defmfun $prev_prime (n)
 784   (unless (and (integerp n) (> n 2))
 785     (merror (intl:gettext "prev_prime: argument must be an integer greater than 2; found: ~M") n))
 786   (if (<= n 11) (return-from $prev_prime (aref *prev_prime_ar* n)))
 787   (if (< n 100000) (return-from $prev_prime (next-prime-det n deltaprimes_prev)))
 788   (next-prime-prob n deltaprimes_prev))
 789
 790
 791 ;;; Find next/prev prime using deterministic test that checks all
 792 ;;; divisors < sqrt(n) and skipping all multiples of 2,3,5,7
 793 ;;; preconditions: 11 < n < 9973*9973
 794 (defun next-prime-det (n deltaprimes)
 795   (incf n (nth (mmod n 210) deltaprimes))
 796   (loop while 1 do
 797        (dolist (p *small-primes*)
 798      (if (= (mmod n p) 0) (return))
 799      (if (>= (* p p) n) (return-from next-prime-det n)))
 800        (incf n (nth (mmod n 210) deltaprimes))))
 801
 802 ;;; Find next/prev prime using probabilistic test and skipping al multiples of
 803 ;;; 2,3,5,7 using deltaprimes list and calculating gcd's with product of
 804 ;;; prime numbers
 805 (defun next-prime-prob (n deltaprimes)
 806   ;; skip all multiples of 2,3,5,7
 807   (incf n (nth (mmod n 210) deltaprimes))
 808   (loop
 809      (and
 810       ;; gcd agaist product of primes in [11..31]
 811       (= (gcd n 955049953) 1)
 812       ;; gcd agaist product of primes in [37..53]
 813       (= (gcd n 162490421) 1)
 814       ;; gcd agaist product of primes in [59..2897]
 815       (= (gcd n bigprimemultiple) 1)
 816       (primep n)
 817       (return-from next-prime-prob n))
 818      ;; skip all multiples of 2,3,5,7"
 819      (incf n (nth (mmod n 210) deltaprimes))))
 820
 821
 822 (defun next-prime (n c)
 823   (when (evenp n) (incf n c))
 824   (loop
 825      (when (primep n)
 826        (return-from next-prime n))
 827      (incf n (* 2 c))))
 828
 829 ;;; return a list of all primes between start and end
 830
 831 (defmfun $primes (start end)
 832   (unless (and (integerp start) (integerp end))
 833     (merror (intl:gettext "primes: arguments must be integers; found: ~M, ~M") start end))
 834   (let ((primes nil))
 835     (cond
 836       ;; take primes from *small-primes* if possible
 837       ((<= start *largest-small-prime*)
 838         (dolist (n *small-primes*)
 839           (when (<= start n end)
 840             (push n primes) ))
 841         (setq start *largest-small-prime*) )
 842       (t
 843         (decf start) )) ; $next_prime returns a value >= argument + 1
 844     ;; search for the rest of primes
 845     (do ((n ($next_prime start) ($next_prime (1+ n))))
 846         ((> n end) (cons '(mlist) (reverse primes)))
 847       (push n primes) )))