## 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))))))))```

Pages: 1 2 3

### 4 Responses to “Fermat’s Divisors Challenges”

1. Rutger said

In Python.

```def sum_divisors_is_square(cube):
sum_divisors = sum([i for i in range(1,cube/2+1) if not cube % i]+[cube])
return sum_divisors ** 0.5 ==  int(sum_divisors**0.5)

for i in range(10):
print i, sum_divisors_is_square(i**3)
```

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.

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

3. Mayur Lokare said

#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

4. Mayur Lokare said

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