fork(1) download
copy 
  1. ; fermat's divisors challenges
  2.  
  3. (define (iroot k n)
  4.   (let ((k-1 (- k 1)))
  5.     (let loop ((u n) (s (+ n 1)))
  6.       (if (<= s u) s
  7.         (loop (quotient (+ (* k-1 u) (quotient n (expt u k-1))) k) u)))))
  8.  
  9. (define (primes n) ; list of primes not exceeding n
  10.   (let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
  11.     (let loop ((i 0) (p 3) (ps (list 2)))
  12.       (cond ((< n (* p p))
  13.               (do ((i i (+ i 1)) (p p (+ p 2))
  14.                    (ps ps (if (vector-ref bits i) (cons p ps) ps)))
  15.                   ((= i len) (reverse ps))))
  16.             ((vector-ref bits i)
  17.               (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
  18.                   ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
  19.                 (vector-set! bits j #f)))
  20.             (else (loop (+ i 1) (+ p 2) ps))))))
  21.  
  22. (define prime? ; strong pseudoprime to prime bases less than 100
  23.   (let* ((ps (primes 100)) (p100 (apply * ps)))
  24.     (lambda (n)
  25.       (define (expm b e m)
  26.         (let loop ((b b) (e e) (x 1))
  27.           (if (zero? e) x
  28.             (loop (modulo (* b b) m) (quotient e 2)
  29.                   (if (odd? e) (modulo (* b x) m) x)))))
  30.       (define (spsp? n a) ; #t if n is a strong pseudoprime base a
  31.         (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
  32.             ((odd? d) (if (= (expm a d n) 1) #t
  33.               (do ((r 0 (+ r 1)))
  34.                   ((or (= (expm a (* d (expt 2 r)) n) (- n 1)) (= r s))
  35.                     (< r s)))))))
  36.       (if (< n 2) #f (if (< 1 (gcd n p100)) (if (member n ps) #t #f)
  37.         (do ((ps ps (cdr ps)))
  38.             ((or (null? ps) (not (spsp? n (car ps)))) (null? ps))))))))
  39.  
  40. (define (next-prime n) ; smallest prime larger than n
  41.   (define (wheel n)
  42.     (vector-ref (vector 1 6 5 4 3 2 1 4 3 2 1 2 1 4
  43.       3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 2) (modulo n 30)))
  44.   (if (< n 2) 2 (if (< n 3) 3 (if (< n 5) 5
  45.     (let loop ((p (+ n (wheel n))))
  46.       (if (prime? p) p (loop (+ p (wheel p)))))))))
  47.  
  48. (define (prev-prime n) ; largest prime smaller than n
  49.   (define (wheel n)
  50.     (vector-ref (vector 1 2 1 2 3 4 5 6 1 2 3 4 1 2
  51.       1 2 3 4 1 2 1 2 3 4 1 2 3 4 5 6) (modulo n 30)))
  52.   (if (<= n 2) #f (if (<= n 3) 2 (if (<= n 5) 3 (if (<= n 7) 5
  53.     (let loop ((p (- n (wheel n))))
  54.       (if (prime? p) p (loop (- p (wheel p))))))))))
  55.  
  56. (define (factors n) ; factors of n in increasing order
  57.   (define (last-pair xs) (if (null? (cdr xs)) xs (last-pair (cdr xs))))
  58.   (define (cycle . xs) (set-cdr! (last-pair xs) xs) xs)
  59.   (define wheel235 (cons 1 (cons 2 (cons 2 (cycle 4 2 4 2 4 6 2 6)))))
  60.   (define (cons< x xs)
  61.     (if (null? xs) (list x)
  62.       (if (< x (car xs)) (cons x xs)
  63.         (cons (car xs) (cons< x (cdr xs))))))
  64.   (let wheel ((n n) (f 2) (fs (list)) (w wheel235))
  65.     (cond ((< n (* f f)) (reverse (cons n fs)))
  66.           ((zero? (modulo n f)) (wheel (quotient n f) f (cons f fs) w))
  67.           ((< f 1000) (wheel n (+ f (car w)) fs (cdr w)))
  68.           (else (let rho ((n n) (t 1) (h 1) (d 1) (c 1) (fs (reverse fs)))
  69.                   (cond ((prime? n) (cons< n fs))
  70.                         ((= d 1)
  71.                           (let* ((h (modulo (+ (* h h) c) n))
  72.                                  (h (modulo (+ (* h h) c) n))
  73.                                  (t (modulo (+ (* t t) c) n)))
  74.                             (rho n t h (gcd (- t h) n) c fs)))
  75.                         ((= d n) (rho n 1 1 1 (+ c 1) fs))
  76.                         ((prime? d)
  77.                           (rho (/ n d) 1 1 1 (+ c 1) (cons< d fs)))
  78.                         (else (rho n 1 1 1 (+ c 1) fs))))))))
  79.  
  80. (define (divisors n) ; divisors of n, including 1 and n
  81.   (let ((divs (list 1)))
  82.     (do ((fs (factors n) (cdr fs))) ((null? fs) (sort < divs))
  83.       (let ((temp (list)))
  84.         (do ((ds divs (cdr ds)))
  85.             ((null? ds) (set! divs (append divs temp)))
  86.           (let ((d (* (car fs) (car ds))))
  87.             (when (not (member d divs))
  88.               (set! temp (cons d temp)))))))))
  89.  
  90. (define (sigma x n . args) ; sum of x'th powers of divisors of n
  91.   (define (uniq-c eql? xs)
  92.     (if (null? xs) xs
  93.       (let loop ((xs (cdr xs)) (prev (car xs)) (k 1) (result '()))
  94.         (cond ((null? xs) (reverse (cons (cons prev k) result)))
  95.               ((eql? (car xs) prev) (loop (cdr xs) prev (+ k 1) result))
  96.               (else (loop (cdr xs) (car xs) 1 (cons (cons prev k) result)))))))
  97.   (define (prod xs) (apply * xs))
  98.   (if (= n 1) 1
  99.     (let ((fs (uniq-c = (if (pair? args) (car args) (factors n)))))
  100.       (if (zero? x)
  101.           (prod (map add1 (map cdr fs)))
  102.           (prod (map (lambda (p a)
  103.                        (/ (- (expt p (* (+ a 1) x)) 1) (- (expt p x) 1)))
  104.                      (map car fs) (map cdr fs)))))))
  105.  
  106. (define (totient n) ; count of positive integers less than n coprime to it
  107.   (let loop ((t n) (p 0) (fs (factors n)))
  108.     (if (null? fs) t
  109.       (let ((f (car fs)))
  110.         (loop (if (= f p) t (* t (/ (- f 1) f))) f (cdr fs))))))
  111.  
  112. (define (moebius n) ; (-1)^k if n has k factors, or 0 if any factors duplicated
  113.   (let loop ((m -1) (f 0) (fs (factors n)))
  114.     (if (null? fs) m
  115.       (if (= f (car fs)) 0
  116.         (loop (- m) (car fs) (cdr fs))))))
  117.  
  118. (define (fermat n) (+ (expt 2 (expt 2 n)) 1)) ; 2^2^n+1
  119.  
  120. (define (mersenne n) (- (expt 2 n) 1)) ; 2^n-1
  121.  
  122. (define (repunit n . bs) ; (b^n-1)/(b-1)
  123.  (let ((b (if (pair? bs) (car bs) 10)))
  124.    (/ (- (expt b n) 1) (- b 1))))
  125.  
  126. (define (factorial n) ; product(1..n)
  127.   (let loop ((n n) (f 1))
  128.     (if (zero? n) f
  129.       (loop (- n 1) (* f n)))))
  130.  
  131. (define (primorial n) ; product(2..p_k..p_n) with p prime
  132.   (let loop ((k 0) (p 2) (prim 1))
  133.     (if (= k n) prim
  134.       (loop (+ k 1) (next-prime p) (* prim p)))))
  135.  
  136. (define (fibonacci n) ; F(n) = F(n-1) + Fb(n-2) with F(1) = F(2) = 1
  137.   (letrec ((square (lambda (x) (* x x))))
  138.     (cond ((= n 0) 0) ((= n 1) 1) ((= n 2) 1)
  139.           ((odd? n) (let ((n2 (+ (quotient n 2) 1)))
  140.                       (+ (square (fibonacci n2))
  141.                          (square (fibonacci (- n2 1))))))
  142.           (else (let ((n2 (quotient n 2)))
  143.                   (* (fibonacci n2) (+ (* (fibonacci (- n2 1)) 2)
  144.                                           (fibonacci n2))))))))
  145.  
  146. (define (lucas n) ; L(n) = L(n-1) + L(n-2) with L(1) = 1 and L(2) = 3
  147.   (/ (fibonacci (+ n n)) (fibonacci n)))
  148.  
  149. (define partition ; number of ways of writing n as sum of integers > 0
  150.   (let ((len 1) (pv (vector 1)))
  151.     (lambda (n)
  152.       (do () ((< n len))
  153.         (let* ((new-len (+ len len))
  154.                (new-pv (make-vector new-len #f)))
  155.           (do ((i 0 (+ i 1))) ((= i len))
  156.             (vector-set! new-pv i (vector-ref pv i)))
  157.           (set! len new-len) (set! pv new-pv)))
  158.       (cond ((negative? n) 0) ((zero? n) 1)
  159.             ((and (< n len) (vector-ref pv n)) => (lambda (x) x))
  160.             (else (let loop ((k 1) (sum 0))
  161.                     (if (< n k) (begin (vector-set! pv n sum) sum)
  162.                       (loop (+ k 1) (+ sum (* (expt -1 (+ k 1))
  163.                         (+ (partition (- n (* k (- (* 3 k) 1) 1/2)))
  164.                            (partition (- n (* k (+ (* 3 k) 1)
  165.                                          1/2))))))))))))))
  166.  
  167. (define random-prime ; random prime with n digits in base b
  168.   (let (;(seed (time-second (current-time))) ; chez
  169.          (seed (current-seconds)) ; racket/chicken
  170.         ;(seed (current-time)) ; guile
  171.         ;(seed (inexact->exact (round (time->seconds (current-time))))) ; gambit
  172.         (rand ; knuth linear congruential method
  173.           (let* ((a 69069) (c 1234567) (m 4294967296))
  174.             (lambda () (set! seed (modulo (+ (* a seed) c) m))
  175.                        (/ seed m))))
  176.         (randint (lambda (lo hi) (+ lo (floor (* (rand) (- hi lo)))))))
  177.     (lambda (n . base)
  178.       (let ((b (if (pair? base) (car base) 2)))
  179.         (let loop ((p (randint 1 b)) (n (- n 1)))
  180.           (if (zero? n) (prev-prime p)
  181.             (loop (+ (* p b) (randint 0 b)) (- n 1))))))))
  182.  
  183. (define (aliquot n . args) ; compute aliquot sequence (default verbose)
  184.   (let ((verbose? (if (pair? args) (car args) #t))) ; try 138, 840 or 3630
  185.     (let loop ((s n) (ss (list n)) (k  0) (fs (factors n)))
  186.       (if verbose? (for-each display `(,n " " ,k " " ,s " " ,fs #\newline)))
  187.       (if (= s 1) (cadr ss)
  188.         (let ((s (- (sigma 1 s fs) s)) (k (+ k 1)))
  189.           (if (member s ss) (member s (reverse ss))
  190.             (loop s (cons s ss) k (factors s))))))))
  191.  
  192. ; fermat's divisors challenges
  193.  
  194. (define (iroot k n)
  195.   (let ((k-1 (- k 1)))
  196.     (let loop ((u n) (s (+ n 1)))
  197.       (if (<= s u) s
  198.         (loop (quotient (+ (* k-1 u) (quotient n (expt u k-1))) k) u)))))
  199.  
  200. (define (primes n) ; list of primes not exceeding n
  201.   (let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
  202.     (let loop ((i 0) (p 3) (ps (list 2)))
  203.       (cond ((< n (* p p))
  204.               (do ((i i (+ i 1)) (p p (+ p 2))
  205.                    (ps ps (if (vector-ref bits i) (cons p ps) ps)))
  206.                   ((= i len) (reverse ps))))
  207.             ((vector-ref bits i)
  208.               (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
  209.                   ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
  210.                 (vector-set! bits j #f)))
  211.             (else (loop (+ i 1) (+ p 2) ps))))))
  212.  
  213. (define prime? ; strong pseudoprime to prime bases less than 100
  214.   (let* ((ps (primes 100)) (p100 (apply * ps)))
  215.     (lambda (n)
  216.       (define (expm b e m)
  217.         (let loop ((b b) (e e) (x 1))
  218.           (if (zero? e) x
  219.             (loop (modulo (* b b) m) (quotient e 2)
  220.                   (if (odd? e) (modulo (* b x) m) x)))))
  221.       (define (spsp? n a) ; #t if n is a strong pseudoprime base a
  222.         (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
  223.             ((odd? d) (if (= (expm a d n) 1) #t
  224.               (do ((r 0 (+ r 1)))
  225.                   ((or (= (expm a (* d (expt 2 r)) n) (- n 1)) (= r s))
  226.                     (< r s)))))))
  227.       (if (< n 2) #f (if (< 1 (gcd n p100)) (if (member n ps) #t #f)
  228.         (do ((ps ps (cdr ps)))
  229.             ((or (null? ps) (not (spsp? n (car ps)))) (null? ps))))))))
  230.  
  231. (define (next-prime n) ; smallest prime larger than n
  232.   (define (wheel n)
  233.     (vector-ref (vector 1 6 5 4 3 2 1 4 3 2 1 2 1 4
  234.       3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 2) (modulo n 30)))
  235.   (if (< n 2) 2 (if (< n 3) 3 (if (< n 5) 5
  236.     (let loop ((p (+ n (wheel n))))
  237.       (if (prime? p) p (loop (+ p (wheel p)))))))))
  238.  
  239. (define (prev-prime n) ; largest prime smaller than n
  240.   (define (wheel n)
  241.     (vector-ref (vector 1 2 1 2 3 4 5 6 1 2 3 4 1 2
  242.       1 2 3 4 1 2 1 2 3 4 1 2 3 4 5 6) (modulo n 30)))
  243.   (if (<= n 2) #f (if (<= n 3) 2 (if (<= n 5) 3 (if (<= n 7) 5
  244.     (let loop ((p (- n (wheel n))))
  245.       (if (prime? p) p (loop (- p (wheel p))))))))))
  246.  
  247. (define (factors n) ; factors of n in increasing order
  248.   (define (last-pair xs) (if (null? (cdr xs)) xs (last-pair (cdr xs))))
  249.   (define (cycle . xs) (set-cdr! (last-pair xs) xs) xs)
  250.   (define wheel235 (cons 1 (cons 2 (cons 2 (cycle 4 2 4 2 4 6 2 6)))))
  251.   (define (cons< x xs)
  252.     (if (null? xs) (list x)
  253.       (if (< x (car xs)) (cons x xs)
  254.         (cons (car xs) (cons< x (cdr xs))))))
  255.   (let wheel ((n n) (f 2) (fs (list)) (w wheel235))
  256.     (cond ((< n (* f f)) (reverse (cons n fs)))
  257.           ((zero? (modulo n f)) (wheel (quotient n f) f (cons f fs) w))
  258.           ((< f 1000) (wheel n (+ f (car w)) fs (cdr w)))
  259.           (else (let rho ((n n) (t 1) (h 1) (d 1) (c 1) (fs (reverse fs)))
  260.                   (cond ((prime? n) (cons< n fs))
  261.                         ((= d 1)
  262.                           (let* ((h (modulo (+ (* h h) c) n))
  263.                                  (h (modulo (+ (* h h) c) n))
  264.                                  (t (modulo (+ (* t t) c) n)))
  265.                             (rho n t h (gcd (- t h) n) c fs)))
  266.                         ((= d n) (rho n 1 1 1 (+ c 1) fs))
  267.                         ((prime? d)
  268.                           (rho (/ n d) 1 1 1 (+ c 1) (cons< d fs)))
  269.                         (else (rho n 1 1 1 (+ c 1) fs))))))))
  270.  
  271. (define (divisors n) ; divisors of n, including 1 and n
  272.   (let ((divs (list 1)))
  273.     (do ((fs (factors n) (cdr fs))) ((null? fs) (sort < divs))
  274.       (let ((temp (list)))
  275.         (do ((ds divs (cdr ds)))
  276.             ((null? ds) (set! divs (append divs temp)))
  277.           (let ((d (* (car fs) (car ds))))
  278.             (when (not (member d divs))
  279.               (set! temp (cons d temp)))))))))
  280.  
  281. (define (sigma x n . args) ; sum of x'th powers of divisors of n
  282.   (define (uniq-c eql? xs)
  283.     (if (null? xs) xs
  284.       (let loop ((xs (cdr xs)) (prev (car xs)) (k 1) (result '()))
  285.         (cond ((null? xs) (reverse (cons (cons prev k) result)))
  286.               ((eql? (car xs) prev) (loop (cdr xs) prev (+ k 1) result))
  287.               (else (loop (cdr xs) (car xs) 1 (cons (cons prev k) result)))))))
  288.   (define (prod xs) (apply * xs))
  289.   (if (= n 1) 1
  290.     (let ((fs (uniq-c = (if (pair? args) (car args) (factors n)))))
  291.       (if (zero? x)
  292.           (prod (map add1 (map cdr fs)))
  293.           (prod (map (lambda (p a)
  294.                        (/ (- (expt p (* (+ a 1) x)) 1) (- (expt p x) 1)))
  295.                      (map car fs) (map cdr fs)))))))
  296.  
  297. (define (totient n) ; count of positive integers less than n coprime to it
  298.   (let loop ((t n) (p 0) (fs (factors n)))
  299.     (if (null? fs) t
  300.       (let ((f (car fs)))
  301.         (loop (if (= f p) t (* t (/ (- f 1) f))) f (cdr fs))))))
  302.  
  303. (define (moebius n) ; (-1)^k if n has k factors, or 0 if any factors duplicated
  304.   (let loop ((m -1) (f 0) (fs (factors n)))
  305.     (if (null? fs) m
  306.       (if (= f (car fs)) 0
  307.         (loop (- m) (car fs) (cdr fs))))))
  308.  
  309. (define (fermat n) (+ (expt 2 (expt 2 n)) 1)) ; 2^2^n+1
  310.  
  311. (define (mersenne n) (- (expt 2 n) 1)) ; 2^n-1
  312.  
  313. (define (repunit n . bs) ; (b^n-1)/(b-1)
  314.  (let ((b (if (pair? bs) (car bs) 10)))
  315.    (/ (- (expt b n) 1) (- b 1))))
  316.  
  317. (define (factorial n) ; product(1..n)
  318.   (let loop ((n n) (f 1))
  319.     (if (zero? n) f
  320.       (loop (- n 1) (* f n)))))
  321.  
  322. (define (primorial n) ; product(2..p_k..p_n) with p prime
  323.   (let loop ((k 0) (p 2) (prim 1))
  324.     (if (= k n) prim
  325.       (loop (+ k 1) (next-prime p) (* prim p)))))
  326.  
  327. (define (fibonacci n) ; F(n) = F(n-1) + Fb(n-2) with F(1) = F(2) = 1
  328.   (letrec ((square (lambda (x) (* x x))))
  329.     (cond ((= n 0) 0) ((= n 1) 1) ((= n 2) 1)
  330.           ((odd? n) (let ((n2 (+ (quotient n 2) 1)))
  331.                       (+ (square (fibonacci n2))
  332.                          (square (fibonacci (- n2 1))))))
  333.           (else (let ((n2 (quotient n 2)))
  334.                   (* (fibonacci n2) (+ (* (fibonacci (- n2 1)) 2)
  335.                                           (fibonacci n2))))))))
  336.  
  337. (define (lucas n) ; L(n) = L(n-1) + L(n-2) with L(1) = 1 and L(2) = 3
  338.   (/ (fibonacci (+ n n)) (fibonacci n)))
  339.  
  340. (define partition ; number of ways of writing n as sum of integers > 0
  341.   (let ((len 1) (pv (vector 1)))
  342.     (lambda (n)
  343.       (do () ((< n len))
  344.         (let* ((new-len (+ len len))
  345.                (new-pv (make-vector new-len #f)))
  346.           (do ((i 0 (+ i 1))) ((= i len))
  347.             (vector-set! new-pv i (vector-ref pv i)))
  348.           (set! len new-len) (set! pv new-pv)))
  349.       (cond ((negative? n) 0) ((zero? n) 1)
  350.             ((and (< n len) (vector-ref pv n)) => (lambda (x) x))
  351.             (else (let loop ((k 1) (sum 0))
  352.                     (if (< n k) (begin (vector-set! pv n sum) sum)
  353.                       (loop (+ k 1) (+ sum (* (expt -1 (+ k 1))
  354.                         (+ (partition (- n (* k (- (* 3 k) 1) 1/2)))
  355.                            (partition (- n (* k (+ (* 3 k) 1)
  356.                                          1/2))))))))))))))
  357.  
  358. (define random-prime ; random prime with n digits in base b
  359.   (let (;(seed (time-second (current-time))) ; chez
  360.          (seed (current-seconds)) ; racket/chicken
  361.         ;(seed (current-time)) ; guile
  362.         ;(seed (inexact->exact (round (time->seconds (current-time))))) ; gambit
  363.         (rand ; knuth linear congruential method
  364.           (let* ((a 69069) (c 1234567) (m 4294967296))
  365.             (lambda () (set! seed (modulo (+ (* a seed) c) m))
  366.                        (/ seed m))))
  367.         (randint (lambda (lo hi) (+ lo (floor (* (rand) (- hi lo)))))))
  368.     (lambda (n . base)
  369.       (let ((b (if (pair? base) (car base) 2)))
  370.         (let loop ((p (randint 1 b)) (n (- n 1)))
  371.           (if (zero? n) (prev-prime p)
  372.             (loop (+ (* p b) (randint 0 b)) (- n 1))))))))
  373.  
  374. (define (aliquot n . args) ; compute aliquot sequence (default verbose)
  375.   (let ((verbose? (if (pair? args) (car args) #t))) ; try 138, 840 or 3630
  376.     (let loop ((s n) (ss (list n)) (k  0) (fs (factors n)))
  377.       (if verbose? (for-each display `(,n " " ,k " " ,s " " ,fs #\newline)))
  378.       (if (= s 1) (cadr ss)
  379.         (let ((s (- (sigma 1 s fs) s)) (k (+ k 1)))
  380.           (if (member s ss) (member s (reverse ss))
  381.             (loop s (cons s ss) k (factors s))))))))
  382.  
  383. (do ((n 1 (+ n 1))) ((< 10 n))
  384.   (let* ((x (sigma 1 (* n n n)))
  385.          (r (iroot 2 x)))
  386.     (when (= x (* r r))
  387.       (display n) (display "^3 = ")
  388.       (display r) (display "^2") (newline))))
  389.  
  390. (do ((n 43097 (+ n 1))) ((< 43099 n))
  391.   (let* ((x (sigma 1 (* n n)))
  392.          (r (iroot 3 x)))
  393.     (when (= x (* r r r))
  394.       (display n) (display "^2 = ")
  395.       (display r) (display "^3") (newline))))
Success #stdin #stdout 0.05s 7272KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
1^3 = 1^2
7^3 = 20^2
43098^2 = 1729.0^3
https://ideone.com/3iqiIH
language:
Scheme (chicken 4.13)
created:
8 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!