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

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

#>
// Assumes endiness of int64_t & double same!
typedef union int64dbl { int64_t i; double d; } punint64dbl;

// 1

static double fp_int_succ( double n )
{ punint64dbl x; x.i = (*(int64_t*)&n) + 1; return x.d; }

static double fp_int_pred( double n )
{ punint64dbl x; x.i = (*(int64_t*)&n) - 1; return x.d; }

// N

static double fp_int_next( double n, int64_t v )
{ punint64dbl x; x.i = (*(int64_t*)&n) + v; return x.d; }

static double fp_int_prev( double n, int64_t v )
{ punint64dbl x; x.i = (*(int64_t*)&n) - v; return x.d; }

//
#define C_rnd_fix() (C_fix(rand()))
<#

(module math-utils

(;export
  degree->radian radian->degree
  fppred fpsucc
  integer-log
  log-with-base log/base
  coprime? pairwise-coprime? fxcoprime?
  prime? primes
  compound-interest simple-interest
  cross-ratio
  chinum *chinum
  sqr square
  cub cube
  avg *average average average*
  least-squares
  trapezoid
  euclidian-distance manhattan-distance
  euclidian-distances manhattan-distances
  slope delta-slope
  @prec to-places
  ~=
  prod summ
  big-pi big-sigma
  binomial binomial/memo *binomial *binomial/memo
  fibonacci/approximate
  fibonacci fibonacci/memo *fibonacci *fibonacci/memo
  factorial/approximate
  factorial factorial/memo *factorial *factorial/memo
  factorial+ factorial+/memo *factorial+ *factorial+/memo
  factorial- factorial-/memo *factorial- *factorial-/memo
  harmonic harmonic/memo *harmonic *harmonic/memo
  fxrnd)

(import scheme)
(import (chicken base))
(import (chicken fixnum))
(import (chicken foreign))
(import (chicken syntax)) ;(check-errors sys)
(import (chicken type))
(import (only (srfi 1) reverse!))
(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))

(: degree->radian (real -> real))
(: radian->degree (real -> real))

(: fppred (float #!optional integer -> float))
(: fpsucc (float #!optional integer -> float))

;NOTE arguments swapped vs srfi-94
(: integer-log (integer integer -> integer (or integer ratnum)))

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

(: prime? (integer -> boolean))
(: primes (integer -> (-> (or integer eof))))

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

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

(: chinum (number -> fixnum))

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

;NOTE kinda violates the *Func/Func/Func* division
(: *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 ((or point number) #!optional (or point number) number number -> number))
(: manhattan-distance ((or point number) #!optional (or point number) number number -> number))

(: euclidian-distances (points -> (list-of number)))
(: manhattan-distances (points -> (list-of number)))

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

(: delta-slope  (number number -> number))
(: slope ((or point number) (or point number) #!optional number number -> number))

(: ~= (number number #!optional number -> boolean))

(: *binomial (integer integer -> integer))
(: binomial (integer integer -> integer))
(: *binomial/memo (integer integer -> integer))
(: binomial/memo (integer integer -> integer))

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

;diverges quickly
(: factorial/approximate (integer -> integer))
(: *factorial (integer -> integer))
(: factorial (integer -> integer))
(: *factorial/memo (integer -> integer))
(: factorial/memo (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))

(: fxrnd (#!optional fixnum -> fixnum))

;(check-errors sys)

(define-syntax check-number
  (syntax-rules ()
    ((check-number ?loc ?obj)
      (begin (##sys#check-number ?obj ?loc) ?obj)) ) )

(define-syntax check-integer
  (syntax-rules ()
    ((check-integer ?loc ?obj)
      (begin (##sys#check-integer ?obj ?loc) ?obj)) ) )

(define-syntax check-real
  (syntax-rules ()
    ((check-real ?loc ?obj)
      (begin (##sys#check-real ?obj ?loc) ?obj)) ) )

;;

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

#; ;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 E (exp 1))
(define SQRT5 (sqrt 5))

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

;;

;pi/180
(define DEG/RAD 5030569068109113/288230376151711744)

(define (degree->radian deg) (* deg DEG/RAD))
(define (radian->degree rad) (/ rad DEG/RAD))

(define fppred
  (case-lambda
    ((n)    ((foreign-lambda double fp_int_pred double) n))
    ((n v)  ((foreign-lambda double fp_int_prev double unsigned-integer64) n v)) ) )

(define fpsucc
  (case-lambda
    ((n)    ((foreign-lambda double fp_int_succ double) n))
    ((n v)  ((foreign-lambda double fp_int_next double unsigned-integer64) n v)) ) )

;Returns two values, the integer-log of <n> in <base>, and the leftmost digit
;in <base>.
(define (integer-log n base)
  (if (< n base)
      (values 0 n)
      (receive (ilog l) (integer-log n (* base base))
        (if (< l base)
            (values (* ilog 2) l)
            (values (+ (* ilog 2) 1) (/ l base))))) )

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

;

;1 ? - 0 is impl detail
(define PRIME-PRE #(0 1 2 3 5 7 11 13 17 23 27))
(define PRIME-PRE-LEN (vector-length PRIME-PRE))
(define PRIME-PRE-MAX (vector-ref PRIME-PRE (- PRIME-PRE-LEN 1)))

;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 (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
  (case-lambda
    ((p1 p2)
      (assume ((p1 point) (p2 point))
        (euclidian-distance (xy-x p1) (xy-y p1) (xy-x p2) (xy-y p2))) )
    ((x1 y1 x2 y2)
      (sqrt (+ (sqr (- x2 x1)) (sqr (- y2 y1)))) ) ) )

(define manhattan-distance
  (case-lambda
    ((p1 p2)
      (assume ((p1 point) (p2 point))
        (manhattan-distance (xy-x p1) (xy-y p1) (xy-x p2) (xy-y p2))) )
    ((x1 y1 x2 y2)
      (+ (abs (- x2 x1)) (abs (- y2 y1))) ) ) )

(define (distances distf ps)
  ;pair-fold
  (if (null? ps) '()
      (let loop ((ps (cdr ps)) (p1 (car ps)) (ds '()))
        (if (null? ps) (reverse! ds)
            (let ((p2 (car ps)))
              (loop (cdr ps) p2 (cons (distf p1 p2) ds)) ) ) ) ) )

(define (euclidian-distances ps) (distances euclidian-distance ps))
(define (manhattan-distances ps) (distances manhattan-distance ps))

(define (@prec p n #!optional (r truncate)) (let ((d (expt 10 p))) (/ (r (* n d)) d)))
(define-syntax-rule (to-places ?prec (?clip ?expr)) (@prec ?prec ?expr ?clip))

(define (delta-slope dx dy)
  (cond ((not (zero? dx)) (/ dy dx))
        ;vertical, horizontal
        ((zero? dy)       (* dx dy))  ;preserve sign! watch for -0.0
        ((negative? dy)   -inf.0)
        (else             +inf.0)) )

(define slope
  (case-lambda
    ((p0 p1)
      (slope (xy-x p0) (xy-y p0) (xy-x p1) (xy-y p1)) )
    ((x0 y0 x1 y1)
      (delta-slope (- x1 x0) (- y1 y0)) ) ) )

;test-support.scm
(define (~= a b #!optional (epsilon 1e-05))
  (cond ((> (abs a) (abs b))  (~= b a epsilon))
        ((zero? a)            (< (abs b) epsilon))
        (else                 (< (abs (/ (- a b) b)) epsilon))) )

;;

;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 #!optional (d 1))
  ;(assert (<= s e) 'range-last-index "bounds mismatch" s e d)
  (let loop ((k s))
    (if (<= e k) (- k d)
        (loop (+ k d))) ) )
(define (range-last-index s e #!optional (d 1))
  ;(assert (<= s e) 'range-last-index "bounds mismatch" s e d)
  ;FIXME (?) s = e => s - d
  (- (+ s (ceiling (/ (- e s) d))) d) )

(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

;lisppad
(define prime?
  (memoize
    (lambda (n)
      (cond ((<= n 1)  #f)
            ((= n 2)   #t)
            ((even? n) #f)
            (else      (do ((d 3 (+ d 2)))
                           ((or (> (square d) n)
                                (zero? (remainder n d)))
                                (> (square d) n)))) ) ) ) )

;1st N primes
(define (primes n)
  (let ((+n+ (the (or integer eof) n))
        (+i+ 2) )
    (lambda ()
      (if (eof-object? +n+) #!eof
          (let loop ((ni +i+))
            (cond ((prime? ni)
                    (set! +n+ (if (= 1 +n+) #!eof (sub1 +n+)))
                    (set! +i+ (add1 ni))
                    ni )
                  (else (loop (add1 ni))) ) ) ) ) ) )

; Binomial

(define (check-binomial-args loc n k)
  (unless (<= 0 (check-integer loc k)) (error loc "must be natural" k))
  (unless (<= k (check-integer loc n)) (error loc "must be greater" n k))
  #;(values n k) )

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

(define (binomial n k)
  (check-binomial-args 'binomial n k)
  (*binomial n k) )

(define *binomial/memo
  (memoize
    (lambda (n k)
      (cond ((zero? k)  1)
            ((= n k)    1)
            (else       (+ (*binomial/memo (- n 1) (- k 1))
                           (*binomial/memo (- n 1) k)))) ) ) )

(define (binomial/memo n k)
  (check-binomial-args 'binomial/memo n k)
  (*binomial/memo n k) )

; 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)
      (if (< n 2) 1
          (* 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)) )

(define **factorial-/memo
  (memoize
    (lambda (s k)
      (if (< k s) 1
          (* k (**factorial-/memo s (- k 1)))) ) ) )

(define **factorial+/memo
  (memoize
    (lambda (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)) )

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

; Random

(define (fxrnd #!optional (lim most-positive-fixnum))
  (fxmod (##core#inline "C_rnd_fix") lim) )

) ;math-utils