;;;; misc-math-utils.scm  -*- Scheme -*-
;;;; Kon Lovett, Oct '21

;;;; 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*
  binomial
  cross-ratio
  square
  cube
  average
  least-squares
  trapezoid
  big-pi big-sigma
  factorial factorial- factorial+
  ;*factorial
  harmonic)

(import scheme
  (chicken base)
  (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))

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

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

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

(: average (number number --> number))

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

(: trapezoid ((number -> number) number number -> (integer -> 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 (number --> number))

;(exact (sequence-ref "https://oeis.org/A002163"))
(define-constant SQRT5 629397181890197/281474976710656)
(define-constant 1/SQRT5 (/ SQRT5))
(define-constant MID1+SQRT5 (/ (+ 1 SQRT5) 2))
(define-constant MID1-SQRT5 (/ (- 1 SQRT5) 2))

(define-constant E 2.7182818284590452353602874713526624977572)   ; e

;from mathh.scm

(define (log-with-base b)
  (let ((lnb (log b)))
    (lambda (x) (/ (log x ) 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)
  (cond ((zero? n)
          0)
        ((negative? n)
          (fibonacci* (- (/ n))))
        (else
          (floor (- (* 1/SQRT5 (expt MID1+SQRT5 n))
                    (* 1/SQRT5 (expt MID1-SQRT5 n))))) ) )

;;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
    ((a b) (/ (+ a b) 2))
    ((ls) (apply average ls))
    (r (/ (apply + r) (length r)))))

(define (least-squares points)
  (do ((points points (cdr points))
       (sigma-x 0 (+ sigma-x (caar points)))
       (sigma-y 0 (+ sigma-y (cdar points)))
       (sigma-x-squared 0 (+ sigma-x-squared (square (caar points))))
       (sigma-xy 0 (+ sigma-xy (* (caar points) (cdar points))))
       (count 0 (add1 count)) )
      ((null? points)
       (let* ((m (/ (- (* count sigma-xy) (* sigma-x sigma-y))
                    (- (* count sigma-x-squared) (* sigma-x sigma-x))))
              (b (/ (- sigma-y (* m sigma-x)) count)))
         (values m b))) ) )

(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 ((x (+ a (* i dx))))
                (l (add1 i) (+ sum (f x))))))
        (l 1 (/ (+ (f a) (f b)) 2)) ) ) ) )

(define-syntax coll-rng
  (syntax-rules ()
    ((coll-rng ?loc ?c ?i ?f ?s ?e)
      (let ((loc ?loc) (c ?c) (i ?i) (f ?f) (s ?s) (e ?e))
        (when (> s e) (error loc "bounds mismatch" s e))
        (let loop ((k s) (r i))
          (if (= k e) r
              (loop (add1 k) (c r (f k))) ) ) ) ) ) )

(define-syntax big-pi
  (syntax-rules ()
    ((big-pi ?f ?s ?e)
      (coll-rng 'big-pi * 1 ?f ?s ?e)) ) )

(define-syntax big-sigma
  (syntax-rules ()
    ((big-sigma ?f ?s ?e)
      (coll-rng 'big-sigma + 0 ?f ?s ?e)) ) )

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

#; ;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- x fall) (big-pi identity (- x (sub1 fall)) (add1 x)))

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

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

) ;misc-math-utils