Changeset 36954 in project


Ignore:
Timestamp:
12/03/18 00:35:33 (8 days ago)
Author:
kon
Message:

fix fprandom, fix fp~=

Location:
release/5/fp-utils/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • release/5/fp-utils/trunk/fp-utils.scm

    r36953 r36954  
    119119;;
    120120
    121 (: fprandom (#!optional (or float fixnum) -> float))
     121(: fprandom (#!optional (or float fixnum) (or float fixnum) -> float))
    122122;
    123 (define (fprandom #!optional (lim most-positive-fixnum) (low 0.0))
    124   (let* (
    125     (low (inexact->exact (fptruncate low)))
    126     (lim
    127       (if (not (flonum? lim))
    128         lim
    129         (let* (
    130           (neg? (fpnegative? lim))
    131           (lim (inexact->exact (fpexpt 10.0 (fpround (fplog10 (fpabs lim))))))
    132           (lim (fxmax most-positive-fixnum lim)) )
    133           (if neg? (fxneg lim) lim)))) )
    134     (if (fx>= low lim)
    135       +nan.0
    136       (fp/ 1.0 (exact->inexact (fxrandom lim low))) ) ) )
     123(define-inline (*fpinv x) (if (fpzero? x) 0.0 (fp/ 1.0 x)))
     124(define-inline (*fpinvfx x) (inexact->exact (*fpinv x)))
     125(define-inline (*fxinvfp x) (*fpinv (exact->inexact x)))
     126(define (fprandom #!optional (lim most-positive-fixnum) (low 0))
     127  (cond
     128    ((and (fixnum? lim) (fixnum? low))
     129      (if (fx>= low lim)
     130        (error 'fprandom "range swapped" lim low)
     131        ;invert maps fx to fp result
     132        (*fxinvfp (fxrandom lim low)) ) )
     133    ((and (flonum? lim) (flonum? low))
     134      (if (fp>= low lim)
     135        (error 'fprandom "range swapped" lim low)
     136        ;invert & swap maps fp to fx parameters
     137        ;result will be in [low .. lim]
     138        (fprandom (*fpinvfx low) (*fpinvfx lim)) ) )
     139    (else
     140      (error 'fprandom "no mixed-mode" lim low) ) ) )
    137141
    138142;;
     
    141145;
    142146(define (fp~= x y #!optional (eps flonum-epsilon))
    143   #; ;enough 4 how i used it but ...
    144   (fp<= (fpabs (fp- x y)) eps)
    145   ;better
     147  ;NOTE minimum/maximum-flonum is smallest/largest positive normal flonum
    146148  (or
    147149    (fp= x y)
    148     (let (
    149       (abs-x (fpabs x))
    150       (abs-y (fpabs y))
    151       (abs-d (fpabs (fp- x y))) )
    152       ;very small relative error?
    153       ;(minimum/maximum-flonum is smallest/largest positive normal flonum)
    154       (if (or (fpzero? abs-x) (fpzero? abs-y) (fp< abs-d minimum-flonum))
    155         (fp< abs-d (fp* eps minimum-flonum))
    156         (fp< (fp/ abs-d (fpmin (fp+ abs-x abs-y) maximum-flonum)) eps) ) ) ) )
     150    (let* (
     151      (d (fp- x y))
     152      (abs-d (fpabs d)) )
     153      (cond
     154        ;either 0 argument or << relative error?
     155        ((or (fpzero? x) (fpzero? y) (fp< abs-d minimum-flonum))
     156          (fp< abs-d (fp* eps minimum-flonum)) )
     157        ;ensure |x| <= |y| so no /0
     158        ((fp> (fpabs x) (fpabs y))
     159          (fp~= y x eps) )
     160        (else
     161          (fp< (fpabs (fp/ d y)) eps) ) ) ) ) )
    157162
    158163(: fp~<= (float float #!optional float --> boolean))
  • release/5/fp-utils/trunk/tests/fp-utils-test.scm

    r36953 r36954  
    4747
    4848(test-group "Fp Approximate ="
     49
     50  ;sweet, sweet repeating expansion
     51  (let ((x (fp+ 0.15 0.15)) (y (fp+ 0.1 0.2)))
     52          (test-assert ".15+.15 <> .1+.2" (not (fp= x y)))
     53          (test-assert ".15+.15 ~= .1+.2" (fp~= x y)) )
     54
     55  ;own epsilon
    4956        (test-assert (fp~= 0.123456 0.123457 5eps))
    5057        (test-assert (fp~<= 0.123456 0.123457 5eps))
     
    7784)
    7885
     86#;
     87(begin
     88  (import (chicken flonum) fp-inlines fp-utils)
     89  (do ((i 0 (add1 i)) )
     90      ((>= i 100000) )
     91    (let ((n (fprandom 0.05 0.002)))
     92      (print n)
     93      (assert (<= 0.002 n))
     94      (assert (<= n 0.05)) ) ) )
     95
    7996;;;
    8097
Note: See TracChangeset for help on using the changeset viewer.