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

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: