;;;; fp-utils.scm  -*- Scheme -*-
;;;; Kon Lovett, Sep '18
;;;; Kon Lovett, May '17

;;;; Issues
;;
;; - all instances of (fl< -0.0 0.0) found ?

(module fp-utils

(;export
  ;
  fpsummation
  ;
  fprandom
  ;
  fpmodulo fpquotient fpremainder
  ;
  fptruncate-with-precision
  fpround-with-precision
  fpceiling-with-precision
  fpfloor-with-precision
  ;
  fpdistance fpdistance*
  ;
  fpmax-and-min
  ;
  fp~= fp~<= fp~>=)

(import scheme)
(import (chicken base))
(import (chicken foreign))
(import (chicken type))
(import (chicken flonum))
(import (chicken fixnum))
(import fx-utils)
(import fp-inlines)

;;;

;;

(: fpsum ((list-of float) --> float))
;
(define (fpsummation ls)
  (if (null? ls)
    0.0
    ;Kahan Summation Formula
    (let sum ((s (car ls)) (c 0.0) (ls (cdr ls)))
      (if (null? ls)
        s
        (let* (
          (y (fp- (car ls) c))
          (t (fp+ s y)) )
          (sum t (fp- (fp- t s) y) (cdr ls)) ) ) ) ) )

;;;

(: C_fmod (float float --> float))
;
(define C_fmod (foreign-lambda double "fmod" double double))

(: C_remainder (float float --> float))
;
(define C_remainder (foreign-lambda double "remainder" double double))

;;

(: fpmodulo (float float --> float))
;
(define (fpmodulo x y) (fptruncate (C_fmod x y)))

(: fpquotient (float float --> float))
;
(define (fpquotient x y) (fptruncate (fp/ x y)))

(: fpremainder (float float --> float))
;
(define (fpremainder x y) (fptruncate (C_remainder x y)))

;;;

;(so no conversion)
(define-constant PRECISION-DEFAULT 4.0)

(define-syntax make-unary-with-precision
  (syntax-rules ()
    ((make-unary-with-precision ?op)
      (lambda (n #!optional (p PRECISION-DEFAULT))
        (if (fpzero? p)
          (?op n)
          (let ((pf (fpprecision-factor p)))
            (fp/ (?op (fp* n pf)) pf) ) ) ) ) ) )

;;

(: fptruncate-with-precision (float #!optional float --> float))
;
(define fptruncate-with-precision (make-unary-with-precision fptruncate))

(: fpround-with-precision (float #!optional float --> float))
;
(define fpround-with-precision (make-unary-with-precision fpround))

(: fpceiling-with-precision (float #!optional float --> float))
;
(define fpceiling-with-precision (make-unary-with-precision fpceiling))

(: fpfloor-with-precision (float #!optional float --> float))
;
(define fpfloor-with-precision (make-unary-with-precision fpfloor))

;;

(: fpdistance (float float float float --> float))
;
(define (fpdistance x1 y1 x2 y2) (fpsqrt (fpdistance* x1 y1 x2 y2)))

(: fpdistance* (float float float float --> float))
;
(define (fpdistance* x1 y1 x2 y2) (fp+ (fpsqr (fp- x1 x2)) (fpsqr (fp- y1 y2))))

;;

(: fpquo-and-rem (float float --> float float))
;
(define (fpquo-and-rem fpn fpd) (values (fpquotient fpn fpd) (fpremainder fpn fpd)))

;;

(: fpmax-and-min (float #!rest float --> float float))
;
(define (fpmax-and-min fp . fps)
  (let loop ((fps fps) (mx fp) (mn fp))
    (if (null? fps)
      (values mx mn)
      (let ((cur (car fps)))
        (loop (cdr fps) (fpmax mx cur) (fpmin mn cur)) ) ) ) )

;;

(: fprandom (#!optional (or float fixnum) (or float fixnum) -> float))
;
(define-inline (*fpinv x) (if (fpzero? x) 0.0 (fp/ 1.0 x)))
(define-inline (*fpinvfx x) (inexact->exact (*fpinv x)))
(define-inline (*fxinvfp x) (*fpinv (exact->inexact x)))
(define (fprandom #!optional (lim most-positive-fixnum) (low 0))
  (cond
    ((and (fixnum? lim) (fixnum? low))
      (if (fx>= low lim)
        (error 'fprandom "range swapped" lim low)
        ;invert maps fx to fp result
        (*fxinvfp (fxrandom lim low)) ) )
    ((and (flonum? lim) (flonum? low))
      (if (fp>= low lim)
        (error 'fprandom "range swapped" lim low)
        ;invert & swap maps fp to fx parameters
        ;result will be in [low .. lim]
        (fprandom (*fpinvfx low) (*fpinvfx lim)) ) )
    (else
      (error 'fprandom "no mixed-mode" lim low) ) ) )

;;

(: fp~= (float float #!optional float --> boolean))
;
(define (fp~= x y #!optional (eps flonum-epsilon))
  ;NOTE minimum/maximum-flonum is smallest/largest positive normal flonum
  (or
    (fp= x y)
    (let* (
      (d (fp- x y))
      (abs-d (fpabs d)) )
      (cond
        ;either 0 argument or << relative error?
        ((or (fpzero? x) (fpzero? y) (fp< abs-d minimum-flonum))
          (fp< abs-d (fp* eps minimum-flonum)) )
        ;ensure |x| <= |y| so no /0
        ((fp> (fpabs x) (fpabs y))
          (fp~= y x eps) )
        (else
          (fp< (fpabs (fp/ d y)) eps) ) ) ) ) )

(: fp~<= (float float #!optional float --> boolean))
;
(define (fp~<= x y #!optional (eps flonum-epsilon)) (or (fp< x y) (fp~= x y eps)))

(: fp~>= (float float #!optional float --> boolean))
;
(define (fp~>= x y #!optional (eps flonum-epsilon)) (or (fp> x y) (fp~= x y eps)))

) ;fp-utils