Merge pull request #309 from uwplse/ival-find-range

Ival find range

Author: pavpanchekha

bench/libraries/jmatjs.fpcore

@@ -2,6 +2,7 @@
 
 (FPCore (z)
  :name "Jmat.Real.gamma, branch z greater than 0.5"
+ :pre (> z 0.5)
  (let ((z* (- z 1)) (g 7))
    (let ((x
           (+
@@ -19,6 +20,7 @@
 
 (FPCore (z)
  :name "Jmat.Real.gamma, branch z less than 0.5"
+ :pre (<= z 0.5)
  (let ((z* (- (- 1 z) 1)) (g 7))
    (let ((x
           (+
@@ -67,6 +69,7 @@
 
 (FPCore (x)
  :name "Jmat.Real.erfi, branch x less than or equal to 0.5"
+ :pre (<= x 0.5)
  (let ((sqrtPI (sqrt PI)))
    (let ((ps (/ 1 sqrtPI)))
      (let ((x* (fabs x)))
@@ -78,6 +81,7 @@
 
 (FPCore (x)
  :name "Jmat.Real.erfi, branch x greater than or equal to 5"
+ :pre (>= x 0.5)
  (let ((sqrtPI (sqrt PI)))
    (let ((ps (/ 1 sqrtPI)))
      (let ((x* (fabs x)))

bench/libraries/mathjs/arithmetic.fpcore

@@ -116,6 +116,7 @@
 
 (FPCore (re im)
  :name "math.sqrt on complex, imaginary part, im greater than 0 branch"
+ :pre (> im 0)
  (* 0.5 (sqrt (* 2.0 (- (sqrt (+ (* re re) (* im im))) re)))))
 
 (FPCore (re im)

bench/regression.fpcore

@@ -56,3 +56,13 @@
 (FPCore (x)
  :name "sqrt E"
  (sqrt (+ (pow x 2) (pow x 2))))
+
+; This used to cause crashes
+
+(FPCore (w l)
+ :name "exp-w crasher"
+ (* (exp (- w)) (pow l (exp w))))
+
+(FPCore (x)
+ :name "expfmod"
+ (* (fmod (exp x) (sqrt (cos x))) (exp (- x))))

src/common.rkt

@@ -8,7 +8,7 @@
          for/append string-prefix call-with-output-files
          take-up-to flip-lists list/true find-duplicates all-partitions
          argmins argmaxs index-of set-disjoint? comparator sample-double
-         random-ranges parse-flag get-seed set-seed!
+         parse-flag get-seed set-seed!
          quasisyntax
          format-time format-bits in-sorted-dict web-resource
          (all-from-out "config.rkt") (all-from-out "debug.rkt"))
@@ -136,30 +136,6 @@
 
 ;; Miscellaneous helper
 
-(define (random-ranges . ranges)
-  (->* () #:rest (cons/c integer? integer?) integer?)
-
-  (define weights
-    (for/list ([(lo hi) (in-dict ranges)])
-      ;; The `max` handles the case lo > hi and similar
-      (max 0 (- hi lo))))
-
-  (define total-weight (apply + weights))
-  (when (= total-weight 0)
-    (error 'random-ranges "Empty ranges"))
-
-  (define num-bits (inexact->exact (ceiling (/ (log total-weight) (log 2)))))
-  (define sample ; Rejection sampling
-    (let loop ()
-      (define sample (if (= num-bits 0) 0 (random-bits num-bits)))
-      (if (< sample total-weight) sample (loop))))
-
-  (let loop ([sample sample] [ranges ranges] [weights weights])
-    ;; The `(car)` is guaranteed to succeed by the construction of `sample`
-    (if (< sample (car weights))
-        (+ (caar ranges) sample)
-        (loop (- sample (car weights)) (cdr ranges) (cdr weights)))))
-
 (define (parse-flag s)
   (match (string-split s ":")
     [(list (app string->symbol category) (app string->symbol flag))

src/config.rkt

@@ -5,14 +5,14 @@
 
 (define all-flags
   #hash([precision . (double fallback)]
-        [setup . (simplify early-exit)]
+        [setup . (simplify early-exit search)]
         [generate . (rr taylor simplify better-rr)]
         [reduce . (regimes avg-error binary-search branch-expressions)]
         [rules . (arithmetic polynomials fractions exponents trigonometry hyperbolic numerics special bools branches)]))
 
 (define default-flags
   #hash([precision . (double fallback)]
-        [setup . (simplify)]
+        [setup . (simplify search)]
         [generate . (rr taylor simplify)]
         [reduce . (regimes avg-error binary-search branch-expressions)]
         [rules . (arithmetic polynomials fractions exponents trigonometry hyperbolic special bools branches)]))
@@ -49,6 +49,9 @@
 ;; Number of iterations of the core loop for improving program accuracy
 (define *num-iterations* (make-parameter 4))
 
+;; The maximum depth for splitting the space when searching for valid areas of points.
+(define *max-find-range-depth* (make-parameter 14))
+
 ;; The maximum number of consecutive skipped points for sampling valid points
 (define *max-skipped-points* (make-parameter 100))
 

src/core/periodicity.rkt

@@ -187,8 +187,9 @@
 			     (let ([context
 				    (prepare-points
 				     program
-                                     `(and ,@(for/list ([(var period) (lp-periods ploc)])
-                                               `(<= 0 ,var ,(* 2 pi var))))
+                                     `(λ ,(program-variables program)
+                                        (and ,@(for/list ([(var period) (lp-periods ploc)])
+                                                 `(<= 0 ,var ,(* 2 pi var)))))
                                      (*output-repr*))])
 			       (parameterize ([*pcontext* context])
 				 (improve-func (make-alt program)))))))

src/core/regimes.rkt

@@ -1,7 +1,7 @@
 #lang racket
 
 (require "../common.rkt" "../alternative.rkt" "../programs.rkt" "../timeline.rkt")
-(require "../syntax/types.rkt" "../interface.rkt")
+(require "../syntax/types.rkt" "../interface.rkt" "../errors.rkt")
 (require "../points.rkt" "../float.rkt") ; For binary search
 
 (module+ test
@@ -86,12 +86,12 @@
                          (critical-subexpression? prog-body expr)))
     expr))
 
-(define (combine-alts best-option repr)
+(define (combine-alts best-option repr sampler)
   (match-define (option splitindices alts pts expr _) best-option)
   (match splitindices
    [(list (si cidx _)) (list-ref alts cidx)]
    [_
-    (define splitpoints (sindices->spoints pts expr alts splitindices repr))
+    (define splitpoints (sindices->spoints pts expr alts splitindices repr sampler))
     (debug #:from 'regimes "Found splitpoints:" splitpoints ", with alts" alts)
 
     (define expr*
@@ -191,7 +191,7 @@
 ;; float form always come from the range [f(idx1), f(idx2)). If the
 ;; float form of a split is f(idx2), or entirely outside that range,
 ;; problems may arise.
-(define (sindices->spoints points expr alts sindices repr)
+(define (sindices->spoints points expr alts sindices repr sampler)
   (define eval-expr
     (eval-prog `(λ ,(program-variables (alt-program (car alts))) ,expr) 'fl repr))
 
@@ -207,7 +207,10 @@
                      [*timeline-disabled* true]
                      [*var-reprs* (dict-set (*var-reprs*) var repr)])
         (define ctx
-          (prepare-points start-prog `(== ,(caadr start-prog) ,v) repr))
+          (prepare-points start-prog
+                          `(λ ,(program-variables start-prog) (== ,(caadr start-prog) ,v))
+                          repr
+                          (λ () (cons v (sampler)))))
         (< (errors-score (errors prog1 ctx repr))
            (errors-score (errors prog2 ctx repr)))))
     (define pt (binary-search-floats pred v1 v2 repr))
@@ -223,20 +226,28 @@
 
     (sp (si-cidx sidx) expr (find-split prog1 prog2 p1 p2)))
 
+  (define (regimes-sidx->spoint sidx)
+    (sp (si-cidx sidx) expr (apply eval-expr (list-ref points (- (si-pidx sidx) 1)))))
+
   (define final-sp (sp (si-cidx (last sindices)) expr +nan.0))
 
+  (define use-binary
+    (and (flag-set? 'reduce 'binary-search)
+         ;; Binary search is only valid if we correctly extracted the branch expression
+         (andmap identity (cons start-prog progs))))
+  (if use-binary
+      (debug #:from 'binary-search "Improving bounds with binary search for" expr "and" alts)
+      (debug #:from 'binary-search "Only using regimes for bounds on" expr "and" alts))
+      
   (append
-   (if (and (flag-set? 'reduce 'binary-search)
-            ;; Binary search is only valid if we correctly extracted the branch expression
-            (andmap identity (cons start-prog progs)))
-       (begin
-         (debug #:from 'binary-search "Improving bounds with binary search for" expr "and" alts)
-         (for/list ([si1 sindices] [si2 (cdr sindices)])
-           (sidx->spoint si1 si2)))
-       (begin
-         (debug #:from 'binary-search "Only using regimes for bounds on" expr "and" alts)
-         (for/list ([sindex (take sindices (sub1 (length sindices)))])
-	   (sp (si-cidx sindex) expr (apply eval-expr (list-ref points (- (si-pidx sindex) 1)))))))
+   (for/list ([si1 sindices] [si2 (cdr sindices)])
+     (cond
+       [use-binary
+        (with-handlers ([exn:fail:user:herbie:sampling?
+                         (lambda (e) (regimes-sidx->spoint si1))])
+          (sidx->spoint si1 si2))]
+       [else
+        (regimes-sidx->spoint si1)]))
    (list final-sp)))
 
 (define (point-with-dim index point val)

src/errors.rkt

@@ -1,16 +1,21 @@
 #lang racket
 (require "config.rkt")
 (provide raise-herbie-error raise-herbie-syntax-error
+         raise-herbie-sampling-error
          herbie-error->string herbie-error-url
          (struct-out exn:fail:user:herbie)
          (struct-out exn:fail:user:herbie:syntax)
+         (struct-out exn:fail:user:herbie:sampling)
          warn warning-log expect-warning)
 
 (struct exn:fail:user:herbie exn:fail:user (url)
         #:extra-constructor-name make-exn:fail:user:herbie)
 
 (struct exn:fail:user:herbie:syntax exn:fail:user:herbie (locations)
-        #:extra-constructor-name make-exn:fail:user:herbie:syntax)
+  #:extra-constructor-name make-exn:fail:user:herbie:syntax)
+
+(struct exn:fail:user:herbie:sampling exn:fail:user (url)
+  #:extra-constructor-name make-exn:fail:user:herbie:sampling)
 
 (define (raise-herbie-error message #:url [url #f] . args)
   (raise (make-exn:fail:user:herbie
@@ -20,6 +25,10 @@
   (raise (make-exn:fail:user:herbie:syntax
           (apply format message args) (current-continuation-marks) url locations)))
 
+(define (raise-herbie-sampling-error message #:url [url #f] . args)
+  (raise (make-exn:fail:user:herbie:sampling
+          (apply format message args) (current-continuation-marks) url)))
+
 (define (herbie-error-url exn)
   (format "https://herbie.uwplse.org/doc/~a/~a"
           *herbie-version* (exn:fail:user:herbie-url exn)))

src/float.rkt

@@ -9,6 +9,7 @@
 (provide 
  get-representation*
  ordinary-value?
+ largest-ordinary-value bound-ordinary-values
  ulp-difference ulps->bits
  midpoint random-generate
  </total <=/total =-or-nan?
@@ -55,6 +56,20 @@
            (ordinary-value? (imag-part x) (get-representation 'binary64)))
       (not (special-value? x repr))))
 
+
+(define (largest-ordinary-value repr)
+  ((representation-repr->bf repr)
+   ((representation-ordinal->repr repr)
+    (- ((representation-repr->ordinal repr) ((representation-bf->repr repr) +inf.bf)) 1))))
+
+(define (bound-ordinary-values repr)
+  (parameterize ([bf-rounding-mode 'nearest])
+    (let loop ([ordinal (bigfloat->ordinal (largest-ordinary-value repr))] [stepsize 1])
+      (define bfval (ordinal->bigfloat ordinal))
+      (if (bfinfinite? ((representation-repr->bf repr) ((representation-bf->repr repr) bfval)))
+          bfval
+          (loop (+ ordinal stepsize) (* stepsize 2))))))
+
 (module+ test
   (define binary64 (get-representation 'binary64))
   (check-true (ordinary-value? 2.5 binary64))

src/interface.rkt

@@ -38,11 +38,19 @@
   1
   null)
 
+(define (shift bits fn)
+  (define shift-val (expt 2 bits))
+  (λ (x) (fn (- x shift-val))))
+
+(define (unshift bits fn)
+  (define shift-val (expt 2 bits))
+  (λ (x) (+ (fn x) shift-val)))
+
 (define-representation (binary64 real)
   bigfloat->flonum
   bf
-  ordinal->flonum
-  flonum->ordinal
+  (shift 63 ordinal->flonum)
+  (unshift 63 flonum->ordinal)
   64
   '(+nan.0 +inf.0 -inf.0))
 
@@ -62,8 +70,23 @@
     [(< x 0) (- (bit-field->single-flonum (- x)))]
     [else (bit-field->single-flonum x)]))
 
+(define (single-flonum-step x n)
+  (ordinal->single-flonum (+ (single-flonum->ordinal x) n)))
+
+(define (bigfloat->single-flonum x)
+  (define loprec (parameterize ([bf-precision 24]) (bf+ 0.bf x)))
+  (define y (real->single-flonum (bigfloat->flonum loprec)))
+  (define x2 (bf y))
+  (match (bf-rounding-mode)
+    ['nearest y]
+    ['up   (if (bf< x2 x) (single-flonum-step y 1) y)]
+    ['down (if (bf> x2 x) (single-flonum-step y -1) y)]
+    ['zero (if (bf< x 0.bf)
+               (if (bf< x2 x) (single-flonum-step y 1) y)
+               (if (bf> x2 x) (single-flonum-step y -1) y))]))
+
 (define-representation (binary32 real)
-  (compose real->single-flonum bigfloat->flonum)
+  bigfloat->single-flonum
   bf
   ordinal->single-flonum
   single-flonum->ordinal