;; ;; Two Scheme-implementations for computing A000203, sigma(n), the sum of the divisors of n. ;; ;; Written by Antti Karttunen over the years. The first version of this stand-alone source file prepared 2017-11-27. ;; See also https://github.com/karttu/IntSeq ;; ;; (declare (usual-integrations)) ;; First follows an implementation of memoization-macro definec. ;; See also http://oeis.org/wiki/Memoization#Scheme ;; Added this 10. July 2002 to avoid allocation catastrophes ;; caused by the careless use of cached integer functions: (define *MAX-CACHE-SIZE-FOR-DEFINEC* 362881) ;; Was: 290512) ;; Was 131072 (define-syntax definec (syntax-rules () ((definec (name arg) e0 ...) (define name (implement-cached-function *MAX-CACHE-SIZE-FOR-DEFINEC* (name arg) e0 ...) ) ;; (define name ...) ) ) ;; syntax-rules ) (define-syntax grow-cache (syntax-rules () ((grow-cache cachename arg) ;; No maxsize specified. (vector-grow cachename (max (1+ arg) (* 2 (vector-length cachename)))) ) ((grow-cache cachename arg 0) ;; Or specified as zero. (vector-grow cachename (max (1+ arg) (* 2 (vector-length cachename)))) ) ((grow-cache cachename arg maxsize) (vector-grow cachename (min maxsize (max (1+ arg) (* 2 (vector-length cachename))))) ) ) ) (define-syntax implement-cached-function (syntax-rules () ((implement-cached-function maxcachesize (funname argname) e0 ...) (letrec ((_cache_ (vector #f #f #f #f #f #f #f #f #f #f #f #f #f #f #f #f)) (funname (lambda (argname) (cond ((null? argname) _cache_) ;; For debugging. ((vector? argname) argname) ;; As well as this: Caches for caches! ((and (not (= 0 maxcachesize)) (>= argname maxcachesize)) e0 ... ) (else (if (>= argname (vector-length _cache_)) (set! _cache_ (grow-cache _cache_ argname maxcachesize)) ) (or (vector-ref _cache_ argname) ((lambda (res) (vector-set! _cache_ argname res) res ) (begin e0 ...) ) ) ) ) ; cond ) ) ) ; letrec-definitions funname ) ; letrec ) ) ) ;; MIT/GNU Scheme has a function vector-grow ;; For some other Schemes you need to implement it separately. ;; Here is one implementation from ;; https://github.com/karttu/IntSeq/blob/master/src/utils/vector-grow.scm ;; ;; (define (vector-grow old-vec new-size) ;; (let ((new-vec (make-vector new-size #f)) ;; (old-size (vector-length old-vec)) ;; ) ;; (let copyloop ((i 0)) ;; (cond ((= i old-size) new-vec) ;; (else ;; (begin ;; (vector-set! new-vec i (vector-ref old-vec i)) ;; (copyloop (+ 1 i)) ;; ) ;; ) ;; ) ;; ) ;; ) ;; ) ;; ;; Comment the following two lines out if you are not using A000203old version of sigma. ;; (Or at least edit an appropriate path there, after installing SLIB !) (load "c:\\program files (x86)\\slib\\mitscheme.init") ;;SLIB Scheme library: http://people.csail.mit.edu/jaffer/SLIB (require 'factor) ;; From SLIB-library. ;; Our version of factor sorts the prime factors into ascending order and caches them also: (definec (ifactor n) (cond ((< n 2) (list)) (else (sort (factor n) <)))) ;; An example: ;; (ifactor 360) --> (2 2 2 3 3 5) ;; (elemcountpairs (ifactor 360)) --> ((2 . 3) (3 . 2) (5 . 1)) (define (elemcountpairs lista) ;; Of numeric elements, already sorted. (let loop ((pairs (list)) (lista lista) (prev #f) ) (cond ((not (pair? lista)) (reverse! pairs)) ((equal? (car lista) prev) (set-cdr! (car pairs) (+ 1 (cdar pairs))) (loop pairs (cdr lista) prev) ) (else ;; A new item? (loop (cons (cons (car lista) 1) pairs) (cdr lista) (car lista)) ) ) ) ) ;; %F A000203 Multiplicative with a(p^e) = (p^(e+1)-1)/(p-1). - David W. Wilson, Aug 01, 2001. (definec (A000203old n) (fold-left (lambda (prod p.e) (* prod (/ (- (expt (car p.e) (+ 1 (cdr p.e))) 1) (- (car p.e) 1)))) 1 (if (= 1 n) (list) (elemcountpairs (ifactor n))) ;; (sort (factor n) <) ) ) ;; The new recursive version that doesn't need any of the factorization machinery, thus no extra libraries: (definec (A000203 n) (if (= 1 n) n (let ((p (A020639 n)) (e (A067029 n))) (* (/ (- (expt p (+ 1 e)) 1) (- p 1)) (A000203 (A028234 n))) ) ) ) ;; A020639 [David W. Wilson] Lpf(n): least prime dividing n (with a(1)=1). ;; This implementation might seem excessively naive, but actually, because of the memoization ;; it works fine when used repeatedly with not too large arguments. ;; It is the function that actually offers the "factorization service" for these implementations. (definec (A020639 n) (if (< n 2) n (let loop ((k 2)) (cond ((zero? (modulo n k)) k) (else (loop (+ 1 k))) ) ) ) ) ;; A028233 [Marc LeBrun] o=1: If n = p_1^e_1 * ... * p_k^e_k, p_1 < ... < p_k primes, then a(n) = p_1^e_1 (definec (A028233 n) (if (< n 2) n (let ((lpf (A020639 n))) (let loop ((m lpf) (n (/ n lpf))) (cond ((not (zero? (modulo n lpf))) m) (else (loop (* m lpf) (/ n lpf))) ) ) ) ) ) ;; A028234 [Marc LeBrun] o=1: If n = p_1^e_1 * ... * p_k^e_k, p_1 < ... < p_k primes, then a(n) = n/p_1^e_1. (define (A028234 n) (/ n (A028233 n))) ;; A032742 [Patrick De Geest] o=1: a(1) = 1; for n > 1, a(n) = largest proper divisor of n. (define (A032742 n) (/ n (A020639 n))) ;; a(1) = 1; for n > 1, a(n) = largest proper divisor of n. ;; A067029 [Reinhard Zumkeller] o=1: Exponent of least prime factor in prime factorization of n, a(1)=0 (definec (A067029 n) (if (< n 2) 0 (let ((mp (A020639 n))) (let loop ((e 0) (n (/ n mp))) (cond ((integer? n) (loop (+ e 1) (/ n mp))) (else e) ) ) ) ) ) ;; ;; That's all for now folks! There are many, many other sequences that can be easily defined ;; by using just those A020639, A028233, A028234 and A067029. I might add them here one day. ;;