possible bug in model for fp.fma

[florianschanda]

Hey,

I am not sure if the attached testcase should be ultimately SAT or UNSAT. However, I am reasonably sure that the model is bogus.

(set-logic QF_FP)
(set-option :produce-models true)
(set-info :source |Florian Schanda|)
(set-info :smt-lib-version 2.5)
(set-info :category crafted)

(declare-const x Float32)
(declare-const y Float32)
(declare-const z Float32)

(define-const x2 Float64 ((_ to_fp 11 53) RNE x))
(define-const y2 Float64 ((_ to_fp 11 53) RNE y))
(define-const z2 Float64 ((_ to_fp 11 53) RNE z))

(define-const r Float32 (fp.fma RNE x y z))
(define-const r2 Float32
  ((_ to_fp 8 24) RNE (fp.add RNE (fp.mul RNE x2 y2) z2)))

(assert (not (= r r2)))
(check-sat)

(get-value (x y z))
(get-value (x2 y2 z2))
(get-value (r r2))
(get-value ((= r r2)))

(exit)

If I run this with the current master branch I get:

sat
((x (fp #b0 #x80 #b10111011101101010010011))
 (y (fp #b1 #xfe #b00000000010100101111000))
 (z (fp #b0 #xfe #b11111000100010011100101)))
((x2 (fp #b0 #b10000000000 #xbbb5260000000))
 (y2 (fp #b1 #b10001111110 #x0052f00000000))
 (z2 (fp #b0 #b10001111110 #xf889ca0000000)))
((r  (fp #b1 #xfe #b10000000000000000000001))
 (r2 (fp #b1 #xfe #b10000000000000000000001)))
(((= r r2) true))

In particular r and r2 seem quite equal in this model, but that is not what I asserted. (And yes, I did intend to use '=' and not 'fp.eq'.)

[wintersteiger]

The model itself (via (get-model)) is actually ok, but something's wrong with the evaluation for get-value.

[florianschanda]

Fair. Probably nothing at all to do with fma then :)

[ikirill]

I can also reproduce this, and in my case the output of (get-model) seems bogus too:

(set-option :produce-models true)
(set-logic QF_FP)
(define-fun zero () (_ FloatingPoint 8 24) (_ +zero 8 24))
(define-fun half () (_ FloatingPoint 8 24) ((_ to_fp 8 24) roundNearestTiesToEven (/ 1 2)))
(define-fun one () (_ FloatingPoint 8 24) ((_ to_fp 8 24) roundNearestTiesToEven (/ 1 1)))
(define-fun four () (_ FloatingPoint 8 24) ((_ to_fp 8 24) roundNearestTiesToEven (/ 4 1)))
(define-fun eps () (_ FloatingPoint 8 24) ((_ to_fp 8 24) roundNearestTiesToEven (/ 1 8388608)))

(declare-fun x () (_ FloatingPoint 8 24))
(declare-fun y () (_ FloatingPoint 8 24))
(declare-fun z () (_ FloatingPoint 8 24))
(declare-fun t () (_ FloatingPoint 8 24))
(define-fun w () (_ FloatingPoint 8 24) (fp.fma roundNearestTiesToEven x (fp.sub roundNearestTiesToEven one t) (fp.mul roundNearestTiesToEven y t)))

(assert
 (and
  (and (fp.isNormal x)
     (fp.isNormal y)
     (fp.isNormal z)
     (fp.isNormal t))
  (and (fp.leq zero t)
     (fp.leq t one))
  (and (fp.leq one x)
     (fp.leq x y)
     (fp.leq y four))
  (fp.eq z w)
  (fp.lt z x)
  ))

(check-sat)
(get-value (z
            w
            (fp.lt zero z)
            (fp.lt z x)
            (fp.eq z w)
            (fp.eq z (fp.fma roundNearestTiesToEven x (fp.sub roundNearestTiesToEven one t) (fp.mul roundNearestTiesToEven y t)))
            (fp.to_real x)
            (fp.to_real y)
            (fp.to_real t)
            (fp.to_real z)
            (fp.to_real w)
            (fp.to_real (fp.fma roundNearestTiesToEven x (fp.sub roundNearestTiesToEven one t) (fp.mul roundNearestTiesToEven y t)))
            ))
(get-model)

Output ((fp.eq z w) was asserted above, but is false):

testfloat $ z3 -smt2 linterp_fma.smt
sat
((z (fp #b0 #x7a #b00000001011011101100110))
 (w (fp #b0 #x80 #b00000100000001011011101))
 ((fp.lt zero z) true)
 ((fp.lt z x) true)
 ((fp.eq z w) false)
 ((fp.eq z (fp.fma roundNearestTiesToEven x (fp.sub roundNearestTiesToEven one t) (fp.mul roundNearestTiesToEven y t))) false)
 ((fp.to_real x) (/ 12113169.0 8388608.0))
 ((fp.to_real y) (/ 8109701.0 2097152.0))
 ((fp.to_real t) (/ 4067395.0 16777216.0))
 ((fp.to_real z) (/ 4217779.0 134217728.0))
 ((fp.to_real w) (/ 8520413.0 4194304.0))
 ((fp.to_real (fp.fma roundNearestTiesToEven x (fp.sub roundNearestTiesToEven one t) (fp.mul roundNearestTiesToEven y t))) (/ 8520413.0 4194304.0)))
(model
  (define-fun z () (_ FloatingPoint 8 24)
    (fp #b0 #x7a #b00000001011011101100110))
  (define-fun t () (_ FloatingPoint 8 24)
    (fp #b0 #x7c #b11110000100000100001100))
  (define-fun y () (_ FloatingPoint 8 24)
    (fp #b0 #x80 #b11101110111110100001010))
  (define-fun x () (_ FloatingPoint 8 24)
    (fp #b0 #x7f #b01110001101010100010001))
)

w gets correctly evaluated in terms of fma(x, 1-t, y*t) both in (get-model) and (get-value), but z is still not equal to w, and the value assigned to z seems bogus.

I also checked that the variable values in get-model match those in get-value. The value for w = fma(x, 1-t, y*t) differs from what I get in Julia for the same values of x, y, t, but it seems only by a single ulp.

julia> z = reinterpret(Float32, (UInt32(0x7a)<<23) | 0b00000001011011101100110)
0.031424902f0

julia> w = reinterpret(Float32, (UInt32(0x80)<<23) | 0b00000100000001011011101)
2.0314248f0

julia> t = reinterpret(Float32, (UInt32(0x7c)<<23) | 0b11110000100000100001100)
0.24243563f0

julia> y = reinterpret(Float32, (UInt32(0x80)<<23) | 0b11101110111110100001010)
3.8670068f0

julia> x = reinterpret(Float32, (UInt32(0x7f)<<23) | 0b01110001101010100010001)
1.4440023f0

julia> x*(1-t)+y*t
2.031425f0

julia> fma(x, 1f0-t, y*t)
2.031425f0

These also match the fp.to_real values above.

I don't get anything like this if I replace fma by straight multiplication-addition.

I'm using z3 on OSX 10.11 installed with brew install z3 --HEAD:

testfloat $ z3 --version
Z3 version 4.5.1 - 64 bit

(Note: I couldn't figure out how to print floating-point values directly in the usual format in smt, I'm not familiar with it at all.)

[ikirill]

A smaller, maybe clearer example. I was wondering whether the bad fma value was somehow random, but it seems to be uniquely, incorrectly, determined. In particular, it behaves as expected when t is defined as a constant with define-fun, unlike when it is a variable (declare-fun) to be found.

(set-option :produce-models true)
(set-logic QF_FP)
(define-fun tval () (_ FloatingPoint 8 24) ((_ to_fp 8 24) roundNearestTiesToEven (/ 4067395 16777216))) ;; 0.24243563
(define-fun one () (_ FloatingPoint 8 24) ((_ to_fp 8 24) roundNearestTiesToEven (/ 1 1))) ;; 1.0
(define-fun xval () (_ FloatingPoint 8 24) ((_ to_fp 8 24) roundNearestTiesToEven (/ 12113169 8388608))) ;; 1.4440023
(define-fun yval () (_ FloatingPoint 8 24) ((_ to_fp 8 24) roundNearestTiesToEven (/ 8109701 2097152))) ;; 3.8670068
(define-fun badzval () (_ FloatingPoint 8 24) ((_ to_fp 8 24) roundNearestTiesToEven (/ 4217779 134217728))) ;; 3.1424902e-2

(declare-fun z () (_ FloatingPoint 8 24))
(declare-fun t () (_ FloatingPoint 8 24)) ;; causes sat
;; causes unsat, and if I remove (fp.lt z xval) then I get the correct value for fma
;; (define-fun t () (_ FloatingPoint 8 24) tval)

(define-fun rhs () (_ FloatingPoint 8 24) (fp.fma roundNearestTiesToEven xval (fp.sub roundNearestTiesToEven one t) (fp.mul roundNearestTiesToEven yval t)))

(assert
 (and
  (fp.eq t tval)
  ;; uncomment one:
  (fp.lt z xval) ;; sat
  ;; (fp.lt z badzval) ;; unsat
  ;; (fp.gt z badzval) ;; unsat
  ;; (fp.eq z badzval) ;; sat
  ;; (fp.eq badzval rhs) ;; sat
  (fp.eq z rhs)))

(check-sat)
(get-value ((fp.eq z rhs)
            (fp.eq z badzval)
            (fp.to_real z)
            (fp.to_real t)
            (fp.to_real rhs)))
(get-model)

[wintersteiger]

Finally found some time to look into this. There's something wrong with macros (define-fun/define-const), so for now just rewriting (define-fun x ... t) into (declare-fun x ...) (= x t) could be a viable workaround until I figure out what's actually wrong.

[ikirill]

That workaround doesn't seem to work for me with current git head. (Did I misunderstand your suggestion?) I got rid of define-funs, but it still seems to fail:

(set-option :produce-models true)
(set-logic QF_FP)

(declare-fun tval () (_ FloatingPoint 8 24))
(assert (fp.eq tval ((_ to_fp 8 24) roundNearestTiesToEven (/ 4067395 16777216)))) ;; 0.24243563
(declare-fun one () (_ FloatingPoint 8 24))
(assert (fp.eq one ((_ to_fp 8 24) roundNearestTiesToEven (/ 1 1)))) ;; 1.0
(declare-fun xval () (_ FloatingPoint 8 24))
(assert (fp.eq xval ((_ to_fp 8 24) roundNearestTiesToEven (/ 12113169 8388608)))) ;; 1.4440023
(declare-fun yval () (_ FloatingPoint 8 24))
(assert (fp.eq yval ((_ to_fp 8 24) roundNearestTiesToEven (/ 8109701 2097152)))) ;; 3.8670068
(declare-fun badzval () (_ FloatingPoint 8 24))
(assert (fp.eq badzval ((_ to_fp 8 24) roundNearestTiesToEven (/ 4217779 134217728)))) ;; 3.1424902e-2

(declare-fun z () (_ FloatingPoint 8 24))
(declare-fun t () (_ FloatingPoint 8 24))

(declare-fun rhs () (_ FloatingPoint 8 24))
(assert (fp.eq rhs (fp.fma roundNearestTiesToEven xval (fp.sub roundNearestTiesToEven one t) (fp.mul roundNearestTiesToEven yval t))))

(assert
 (and
  (fp.eq t tval)
  (fp.eq z rhs)
  (fp.lt z xval) ;; sat, but should be unsat
  ))

(check-sat)
(get-value ((fp.eq z rhs)
            (fp.eq z badzval)
            (fp.to_real z)
            (fp.to_real t)
            (fp.to_real rhs)))
(get-model)

Output:

testfloat $ ~/software/z3/build/z3 -smt2 test_t7_3.smt
sat
(((fp.eq z rhs) true)
 ((fp.eq z badzval) true)
 ((fp.to_real z) (/ 4217779.0 134217728.0))
 ((fp.to_real t) (/ 4067395.0 16777216.0))
 ((fp.to_real rhs) (/ 4217779.0 134217728.0)))
(model
  (define-fun t () (_ FloatingPoint 8 24)
    (fp #b0 #x7c #b11110000100000100001100))
  (define-fun tval () (_ FloatingPoint 8 24)
    (fp #b0 #x7c #b11110000100000100001100))
  (define-fun one () (_ FloatingPoint 8 24)
    (fp #b0 #x7f #b00000000000000000000000))
  (define-fun badzval () (_ FloatingPoint 8 24)
    (fp #b0 #x7a #b00000001011011101100110))
  (define-fun rhs () (_ FloatingPoint 8 24)
    (fp #b0 #x7a #b00000001011011101100110))
  (define-fun z () (_ FloatingPoint 8 24)
    (fp #b0 #x7a #b00000001011011101100110))
  (define-fun yval () (_ FloatingPoint 8 24)
    (fp #b0 #x80 #b11101110111110100001010))
  (define-fun xval () (_ FloatingPoint 8 24)
    (fp #b0 #x7f #b01110001101010100010001))
)

[wintersteiger]

@ikirill Yes, I was wrong. There's a proper bug in fp.fma and I have an input that triggers it, but I haven't had time to dig deep enough yet.

[ikirill]

(Edit: call this case 1) Interesting: my testcase above now produces unsat with the commit 0f158330, but I still have this one (it asserts that linearly interpolating between x and y might produce a value less than x—could be either sat/unsat I guess, but the model is wrong):

(set-option :produce-models true)
(set-logic QF_FP)
(define-fun zero () (_ FloatingPoint 8 24) (_ +zero 8 24))
(define-fun half () (_ FloatingPoint 8 24) ((_ to_fp 8 24) roundNearestTiesToEven (/ 1 2)))
(define-fun one () (_ FloatingPoint 8 24) ((_ to_fp 8 24) roundNearestTiesToEven (/ 1 1)))
(define-fun four () (_ FloatingPoint 8 24) ((_ to_fp 8 24) roundNearestTiesToEven (/ 4 1)))
(define-fun eps () (_ FloatingPoint 8 24) ((_ to_fp 8 24) roundNearestTiesToEven (/ 1 8388608)))

(declare-fun x () (_ FloatingPoint 8 24))
(declare-fun y () (_ FloatingPoint 8 24))
(declare-fun z () (_ FloatingPoint 8 24))
(declare-fun t () (_ FloatingPoint 8 24))
(define-fun w () (_ FloatingPoint 8 24) (fp.fma roundNearestTiesToEven x (fp.sub roundNearestTiesToEven one t) (fp.mul roundNearestTiesToEven y t)))

(assert
 (and
  (and (fp.isNormal x) (fp.isNormal y) (fp.isNormal z) (fp.isNormal t))
  (and (fp.leq zero t) (fp.leq t one))
  (and (fp.leq one x) (fp.leq x y) (fp.leq y four))
  (fp.eq z w)
  (fp.lt z x)))

(check-sat)
(get-value (z
            w
            (fp.lt zero z)
            (fp.lt z x)
            (fp.eq z w)
            (fp.eq z (fp.fma roundNearestTiesToEven x (fp.sub roundNearestTiesToEven one t) (fp.mul roundNearestTiesToEven y t)))
            (fp.to_real x)
            (fp.to_real y)
            (fp.to_real t)
            (fp.to_real z)
            (fp.to_real w)
            (fp.to_real (fp.fma roundNearestTiesToEven x (fp.sub roundNearestTiesToEven one t) (fp.mul roundNearestTiesToEven y t)))
            ))
(get-model)

with the output

sat
((z (fp #b0 #x7f #b01101000000110101011110))
 (w (fp #b0 #x80 #b01101000000110101011110))
 ((fp.lt zero z) true)
 ((fp.lt z x) true)
 ((fp.eq z w) false)
 ((fp.eq z (fp.fma roundNearestTiesToEven x (fp.sub roundNearestTiesToEven one t) (fp.mul roundNearestTiesToEven y t))) false)
 ((fp.to_real x) (/ 11799901.0 4194304.0))
 ((fp.to_real y) (/ 3795615.0 1048576.0))
 ((fp.to_real t) (/ 8886119.0 35184372088832.0))
 ((fp.to_real z) (/ 5899951.0 4194304.0))
 ((fp.to_real w) (/ 5899951.0 2097152.0))
 ((fp.to_real (fp.fma roundNearestTiesToEven x (fp.sub roundNearestTiesToEven one t) (fp.mul roundNearestTiesToEven y t))) (/ 5899951.0 2097152.0)))
(model
  (define-fun z () (_ FloatingPoint 8 24)
    (fp #b0 #x7f #b01101000000110101011110))
  (define-fun t () (_ FloatingPoint 8 24)
    (fp #b0 #x69 #b00001111001011101100111))
  (define-fun y () (_ FloatingPoint 8 24)
    (fp #b0 #x80 #b11001111010101001111100))
  (define-fun x () (_ FloatingPoint 8 24)
    (fp #b0 #x80 #b01101000000110101011101))
)

where it asserts (fp.eq z w), but then in the model it's false.

Checking in Julia shows z is off by just a factor of 2 from w, but w is correct (in single precision—I made an error pasting the values because double is the default, hence the edits).

julia> x, y, z, t, w = 11799901.0f0 /4194304.0f0, 3795615.0f0 /1048576.0f0, 5899951.0f0/ 4194304.0f0, 8886119.0f0 / 35184372088832.0f0, 5899951.0f0 /2097152.0f0;
julia> fma(x, 1-t, y*t)
2.8133159f0
julia> w
2.8133159f0
julia> z
1.4066579f0
julia> 2z
2.8133159f0

[ikirill]

(Edit: call it case 2). Testcase with round numbers, in case it helps (edit: even simpler now):

(set-option :produce-models true)
(set-logic QF_FP)
(define-fun zero () (_ FloatingPoint 8 24) (_ +zero 8 24))
(define-fun one () (_ FloatingPoint 8 24) ((_ to_fp 8 24) roundNearestTiesToEven (/ 1 1)))

(declare-fun t () (_ FloatingPoint 8 24))
(declare-fun z () (_ FloatingPoint 8 24))
(define-fun w () (_ FloatingPoint 8 24) (fp.fma roundNearestTiesToEven one t one))

(assert
 (and
  (fp.eq z w)
  (fp.leq zero t)
  (fp.lt z one)))

(check-sat)
(get-value ((fp.to_real t)
            (fp.to_real z)
            (fp.to_real w)
            (fp.lt z one)
            (fp.lt w one)))
(get-model)

sat
(((fp.to_real t) 0.0)
 ((fp.to_real z) (/ 1.0 2.0))
 ((fp.to_real w) 1.0)
 ((fp.lt z one) true)
 ((fp.lt w one) false))
(model
  (define-fun z () (_ FloatingPoint 8 24)
    (fp #b0 #x7e #b00000000000000000000000))
  (define-fun t () (_ FloatingPoint 8 24)
    (_ +zero 8 24))
)

[wintersteiger]

There are at least two different bugs involved here; could you check your examples again, but without using define-fun?

[ikirill]

The second testcase (https://github.com/Z3Prover/z3/issues/872#issuecomment-317286855) (I think this is what you mean with define-fun, right?)

(set-option :produce-models true)
(set-logic QF_FP)
(declare-fun zero () (_ FloatingPoint 8 24))
(assert (= zero (_ +zero 8 24)))
(declare-fun one () (_ FloatingPoint 8 24))
(assert (= one ((_ to_fp 8 24) RNE 1.0 0)))

(declare-fun t () (_ FloatingPoint 8 24))
(declare-fun z () (_ FloatingPoint 8 24))
(declare-fun w () (_ FloatingPoint 8 24))
(assert (= w (fp.fma RNE one t one)))

(assert
 (and
  (fp.eq z w)
  (fp.leq zero t)
  (fp.lt z one)))

(check-sat)
(get-value ((fp.to_real t)
            (fp.to_real z)
            (fp.to_real w)
            (fp.lt z one)
            (fp.lt w one)))
(get-model)

Output:

sat
(((fp.to_real t) (/ 1.0 16384.0))
 ((fp.to_real z) (/ 16385.0 32768.0))
 ((fp.to_real w) (/ 16385.0 32768.0))
 ((fp.lt z one) true)
 ((fp.lt w one) true))
(model
  (define-fun one () (_ FloatingPoint 8 24)
    (fp #b0 #x7f #b00000000000000000000000))
  (define-fun w () (_ FloatingPoint 8 24)
    (fp #b0 #x7e #b00000000000001000000000))
  (define-fun t () (_ FloatingPoint 8 24)
    (fp #b0 #x71 #b00000000000000000000000))
  (define-fun z () (_ FloatingPoint 8 24)
    (fp #b0 #x7e #b00000000000001000000000))
  (define-fun zero () (_ FloatingPoint 8 24)
    (_ +zero 8 24))
)

Here fma(1, t, 1) is off by a factor of 2:

julia> t, z, w = Float32(1.0/16384.0), Float32(16385.0/32768.0), Float32(16385.0/32768.0)
(6.1035156f-5,0.5000305f0,0.5000305f0)
julia> fma(1f0, t, 1f0)
1.000061f0
julia> 2z
1.000061f0

[ikirill]

Same thing with the first testcase (https://github.com/Z3Prover/z3/issues/872#issuecomment-317284153): (fp.eq z w) is now good, but the value of z is off by a factor of 2.

[wintersteiger]

The second testcase

Did you at any point copy/paste any of your test cases? If so, please stop, neither github nor my browser does do a diff for me; just number them and refer to them by number if you need to refer to a test case a second time.

[ikirill]

I inserted links/numbers into my earlier comments. Is it clearer now? I don't know how to post github comments in a way that allows diffs.

[wintersteiger]

@ikirill I'm not aware of any such features either, that's why I asked for it to be done manually. Could you clarify which of your test cases you expect to be satisifable? "Off by factor of x" really doesn't make sense for unsatisfiable problems.

[ikirill]

@wintersteiger The only change I made was replacing define-fun by declare-fun + (assert (= ...)), so I expect satisfiability not to change.

I don't know if case 1 is sat/unsat, but case 2 asserts that fma(1, t, 1) < 1 and so it should be unsat.

I didn't think satisfiability mattered for what I wrote, but maybe that's wrong, and you certainly understand this better than me. All I saw was that when I asserted z = fma(1, t, 1), and the model gave me the values of z and t, I plugged them into a correct fma implementation and saw that z's value was 0.5*fma(1, t, 1), not the correct value fma(1,t,1). I interpreted this to mean that z3 "evaluated" fma(1,t,1) incorrectly, producing a value off by a factor of two—presumably some problem with computing the final exponent, because the significand is correct.

[wintersteiger]

Two fixes added. Thanks again for the test cases, they revealed two bugs in the floats; one in the symbolic encoding of fp.fma and the second was a corner case in the concrete float (mpf) rounder, which is also used for all other fp.* operations, so there's a chance I broke one thing while fixing the other. I'm not getting any problems reported on the SMTLIB benchmarks with or without model_validate=true, but please re-run any benchmarks or test-cases on your end. I'm keeping the issue open for now in case something new comes up.

[wintersteiger]

@ikirill sounds good, I think I read your example the right way round then. The "evalautor" (mpfs) was indeed off. The exponent was fine, but the sticky bit in the normalization shift in the rounder was misaligned for some combinations of inputs.

[ikirill]

Okay, I think I have one more (case 3). This one asserts that fma(x, y, z) != Float32(fma(Float128(x), Float128(y), Float128(z)). I suspect it might be sat—I couldn't find a reference on how much precision is necessary for double rounding to be harmless for fma, and z3 gave me a valid counterexample for Float64, but the model it gives using Float128 seems incorrect.

(set-option :produce-models true)
(set-option :pp.fp_real_literals true)
(set-logic QF_FP)

(declare-const zero Float32)
(assert (= zero (_ +zero 8 24)))
(declare-const one Float32)
(assert (= one ((_ to_fp 8 24) RNE 1.0)))
;; (declare-const eps Float32)
;; (assert (= eps ((_ to_fp 8 24) RNE 1.0 -23)))
;; (declare-const eps2 Float32)
;; (assert (= eps2 ((_ to_fp 8 24) RNE 1.0 -46)))

(declare-const low Float32)
(assert (= low ((_ to_fp 8 24) RNE 1.0)))
(declare-const high Float32)
(assert (= high ((_ to_fp 8 24) RNE 2.0)))

(declare-const x Float32)
(declare-const y Float32)
(declare-const z Float32)

;; Different constraints, so that the counterexamples aren’t too extreme
;; (assert (and (fp.leq zero x) (fp.leq x y) (fp.leq y one)))
(assert (and (fp.leq low x) (fp.leq x y) (fp.leq y high)))
;; (assert (and (fp.leq eps z) (fp.leq z one)))
(assert (and (fp.leq zero z) (fp.leq z one)))
(assert (and (fp.isNormal x) (fp.isNormal y) (fp.isNormal z)))

(declare-const sfma Float32)
(assert (= sfma (fp.fma RNE x y z)))

(declare-const dfma Float32)
(assert (= dfma ((_ to_fp 8 24) RNE (fp.fma RNE ((_ to_fp 15 113) RNE x) ((_ to_fp 15 113) RNE y) ((_ to_fp 15 113) RNE z)))))
;; (assert (= dfma ((_ to_fp 8 24) RNE (fp.fma RNE ((_ to_fp 11 53) RNE x) ((_ to_fp 11 53) RNE y) ((_ to_fp 11 53) RNE z)))))

(declare-const dma Float32)
(assert (= dma ((_ to_fp 8 24) RNE (fp.add RNE (fp.mul RNE ((_ to_fp 11 53) RNE x) ((_ to_fp 11 53) RNE y)) ((_ to_fp 11 53) RNE z)))))

(assert (not (= sfma dfma)))
;; (assert (not (= sfma dma)))
;; (assert (or (not (= sfma dfma)) (not (= sfma dma))))

(check-sat)
(get-value ((fp.to_real x)
            (fp.to_real y)
            (fp.to_real z)
            (fp.to_real sfma)
            (fp.to_real dfma)
            (fp.to_real dma)))
(get-model)

sat
(((fp.to_real x) (/ 9943857.0 8388608.0))
 ((fp.to_real y) (/ 16265305.0 8388608.0))
 ((fp.to_real z) (/ 12251097.0 281474976710656.0))
 ((fp.to_real sfma) (/ 18829.0 8192.0))
 ((fp.to_real dfma) (/ 9640447.0 4194304.0))
 ((fp.to_real dma) (/ 9640447.0 4194304.0)))
(model
  (define-fun y () Float32
    ((_ to_fp 8 24) RTZ 1.93897545337677001953125 0))
  (define-fun zero () Float32
    (_ +zero 8 24))
  (define-fun sfma () Float32
    ((_ to_fp 8 24) RTZ 1.14923095703125 1))
  (define-fun dma () Float32
    ((_ to_fp 8 24) RTZ 1.14923083782196044921875 1))
  (define-fun low () Float32
    ((_ to_fp 8 24) RTZ 1.0 0))
  (define-fun one () Float32
    ((_ to_fp 8 24) RTZ 1.0 0))
  (define-fun z () Float32
    ((_ to_fp 8 24) RTZ 1.46044456958770751953125 -25))
  (define-fun dfma () Float32
    ((_ to_fp 8 24) RTZ 1.14923083782196044921875 1))
  (define-fun high () Float32
    ((_ to_fp 8 24) RTZ 1.0 1))
  (define-fun x () Float32
    ((_ to_fp 8 24) RTZ 1.18540012836456298828125 0))
)

I checked using MPFR (i.e.: Float32(fma(BigFloat(x), BigFloat(y), BigFloat(z)))) and the value it gives for fma(x,y,z) (sfma = 2.298462f0) is one ulp more than what MPFR gives (dfma = 2.2984617f0—it also agrees with double rounding through Float64).

[florianschanda]

Great - I'll rerun all my benchmarks and let you know :)

[ikirill]

(Case 4, smaller+clearer, values taken from the output of case 3.)

(set-logic QF_FP)
(set-info :status unsat)
(declare-const x Float32)
(declare-const y Float32)
(declare-const z Float32)
(declare-const r Float32)
(assert (= x (fp #b0 #b01111111 #b11011101101101101000101))) ; 1.8660666f0
(assert (= y (fp #b0 #b01111111 #b11101100000111000100001))) ; 1.9223062f0
(assert (= z (fp #b0 #b01100111 #b10100111110110000110101))) ; 9.868419f-8
(assert (= r (fp #b0 #b10000000 #b11001011001001111100011))) ; 3.5871513f0 ; right
;; (assert (= r (fp #b0 #b10000000 #b11001011001001111100100))) ; wrong
(assert (not (= (fp.fma RNE x y z) r)))
(check-sat)
(exit)

sat
(error "line 15 column 10: check annotation that says unsat")

[ikirill]

(Case 5, similar to case 4, numbers from output of case 3) Interesting: both this and case 4 have z < eps(Float32), and I couldn't find a counterexample with z > eps(Float32).

(set-logic QF_FP)
(set-info :status unsat)
(declare-const x Float32)
(declare-const y Float32)
(declare-const z Float32)
(declare-const r Float32)
(assert (= x (fp #b0 #b01111111 #b00101111011101100110001))) ; 1.1854001f0
(assert (= y (fp #b0 #b01111111 #b11110000011000001011001))) ; 1.9389755f0
(assert (= z (fp #b0 #b01100110 #b01110101110111111011001))) ; 4.352464f-8
(assert (= r (fp #b0 #b10000000 #b00100110001100111111111))) ; 2.2984617f0
;; (assert (= r (fp #b0 #b10000000 #b00100110001101000000000))) ; wrong
(assert (not (= (fp.fma RNE x y z) r)))
(check-sat)
(get-model)
(exit)

sat
(error "line 15 column 10: check annotation that says unsat")

[wintersteiger]

@Ikirill: regarding "right" and "wrong", who is right or wrong? Z3, you, or me?

[wintersteiger]

@ikirill: case 3: What are all the commented lines? Are those supposed to be alternatives that could/should/must be sat/unsat/right/wrong as well or are they intended to be ignored?

[ikirill]

@wintersteiger I understand your point about right/wrong, that crossed my mind too. I used Float32, Float64, MPFR, and exact rational arithmetic in Julia to check. So for case 4 I get

julia> x = Float32(Rational(1.8660666f0) * Rational(1.9223062f0) + Rational(9.868419f-8))
3.5871513f0
julia> println(bits(x))
01000000011001011001001111100011

and case 5:

julia> Float32(Rational{BigInt}(1.1854001f0) * Rational{BigInt}(1.9389755f0) + Rational{BigInt}(4.352464f-8))
2.2984617f0

and so I marked this "right". "Wrong" is what z3 gave me. (If there's a bug in Julia or openlibm, then that's also interesting but doesn't belong here.)

Oh right: ignore the commented lines, I had them in my files because I was figuring out how I could modify the testcase and what the differences were, but they aren't for running. I didn't want to drop the unused code because it would be so much more unclear what the motivation was. I don't know if case 3 is sat/unsat, but the values in the model seem incorrect.

[ikirill]

Also, the bit patterns in the significand of case 5 look a little bit special

julia> bits(Float64(Rational{BigInt}(1.1854001f0) * Rational{BigInt}(1.9389755f0) + Rational{BigInt}(4.352464f-8)))
"0100000000000010011000110011111111101111111111111111111111101000"
julia> bits(Float32(Rational{BigInt}(1.1854001f0) * Rational{BigInt}(1.9389755f0) + Rational{BigInt}(4.352464f-8)))
"01000000000100110001100111111111"

Compare:

0010011000110011111111101111111111111111111111101000
00100110001100111111111

(So if there were an extra one at the end "001001100011001111111111", as if rounding 011111... to 10000..., then that would give z3's answer with RNE—you can't for example round in two steps, from 53 to 25 to 24.)

[wintersteiger]

Yes, the core of the issue (I think) is that the rounder circuits keep an extra bit of precision, so they can "cheat" when they shouldn't be able to. I'm cleaning it up at the moment, but may take a while as it's a bit tricky to fix without breaking other things.

[ikirill]

(Case 6) In this one, the result is really off, not just a single ulp.

(set-option :pp.fp_real_literals true)
(set-logic QF_FP)
(set-info :status unsat)
(declare-const x Float32)
(declare-const y Float32)
(declare-const z Float32)
(declare-const r Float32)
(declare-const w Float32)
(assert (= x (fp #b0 #b01111111 #b00000000101111111111111))) ; 1.0029296f0
(assert (= y (fp #b0 #b01111111 #b00000001001110001010111))) ; 1.0047711f0
(assert (= z (fp #b1 #b01111111 #b00000001111110011001011))) ; -1.0077146f0
(assert (= r (fp #b0 #b01100110 #b00000000000111010100100))) ; 2.9815638f-8
(assert (= w (fp.fma RNE x y z)))
(assert (not (= w r)))
(check-sat)
(get-model)
(exit)

sat
(error "line 15 column 10: check annotation that says unsat")
(model
  (define-fun w () Float32
    ((_ to_fp 8 24) RTZ 1.5 -25))
  (define-fun z () Float32
    ((_ to_fp 8 24) RTZ -1.00771462917327880859375 0))
  (define-fun y () Float32
    ((_ to_fp 8 24) RTZ 1.00477111339569091796875 0))
  (define-fun x () Float32
    ((_ to_fp 8 24) RTZ 1.00292956829071044921875 0))
  (define-fun r () Float32
    ((_ to_fp 8 24) RTZ 1.000446796417236328125 -25))
)

The exponent is good, but the significand is 1.5 instead of 1.0004...

[ikirill]

(Case 7) I discovered case 6 by trying to prove that 2MultFMA is correct.

(set-info :smt-lib-version 2.6)
(set-option :pp.fp_real_literals true)
(set-logic QF_FP)
(set-info :status unsat) ; 2MultFMA is correct

(declare-const x Float32)
(declare-const y Float32)
(assert (and (fp.isNormal x) (fp.isNormal y)))

;; Might be sat due to over-/under-flows if you remove these bounds.
(assert (and (fp.leq ((_ to_fp 8 24) RNE 1.0) x)
           (fp.leq x y)
           (fp.leq y ((_ to_fp 8 24) RNE 2.0))))

(declare-const xy32 Float32)
(assert (= xy32 (fp.mul RNE x y)))

(declare-const xy64 Float64)
(assert (= xy64 (fp.mul RNE ((_ to_fp 11 53) RNE x) ((_ to_fp 11 53) RNE y))))
(assert (fp.isNormal ((_ to_fp 8 24) RNE xy64)))

(declare-const hi Float32)
(declare-const lo Float32)
(assert (= hi (fp.mul RNE x y)))
(assert (= lo (fp.fma RNE x y (fp.neg hi))))
(declare-const xy Float64)
(assert (= xy (fp.add RNE ((_ to_fp 11 53) RNE hi) ((_ to_fp 11 53) RNE lo))))

(assert (not (= xy64 xy)))

(check-sat)
(get-model)

[ikirill]

When writing these tests, what's a good way to express a constraint like exponent(x) >= ...?

[wintersteiger]

There are no accessor functions for the components of SMT floats because they are not assumed to be represented by bit-vectors. Z3 has to_ieee_bv, but that's an extension with whatever semantics I thought suitable at the time I added it. I guess the best way to do that is to simply build the smallest float you consider in range, e.g. (fp.leq ... (fp #b0 #bsmallest_exponent #b0...0)) (+- sign).

From now on I will not consider any more new test cases for this issue if they don't come with an explanation of what is going wrong and a suggested fix. (I have my own random test generator...)

[wintersteiger]

@ikirill For case 4, are you absolutely positive that the "right" result is the correct one? Is there any documentation on what Julia does for FMA when the hardware doesn't support it? I'm using the pre-compiled binaries for Windows and my CPU doesn't support any FMA instructions, but it runs fine, so obviously it's doing something in software.

I've gone through these arguments many times yesterday and today and every time I reach the same conclusion: the rounding bits are such that we have to round up to ...100 instead of ...011.

[ikirill]

Do you see the same also for case 6? That one isn't a single ulp.

Julia's fma falls back on openlibm (you can check whether it does with code_native(fma, (Float32,Float32,Float32))), as here: https://github.com/JuliaLang/julia/blob/85a2555e539542c048fea99cd326d3fd8fb6da68/base/floatfuncs.jl#L263 And fmaf in openlibm is here: https://github.com/JuliaLang/openlibm/blob/1581174c85f7b645b15ba1ac1c3a98fb601f0fe7/src/s_fmaf.c#L43 Both Base.fma_llvm (native) and Base.fma_libm produce the same result for me on OSX, as well as C++'s fma, so either my native fma is wrong or something went wrong setting up the test.

(Incidentally, at this point, I believe there is still a problem where julia doesn't use native fma out of the box and you need to rebuild its system image. If your cpu is literally missing the FMA cpu flag, not just the compiler not being aware of it, you could also try C++11's std::fma, whatever that fallback is.)

I don't really understand what you mean. Based on my calculations (BigFloat, which is MPFR), the exact answer has the significand (which fits into Float64)

1100101100100111110001101111111111111111111111110000

and so there is no tie and 0.01111....1110000 should be rounded to zero giving

11001011001001111100011

because in binary 0.0111.... < 0.1 = 1/2. Did I do this wrong?

Thus I get

julia> fma(1.86606657505035400390625, 1.92230618000030517578125, 9.868418970881975837983191013336181640625e-8)
3.5871514081954885
julia> Float32(3.5871514081954885)
3.5871513f0
julia> println(bits(3.5871513f0))
01000000011001011001001111100011

I guess the reason I trust the software fma is because I find the extended precision arithmetic there so much easier to follow than all the bit manipulations—I can read openlibm's fma, but I don't understand what mk_fma is doing.

[ikirill]

I have my own random test generator

It seems to me like if there are tricky corner cases in rounding, a random test generator might not be able to catch them if the probability of generating one is sufficiently low. All of the test cases I have I generated by trying to prove properties of floating-point arithmetic that I suspect to be true, to which z3 gave me counterexamples. That's why I think there might still be a problem.

[ikirill]

Another thing I noticed: if I drop "f0" and treat the decimal constants as Float64's, then I get z3's answer:

julia> Float32(fma(Float64(1.8660666f0), Float64(1.9223062f0), Float64(9.868419f-8)))
3.5871513f0
julia> Float32(fma(1.8660666, 1.9223062, 9.868419e-8))
3.5871515f0
julia> FMA.smtbv(Float32(fma(1.8660666, 1.9223062, 9.868419e-8)))
"#(fp #b0 #b10000000 #b11001011001001111100100)"

but I don't think this matters, the choice of numbers here must be really sensitive.

[wintersteiger]

Case 6 is completely different, probably due to a bug I introduced while added the previous fix. I haven't look at it at all yet.

The precise result of the multiplication in case 4 is (2 sbits long, just bit-vector multiplication of the normalized significands): 111001011001001111100011000101100000100111100101

The aligned third input is (bunch of zeros because its exponent of -24 is aligned with the 0-exponent of the other numbers): 0000000000000000000000000110100111110110000110101

Adding these two we get: 111001011001001111100011111010011111011000011010

Extracting (24 + 4) digits with sticky flag produces 1.11001011001001111100011,111

(dot indicating the uppermost bit and comma separating the three extra rounding bits). Clearly, ... 1,111 must be rounded up for RNE.

So, at least one of the tools is off by at least one bit. I'm essentially rewriting the concrete fp.fma now (the symbolic encoding should be less buggy, but we can't compare test results between them at the moment). This will take a considerable amount of time which I currently don't have, so I expect quite some time will pass until this is fixed.

[ikirill]

I'm confused, it looks to me like adding the bitvectors directly (numbers below copy-pasted from your comment) there are no carry bits, so I don't get your result (again, apologies if I'm being slow here):

  111001011001001111100011000101100000100111100101 (copy-pasted)
+ 0000000000000000000000000110100111110110000110101 (copy-pasted)
= 1110010110010011111000110111111111111111111111111 (computed by eye, matches what I had above)
≠ 111001011001001111100011111010011111011000011010 (copy-pasted)

[ikirill]

Did you shift aligned third input to the left by one there?

julia> bits(0b1110010110010011111000110001011000001001111001010 + 0b0000000000000000000000000110100111110110000110101)
"0000000000000001110010110010011111000110111111111111111111111111"
julia> bits(0b1110010110010011111000110001011000001001111001010 + (0b0000000000000000000000000110100111110110000110101 << 1))
"0000000000000001110010110010011111000111110100111110110000110100"

[wintersteiger]

Yeah, I added the leading zeros just for the comment here, and there's one too many now. There are carries, right around the comma where they make the difference.

111001011001001111100011000101100000100111100101 = A*B
000000000000000000000000110100111110110000110101 = aligned C
111001011001001111100011111010011111011000011010 = A*B+C

It's some sort of alignment problem (possibly due to exponent differences, etc). But, so far it doesn't make sense to me.

[ikirill]

I did it "by hand" (only using frexp and bits), moving the decimal point by +2 and -23, and I got no carries:

   11.100101100100111110001100010110000010011110010100000
 +  0.0000000000000000000000011010011111011000011010100000000000000000000000000000

julia> frexp(Float64(1.8660666f0) * Float64(1.9223062f0))
(0.8967878273778247,2)
julia> frexp(Float64(9.868419f-8))
(0.8278229832649231,-23)

[wintersteiger]

The exponent difference is -24, not -23, isn't it?

[ikirill]

If frexp says the number has exponent +2 or -23, there should be -2 or +23 zeros directly after point, shouldn't there? I wasn't doing it with an exponent difference, just directly numbers. I've pasted them into Wolfram Alpha now, to check again: https://www.wolframalpha.com/input/?i=binary+0.0000000000000000000000011010011111011000011010100001010011111011010100111101+in+decimal and it's not off by a factor of 2, so 23 zeros is right. (The other one: https://www.wolframalpha.com/input/?i=binary+11.100101100100111110001100010110000010011110010100000+in+decimal)

[ikirill]

Sorry, I typed in the number wrong, had to edit my comment, some of the trailing bits of the second are zeros now.

[wintersteiger]

I'm just looking at the input numbers: the first two have exponent 0 and the third has exponent -24. I'm still working out what that means for the exponent difference after the multiplication, it obviously isn't always right. Because the numbers are normalized to 1.xxx, I might well have to add/sub at least 1 to the exponent or the exponent difference.

[wintersteiger]

Fixed! All numbers were aligned properly, the problem was a rather subtle point about one of the sticky flags. All my tests and up until case 5 here pass now. Cases 6/7 are too far off to be a subtle issue, but they will have to wait until later.

[ikirill]

(Case 8) After 175f042db88bcbd8f2ce21316fb5c321355d601a everything above passes for me, but running Case 7 produces this counterexample. The main reason Case 7 is supposed to be unsat is because rnd(x*y) + fma(x,y,-rnd(x*y)) should equal x*y in exact arithmetic.

(set-option :pp.fp_real_literals true)
(set-logic QF_FP)
(set-info :status unsat)
(declare-const x Float32)
(declare-const y Float32)
(declare-const z Float32)
(declare-const r Float32)
(declare-const w Float32)
;; julia> fma(1.5001205f0, 1.5312883f0, -2.297117f0)
;; -4.381691f-8
(assert (= x (fp #b0 #b01111111 #b10000000000001111110011))) ; 1.5001205f0
(assert (= y (fp #b0 #b01111111 #b10001000000000101000001))) ; 1.5312883f0
(assert (= z (fp #b1 #b10000000 #b00100110000001111110111))) ; -2.297117f0
(assert (= r (fp #b1 #b01100110 #b01111000011000100110100))) ; -4.381691f-8
(assert (= w (fp.fma RNE x y z)))
(assert (not (= w r)))
(check-sat)
(get-model)

Explanation:

1.10000000000001111110011 * 1.10001000000000101000001:
julia> bits(UInt64(0b110000000000001111110011) * UInt64(0b110001000000000101000001))[17:end]
"100100110000001111110110110100001111001110110011"

    10.010011000000111111011011010000111100111011001100000
  - 10.0100110000001111110111
= -  0.000000000000000000000000101111000011000100110100
≠ -  0.000000000000000000000000101111000011000100111000

[wintersteiger]

That's fully expected, 175f042 fixed just one of many problems.

[wintersteiger]

As far as I can tell Z3 is right about the counter-example to case 7. The latest version finds these values:

(model
  (define-fun lo () Float32
    (fp #b1 #x66 #b10101101110011010001000))
  (define-fun hi () Float32
    (fp #b0 #x80 #b11011001101101110111011))
  (define-fun y () Float32
    (fp #b0 #x7f #b11111100010010000011101))
  (define-fun xy32 () Float32
    (fp #b0 #x80 #b11011001101101110111011))
  (define-fun x () Float32
    (fp #b0 #x7f #b11011101001011101000001))
  (define-fun xy64 () Float64
    (fp #b0 #b10000000000 #xd9b775948cba0))
  (define-fun xy () Float64
    (fp #b0 #b10000000000 #xd9b775948cbc0))
)

Those numbers are exactly what my hardware finds. The only thing I didn't check rigorously were the casts, but those are all from small to large.

In reals:

(model
  (define-fun lo () Float32
    ((_ to_fp 8 24) RTZ -1.67891025543212890625 -25))
  (define-fun hi () Float32
    ((_ to_fp 8 24) RTZ 1.85045564174652099609375 1))
  (define-fun y () Float32
    ((_ to_fp 8 24) RTZ 1.98547708988189697265625 0))
  (define-fun xy32 () Float32
    ((_ to_fp 8 24) RTZ 1.85045564174652099609375 1))
  (define-fun x () Float32
    ((_ to_fp 8 24) RTZ 1.86399090290069580078125 0))
  (define-fun xy64 () Float64
    ((_ to_fp 11 53) RTZ 1.85045561672880154446829692460596561431884765625 1))
  (define-fun xy () Float64
    ((_ to_fp 11 53) RTZ 1.8504556167288086498956545256078243255615234375 1))
)

[wintersteiger]

Nevermind, there's another precision issue to sort out first.

[ikirill]

There is a typo in https://github.com/wintersteiger/z3/commit/52cf80d63722ec6c322daa01c01fd7755306d9f8#diff-3a5e58244e798960480d3ee062492ba3R1670 causing an assertion failure, but after replacing 4 → 5, Z3 now fails to find a counterexample to Case 7 for me.