Question about 8000 sampling points

[zuoyanzhang]

I passed herbie test ((35000000.0 + ((0.401 * (1000.0 / x)) * (1000.0 / x))) * (x - (1000.0 * 4.27e-5))) - ((1.3806503e-23 * 1000.0) * 300.0) expression, the maximum error given is 6.6bits, the test interval is [0.01, 100000], the corresponding input x to find this point is 13809681318699300000, I speculate that x here should be decimal, convert it to hexadecimal The hexadecimal number is 0xbfa5d7fe166b74a0. When converting the hexadecimal number to the corresponding floating-point number, it is -0.042663517224107. Why is the corresponding negative number here? My test range is obviously a positive number. I would like to know how each point of herbie's 8000 sampling points is converted into a corresponding floating point number, or whether herbie can give the input corresponding to the maximum error in a given interval. Maybe I didn't understand what 13809681318699300000 is here, I hope to get your reply, thank you

[elmisback]

That calculation looks correct to me. I wrote the code that's reporting those values, but not the code that does the actual translation from the internal sample point representation to the ordinals you're seeing. I'll have a look at the data from that expression to see if I can figure out what's going on, probably ask someone more knowledgeable, and get back to you.

Regarding your question about the maximum error input, I doubt Herbie will necessarily find that input point, since it's using a uniform sample.

[elmisback]

Sorry, I just forgot how these are laid out. Basically, these values are used as a positive, monotonically increasing index over the floats. This is different from what you get if you just directly convert from a signed integer.

Checking the reported values for the ticks which are used to lay out the plot for a simple function that covers the whole range like x - 42, we get these values:

"ticks_by_varidx":[[["-2e+308",5442597155835700],["-e+206",1534620040627248120],["-e+103",3075435233412954392],["-1",4616189618054758400],["-e-102",6142370034567588547],["-e-205",7683126141587503524],["0",9223372036854775808],["e-206",10748314587012925997],["e-103",12289385714708364181],["1",13830554455654793216],["e+102",15356310588812194163],["e+205",16897476253732898419],["e+308",18440533523121328391]]]

These values are approximately matched by a javascript function like this. (The part you'd care about is the last couple lines):

function herbie_ordinal_to_real(ordinal) {
  // (approximate) -- using mathjs
  function real_from_signed_int (signed) {
    const buffer = new ArrayBuffer(8)
    const view = new DataView(buffer)
    view.setBigInt64(0, signed)
    return view.getFloat64(0)
  }
  function to_signed_int (float64) {
    const buffer = new ArrayBuffer(8)
    const view = new DataView(buffer)
    view.setFloat64(0, float64)
    return view.getBigInt64(0)
  }
  
  const straight_conversion = real_from_signed_int(BigInt(ordinal))  // conversion you were doing
  const zero = math.bignumber(to_signed_int(-0.0).toString())
  // Flip the sign, then the negative half of the range is also reflected
  return - straight_conversion < 0 ? - real_from_signed_int(math.subtract(zero, math.bignumber(ordinal))) : - straight_conversion
}

[zuoyanzhang]

Thank you very much for your patient answer, first of all the herbie tool is doing very well.
Through the conversion code you gave, I already know how the conversion is done internally. Let me talk about a problem I encountered, for the carbonGas expression ((35000000.0 + ((0.401 * (1000.0 / x)) * (1000.0 / x))) * (x - (1000.0 * 4.27e-5))) - ((1.3806503e-23 * 1000.0) * 300.0), the interval is set to [0.01,10], the maximum error obtained by using the herbie tool is 8.6bits, and the corresponding sampling point is 13809684655913878000, which is converted to float by the conversion code you gave The number of points is 0.04268667380199531, the ULP error I got by comparing it with the mpfr high-precision library is 7.063657e-01, and the bits error is 0.0.
In the same way, I verified the (exp(x)-1)/x expression, and the result obtained by the maximum error point corresponding to this expression is correct. So I don't know if herbie has a bug with the carbonGas expression or it's my problem.

[elmisback]

Maybe @pavpanchekha might have more insight into this question, since I'm not sure how to do the MPFR comparison myself. One possibility I can think of is that the point being stored in the JSON might be slightly different from the actually sampled and tested point due to some kind of truncation when storing it as a number rather than a string. I didn't really worry about this issue when writing the code since I assumed the points were just going to be used for plotting.

For Pavel--the point in question is at index 4907 in the arrays of the points.json of this run https://herbie.uwplse.org/demo/8d8d99731736e6050c5de929d6627859c88ae00a.1.6/graph.html.

[pavpanchekha]

Hi all, I'll definitely take a look when I have a moment.

[pavpanchekha]

Ok, so here's how Herbie is getting the answer it's getting. This is a transcription of a live REPL session launched in src/points.rkt in Herbie's codebase. I'm doing that in case you want to chase down how various things are defined:

$ cd src/
$ racket -i

; Load stuff we'll use and 
> (require "load-plugin.rkt" "programs.rkt" "float.rkt" math/bigfloat math/flonum)
> (define ctx (make-debug-context '(x)))

; Define program and input point
> (define prog '(lambda (x) (-.f64 (*.f64 (+.f64 35000000.0 (*.f64 (*.f64 0.401 (/.f64 1000.0 x)) (/.f64 1000.0 x))) (-.f64 x (*.f64 1000.0 4.27e-5))) (*.f64 (*.f64 1.3806503e-23 1000.0) 300.0))))
> (define x 0.04268667380199531)

; Compile the program in float and bigfloat modes
> (define fn.fl (eval-prog prog 'fl ctx))
> (define fn.bf (eval-prog prog 'bf ctx))

; Evaluate the program in float and bigfloat modes
> (define y.fl (fn.fl x))
> y.fl
-3399.105886487708
; Use 10,000 bits to make sure rounding isn't an issue
> (define y.bf (bigfloat->flonum (parameterize ([bf-precision 10000]) (fn.bf x))))
> y.bf
-3399.105886487528

; Compute error between float and bigfloat modes
> (flonums-between y.fl y.bf)
396
> (ulps->bits (+ 1 (abs (flonums-between y.fl y.bf))))
8.632995197142959

I know you said you also used MPFR to confirm the correct answer, and got something different, and I'm not sure why or how. Can you explain further?

[zuoyanzhang]

Good to know how the given point error of the herbie tool is calculated. Sorry, I double-checked and it was my mpfr program that was written wrong. Thank you for your patience in answering.

[pavpanchekha]

Always happy to have such diligent users!

[zuoyanzhang]

@pavpanchekha Dear Professor Pavel Panchekha, I have a question to confirm. The specific experimental process is as follows: For the benchmark NMSEproblem331, the expression is (1.0 / (x + 1.0)) - (1.0 / x), the test interval I chose is [0.01, 100000], and the corresponding fpcore program is:
(FPCore (x)
:name "1.0/(x+1) - (1.0/x)"
:pre (and (<= 0.01 x 100000))
(- (/ 1.0 (+ x 1)) (/ 1.0 x)))
I execute "herbie report ./NMSEproblem331.fpcore ./xxx" in the terminal. In the generated points.json file, the maximum error reported is 16.6bits. I found three inputs that caused 16.6bits, which are [13904968337258583276], [13904869120057830665], [13904670248646492365], converted to floating point numbers are 99824.49691209197, 98380.6966175437, 95486.73669838905 respectively. Calculate the error of these three points respectively, the results are 14.962bits, 14.914bits, 14.996bits, and the three points reported by herbie are 16.6bits. Where did I go wrong?