;;;; fp-utils.scm
;;;; Kon Lovett, May '17

;;;; Issues
;;;;

(module fp-utils

(;export
  ;
  fprandom
  ;
  fpzero? fppositive? fpcardinal? fpnegative? fpeven? fpodd?
  ;
  fpclosed-right? fpclosed? fpclosed-left? fpclosedr? fpclosedl?
  ;
  fpadd1 fpsub1
  ;
  fpmodulo
  fpquotient fpremainder
  ;
  fpfraction
  ;
  fptruncate/precision fpround/precision fpceiling/precision fpfloor/precision
  ;
  fp~= fp~<= fp~>=
  ;
  fpsqr fpcub
  ;
  fpdegree->radian fpradian->degree
  ;
  fpdistance fpdistance*
  ;
  fpmax-and-min
  ;
  fpprecision-factor)

(import scheme)

(import chicken foreign extras)

(declare
  (bound-to-procedure
    ##sys#flonum-fraction
    ##sys#check-inexact) )

;;;

(define C_fmod
  (foreign-lambda double "fmod" double double))

(define C_remainder
  (foreign-lambda double "remainder" double double))

(define (*fpeven? n)
  (fp= 0.0 (##sys#flonum-fraction (fp/ n 2.0))) )

;;;

;;

(define (fprandom #!optional (x most-positive-fixnum))
  (fp/ 1.0 (exact->inexact (random x))) )

;;

(define (fpzero? n)
  (fp= 0.0 n) )

(define (fppositive? n)
  (fp< 0.0 n) )

(define (fpcardinal? n)
  (fp<= 0.0 n) )

(define (fpnegative? n)
  (fp> 0.0 n) )

(define (fpeven? n)
  (and
    (fpinteger? n)
    (*fpeven? n)) )

(define (fpodd? n)
  (and
    (fpinteger? n)
    (not (*fpeven? n))) )

;;

(define (fpclosed-right? l x h)
  (and (fp< l x) (fp<= x h)) )

(define (fpclosed? l x h)
  (and (fp<= l x) (fp<= x h)) )

(define (fpclosed-left? l x h)
  (and (fp<= l x) (fp< x h)) )

(define fpclosedr? fpclosed-right?)
(define fpclosedl? fpclosed-left?)

;;

(define (fpadd1 n)
  (fp+ n 1.0) )

(define (fpsub1 n)
  (fp- n 1.0) )

;;

(define (fpmodulo x y)
  (##sys#check-inexact x 'fpmodulo)
  (##sys#check-inexact y 'fpmodulo)
  (fptruncate (C_fmod x y)) )

(define (fpquotient x y)
  (fptruncate (fp/ x y)) )

(define (fpremainder x y)
  (##sys#check-inexact x 'fpremainder)
  (##sys#check-inexact y 'fpremainder)
  (fptruncate (C_remainder x y)) )

;;

(define (fpfraction n)
	(##sys#flonum-fraction n) )

;;;

;;

(define (fp~= x y #!optional (eps flonum-epsilon))
  (let ((diff (fp- x y)))
    (or
      (fpzero? diff)
      (fp<= (fpabs diff) eps) ) ) )

(define (fp~<= x y #!optional (eps flonum-epsilon))
  (or
    (fp< x y)
    (fp~= x y eps) ) )

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

;;;

;;

(define (fpsqr n)
  (fp* n n) )

(define (fpcub n)
  (fp* n (fp* n n)) )

;;;

(define-constant PRECISION-DEFAULT 4.0)

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

;;

(define fptruncate/precision (make-unary/precision fptruncate))
(define fpround/precision (make-unary/precision fpround))
(define fpceiling/precision (make-unary/precision fpceiling))
(define fpfloor/precision (make-unary/precision fpfloor))

;;

(define-constant DEGREE 0.0174532925199432957692369076848861271344) ;pi/180

(define (fpdegree->radian deg)
  (fp* deg DEGREE) )

(define (fpradian->degree rad)
  (fp/ rad DEGREE) )

;;

(define (fpdistance x1 y1 x2 y2)
  (fpsqrt (fpdistance* x1 y1 x2 y2)) )

(define (fpdistance* x1 y1 x2 y2)
  (fp+ (fpsqr (fp- x1 x2)) (fpsqr (fp- y1 y2))) )

;;

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

;;

(define (fpprecision-factor p #!optional (base 10.0))
  (fpexpt base (exact->inexact p)) )

) ;fp-utils