math-r7rs/math.body.scm

;; Jason K. MacDuffie
;; I dedicate this to the public domain, for the benefit of all
;; due to the triviality of the library.

(define pi 3.141592653589793)
(define e 2.718281828459045)
(define i 0+i)
(define (deg->rad d)
  (/ (* d pi) 180))
(define (rad->deg r)
  (/ (* r 180) pi))
(define (dsin d)
  (sin (deg->rad d)))
(define (dcos d)
  (cos (deg->rad d)))
(define (dtan d)
  (tan (deg->rad d)))
(define (dasin a)
  (rad->deg (asin a)))
(define (dacos a)
  (rad->deg (acos a)))
(define (datan a)
  (rad->deg (atan a)))

(define (integrate f minval maxval dx)
  ;; simple integration
  (let loop ((x minval)
             (sum 0.0))
    (if (<= x maxval)
        (loop (+ x dx)
              (+ sum (* dx (f x))))
        sum)))

(define (factorial x)
  ;; Technically this is a continuation of factorial
  ;; But it is convenient sometimes
  (if (and (exact? x) (integer? x) (not (negative? x)))
      (let loop ((i 1)
                 (p 1))
        (if (<= i x)
            (loop (+ i 1) (* p i))
            p))
      (gamma (+ x 1))))

(define (stirling-gamma z)
  (* (sqrt (* 2 pi (/ z)))
     (expt (* (/ e) (+ z (/ (- (* 12 z) (/ 1/10 z)))))
           z)))

(define (gamma z)
  (define a (- (real-part z) (floor (real-part z))))
  (define b (imag-part z))
  (cond
   ((and (integer? z)
         (positive? z))
    (let loop ((i 1)
               (p 1.0))
      (if (< i z)
          (loop (+ i 1) (* i p))
          p)))
   ((< (real-part z) 150)
    (let loop ((i z)
               (p (stirling-gamma (+ 150 a (* 0+i b)))))
      (if (< (real-part i) 150)
          (loop (+ i 1)
                (/ p i))
          p)))
   ((< (real-part z) 151)
    (stirling-gamma z))
   (else
    (let loop ((i z)
               (p (stirling-gamma (+ 150 a (* 0+i b)))))
      (if (> (real-part i) 151)
          (loop (- i 1)
                (* p (- i 1)))
          p)))))

(define (inverse-gamma x)
  (if (<= x 0.8856032)
      (error "inverse-gamma" "x must be greater than 0.8856032"))
  (if (>= x 2e300)
      (error "inverse-gamma" "x must be less than 2.0e300"))
  (let ((upper-bound (let loop ((u 1.46164))
                       (if (< (gamma u) x)
                           (loop (+ u 1))
                           u))))
    (let loop ((frac 0.5)
               (prev 0.0)
               (r (- upper-bound 0.5)))
      (if (not (= prev r))
          (loop (* frac 0.5)
                r
                (if (< (gamma r) x)
                    (+ r frac)
                    (- r frac)))
          (if (= x (gamma (round r)))
              (round r)
              r)))))

(define (inverse-factorial x)
  (let ((y (- (inverse-gamma x) 1)))
    (if (and (exact? x) (integer? y))
        (exact y)
        y)))

;; Inverse of the function (* n (log n))
(define (anti-nlogn x)
  (define (nlogn n)
    (* n (log n)))
  (let ((stopat
         (cond
          ((< x 2.75)
           (error "value must be at least 2.75"))
          ((< x 10.0)
           2000)
          (else
           100))))
    (let loop ((guess x)
               (i 0))
      (if (or (= i stopat)
              (= (nlogn guess) x))
          guess
          (loop (/ x (log guess))
                (+ i 1))))))

;; Some statistics functions
(define (mean l)
  (/ (apply + l) (length l)))

(define (variance l . sample-option)
  (define sample? (and (not (null? sample-option))
                       (eq? (car sample-option) 'sample)))
  (define a (square (apply + l)))
  (define b (apply + (map square l)))
  (define n (length l))
  (define ss (- b (/ a n)))
  (/ ss (if sample? (- n 1) n)))

(define (deviation l . sample-option)
  (define sample? (and (not (null? sample-option))
                       (eq? (car sample-option) 'sample)))
  (sqrt (apply variance
               (list l (if sample? 'sample 'population)))))

(define (z-scores l)
  (define sd (deviation l))
  (define m (mean l))
  (map (lambda (x)
         (/ (- x m) sd))
       l))

(define (normal-pmf mu sigsq)
  (define C (/ (sqrt (* 2 sigsq pi))))
  (lambda (x)
    (* C (exp (- (/ (square (- x mu))
                    (* 2 sigsq)))))))

(define (normal-cdf mu sigsq)
  ;; For simplicity, make internal mu 0
  (define f (normal-pmf 0 sigsq))
  (define (simple-cdf x)
    (if (< x -10)
        0.0
        (integrate f (* -4 sigsq) x (* 5.0e-5 sigsq))))
  (lambda (x)
    (if (< x mu)
        (simple-cdf (- x mu))
        (- 1 (simple-cdf (- mu x))))))

(define (primes upto-minus-one)
  ;; Sieve of Eratosthenes
  (cond
   ((> upto-minus-one 10000000)
    (error "primes"
      "Cannot calculate for values greater than 10000000"
      upto-minus-one))
   ((< upto-minus-one 2)
    '())
   (else
    (let ((upto (+ upto-minus-one 1)))
      (define stop-at (exact (floor (sqrt upto))))
      (define sieve-vector (make-vector upto #t))
      (define (believed-prime? x)
        (vector-ref sieve-vector x))
      (vector-set! sieve-vector 0 #f)
      (vector-set! sieve-vector 1 #f)
      (let loop1 ((i 2))
        (if (<= i stop-at)
            (if (believed-prime? i)
                (let loop2 ((j (* i 2)))
                  (if (< j upto)
                      (begin
                        (vector-set! sieve-vector j #f)
                        (loop2 (+ j i)))
                      (loop1 (+ i 1))))
                (loop1 (+ i 1)))
            ;; now return a list of primes
            (let make-primes ((out '())
                              (j (- upto 1)))
              (if (>= j 0)
                  (if (believed-prime? j)
                      (make-primes (cons j out) (- j 1))
                      (make-primes out (- j 1)))
                  out))))))))

(define (prime? n)
  (define upto
    (floor (exact (sqrt n))))
  (if (< n 2)
      #f
      (let loop ((next (primes upto)))
        (cond
         ((null? next)
          #t)
         ((= 0 (modulo n (car next)))
          #f)
         (else
          (loop (cdr next)))))))

(define (approx-equal a b eps)
  (< (abs (- (/ a b) 1)) eps))

(define (~= x y . rest)
  (and (approx-equal x y 1e-13)
       (or (null? rest)
           (apply ~= x rest))))

(define (percent-scale l)
  (let ((total (apply + l)))
    (map (lambda (i) (* 100 (/ i total))) l)))