## Fermat’s Divisors Challenges

### September 8, 2015

(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))))))))

In Python.

Finding divisors could be more efficient (iterate up through sqrt(n), if i is a divisor, then so is n/i). Also, Python’s range() function doesnt like big numbers, so should be transformed in a while loop.

[…] makes solving the Fermat Exponent problem […]

#include

int main()

{

int i,sum,sq=0,l=0,n;

printf(“Enter the Range\n”);

scanf(“%d”,&n);

for(i=1;i<n;i++)

{

sum=1+i+i*i+i*i*i;

while(sq %d\n”,i,l);

l=1;

}

}

O/P:

Enter the Range

1000

1 -> 2

7 -> 20

EX 2.

#include

int main()

{

int i,sum,cu=0,l=0,n;

printf(“Enter the Range\n”);

scanf(“%d”,&n);

for(i=1;i<n;i++)

{

sum=1+i+i*i;

while(cu %d\n”,l,i);

l=1;

cu=0;

}

}

O/P:

Enter the Range

1000

7 -> 18