PrimeETs.ss

;; LIBRARY

;; base-x log of a number n
(define logx (lambda (x n) (/ (log n) (log x))))

(define log2 (lambda (n) (logx 2 n)))

;; lists integers from i to j, inclusive
(define range
  (lambda (i j)
    (if (> i j)
      '()
      (cons i (range (+ i 1) j)))))

;; lists odd or even natural numbers <= n, matching the parity of n
(define odd-ls
  (lambda (n)
    (if (< n 1)
      (list)
      (cons n (odd-ls (- n 2))))))

(define mean
  (lambda (ls)
    (if (null? ls)
      0
      (/ (apply + ls) (length ls)))))

(define prime?
  (lambda (n)
    (letrec ((loop
               (lambda (n test max)
                 (cond
                   [(> test max) #t]
                   [(= (modulo n test) 0) #f]
                   [else (loop n (+ test 2) max)]))))
      (cond
        [(= n 2) #t]
        [(or (< n 2) (= (modulo n 2) 0)) #f]
        [else (loop n 3 (floor (sqrt n)))]))))

;; inverse of a prime counting function by Minac (Ribenboim 1995, p. 181)
;; returns the nth prime
(define pi-inv
  (lambda (n)
    (letrec ((loop
               (lambda (j acc n)
               (if (zero? n)
                 (- j 1)
                 (loop
                   (+ j 1)
                   (* j acc)
                   (- n (floor (- (/ (+ acc 1) j)
                                  (floor (/ acc j))))))))))
      (loop 2 1 n))))

(define remove-multiples
  (lambda (x ls)
    (cond
      [(null? ls) '()]
      [(= (modulo (car ls) x) 0) (remove-multiples x (cdr ls))]
      [else (cons (car ls) (remove-multiples x (cdr ls)))])))

(define n-primes
  (lambda (n)
    (letrec ((loop
               (lambda (ls)
                 (if (null? ls)
                   '()
                   (cons (car ls)
                     (loop (remove-multiples (car ls) (cdr ls))))))))
      (cond
        [(= n 0) '()]
        [(= n 1) '(2)]
        [(cons 2 (loop (cdr (reverse (odd-ls (pi-inv n))))))]))))

(define n-composites
  (lambda (n)
    (letrec ((inner
               (lambda (n int comps)
                 (cond
                   [(= n 0) (cddr (reverse comps))]
                   [(prime? int) (inner n (+ int 1) comps)]
                   [else (inner (- n 1) (+ int 1) (cons int comps))]))))
      (inner (+ n 2) 0 '()))))

;; MAIN FUNCTIONS

(define patent-val-div
  (lambda (et basis)
    (map (lambda (x)
           (inexact->exact (round (* et (logx (apply min basis) x)))))
      basis)))

(define remove-torsion
  (lambda (ls)
    (cond
      [(null? ls) '()]
      [(> (apply gcd (car ls)) 1) (remove-torsion (cdr ls))]
      [else (cons (car ls) (remove-torsion (cdr ls)))])))

;; TOP tuning of basis given val
(define top-val
  (lambda (val basis)
    (let* ((log2-primes (map log2 basis))
           (exact-val (map / val log2-primes))
           (big (apply max exact-val))
           (small (apply min exact-val))
           (mid (/ (+ big small) 2))
           (tuned-val (map (lambda (x) (/ x mid)) val)))
      (map (lambda (x) (* x 1200))
        tuned-val))))

(define top-damage
  (lambda (val basis)
    (let*
      ((weights (map log2 basis))
       (ji (map (lambda (x) (* x 1200)) weights))
       (errors (map abs (map - ji (top-val val basis)))))
       (apply max (map / errors weights)))))

(define best-val-err
  (lambda (et chord errorfunction)
    (apply min
      (map (lambda (val)
             (errorfunction val chord))
        (remove-torsion
          (map (lambda (exact-et)
                 (patent-val-div exact-et chord))
            (map (lambda (x) (+ et x))
              '(-0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4))))))))

(define range-prime-comp
  (lambda (min max)
    (letrec ((inner
               (lambda (min primes comps)
                 (cond
                   [(> min max) (list (reverse primes) (reverse comps))]
                   [(prime? min) (inner (+ min 1) (cons min primes) comps)]
                   [else (inner (+ min 1) primes (cons min comps))]))))
      (inner min '() '()))))

(define avg-rel-err
  (lambda (basis divs)
    (letrec ((step-size (lambda (div)
                          (/ (* (log2 (car basis)) 1200) div))))
      (mean (map / (map (lambda (x)
                          (best-val-err x basis top-damage))
                     divs)
                   (map (lambda (x)
                          (step-size x))
                     divs))))))

;; EXAMPLES

;; (avg-rel-err '(2 3 5 7) (cddr (n-primes 1002)))
;; 0.11286850458406111

;; (avg-rel-err '(2 3 5 7) (n-composites 1000))
;; 0.11777593026693876

;; (avg-rel-err '(2 3 5 7) (car (range-prime-comp 5 100))) ;; primes
;; 0.11412166292698245

;; (avg-rel-err '(2 3 5 7) (cadr (range-prime-comp 5 100))) ;; comps
;; 0.11763126363512029