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

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

(module math-utils

(;export
  log-with-base log/base
  coprime? pairwise-coprime? fxcoprime?
  compound-interest simple-interest
  fibonacci fibonacci*
  binomial
  cross-ratio
  square
  cube
  average
  least-squares
  trapezoid
  big-pi big-sigma
  factorial factorial- factorial+
  ;*factorial
  harmonic
  euclidian-distance
  manhattan-distance
  with-places
  ;DEPRECATED
  number->sign)

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

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

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

(: number->sign (deprecated signum))
(define number->sign signum)

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

(: fxcoprime? (fixnum fixnum --> boolean))
(: coprime? (integer #!rest integer --> boolean))
(: pairwise-coprime? (integer #!rest integer --> boolean))

(: 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)) ) ) )

;mathh-constants
(define-constant SQRT5      2.2360679774997896964091736687312762354406)   ; sqrt(5)

(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)

#; ;1/2 the speed of `coprime?'
(define (fxcoprime? a b)
  (or (fx= a b)
      (if (fx< a b)
          (fxcoprime? b a)
          (or (and (fx= 1 b) (or (fx= 2 a) (fx= 3 a)))
              (let ((2b (fx+ b b)))
                (if (and (fx< b a) (fx< a 2b))
                    (fxcoprime? b (fx- 2b a))
                    (let ((3b (fx+ b 2b)))
                      (cond ((fx> a 3b)
                              (fxcoprime? (fx- a 2b) b))
                            ((and (fx< 2b a) (fx< a 3b))
                              (fxcoprime? b (fx- a 2b)))
                            (else
                              (and (not (fx= 0 b))
                                   (not (fx= 0 (fxrem a b)))) ) ) ) ) ) ) ) ) )
(define (fxcoprime? a b)
  (fx= 1 (fxgcd a b)) )

(define-inline (any func . args)
  (let loop ((args args))
    (and (not (null? args))
         (or (func (car args))
             (loop (cdr args)))) ) )

(define (andmap func . args)
  (or (apply any null? args)
      (and (apply func (map car args))
           (apply andmap func (map cdr args)))) )

;math egg
(define (coprime? a . bs)
  (= 1 (apply gcd (cons a bs))) )

;math egg
(define (pairwise-coprime? a . bs)
  (or (null? bs)
      (and (andmap (cut coprime? a <>) bs)
           (apply pairwise-coprime? bs))) )

;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)))

(define (euclidian-distance x1 y1 x2 y2)
  (sqrt (+ (square (- x2 x1)) (square (- y2 y1)))) )

(define (manhattan-distance x1 y1 x2 y2)
  (+ (abs (- x2 x1)) (abs (- y2 y1))) )

(define-syntax with-places
  (syntax-rules ()
    ((with-places ?prec (?clip ?expr))
      (let ((f (expt 10 ?prec))) (/ (?clip (* ?expr f)) f)) ) ) )

) ;math-utils