;;;; misc-math-utils.scm  -*- Scheme -*-
;;;; Kon Lovett, Oct '21
;;;; some from utilities.lisp

;;;; Issues
;;
;; - Recast as functor
;;
;; - Use error correcting & pairwise - see fp-utils

(module misc-math-utils

(;export
  log-with-base log/base
  compound-interest simple-interest
  fibonacci fibonacci*
  binomial
  cross-ratio
  square
  cube
  average
  least-squares
  trapezoid
  big-pi big-sigma
  factorial factorial- factorial+
  ;*factorial
  harmonic)

(import scheme
  (chicken base)
  (chicken fixnum)
  (chicken type))

(define-type real (or integer ratnum float))

(define-type point (pair number number))
(define-type points (list-of point))

;FIXME cannot state purity w/ `-->' for argument or result type
(: log-with-base (real --> (real -> real)))
(: log/base (real --> (real -> real)))

(: simple-interest (real real #!optional real --> real))
(: compound-interest (real real real #!optional real --> real))

(: fibonacci* (real --> real))
(: fibonacci (integer --> integer))

(: binomial (integer integer --> integer))

(: cross-ratio (number number number number --> number))

(: square (number --> number))
(: cube (number --> number))

;close
(: average ((or (list-of number) number) #!rest number --> number))

(: least-squares (points --> number number))

(: trapezoid ((number -> number) number number -> (fixnum -> number)))

(: big-pi ((number -> number) number number --> number))
(: big-sigma ((number -> number) number number --> number))

(: factorial (integer --> integer))
#; ;diverges quickly
(: *factorial (integer --> integer))
(: factorial- (integer integer --> integer))
(: factorial+ (integer integer --> integer))

(: harmonic (integer --> number))

;;

(define-inline (reduce* f s l)
  ;(foldl (cut f <> <>) s l)
  (let loop ((l l) (s s) (len 0))
    (if (null? l) (values s len)
        (loop (cdr l) (f s (car l)) (fx+ 1 len)) ) ) )

(include-relative "mathh-constants")

(define-constant 1/SQRT5 (/ SQRT5))
(define-constant MID1+SQRT5 (/ (+ 1 SQRT5) 2))
(define-constant MID1-SQRT5 (/ (- 1 SQRT5) 2))

;from mathh.scm

(define (log-with-base b)
  (let ((lnb (log b)))
    (lambda (a) (/ (log a ) lnb)) ) )

(define log/base log-with-base)

;rate is per time period units, freq & perd must be in same time units
;
(define (simple-interest rate time #!optional (prin 1))
  (* prin (+ 1 (* rate time))) )

(define (compound-interest rate freq time #!optional (prin 1))
  (* prin (expt (+ 1 (/ rate freq)) (* freq time))) )

;approximate!
(define (fibonacci* n)
  (assert (real? n) 'fibonacci* "bad argument type - not a real" n)
  (cond ((zero? n)
          0)
        ;assumes successor operation is in the "direction" of the parameter
        ((negative? n)
          (- (fibonacci* (- n))) )
        (else
          (- (* 1/SQRT5 (expt MID1+SQRT5 n))
             (* 1/SQRT5 (expt MID1-SQRT5 n))) ) ) )

(define (fibonacci n)
  (assert (integer? n) 'fibonacci "bad argument type - not an integer" n)
  (cond ((zero? n)
          0)
        ;assumes successor operation is in the "direction" of the parameter
        ((negative? n)
          (- (fibonacci (- n))) )
        (else
          (let loop ((a 1) (b 1) (n n))
            (if (<= n 2) b
                (loop b (+ a b) (- n 1)) ) ) ) ) )

;; from utilities.lisp

(define (binomial n k)
  (cond ((zero? k)
          1)
        ((= n k)
          1)
        (else
          (+ (binomial (sub1 n) (sub1 k))
             (binomial (sub1 n) k)))) )

(define (cross-ratio a b c d)
  (/ (* (- a c) (- b d))
     (* (- a d) (- b c))) )

(define (square n) (* n n))

(define (cube n) (* n n n))

(define average
  (case-lambda
    ((ls) (let-values (((s l) (reduce* + 0 ls))) (/ s l)))
    ((a b) (/ (+ a b) 2))
    (r (average r)) ) )

;returns b e such that y ~= b * x + e
;FIXME take f s e i, f(n) where s <= n <= e & n' = n+1
(define (least-squares points)
  (do ((points points (cdr points))
       (sigma-a 0 (+ sigma-a (caar points)))
       (sigma-b 0 (+ sigma-b (cdar points)))
       (sigma-a-squared 0 (+ sigma-a-squared (square (caar points))))
       (sigma-xy 0 (+ sigma-xy (* (caar points) (cdar points))))
       (count 0 (fx+ count 1)) )
      ((null? points)
       (let* ((m (/ (- (* count sigma-xy) (* sigma-a sigma-b))
                    (- (* count sigma-a-squared) (* sigma-a sigma-a))))
              (b (/ (- sigma-b (* m sigma-a)) count)))
         (values m b))) ) )

;trapezoid rule for area. returns function taking "precision"
(define (trapezoid f a b)
  (let ((b-a (- b a)))
    (lambda (n)
      (let ((dx (/ b-a n)))
        (define (l i sum)
          (if (>= i n) (* sum dx)
              (let ((a (+ a (* i dx))))
                (l (add1 i) (+ sum (f a))))))
        (l 1 (/ (+ (f a) (f b)) 2)) ) ) ) )

(define-syntax coll-rng
  (syntax-rules ()
    ((coll-rng ?f ?i ?g ?s ?e)
      (let ((f ?f) (i ?i) (g ?g) (s ?s) (e ?e))
        (let loop ((k s) (r i))
          (if (= k e) r
              (loop (add1 k) (f r (g k))) ) ) ) ) ) )

(define-syntax big-pi
  (syntax-rules ()
    ((big-pi ?g ?s ?e)
      (let ((s ?s) (e ?e))
        (when (> s e) (error 'big-pi "bounds mismatch" s e))
        (coll-rng * 1 ?g s e)) ) ) )

(define-syntax big-sigma
  (syntax-rules ()
    ((big-sigma ?g ?s ?e)
      (let ((s ?s) (e ?e))
        (when (> s e) (error 'big-sigma "bounds mismatch" s e))
        (coll-rng + 0 ?g s e)) ) ) )

(define (factorial a) (big-pi identity 1 (add1 a)))

#; ;naive - w/o "large" float diverges at 18! (= 1..17 !)
(define (*factorial n) (inexact->exact (round (expt E (big-sigma log 1 (add1 n))))))

(define (factorial- a fall) (big-pi identity (- a (sub1 fall)) (add1 a)))

(define (factorial+ a rise) (big-pi identity a (+ a rise)))

(define (harmonic n) (big-sigma / 1 (+ n 1 1)))

) ;misc-math-utils