Re: Calculation discrepancy between Perl versions

FWIW: I get the same results on I686 Windows, as you get on I686 linux, which suggests to me that the difference if down to differences in the floating point hardware rather either the way Perl is built, or the underlying CRT math functions.

I thought for a while that it might be down to the IEEE rounding mode is use, but I tried it with all 4 modes and whilst the results do vary, the differences are far less than you are seeing:

C:\test>junk54 MaunaLoa 636 1
0.011256        -43.244877      41833.442623
Roots: 2093.870 or 1747.905

C:\test>junk54 MaunaLoa 636 2
0.011256        -43.245005      41833.967213
Roots: 2093.804 or 1747.982

C:\test>junk54 MaunaLoa 636 3
0.011256        -43.245005      41833.442623
Roots: 2093.939 or 1747.847

C:\test>junk54 MaunaLoa 636 4
0.011256        -43.245005      41833.967213
Roots: 2093.804 or 1747.982
[download]

If I add use bignum; I get the same result regardless of rounding mode (as you'd expect with infinite precision):

C:\test>junk54 MaunaLoa 636 0
0.011250        -43.220000      41809.395000
Roots: 2093.969 or 1747.809
[download]

But the result is still not far from the standard double precision results, and far away from your results on those other platforms.

So then I decided to "check the math".

I fed your sample data into wolfram/alpha's quadratic fit and got a = 0.01125 b=-43.33 c=41809.39500009608 (your results are close enough!)

I then walked manually through the calculation:

a = 0.01125 
b=-43.33 
c=41809.39500009608 - 636 = 41173.39500009608

b^2 = 1867.9684

4ac = 4 * 0.01125 * 41173.39500009608 = 1852.8027750043236 
[1852.8027750043236]

2a =  0.0225

sqrt( b^2 - 4ac) = sqrt( 1867.9684 - 1852.8027750043236 ) = 3.89430674
+13438813500081381515803
[3.8943067413438813500081381515802609939934472412114543853497515883707
+998992565436459214366842601252101413353364558519950923345771453315158
+871886529115406835628748856351213011743063789464283329454353106415694
+831827916183018713201714487448856249092857883643766930086114150302755
+812386977175989440678382944531781851160460991847426853996720283630521
+154852957729948705070366158701306031797881524733027828871220198125664
+64437007666673563847997186038400782988968755635538]

x' = -b + sqrt(b^2 - 4ac) = 43.22 + 3.8943067413438813500081381515803 
+= 47.11430674134388135000813815158
[47.114306741343881350008138151580260993993447241211454385349751588370
+799899256543645921436684260125210141335336455851995092334577145331515
+887188652911540683562874885635121301174306378946428332945435310641569
+4831827916183018713201714...]

y' = -b - sqrt(b^2 - 4ac) = 43.22 - 3.8943067413438813500081381515803 
+= 39.32569325865611864999186184842
[-39.32569325865611864999186184841973900600655275878854561465024841162
+920010074345635407856331573987478985866466354414800490766542285466848
+411281134708845931643712511436487869882569362105357166705456468935843
+0516817208381698128679..]

x = x' / 2a = 47.11430674134388135000813815158 / 0.0225 = 2093.9691885
+041725044448061400702
[2093.9691885041725044448061400702338219552643218316201949044334039275
+911066336241620409527415226722315618371260647045331152148700953480673
+727639401294018081583499949171165022744136168420634814642415693618475
+3258590184970230539200761...]

y = y' / 2a = 39.32569325865611864999186184842 / 0.0225 = 1747.8085892
+736052733329716377076
[1747.8085892736052733329716377075439558225134559461575828733443738501
+866711441536157368250362551055462159406517130732446625629076824297104
+050138376483759696194277828606612755033641609357142963135362084159302
+451918759280754723857...]
[download]

The small results are done using calc.exe, the extended numbers in brackets are produced using wolfram alpha.

All of which suggests that your I686 results are correct, and the "bug" is on the other platforms.

Perhaps they're using single precision? Or at least, perhaps lookup tables for sqrt() &| pow() that are only single precision?

Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.

"Science is about questioning the status quo. Questioning authority".

In the absence of evidence, opinion is indistinguishable from prejudice.

RIP an inspiration; A true Folk's Guy

Comment on Re: Calculation discrepancy between Perl versions Select or Download Code

Replies are listed 'Best First'.
Re^2: Calculation discrepancy between Perl versions by oko1 (Deacon) on Oct 05, 2010 at 16:17 UTC
Thanks, all, for the very useful responses (I'm responding to BrowserUk's reply specifically because it's so detailed and helpful in so many aspects) - this was very useful for both confirmation and more direction in deciding where to look for the error. I'm not all that familiar with the guts of Perl, but it's looking likely that there's a major compilation difference responsible here: either single-precision lookup tables, or perhaps a radically different math lib. Annoying that something like that could affect a Perl program... but on the other hand, good to learn that it can. Again, thank you all very much. -- "Language shapes the way we think, and determines what we can think about." -- B. L. Whorf	[reply]
Re^3: Calculation discrepancy between Perl versions by BrowserUk (Patriarch) on Oct 05, 2010 at 17:17 UTC
either single-precision lookup tables, or perhaps a radically different math lib. Annoying that something like that could affect a Perl program... I'm very unsure of my ground here as I'm not familiar with the other platforms, but it might not be software--Perl or the libs--but the floating point hardware. If their FPU's are only single precision, that might account for the results. That said, my best efforts to perform the calculation in single precision doesn't get close to the inaccuracies you;re seeing: #include <stdio.h> float sqrtf( float n ) { float g = n / (float)2.0; while( ( ( gg ) - n ) > 0.000001 ) { g = ( g + n/g ) / (float)2.0; } return g; } void main( void ) { float a = (float)0.01125, b = (float)-43.22, c = (float)41173.3950 +0009608; float two = 2.0, four = 4.0; float b2 = b b; float ac4 = four * a * c; float sqrtbit = sqrtf( b2 - ac4 ); float a2 = two * a; float r1 = ( -b + sqrtbit ) / a2; float r2 = ( -b - sqrtbit ) / a2; printf( "b2:%f ac4:%f sqrt:%f a2:%f\n", b2, ac4, sqrtbit, a2 ); printf( "roots: %f , %f\n", r1, r2 ); } [download] Produces: `C:\test>quads b2:1867.968506 ac4:1852.802856 sqrt:3.894310 a2:0.022500 roots: 2093.969238 , 1747.808472` [download] Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error. "Science is about questioning the status quo. Questioning authority". In the absence of evidence, opinion is indistinguishable from prejudice. RIP an inspiration; A true Folk's Guy	[reply] [d/l] [select]