factor

(define start
    ;;; This program is a tribute to the test program 'factor.nj', which severely penetrated my implementation of the garbage collector for the Ninja VM.
    (lambda ()
        ;(small-primes 100)
        (display "upper limit: ")
        (let ((lim (read)))
            (compute 1 lim)
        )
    )
)
(define small-primes-result '())
(define (small-primes limit . rest) ; get a list of small primes up to limit
    (if (< limit 0)
        (begin
            (set! small-primes-result '()) ;flush primes with a negative arg
            ()
        )
        (if (null? small-primes-result)
            (let (( i (if (null? rest) 0 (car rest)) ))
                (if (>= i limit)
                    ()
                    (if (= i 0)
                        (begin
                            (set! small-primes-result (cons 2 (cons 3 (small-primes limit 5))))
                            small-primes-result
                        )
                        (if (is-prime i)
                            (cons i (small-primes limit (+ i 2)))
                            (small-primes limit (+ i 2))
                        )
                    )
                )
            )
            small-primes-result
        )
    )
)
(define (is-prime n . rest)
    (let (( i (if (null? rest) 0 (car rest)) ))
        (if (< i 3)
            (is-prime n 3) ;set start value -> i=3
            (if (<= (square i) n)
                (if (= (modulo n i) 0)
                    #f
                    (is-prime n (+ i 2))
                )
                #t
            )
        )
    )
)
(define (compute n limit)
    (if (= n limit)
        #t
        (begin
            (display "10^")
            (display n)
            (display "+1 = ")
            (let ((m (+ (expt 10 n) 1)))
                (fmt m)
                (let find-fac ((factors (factorize m)))
                    (fmt "result: " 0)
                    (fmt factors)
                    ; check if all primes -> if yes, check product. if no, factorize remaining composites
                    (fmt "check: product = " 0)
                    (let ((prod (reduce factors (lambda (x y) (* x y)) 1)))
                        (fmt prod)
                        (if (not (= m prod))
                            #f
                            (compute (+ n 1) limit)
                        )
                    )
                )
            )
        )
    )
)
(define (all-primes l)
    (if (null? l)
        #t
        (if (is-prime (car l))
            (all-primes (cdr l))
            #f
        )
    )
)
(define (factorize n) ;factorize n using trial division
    (if (< n 2)
        ()
        (let ((p (small-prime-factor n)))
            (if (> p -1)
                (if (< (square p) n)
                    (begin
                        (fmt "detected small prime factor " 0)
                        (fmt p)
                        (cons p (factorize (/ n p)))
                    )
                    (begin
                        (fmt "the remaining factor " 0)
                        (fmt n 0)
                        (fmt " is a small prime")
                        (cons n ())
                    )
                )
                (begin
                    (fmt "the remaining factor " 0)
                    (fmt n 0)
                    (fmt " doesn't have any prime factors < " 0)
                    (fmt SMALL_PRIMES_LIM)
                    (if #t  ;;; disable composite testing for now
                        (cons n ())
                        (if (is-composite n)
                            (begin
                                (fmt "but is definitely composite")
                                (let ((f1 (find-factor n)))
                                    (if (= f1 -1)
                                        (begin
                                            (fmt "cannot factorize " 0)
                                            (fmt n 0)
                                            (fmt ", giving up")
                                            (cons n ())
                                        )
                                        (let ((f2 (/ n f1)))
                                            (fmt "this number can be split into " 0)
                                            (fmt f1 0)
                                            (fmt " and " 0)
                                            (fmt f2)
                                            (cons (factorize f1) (cons (factorize f2)))
                                        )
                                    )
                                )
                            )
                            (begin
                                (fmt "and is very probably prime")
                                ;if prove-prime(x) then
                                ;   primality has been proven
                                (cons n ())
                            )
                        )
                    )
                )
            )
        )
    )
)
(define SMALL_PRIMES_LIM 2000)
(define (small-prime-factor n)
    ; find a small prime factor of n (using prime list plist)
    ; -1: not found
    (let cont ((plist (small-primes SMALL_PRIMES_LIM)))
        (if (< (length plist) 1)
            -1
            (if (= (modulo n (car plist)) 0)
                (car plist)
                (cont (cdr plist))
            )
        )
    )
)
(define (is-composite n)     ;Rabin-Miller test
    (let cont1 ((q (- n 1)) (t 0))
        (if (= (modulo q 2) 0)
            (cont1 (/ q 2) (+ t 1))
            'todo ;...
        )
    )
)
;procedure pollard(n):
;   x<-2, y<-2, d<-1
;   while d=1 do
;       x<-g(x)=(x^2-1) mod n
;       y<-g(g(y))=((y^2-1) mod n)^2-1 mod n)
;       d<-gcd(|x-y|,n)
;   if d=n then
;       return FAIL
;   else
;       return n

(define (find-factor n)      ;Pollard rho method
 'todo
)

(start)

;module tests
(define (test-is-prime)
    (fmt (is-prime 71 0))
)
(define (test-small-primes)
    (fmt (small-primes 100))
)
(define (test-factorize)
    (assert-eq (factorize 144) '(2 2 2 2 3 3))
)
(define (test-compute)
    (compute 1 40)
)