; fermat's divisors challenges
(define (iroot k n)
(let ((k-1 (- k 1)))
(let loop ((u n) (s (+ n 1)))
(if (<= s u) s
(loop (quotient (+ (* k-1 u) (quotient n (expt u k-1))) k) u)))))
(define (primes n) ; list of primes not exceeding n
(let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
(let loop ((i 0) (p 3) (ps (list 2)))
(cond ((< n (* p p))
(do ((i i (+ i 1)) (p p (+ p 2))
(ps ps (if (vector-ref bits i) (cons p ps) ps)))
((= i len) (reverse ps))))
((vector-ref bits i)
(do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
(vector-set! bits j #f)))
(else (loop (+ i 1) (+ p 2) ps))))))
(define prime? ; strong pseudoprime to prime bases less than 100
(let* ((ps (primes 100)) (p100 (apply * ps)))
(lambda (n)
(define (expm b e m)
(let loop ((b b) (e e) (x 1))
(if (zero? e) x
(loop (modulo (* b b) m) (quotient e 2)
(if (odd? e) (modulo (* b x) m) x)))))
(define (spsp? n a) ; #t if n is a strong pseudoprime base a
(do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
((odd? d) (if (= (expm a d n) 1) #t
(do ((r 0 (+ r 1)))
((or (= (expm a (* d (expt 2 r)) n) (- n 1)) (= r s))
(< r s)))))))
(if (< n 2) #f (if (< 1 (gcd n p100)) (if (member n ps) #t #f)
(do ((ps ps (cdr ps)))
((or (null? ps) (not (spsp? n (car ps)))) (null? ps))))))))
(define (next-prime n) ; smallest prime larger than n
(define (wheel n)
(vector-ref (vector 1 6 5 4 3 2 1 4 3 2 1 2 1 4
3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 2) (modulo n 30)))
(if (< n 2) 2 (if (< n 3) 3 (if (< n 5) 5
(let loop ((p (+ n (wheel n))))
(if (prime? p) p (loop (+ p (wheel p)))))))))
(define (prev-prime n) ; largest prime smaller than n
(define (wheel n)
(vector-ref (vector 1 2 1 2 3 4 5 6 1 2 3 4 1 2
1 2 3 4 1 2 1 2 3 4 1 2 3 4 5 6) (modulo n 30)))
(if (<= n 2) #f (if (<= n 3) 2 (if (<= n 5) 3 (if (<= n 7) 5
(let loop ((p (- n (wheel n))))
(if (prime? p) p (loop (- p (wheel p))))))))))
(define (factors n) ; factors of n in increasing order
(define (last-pair xs) (if (null? (cdr xs)) xs (last-pair (cdr xs))))
(define (cycle . xs) (set-cdr! (last-pair xs) xs) xs)
(define wheel235 (cons 1 (cons 2 (cons 2 (cycle 4 2 4 2 4 6 2 6)))))
(define (cons< x xs)
(if (null? xs) (list x)
(if (< x (car xs)) (cons x xs)
(cons (car xs) (cons< x (cdr xs))))))
(let wheel ((n n) (f 2) (fs (list)) (w wheel235))
(cond ((< n (* f f)) (reverse (cons n fs)))
((zero? (modulo n f)) (wheel (quotient n f) f (cons f fs) w))
((< f 1000) (wheel n (+ f (car w)) fs (cdr w)))
(else (let rho ((n n) (t 1) (h 1) (d 1) (c 1) (fs (reverse fs)))
(cond ((prime? n) (cons< n fs))
((= d 1)
(let* ((h (modulo (+ (* h h) c) n))
(h (modulo (+ (* h h) c) n))
(t (modulo (+ (* t t) c) n)))
(rho n t h (gcd (- t h) n) c fs)))
((= d n) (rho n 1 1 1 (+ c 1) fs))
((prime? d)
(rho (/ n d) 1 1 1 (+ c 1) (cons< d fs)))
(else (rho n 1 1 1 (+ c 1) fs))))))))
(define (divisors n) ; divisors of n, including 1 and n
(let ((divs (list 1)))
(do ((fs (factors n) (cdr fs))) ((null? fs) (sort < divs))
(let ((temp (list)))
(do ((ds divs (cdr ds)))
((null? ds) (set! divs (append divs temp)))
(let ((d (* (car fs) (car ds))))
(when (not (member d divs))
(set! temp (cons d temp)))))))))
(define (sigma x n . args) ; sum of x'th powers of divisors of n
(define (uniq-c eql? xs)
(if (null? xs) xs
(let loop ((xs (cdr xs)) (prev (car xs)) (k 1) (result '()))
(cond ((null? xs) (reverse (cons (cons prev k) result)))
((eql? (car xs) prev) (loop (cdr xs) prev (+ k 1) result))
(else (loop (cdr xs) (car xs) 1 (cons (cons prev k) result)))))))
(define (prod xs) (apply * xs))
(if (= n 1) 1
(let ((fs (uniq-c = (if (pair? args) (car args) (factors n)))))
(if (zero? x)
(prod (map add1 (map cdr fs)))
(prod (map (lambda (p a)
(/ (- (expt p (* (+ a 1) x)) 1) (- (expt p x) 1)))
(map car fs) (map cdr fs)))))))
(define (totient n) ; count of positive integers less than n coprime to it
(let loop ((t n) (p 0) (fs (factors n)))
(if (null? fs) t
(let ((f (car fs)))
(loop (if (= f p) t (* t (/ (- f 1) f))) f (cdr fs))))))
(define (moebius n) ; (-1)^k if n has k factors, or 0 if any factors duplicated
(let loop ((m -1) (f 0) (fs (factors n)))
(if (null? fs) m
(if (= f (car fs)) 0
(loop (- m) (car fs) (cdr fs))))))
(define (fermat n) (+ (expt 2 (expt 2 n)) 1)) ; 2^2^n+1
(define (mersenne n) (- (expt 2 n) 1)) ; 2^n-1
(define (repunit n . bs) ; (b^n-1)/(b-1)
(let ((b (if (pair? bs) (car bs) 10)))
(/ (- (expt b n) 1) (- b 1))))
(define (factorial n) ; product(1..n)
(let loop ((n n) (f 1))
(if (zero? n) f
(loop (- n 1) (* f n)))))
(define (primorial n) ; product(2..p_k..p_n) with p prime
(let loop ((k 0) (p 2) (prim 1))
(if (= k n) prim
(loop (+ k 1) (next-prime p) (* prim p)))))
(define (fibonacci n) ; F(n) = F(n-1) + Fb(n-2) with F(1) = F(2) = 1
(letrec ((square (lambda (x) (* x x))))
(cond ((= n 0) 0) ((= n 1) 1) ((= n 2) 1)
((odd? n) (let ((n2 (+ (quotient n 2) 1)))
(+ (square (fibonacci n2))
(square (fibonacci (- n2 1))))))
(else (let ((n2 (quotient n 2)))
(* (fibonacci n2) (+ (* (fibonacci (- n2 1)) 2)
(fibonacci n2))))))))
(define (lucas n) ; L(n) = L(n-1) + L(n-2) with L(1) = 1 and L(2) = 3
(/ (fibonacci (+ n n)) (fibonacci n)))
(define partition ; number of ways of writing n as sum of integers > 0
(let ((len 1) (pv (vector 1)))
(lambda (n)
(do () ((< n len))
(let* ((new-len (+ len len))
(new-pv (make-vector new-len #f)))
(do ((i 0 (+ i 1))) ((= i len))
(vector-set! new-pv i (vector-ref pv i)))
(set! len new-len) (set! pv new-pv)))
(cond ((negative? n) 0) ((zero? n) 1)
((and (< n len) (vector-ref pv n)) => (lambda (x) x))
(else (let loop ((k 1) (sum 0))
(if (< n k) (begin (vector-set! pv n sum) sum)
(loop (+ k 1) (+ sum (* (expt -1 (+ k 1))
(+ (partition (- n (* k (- (* 3 k) 1) 1/2)))
(partition (- n (* k (+ (* 3 k) 1)
1/2))))))))))))))
(define random-prime ; random prime with n digits in base b
(let (;(seed (time-second (current-time))) ; chez
(seed (current-seconds)) ; racket/chicken
;(seed (current-time)) ; guile
;(seed (inexact->exact (round (time->seconds (current-time))))) ; gambit
(rand ; knuth linear congruential method
(let* ((a 69069) (c 1234567) (m 4294967296))
(lambda () (set! seed (modulo (+ (* a seed) c) m))
(/ seed m))))
(randint (lambda (lo hi) (+ lo (floor (* (rand) (- hi lo)))))))
(lambda (n . base)
(let ((b (if (pair? base) (car base) 2)))
(let loop ((p (randint 1 b)) (n (- n 1)))
(if (zero? n) (prev-prime p)
(loop (+ (* p b) (randint 0 b)) (- n 1))))))))
(define (aliquot n . args) ; compute aliquot sequence (default verbose)
(let ((verbose? (if (pair? args) (car args) #t))) ; try 138, 840 or 3630
(let loop ((s n) (ss (list n)) (k 0) (fs (factors n)))
(if verbose? (for-each display `(,n " " ,k " " ,s " " ,fs #\newline)))
(if (= s 1) (cadr ss)
(let ((s (- (sigma 1 s fs) s)) (k (+ k 1)))
(if (member s ss) (member s (reverse ss))
(loop s (cons s ss) k (factors s))))))))
; fermat's divisors challenges
(define (iroot k n)
(let ((k-1 (- k 1)))
(let loop ((u n) (s (+ n 1)))
(if (<= s u) s
(loop (quotient (+ (* k-1 u) (quotient n (expt u k-1))) k) u)))))
(define (primes n) ; list of primes not exceeding n
(let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
(let loop ((i 0) (p 3) (ps (list 2)))
(cond ((< n (* p p))
(do ((i i (+ i 1)) (p p (+ p 2))
(ps ps (if (vector-ref bits i) (cons p ps) ps)))
((= i len) (reverse ps))))
((vector-ref bits i)
(do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
(vector-set! bits j #f)))
(else (loop (+ i 1) (+ p 2) ps))))))
(define prime? ; strong pseudoprime to prime bases less than 100
(let* ((ps (primes 100)) (p100 (apply * ps)))
(lambda (n)
(define (expm b e m)
(let loop ((b b) (e e) (x 1))
(if (zero? e) x
(loop (modulo (* b b) m) (quotient e 2)
(if (odd? e) (modulo (* b x) m) x)))))
(define (spsp? n a) ; #t if n is a strong pseudoprime base a
(do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
((odd? d) (if (= (expm a d n) 1) #t
(do ((r 0 (+ r 1)))
((or (= (expm a (* d (expt 2 r)) n) (- n 1)) (= r s))
(< r s)))))))
(if (< n 2) #f (if (< 1 (gcd n p100)) (if (member n ps) #t #f)
(do ((ps ps (cdr ps)))
((or (null? ps) (not (spsp? n (car ps)))) (null? ps))))))))
(define (next-prime n) ; smallest prime larger than n
(define (wheel n)
(vector-ref (vector 1 6 5 4 3 2 1 4 3 2 1 2 1 4
3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 2) (modulo n 30)))
(if (< n 2) 2 (if (< n 3) 3 (if (< n 5) 5
(let loop ((p (+ n (wheel n))))
(if (prime? p) p (loop (+ p (wheel p)))))))))
(define (prev-prime n) ; largest prime smaller than n
(define (wheel n)
(vector-ref (vector 1 2 1 2 3 4 5 6 1 2 3 4 1 2
1 2 3 4 1 2 1 2 3 4 1 2 3 4 5 6) (modulo n 30)))
(if (<= n 2) #f (if (<= n 3) 2 (if (<= n 5) 3 (if (<= n 7) 5
(let loop ((p (- n (wheel n))))
(if (prime? p) p (loop (- p (wheel p))))))))))
(define (factors n) ; factors of n in increasing order
(define (last-pair xs) (if (null? (cdr xs)) xs (last-pair (cdr xs))))
(define (cycle . xs) (set-cdr! (last-pair xs) xs) xs)
(define wheel235 (cons 1 (cons 2 (cons 2 (cycle 4 2 4 2 4 6 2 6)))))
(define (cons< x xs)
(if (null? xs) (list x)
(if (< x (car xs)) (cons x xs)
(cons (car xs) (cons< x (cdr xs))))))
(let wheel ((n n) (f 2) (fs (list)) (w wheel235))
(cond ((< n (* f f)) (reverse (cons n fs)))
((zero? (modulo n f)) (wheel (quotient n f) f (cons f fs) w))
((< f 1000) (wheel n (+ f (car w)) fs (cdr w)))
(else (let rho ((n n) (t 1) (h 1) (d 1) (c 1) (fs (reverse fs)))
(cond ((prime? n) (cons< n fs))
((= d 1)
(let* ((h (modulo (+ (* h h) c) n))
(h (modulo (+ (* h h) c) n))
(t (modulo (+ (* t t) c) n)))
(rho n t h (gcd (- t h) n) c fs)))
((= d n) (rho n 1 1 1 (+ c 1) fs))
((prime? d)
(rho (/ n d) 1 1 1 (+ c 1) (cons< d fs)))
(else (rho n 1 1 1 (+ c 1) fs))))))))
(define (divisors n) ; divisors of n, including 1 and n
(let ((divs (list 1)))
(do ((fs (factors n) (cdr fs))) ((null? fs) (sort < divs))
(let ((temp (list)))
(do ((ds divs (cdr ds)))
((null? ds) (set! divs (append divs temp)))
(let ((d (* (car fs) (car ds))))
(when (not (member d divs))
(set! temp (cons d temp)))))))))
(define (sigma x n . args) ; sum of x'th powers of divisors of n
(define (uniq-c eql? xs)
(if (null? xs) xs
(let loop ((xs (cdr xs)) (prev (car xs)) (k 1) (result '()))
(cond ((null? xs) (reverse (cons (cons prev k) result)))
((eql? (car xs) prev) (loop (cdr xs) prev (+ k 1) result))
(else (loop (cdr xs) (car xs) 1 (cons (cons prev k) result)))))))
(define (prod xs) (apply * xs))
(if (= n 1) 1
(let ((fs (uniq-c = (if (pair? args) (car args) (factors n)))))
(if (zero? x)
(prod (map add1 (map cdr fs)))
(prod (map (lambda (p a)
(/ (- (expt p (* (+ a 1) x)) 1) (- (expt p x) 1)))
(map car fs) (map cdr fs)))))))
(define (totient n) ; count of positive integers less than n coprime to it
(let loop ((t n) (p 0) (fs (factors n)))
(if (null? fs) t
(let ((f (car fs)))
(loop (if (= f p) t (* t (/ (- f 1) f))) f (cdr fs))))))
(define (moebius n) ; (-1)^k if n has k factors, or 0 if any factors duplicated
(let loop ((m -1) (f 0) (fs (factors n)))
(if (null? fs) m
(if (= f (car fs)) 0
(loop (- m) (car fs) (cdr fs))))))
(define (fermat n) (+ (expt 2 (expt 2 n)) 1)) ; 2^2^n+1
(define (mersenne n) (- (expt 2 n) 1)) ; 2^n-1
(define (repunit n . bs) ; (b^n-1)/(b-1)
(let ((b (if (pair? bs) (car bs) 10)))
(/ (- (expt b n) 1) (- b 1))))
(define (factorial n) ; product(1..n)
(let loop ((n n) (f 1))
(if (zero? n) f
(loop (- n 1) (* f n)))))
(define (primorial n) ; product(2..p_k..p_n) with p prime
(let loop ((k 0) (p 2) (prim 1))
(if (= k n) prim
(loop (+ k 1) (next-prime p) (* prim p)))))
(define (fibonacci n) ; F(n) = F(n-1) + Fb(n-2) with F(1) = F(2) = 1
(letrec ((square (lambda (x) (* x x))))
(cond ((= n 0) 0) ((= n 1) 1) ((= n 2) 1)
((odd? n) (let ((n2 (+ (quotient n 2) 1)))
(+ (square (fibonacci n2))
(square (fibonacci (- n2 1))))))
(else (let ((n2 (quotient n 2)))
(* (fibonacci n2) (+ (* (fibonacci (- n2 1)) 2)
(fibonacci n2))))))))
(define (lucas n) ; L(n) = L(n-1) + L(n-2) with L(1) = 1 and L(2) = 3
(/ (fibonacci (+ n n)) (fibonacci n)))
(define partition ; number of ways of writing n as sum of integers > 0
(let ((len 1) (pv (vector 1)))
(lambda (n)
(do () ((< n len))
(let* ((new-len (+ len len))
(new-pv (make-vector new-len #f)))
(do ((i 0 (+ i 1))) ((= i len))
(vector-set! new-pv i (vector-ref pv i)))
(set! len new-len) (set! pv new-pv)))
(cond ((negative? n) 0) ((zero? n) 1)
((and (< n len) (vector-ref pv n)) => (lambda (x) x))
(else (let loop ((k 1) (sum 0))
(if (< n k) (begin (vector-set! pv n sum) sum)
(loop (+ k 1) (+ sum (* (expt -1 (+ k 1))
(+ (partition (- n (* k (- (* 3 k) 1) 1/2)))
(partition (- n (* k (+ (* 3 k) 1)
1/2))))))))))))))
(define random-prime ; random prime with n digits in base b
(let (;(seed (time-second (current-time))) ; chez
(seed (current-seconds)) ; racket/chicken
;(seed (current-time)) ; guile
;(seed (inexact->exact (round (time->seconds (current-time))))) ; gambit
(rand ; knuth linear congruential method
(let* ((a 69069) (c 1234567) (m 4294967296))
(lambda () (set! seed (modulo (+ (* a seed) c) m))
(/ seed m))))
(randint (lambda (lo hi) (+ lo (floor (* (rand) (- hi lo)))))))
(lambda (n . base)
(let ((b (if (pair? base) (car base) 2)))
(let loop ((p (randint 1 b)) (n (- n 1)))
(if (zero? n) (prev-prime p)
(loop (+ (* p b) (randint 0 b)) (- n 1))))))))
(define (aliquot n . args) ; compute aliquot sequence (default verbose)
(let ((verbose? (if (pair? args) (car args) #t))) ; try 138, 840 or 3630
(let loop ((s n) (ss (list n)) (k 0) (fs (factors n)))
(if verbose? (for-each display `(,n " " ,k " " ,s " " ,fs #\newline)))
(if (= s 1) (cadr ss)
(let ((s (- (sigma 1 s fs) s)) (k (+ k 1)))
(if (member s ss) (member s (reverse ss))
(loop s (cons s ss) k (factors s))))))))
(do ((n 1 (+ n 1))) ((< 10 n))
(let* ((x (sigma 1 (* n n n)))
(r (iroot 2 x)))
(when (= x (* r r))
(display n) (display "^3 = ")
(display r) (display "^2") (newline))))
(do ((n 43097 (+ n 1))) ((< 43099 n))
(let* ((x (sigma 1 (* n n)))
(r (iroot 3 x)))
(when (= x (* r r r))
(display n) (display "^2 = ")
(display r) (display "^3") (newline))))