;;;; 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
  binomial
  cross-ratio
  chinum *chinum
  sqr square
  cub cube
  avg *average average average*
  least-squares
  trapezoid
  euclidian-distance
  manhattan-distance
  @prec to-places
  prod summ
  big-pi big-sigma
  fibonacci/approximate
  fibonacci *fibonacci
  fibonacci/memo *fibonacci/memo
  ;factorial/approximate
  factorial *factorial
  factorial/memo *factorial/memo
  factorial- factorial+ *factorial- *factorial+
  factorial-/memo factorial+/memo *factorial-/memo *factorial+/memo
  harmonic *harmonic
  harmonic/memo *harmonic/memo)

(import scheme)
(import (chicken base))
(import (chicken fixnum))
(import (chicken type))
(import (only miscmacros define-syntax-rule))
(import memoize)

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

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

(: binomial (integer integer -> integer))

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

(: chinum (number -> fixnum))

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

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

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

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

(: euclidian-distance (number number number number -> number))
(: manhattan-distance (number number number number -> number))

(: @prec (fixnum number #!optional (number -> number) -> number))

(: fibonacci/approximate (integer -> real))
(: *fibonacci (integer -> integer))
(: fibonacci (integer -> integer))
(: *fibonacci/memo (integer -> integer))
(: fibonacci/memo (integer -> integer))

(: *factorial (integer -> integer))
(: factorial (integer -> integer))
(: *factorial/memo (integer -> integer))
(: factorial/memo (integer -> integer))
#; ;diverges quickly
(: factorial/approximate (integer -> integer))

;real (not integer) due to use as Pochhammer function
(: *factorial- (real integer -> real))
(: *factorial+ (real integer -> real))
(: factorial- (real integer -> real))
(: factorial+ (real integer -> real))
(: *factorial-/memo (real integer -> real))
(: *factorial+/memo (real integer -> real))
(: factorial-/memo (real integer -> real))
(: factorial+/memo (real integer -> real))

(: *harmonic (integer -> real))
(: harmonic (integer -> real))
(: *harmonic/memo (integer -> real))
(: harmonic/memo (integer -> real))

;;

(define-inline (xy x y) (cons x y))
(define-inline (xy-x p) (car p))
(define-inline (xy-y p) (cdr p))

(define-inline (check-integer loc obj)
  (assert (integer? obj) loc "bad argument type - not an integer" obj)
  obj )

(define-inline (check-real loc obj)
  (assert (real? obj) loc "bad argument type - not an real" obj)
  obj )

(define-inline (check-number loc obj)
  (assert (number? obj) loc "bad argument type - not an number" obj)
  obj )

#; ;FIXME -inline vs -syntax-rule
(define-syntax-rule (check-integer-range loc s e)
  (when (> (check-integer loc s) (check-integer loc e))
    (error loc "bounds mismatch" s e)) )

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

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

;;

#;(define-constant E (exp 1))
(define-constant SQRT5 (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?', d'uh
(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 (coprime? a . bs) (= 1 (apply gcd a bs))) ;math egg

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

(define (binomial n k)
  (cond ((zero? (check-integer 'binomial k))
          1)
        ((= (check-integer 'binomial n) k)
          1)
        (else
          (+ (binomial (- n 1) (- k 1))
             (binomial (- n 1) k)))) )

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

;FIXME handle complex (a'la `signum'?)
(define-syntax-rule (*chinum n) (if (even? n) 1 -1)) ;(expt -1 n)
(define (chinum n) (*chinum n))

(define-syntax-rule (sqr x) (* x x))
(define-syntax-rule (cub x) (* x x x))

(define (square n) (sqr n))
(define (cube n) (cub n))

(define-syntax-rule (avg ?n0 ...) (/ (+ ?n0 ...) (length ?n0 ...)))

(define (*average ls)
  (receive (s l) (reduce* + 0 ls)
    (if (fx= 0 l) 0
        (/ s l)) ) )

(define average
  (case-lambda
    (()     0)
    ((a b)  (/ (+ a b) 2))
    (r      (*average r)) ) )

(define (average* . nums) (*average (flatten nums)))

;returns slope+intercept m b such that y ~= m * x + b
;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-x 0 (+ sigma-x (caar points))) ;(point-x (car points))
       (sigma-y 0 (+ sigma-y (cdar points))) ;(point-y (car points))
       (sigma-x-squared 0 (+ sigma-x-squared (sqr (caar points))))
       (sigma-xy 0 (+ sigma-xy (* (caar points) (cdar points))))
       (count 0 (fx+ count 1)) )
      ((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))) ) )

;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 (fx>= i n) (* sum dx)
              (let ((a (+ a (* i dx))))
                (l (fx+ i 1) (+ sum (f a))))))
        (l 1 (/ (+ (f a) (f b)) 2)) ) ) ) )

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

(define (@prec p n #!optional (r truncate))
  (let ((f (expt 10 p))) (/ (r (* n f)) f)) )

(define-syntax-rule (to-places ?prec (?clip ?expr)) (@prec ?prec ?expr ?clip))

;;

;NOTE Limit (end) is exclusive.

(: range-last-index (real real -> real))
(: **factorial-/memo (real real -> real))
(: **factorial+/memo (real real -> real))

;(define (range-last-index s e) (- (+ s (- e s)) 1))
;(define (range-last-index s e) (- (+ s (floor (- e s))) 1))
;(define (range-last-index s e) (+ s (- (floor e) (floor s))))
;(define (range-last-index s e) (+ s (floor (- e s))))
(define (range-last-index s e)
  ;(assert (< s e))
  (let loop ((k s))
    (if (<= e k) #;k (- k 1)
        (loop (+ k 1))) ) )

(define-syntax-rule (range-collect ?f ?v ?s ?d ?e ?i ?arg1 ...)
  (let ((f ?f) (e ?e) (d ?d))
    (let loop ((?v ?s) (a ?i))
      (if (<= e ?v) a
          (loop (+ ?v d) (f a ?arg1 ...))) ) ) )

; rep - : (summ x (...) (- ...))
; rep / : (prod x (...) (/ ...))

(define-syntax prod
  (syntax-rules ()
    ((prod ?v () ?arg1 ...)         (prod ?v (1) ?arg1 ...))
    ((prod ?v (?s) ?arg1 ...)       (prod ?v (?s +inf.0) ?arg1 ...))
    ((prod ?v (?s ?e) ?arg1 ...)    (prod ?v (?s ?e 1) ?arg1 ...))
    ((prod ?v (?s ?e ?i) ?arg1 ...) (range-collect * ?v ?s ?i ?e 1 ?arg1 ...) ) ) )

(define-syntax summ
  (syntax-rules ()
    ((summ ?v () ?arg1 ...)         (summ ?v (1) ?arg1 ...))
    ((summ ?v (?s) ?arg1 ...)       (summ ?v (?s +inf.0) ?arg1 ...))
    ((summ ?v (?s ?e) ?arg1 ...)    (summ ?v (?s ?e 1) ?arg1 ...))
    ((summ ?v (?s ?e ?i) ?arg1 ...) (range-collect + ?v ?s ?i ?e 0 ?arg1 ...) ) ) )

(define-syntax coll-rng
  (syntax-rules (identity)
    ;short-circuit common pattern
    ((coll-rng ?f ?i identity ?s ?e)
      (range-collect ?f k ?s 1 ?e ?i k) )
    ;has next-value generator (or unrecognized "identity")
    ((coll-rng ?f ?i ?g ?s ?e)
      (let ((g ?g))
        (range-collect ?f k ?s 1 ?e ?i (g k)) ) ) ) )

(define-syntax-rule (big-pi ?g ?s ?e)     (coll-rng * 1 ?g ?s ?e))
(define-syntax-rule (big-sigma ?g ?s ?e)  (coll-rng + 0 ?g ?s ?e))

;; Series

; Fibonacci

;approximate! error (f* -f) grows > 1 @ n=71
(define (fibonacci/approximate n)
  (cond ((zero? (check-integer 'fibonacci/approximate n))
          0)
        ;assumes successor operation is in the "direction" of the parameter
        ((negative? n)
          (- (fibonacci/approximate (- n))) )
        (else
          (- (* 1/SQRT5 (expt MID1+SQRT5 n))
             (* 1/SQRT5 (expt MID1-SQRT5 n))) ) ) )

(define (*fibonacci 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)) ) ) ) ) )

(define (fibonacci n) (*fibonacci (check-integer 'fibonacci n)))

(define *fibonacci/memo
  (memoize
    (lambda (n)
      (cond ((zero? n)  0)
            ((= 1 n)    1)
            (else       (+ (*fibonacci/memo (- n 1))
                           (*fibonacci/memo (- n 2)))) ) ) ) )

(define (fibonacci/memo n) (*fibonacci/memo (check-integer 'fibonacci/memo n)))

; Factorial

#; ;w/o "large" float, or rational expt, diverges at 18! (= 1..17 !),
(define (factorial/approximate n)
  (check-integer 'factorial/approximate n)
  (inexact->exact (round (expt E (big-sigma log 1 (+ n 1))))) )

(define (*factorial n) (big-pi identity 1 (+ n 1)))
(define (factorial n) (*factorial (check-integer 'factorial n)))

(define *factorial/memo
  (memoize
    (lambda (n)
      (cond ((or (zero? n) (= 1 n))  1)
            (else                    (* n (*factorial/memo (- n 1)))) ) ) ) )

(define (factorial/memo n) (*factorial/memo (check-integer 'factorial/memo n)))

; Factorial- Factorial+

(define (*factorial- a fall) (big-pi identity (- a (- fall 1)) (+ a 1)))    ;order!
(define (*factorial+ a rise) (big-pi identity a                (+ a rise)))

(define (factorial- a fall)
  (*factorial- (check-real 'factorial- a) (check-integer 'factorial- fall)) )

(define (factorial+ a rise)
  (*factorial+ (check-real 'factorial+ a) (check-integer 'factorial+ rise)) )

(cond-expand
  (factorial±/memo
(define factorial-/memo factorial-)
(define factorial+/memo factorial+)
(define *factorial-/memo *factorial-)
(define *factorial+/memo *factorial+)
  )
  (else
(define **factorial-/memo
  (memoize
    (lambda (s k)
(print "**-! " s " " k)
      (if (<= k s) 1
          (* k (**factorial-/memo s (- k 1)))) ) ) )

(define **factorial+/memo
  (memoize
    (lambda (s k)
(print "**+! " s " " k)
      (if (<= k s) 1
          (* k (**factorial+/memo s (- k 1)))) ) ) )

(define (*factorial-/memo a fall)
  (let ((s (- a (- fall 1)))) ;order!
    (**factorial-/memo s (range-last-index s (+ a 1))) ) )

(define (*factorial+/memo a rise)
  (**factorial+/memo a (range-last-index a (+ a rise))) )

(define (factorial-/memo a fall)
  (*factorial-/memo (check-real 'factorial-/memo a) (check-real 'factorial-/memo fall)) )

(define (factorial+/memo a rise)
  (*factorial+/memo (check-real 'factorial+/memo a) (check-real 'factorial+/memo rise)) )
) ) ;cond-expand factorial±/memo

; Harmonic

; 1/1 + 1/2 + ... + 1/n + 1/n+1
(define (*harmonic n) (big-sigma / 1 (+ n 1)))
(define (harmonic n) (*harmonic (check-integer 'harmonic n)))

(define **harmonic/memo
  (memoize
    (lambda (n)
      (let ((n (- n 1)))
        (if (zero? n) 0
            (+ (/ n) (**harmonic/memo n))) ) ) ) )

(define (*harmonic/memo n) (**harmonic/memo (+ n 1)))
(define (harmonic/memo n) (*harmonic/memo (check-integer 'harmonic/memo n)))

) ;math-utils