Cammy

Fork

(module fp (nat? fp?
            build-fp fp->number
            fp-sign° fp-pos° fp-neg°
            fp-<° fp-<=°
            fp-zero+ fp-zero- fp-inf+ fp-inf-
            fp-+° fp-*°)
  (import scheme (chicken base)
          mini-kanren
          (only matchable match-lambda)
          (only mathh frexp)
          rels
          (chicken pretty-print))

  ; An implementation of relational floating-point arithmetic, as seen in:
  ; https://www.cs.toronto.edu/~lczhang/sandre_float2021.pdf
  ; Zeroes and infinities are given their own distinct encodings. Otherwise, we
  ; follow the paper relatively closely.

  ; Like (not (var? x)) but checks entire lists.
  (define nat? (match-lambda [(or (0 . n) (1 . n)) (nat? n)] [() #t] [_ #f]))

  ; Whether a value is a floating-point number.
  ; In particular, this will not ratify logic variables that are not ground.
  ; Float -> Bool
  (define fp?
    (match-lambda
      [((or 'pos 'neg) (or 'zero 'inf)) #t]
      [((or 'pos 'neg) (? nat?) (? nat?)) #t]
      [_ #f]))

  ; Overall precision.
  (define exponent-length 3)
  (define mantissa-length 4)
  (define mantissa-factor (expt 2.0 mantissa-length))
  (define exponent-bias (- (expt 2 (- exponent-length 1)) 1))
  (define exponent-max (- (expt 2 exponent-length) 1))
  
  ; Convert an inexact number into a relational term.
  ; float -> Float
  (define build-fp
    (match-lambda
      [+inf.0 '(pos inf)] [-inf.0 '(neg inf)]
      [+0.0 '(pos zero)] [-0.0 '(neg zero)]
      ; XXX doesn't handle denorms
      ; XXX doesn't overflow to infinity
      ; XXX rounding?
      [f (let-values (((m e) (frexp f)) ((s) (if (<= 0.0 f) 'pos 'neg)))
           (list s
             ; NB: frexp's image is in [0.5, 1.0) but we are working with [1.0, 2.0)
             ;     So, shift it over by one, and remove it from the exponent
             (build-num (+ -1 e exponent-bias))
             (build-num (abs (inexact->exact (round (* m mantissa-factor)))))))]))

  ; Convert a ground floating-point number into a Scheme inexact number.
  ; Float -> float
  (define fp->number
    (match-lambda
      ['(pos inf) +inf.0] ['(neg inf) -inf.0]
      ['(pos zero) 0.0] ['(neg zero) -0.0]
      [(s e m) (let
                 ((sign (if (eqv? s 'neg) -1.0 1.0))
                  ; NB: shift back for frexp/native range
                  (ex (expt 2.0 (- (little-endian->number e) exponent-bias -1)))
                  (mant (/ (little-endian->number m) mantissa-factor)))
                 (* sign ex mant))]
      [_ +nan.0]))

  ; The sign of a floating-point number, either 'pos or 'neg.
  (define (fp-sign° f s) (caro f s))
  (define (fp-pos° f) (fp-sign° f 'pos))
  (define (fp-neg° f) (fp-sign° f 'neg))

  ; The involutive negation of a floating-point number.
  (define (fp-negate° f1 f2)
    (conde
      ((fp-pos° f1) (fp-neg° f2))
      ((fp-neg° f1) (fp-pos° f2))))

  ; Signed zeroes.
  (define fp-zero+ '(pos zero))
  (define fp-zero- '(neg zero))
  (define (fp-zero° f) (conde ((== f fp-zero+)) ((== f fp-zero-))))

  ; Signed infinities.
  (define fp-inf+ '(pos inf))
  (define fp-inf- '(neg inf))
  (define (fp-inf° f) (conde ((== f fp-inf+)) ((== f fp-inf-))))

  ; Natural numbers which are exponents.
  (define (exponent° e)
    (fresh (n) (<=o n (build-num exponent-length)) (lengtho e n)))

  ; Natural numbers which are mantissas.
  (define (mantissa° m) (lengtho m (build-num mantissa-length)))
  (define (mantissa-short° m) (lengtho m (build-num (- mantissa-length 1))))

  ; Decompose nonzero finite floating-point numbers.
  (define (fp-decomp° f s e m)
    (fresh () (== f `(,s ,e ,m))
           (conde ((== s 'pos)) ((== s 'neg))) (exponent° e) (mantissa° m)))

  (define (fp-finite° f)
    (conde
      ((fp-zero° f))
      ((fresh (s e m) (fp-decomp° f s e m)))))

  ; The less-than relation.
  ; See p3-4.
  (define (fp-<° x y)
    (conde
      ((== x fp-inf-) (fp-finite° y))
      ((== y fp-inf+) (fp-finite° x))
      ((== x fp-zero+) (fp-pos° y))
      ((fp-neg° x) (== y fp-zero-))
      ((fp-neg° x) (fp-pos° y))
      ((fresh (e1 e2 m1 m2) (== x `(pos e1 m1)) (== y `(pos e2 m2))
              (conde ((<o e1 e2)) ((== e1 e2) (<o m1 m2)))))
      ((fresh (e1 e2 m1 m2) (== x `(neg e1 m1)) (== y `(neg e2 m2))
              (conde ((<o e1 e2)) ((== e1 e2) (<o m1 m2)))))))

  ; The less-than-or-equal relation.
  ; See p4.
  (define (fp-<=° x y) (conde ((== x y)) ((fp-<° x y))))

  ; Shift a mantissa.
  (define (mantissa-shift° mantissa diff shifted)
    (fresh (prefix)
           (appendo prefix shifted mantissa)
           (lengtho prefix diff)))

  ; Trim a mantissa, saving the extra bits.
  (define (drop-leastsig-bit° m rm bits)
    (fresh () (mantissa° rm) (appendo bits rm m)))

  ; Overflow check.
  ; Since we can't confuse infinite and finite numbers, we don't need to involve the mantissa.
  (define (fp-overflow° exponent) (<o (build-num exponent-max) exponent))

  ; Helper for addition.
  ; Requires that e1 <= e2.
  ; See p4-5.
  ; NB: Because of how we encode infinities, we pass a possibly-infinite return
  ; value rather than a return exponent and mantissa.
  (define (fp-samesignadder° e1 m1 e2 m2 r)
    (fresh (re rm expo-diff shifted-m1 msum bits)
           (conde
             ((fp-inf° r) (fp-overflow° re))
             ((fresh (s) (fp-decomp° r s re rm))))
           ; NB: This precondition is always required, so we compute it here
           ; instead of in callers.
           (+o expo-diff e1 e2)
           (mantissa-shift° m1 expo-diff shifted-m1)
           (conde
             ((nullo bits) (== e2 re))
             ((fresh (x) (caro bits x) (succ° e2 re))))
           (drop-leastsig-bit° msum rm bits)
           (+o shifted-m1 m2 msum)))

  ; Helper for addition at infinities.
  (define (fp-infadd° inf f1 f2 r)
    (fresh () (== r inf)
           (conde
             ((== f1 inf) (fp-finite° f2))
             ((fp-finite° f1) (== f2 inf))
             ((== f1 inf) (== f2 inf)))))

  ; Addition.
  (define (fp-+° f1 f2 r)
    (conde
      ((== f2 r) (fp-zero° f1))
      ((== f1 r) (fp-zero° f2))
      ; XXX more forgiving than IEEE 754! Allows ∞ - ∞ = 0
      ((fp-zero° r) (fp-negate° f1 f2))
      ((fp-infadd° fp-inf+ f1 f2 r))
      ((fp-infadd° fp-inf- f1 f2 r))
      ((fresh (s1 e1 m1 s2 e2 m2)
              ; (project (f1 f2 r) (begin (pp `(fp-+? ,f1 ,f2 ,r)) succeed))
              (fp-decomp° f1 s1 e1 m1)
              (fp-decomp° f2 s2 e2 m2)
              (conde
                ((== s1 s2) (fp-sign° r s1)
                 ; NB: fp-swap inlined here, other branches recurse
                 (conde
                   ((<o e1 e2) (fp-samesignadder° e1 m1 e2 m2 r))
                   ((<=o e2 e1) (fp-samesignadder° e2 m2 e1 m1 r))))
                ((sign-not° s1 s2)
                 (conde
                   ((fp-sign° r s1)
                    ; f1 + f2 = r <-> r + -f2 = f1
                    (fresh (t) (fp-negate° f2 t) (fp-+° r t f1)))
                   ((fp-sign° r s2)
                    ; f1 + f2 = r <-> r + -f1 = f2
                    (fresh (t) (fp-negate° f1 t) (fp-+° r t f2))))))
              ; (project (f1 f2 r) (begin (pp `(fp-+° ,f1 ,f2 ,r)) succeed))
              ))))

  ; NOT for sign bits.
  (define (sign-not° s1 s2)
    (conde ((== s1 'pos) (== s2 'neg)) ((== s1 'neg) (== s2 'pos))))

  ; XOR for sign bits.
  ; For refutational completeness, this is given as a table.
  (define (sign-xor° s1 s2 s3)
    (conde
      ((== s1 'pos) (== s2 'pos) (== s3 'pos))
      ((== s1 'pos) (== s2 'neg) (== s3 'neg))
      ((== s1 'neg) (== s2 'pos) (== s3 'neg))
      ((== s1 'neg) (== s2 'neg) (== s3 'pos))))

  ; Add exponents for multiplication.
  (define (compute-exp° e1 e2 bits re)
    (fresh (pre-re esum)
           (conde
             ((mantissa° bits) (succ° pre-re re))
             ((mantissa-short° bits) (== pre-re re)))
           (+o e1 e2 esum)
           (+o (build-num exponent-bias) pre-re esum)))

  ; Multiplication.
  ; See p6.
  (define (fp-*° f1 f2 r)
    (fresh (s1 s2 rs)
           ; (project (f1 f2 r) (begin (pp `(fp-*? ,f1 ,f2 ,r)) succeed))
           (sign-xor° s1 s2 rs)
           (fp-sign° f1 s1) (fp-sign° f2 s2) (fp-sign° r rs)
           ; (project (f1 f2 r) (begin (pp `(fp-* signs ,f1 ,f2 ,r)) succeed))
           (conde
             ; XXX wrong WRT IEEE 754!
             ((fp-zero° r) (conde ((fp-zero° f1)) ((fp-zero° f2))))
             ((fp-inf° r) (conde ((fp-inf° f1)) ((fp-inf° f2))))
             ((fresh (e1 m1 e2 m2 m1m2 re rm lsb)
                     (fp-decomp° f1 s1 e1 m1)
                     (fp-decomp° f2 s2 e2 m2)
                     (*o m1 m2 m1m2)
                     ; (project (m1 m2 m1m2) (begin (pp `(fp-* mantissas ,m1 ,m2 ,m1m2)) succeed))
                     (drop-leastsig-bit° m1m2 rm lsb)
                     (compute-exp° e1 e2 lsb re)
                     ; (project (e1 e2 re) (begin (pp `(fp-* exponents ,e1 ,e2 ,re)) succeed))
                     (conde
                       ((fp-inf° r) (fp-overflow° re))
                       ((fp-decomp° r rs re rm))))))
           ; (project (f1 f2 r) (begin (pp `(fp-*° ,f1 ,f2 ,r)) succeed))
           ))
)