;;;; math-utils-test.scm  -*- Scheme -*-
;;;; Kon Lovett, Oct '22

(import (test-utils gloss) (only (chicken format) format))

(import test)

(test-begin "Math Utils")

;;;

(import math-utils)
(import (chicken fixnum) (chicken random))

(cond-expand
  (compiling
    (glossln)
    (gloss "*****")
    (gloss "* Compiler warnings are expected")
    (gloss "*****") )
  (else) )

;;

(define ((least-squares-error m e) pt)
  (- (cdr pt) (+ (* m (car pt)) e)) )

(define (average-least-squares-error m e pts)
  (apply average (map (least-squares-error m e) pts)) )

(define (fibonacci/naive n)
  (case n
    ((0) 0)
    ((1) 1)
    (else (+ (fibonacci/naive (- n 1)) (fibonacci/naive (- n 2)))) ) )

#; ;UNUSED
(define (approx-points b e pts)
  (map (lambda (c)
         (cons (exact->inexact (car c))
               (exact->inexact (round (+ (* b (car c)) e)))))
       pts))

;;

(define some-pts '((0 . 0) (1 . 2) (2 . 4) (3 . 8) (4 . 16)))
(define data-pts '((1.47 . 52.21) (1.50 . 53.12) (1.52 . 54.48)
                   (1.55 . 55.84) (1.57 . 57.20) (1.60 . 58.57)
                   (1.63 . 59.93) (1.65 . 61.29) (1.68 . 63.11)
                   (1.70 . 64.47) (1.73 . 66.28) (1.75 . 68.10)
                   (1.78 . 69.92) (1.80 . 72.19) (1.83 . 74.46)))

(test-group "trapezoid"
   (define area-3ln (trapezoid (lambda (x) (* 3 (log x))) 3 4))
   (define data-pts/exact
           (map (lambda (x) (cons (inexact->exact (car x)) (inexact->exact (cdr x))))
                data-pts))

   (test-assert (procedure? area-3ln))
   (test 3.7478 (area-3ln 10))
   (test 3.7480 (area-3ln 10000))
)

(test-group "average"
  (test-error (*average 23 '(4 5 4 5)))
  (test 9/2 (*average '(4 5 4 5)))

  (test-error (average '(4 5 4 5)))
  (test 1 (average 1))
  (test 2 (average 1 3))
  (test 3 (average 1 3 5))
  (test 9/2 (average 4 5))
  (test 9/2 (average 4 5 4 5))

  (test 183/7 (average* 69 '((1 3 5) (4 5)) 96))
)

(test-group "binomial"
  (test 70 (binomial 8 4))
)

(test-group "Least Squares (error 0)"
  (receive (m e) (least-squares some-pts)
    (test 0 (average-least-squares-error m e some-pts)) )
  (receive (m e) (least-squares data-pts)
    (test 0.0 (average-least-squares-error m e data-pts)) )
)

;https://en.wikipedia.org/wiki/Cross-ratio
(test-group "Cross Ratio (identity)"
  (let* ((a 1) (b 2) (c 3) (d 4)
         (l (cross-ratio a b c d)))
    (test 4/3 l)
    (test (/ 1 l)       (cross-ratio a b d c))
    (test (/ 1 (- 1 l)) (cross-ratio a c d b))
    (test (- 1 l)       (cross-ratio a c b d))
    (test (/ l (- l 1)) (cross-ratio a d c b))
    (test (/ (- l 1) l) (cross-ratio a d b c)) )
)

(test-group "interest"
  (test 1758.0 (simple-interest 0.043 4 1500))
  (test 1887.0 (simple-interest 0.043 6 1500))

  (test 1780.42 (compound-interest 0.043 6 4 1500))
  (test 1938.84 (compound-interest 0.043 4 6 1500))
)

(test-group "log-with-base"
  (define log256 (log-with-base 256))
  (test-assert (procedure? log256))
  (test "log/base 256" 9 (inexact->exact (round (log256 995972255318008407568))))
)

(test-group "coprime"
  (test-assert (coprime? 2 1))
  (test-assert (coprime? 3 1))
  (test-assert (coprime? 81 19))
  (test-assert (coprime? 1 1))
  (test-assert (not (coprime? 9 3)))
  (test-assert (coprime? (* 3 7) (* 5 19)))
  (test-assert (not (coprime? (* 3 7 5) (* 5 19 2))))
  (test-assert (coprime? 6 10 15))
)

(test-group "pairwise coprime"
  (test-assert (pairwise-coprime? 10 7 33 13))
  (test-assert (not (pairwise-coprime? 10 7 33 14)))
  (test-assert (not (pairwise-coprime? 6 10 15)))
)

(test-group "*coprime"
  (test-assert (fxcoprime? 2 1))
  (test-assert (fxcoprime? 3 1))
  (test-assert (fxcoprime? 81 19))
  (test-assert (fxcoprime? 1 1))
  (test-assert (not (fxcoprime? 9 3)))
  (test-assert (fxcoprime? (* 3 7) (* 5 19)))
  (test-assert (not (fxcoprime? (* 3 7 5) (* 5 19 2))))
)

(test-group "to-places"
  (test 3.141 (to-places 3 (truncate (acos -1))))
  (test 3.0 (to-places 0 (truncate (acos -1))))
  (test 3.14159 (to-places 5 (truncate (acos -1))))
)

;; Series

(test-group "fibonacci"
  (test 55 (fibonacci/naive 10))

  (test-group "approx fib"
    (test-error (fibonacci/approximate 'a))
    (test (fibonacci/naive 10) (inexact->exact (floor (fibonacci/approximate 10))))
  )

  (test-group "true fib"
    (test-error (*fibonacci 'a))
    (test-error (fibonacci 'a))
    (test (fibonacci/naive 10) (*fibonacci 10))
    (test (fibonacci/naive 10) (fibonacci 10))
    (test-error (*fibonacci/memo 'a))
    (test-error (fibonacci/memo 'a))
    (test (fibonacci/naive 10) (*fibonacci/memo 10))
    (test (fibonacci/naive 10) (fibonacci/memo 10))
  )
)

(test-group "factorial"
  (test-error (factorial 'a))
  (test 1 (factorial 0))
  (test 1 (factorial 1))
  (test 2 (factorial 2))
  (test 6 (factorial 3))
  (test 720 (factorial 6))

  (test-error (*factorial/memo 'a))
  (test-error (factorial/memo 'a))
  (do ((i 0 (add1 i)))
      ((< 6 i))
    (test (factorial i) (*factorial/memo i))
    (test (factorial i) (factorial/memo i)) )
)

(test-group "factorialą"
  (define (fxrand #!optional (rng most-positive-fixnum)) (pseudo-random-integer rng))

  (define *tr* 5 #;(fxrand))

  (test-error (*factorial- 'a 4))
  (test-error (*factorial+ 'a 4))
  (test-error (factorial- 'a 4))
  (test-error (factorial+ 'a 4))

  ;https://en.wikipedia.org/wiki/Falling_and_rising_factorials
  (test 1 (*factorial- *tr* 0))
  (test 1 (*factorial+ *tr* 0))
  (test *tr* (*factorial- *tr* 1))
  (test *tr* (*factorial+ *tr* 1))
  (test (- (square *tr*) *tr*) (factorial- *tr* 2))
  (test (+ (square *tr*) *tr*) (factorial+ *tr* 2))
  ;quit while ahead
  (test (+ (- (expt *tr* 4) (* 6 (cube *tr*))) (- (* 11 (square *tr*)) (* 6 *tr*)))
        (*factorial- *tr* 4))
  (test (+ (expt *tr* 4) (* 6 (cube *tr*)) (* 11 (square *tr*)) (* 6 *tr*))
        (*factorial+ *tr* 4))

  (test-error (*factorial-/memo 'a 4))
  (test-error (*factorial+/memo 'a 4))
  (test-error (factorial-/memo 'a 4))
  (test-error (factorial+/memo 'a 4))
  (do ((i 0 (+ i 1)))
      ((< 5 i))
    (gloss 'factorial- (- (*factorial- *tr* i) (*factorial-/memo *tr* i)))
    ;#;
    (test (*factorial- *tr* i) (*factorial-/memo *tr* i))
    (gloss 'factorial+ (- (*factorial+ *tr* i) (*factorial+/memo *tr* i)))
    ;#;
    (test (*factorial+ *tr* i) (*factorial+/memo *tr* i)) )
)

(test-group "harmonic"
  ;lush
  (define (harmonic/baseline n)
    (let loop ((i 1) (z 0))
      (if (< n i) z
          (loop (+ i 1) (+ z (/ i)))) ) )

  (test-error (harmonic 'q))
  (test 25/12 (harmonic/baseline 4))
  (test (harmonic/baseline 4) (harmonic 4))
  (test-error (*harmonic/memo 'q))
  (test (*harmonic 4) (*harmonic/memo 4))
  (test-error (harmonic/memo 'q))
  (test (harmonic 4) (harmonic/memo 4))
)

;;;

(test-end "Math Utils")

(test-exit)