[ create a new paste ] login | about

Project: programmingpraxis
Link: http://programmingpraxis.codepad.org/S6l8kId7    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Aug 6:
; lenstra's algorithm

(define sort #f)
(define merge #f)
(let ()
  (define dosort
    (lambda (pred? ls n)
      (if (= n 1)
          (list (car ls))
          (let ((i (quotient n 2)))
            (domerge pred?
                     (dosort pred? ls i)
                     (dosort pred? (list-tail ls i) (- n i)))))))
  (define domerge
    (lambda (pred? l1 l2)
      (cond
        ((null? l1) l2)
        ((null? l2) l1)
        ((pred? (car l2) (car l1))
         (cons (car l2) (domerge pred? l1 (cdr l2))))
        (else (cons (car l1) (domerge pred? (cdr l1) l2))))))
  (set! sort
    (lambda (pred? l)
      (if (null? l) l (dosort pred? l (length l)))))
  (set! merge
    (lambda (pred? l1 l2)
      (domerge pred? l1 l2))))

(define rand
  (let* ((a 3141592653) (c 2718281829)
         (m (expt 2 35)) (x 5772156649)
         (next (lambda ()
                 (let ((x-prime (modulo (+ (* a x) c) m)))
                   (set! x x-prime) x-prime)))
         (k 103)
         (v (list->vector (reverse
              (let loop ((i k) (vs (list x)))
                (if (= i 1) vs
                  (loop (- i 1) (cons (next) vs)))))))
         (y (next))
         (init (lambda (s)
                 (set! x s) (vector-set! v 0 x)
                 (do ((i 1 (+ i 1))) ((= i k))
                   (vector-set! v i (next))))))
    (lambda seed
      (cond ((null? seed)
              (let* ((j (quotient (* k y) m))
                     (q (vector-ref v j)))
                (set! y q)
                (vector-set! v j (next)) (/ y m)))
            ((eq? (car seed) 'get) (list a c m x y k v))
            ((eq? (car seed) 'set)
              (let ((state (cadr seed)))
                (set! a (list-ref state 0))
                (set! c (list-ref state 1))
                (set! m (list-ref state 2))
                (set! x (list-ref state 3))
                (set! y (list-ref state 4))
                (set! k (list-ref state 5))
                (set! v (list-ref state 6))))
            (else (init (modulo (numerator
                    (inexact->exact (car seed))) m))
                  (rand))))))

(define (randint . args)
  (cond ((null? (cdr args))
          (inexact->exact (floor (* (rand) (car args)))))
        ((< (car args) (cadr args))
          (+ (inexact->exact (floor (* (rand) (- (cadr args) (car args))))) (car args)))
        (else (+ (inexact->exact (ceiling (* (rand) (- (cadr args) (car args))))) (car args)))))

(define (primes n)
  (let* ((max-index (quotient (- n 3) 2))
         (v (make-vector (+ 1 max-index) #t)))
    (let loop ((i 0) (ps '(2)))
      (let ((p (+ i i 3)) (startj (+ (* 2 i i) (* 6 i) 3)))
        (cond ((>= (* p p) n)
               (let loop ((j i) (ps ps))
                  (cond ((> j max-index) (reverse ps))
                        ((vector-ref v j)
                          (loop (+ j 1) (cons (+ j j 3) ps)))
                        (else (loop (+ j 1) ps)))))
              ((vector-ref v i)
                (let loop ((j startj))
                  (if (<= j max-index)
                      (begin (vector-set! v j #f)
                             (loop (+ j p)))))
                      (loop (+ 1 i) (cons p ps)))
              (else (loop (+ 1 i) ps)))))))

(define (expm b e m)
  (define (m* x y) (modulo (* x y) m))
  (cond ((zero? e) 1)
        ((even? e) (expm (m* b b) (/ e 2) m))
        (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))

(define (check? a n)
  (let loop ((r 0) (s (- n 1)))
    (if (even? s) (loop (+ r 1) (/ s 2))
      (if (= (expm a s n) 1) #t
        (let loop ((j 0) (s s))
          (cond ((= j r) #f)
                ((= (expm a s n) (- n 1)) #t)
                (else (loop (+ j 1) (* s 2)))))))))

(define (prime? n)
  (cond ((< n 2) #f) ((= n 2) #t) ((even? n) #f)
        (else (let loop ((k 50))
                (cond ((zero? k) #t)
                      ((not (check? (randint 1 n) n)) #f)
                      (else (loop (- k 1))))))))

(define (lenstra n limit)
  (call-with-current-continuation (lambda (return)
    (let* ((x (randint n)) (y (randint n)) (a (randint n))
           (b (- (* y y) (* x x x) (* a x)))
           (d (+ (* 4 a a a) (* 27 b b))) (g (gcd d n))
           (ps (primes limit)))
      (define (inverse x)
        (let loop ((x (modulo x n)) (a 1))
          (cond ((zero? x) (return #f))
                ((= x 1) a)
                (else (let ((q (- (quotient n x))))
                        (loop (+ n (* q x)) (modulo (* q a) n)))))))
      (define (add p1 p2)
        (let* ((x1 (car p1)) (y1 (cdr p1))
               (x2 (car p2)) (y2 (cdr p2))
               (g (gcd (- x1 x2) n)))
          (if (zero? g) (return #f)) (if (< 1 g n) (return g))
          (let* ((slope (* (inverse (- x1 x2)) (- y1 y2)))
                 (newx (modulo (- (* slope slope) x1 x2) n))
                 (newy (modulo (- (* slope (- x1 newx)) y1) n)))
            (cons newx newy))))
      (define (double p)
        (let* ((x (car p)) (y (cdr p)) (g (gcd (+ y y) n)))
          (if (zero? g) (return #f)) (if (< 1 g n) (return g))
          (let* ((slope (* (inverse (+ y y)) (+ (* 3 x x) b)))
                 (newx (modulo (- (* slope slope) x x) n))
                 (newy (modulo (- (* slope (- x newx)) y) n)))
            (cons newx newy))))
      (define (mult k p)
        (cond ((= k 1) p)
              ((even? k) (mult (quotient k 2) (double p)))
              (else (add (mult (quotient k 2) (double p)) p))))
      (if (zero? g) (return #f)) (if (< 1 g n) (return g))
      (let loop ((p (car ps)) (q (car ps)) (ps (cdr ps)) (xy (cons x y)))
        (cond ((null? ps) (return #f))
              ((< limit q) (loop (car ps) (car ps) (cdr ps) xy))
              (else (loop p (* p q) ps (mult p xy)))))))))

(define (factor n limit curves)
  (let loop ((curves curves))
    (if (zero? curves) #f
      (let ((f (lenstra n limit)))
        (if f f (loop (- curves 1)))))))

(define (trial-factors n ps)
  (let loop ((ps ps) (facts '()) (cofact n))
    (cond ((or (null? ps) (= cofact 1))
            (values (reverse facts) cofact))
          ((< cofact (* (car ps) (car ps)))
            (values (reverse (cons cofact facts)) 1))
          ((zero? (modulo cofact (car ps)))
            (loop ps (cons (car ps) facts)
                  (quotient cofact (car ps))))
          (else (loop (cdr ps) facts cofact)))))

(define (factors n . args) ; factors n limit curves trial-limit
  (let ((limit (if (< 1 (length args)) (car args) 10000))
        (curves (if (< 2 (length args)) (cadr args) 1000))
        (trial-limit (if (< 3 (length args)) (caddr args) #e1e6)))
    (call-with-values
      (lambda () (trial-factors n (primes trial-limit)))
      (lambda (facts cofact)
        (let fact ((n cofact) (facts facts))
          (cond ((= n 1) (sort < facts))
                ((prime? n) (sort < (cons n facts)))
                (else (let ((f (factor n limit curves)))
                        (cond ((not f)
                                (display "failed with factors")
                                (for-each (lambda (f)
                                            (display " ") (display f))
                                          (sort < facts))
                                (display " and co-factor ")
                                (display n) (newline))
                             ((prime? f) (fact (/ n f) (cons f facts)))
                             (else (let ((fs (fact (/ n f) (list f))))
                                     (fact (apply / n fs) (append fs facts)))))))))))))

(display (factors (/ (- (expt 10 47) 1) 9)))


Output:
1
(35121409 316362908763458525001406154038726382279)


Create a new paste based on this one


Comments: