Page 1 of 1
Same Calculation, different results
Posted: Wed Apr 22, 2020 3:10 pm
by CalamityJames
I was interested in the PureBasic routines for the Haversine (a calculation of the distances between two points on a curved surface) found here:
https://www.purebasic.fr/english/viewtopic.php?t=47853. This question is not about the formula but about two ways of calculating it. The code below shows two functionally identical procedures, the first is the way any sane person would write it (and which is taken from the link), and the second which avoids any intermediate variables. They give slightly different results. The difference is only 145m in 540km so it's not a big deal. I also know that this does not indicate any fault with PureBasic, but is to do with the number of bits of accuracy used. I would appreciate any additional clarification anyone can offer. (The places in the formula are the The Royal Observatory in Greenwich and the Royal Mile in Edinburgh.)
Code: Select all
EnableExplicit
#Earth_Radius = 6371.000785
Procedure.d Haversine(PlanetRadius.d, Lat1.d, Lon1.d, Lat2.d, Lon2.d)
Protected Sin_dLat.d, Sin_dLon.d, a.d
Sin_dLat = Sin(Radian(Lat2 - Lat1) / 2)
Sin_dLon = Sin(Radian(Lon2 - Lon1) / 2)
a = Sin_dLat * Sin_dLat + Sin_dLon * Sin_dLon * Cos(Radian(Lat1)) * Cos(Radian(Lat2))
ProcedureReturn PlanetRadius * 2 * ATan2(Sqr(1 - a), Sqr(a))
EndProcedure
Procedure.d Haversine2(PlanetRadius.d, Lat1.d, Lon1.d, Lat2.d, Lon2.d)
ProcedureReturn ATan2(Sqr(1 - Sin(Radian(Lat2 - Lat1) / 2) * Sin(Radian(Lat2 - Lat1) / 2) + Sin(Radian(Lon2 - Lon1) / 2) * Sin(Radian(Lon2 - Lon1) / 2) * Cos(Radian(Lat1)) * Cos(Radian(Lat2))), Sqr(Sin(Radian(Lat2 - Lat1) / 2) * Sin(Radian(Lat2 - Lat1) / 2) + Sin(Radian(Lon2 - Lon1) / 2) * Sin(Radian(Lon2 - Lon1) / 2) * Cos(Radian(Lat1)) * Cos(Radian(Lat2)))) * PlanetRadius * 2
EndProcedure
Debug Haversine(#Earth_Radius, 51.476852, -0.000500, 55.950558, -3.185556)
Debug Haversine2(#Earth_Radius, 51.476852, -0.000500, 55.950558, -3.185556)
Debug "Difference = " + StrD(Haversine(#Earth_Radius, 51.476852, -0.000500, 55.950558, -3.185556) - Haversine2(#Earth_Radius, 51.476852, -0.000500, 55.950558, -3.185556), 17)
Re: Same Calculation, different results
Posted: Wed Apr 22, 2020 4:49 pm
by Josh
I haven't checked to see if your two calculations are actually identical. I don't want to say that PureBasic can't do math, but it definitely has certain peculiarities where you never really know what the result will be. With only two results for one and the same calculation you are still one of the lucky ones
Re: Same Calculation, different results
Posted: Wed Apr 22, 2020 5:51 pm
by PeDe
I think the order of the calculations is changing. If you put brackets, the result is identical.
Code: Select all
Procedure.d Haversine2(PlanetRadius.d, Lat1.d, Lon1.d, Lat2.d, Lon2.d)
ProcedureReturn ATan2(Sqr(1 - ((Sin(Radian(Lat2 - Lat1) / 2)) * (Sin(Radian(Lat2 - Lat1) / 2)) + (Sin(Radian(Lon2 - Lon1) / 2)) * (Sin(Radian(Lon2 - Lon1) / 2)) * Cos(Radian(Lat1)) * Cos(Radian(Lat2)))), Sqr((Sin(Radian(Lat2 - Lat1) / 2)) * (Sin(Radian(Lat2 - Lat1) / 2)) + (Sin(Radian(Lon2 - Lon1) / 2)) * (Sin(Radian(Lon2 - Lon1) / 2)) * Cos(Radian(Lat1)) * Cos(Radian(Lat2)))) * PlanetRadius * 2
EndProcedure
[18:51:02] Waiting for executable to start...
[18:51:02] Executable type: Windows - x86 (32bit, Unicode, Thread, Purifier)
[18:51:02] Executable started.
[18:51:02] [Debug] 539.66501330444737
[18:51:02] [Debug] 539.66501330444737
[18:51:02] [Debug] Difference = 0.00000000000000000
[18:51:02] The Program execution has finished.
Re: Same Calculation, different results
Posted: Wed Apr 22, 2020 5:59 pm
by skywalk
It would be better if you compared your floating point results with another app?
It is not clear whether some of the implied types in the shorter procedure default to FLOAT instead of DOUBLE.
That is why we need a known comparison to make conclusions.
Re: Same Calculation, different results
Posted: Wed Apr 22, 2020 6:13 pm
by wilbert
For most floating point operations, PureBasic uses the fpu which is 80 bits internally.
Every time an intermediate result needs to be stored, it is stored as a 64 bit (double) or 32 bit (float) value so accuracy is lost.
The most accurate calculation for this formula would be to use assembly code and only store the final result.
Re: Same Calculation, different results
Posted: Wed Apr 22, 2020 6:41 pm
by skywalk
Yes, I am curious why 32-bit would still be used for an all DOUBLE calculation?
I thought Fred updated all the math/trig functions to DOUBLE?
Re: Same Calculation, different results
Posted: Wed Apr 22, 2020 7:20 pm
by Little John
wilbert wrote:For most floating point operations, PureBasic uses the fpu which is 80 bits internally.
Every time an intermediate result needs to be stored, it is stored as a 64 bit (double) or 32 bit (float) value so accuracy is lost.
It's a pity that PureBasic does not have 80 bit numeric variables (e.g. PowerBasic has them since many years).
Re: Same Calculation, different results
Posted: Thu Apr 23, 2020 6:28 am
by Helle
PureBasic don´t use FPU-instructions, it uses the msvcrt.dll (routines without FPU).
Re: Same Calculation, different results
Posted: Fri Apr 24, 2020 11:58 am
by CalamityJames
Thanks to those who replied, particulary to PeDe whose observation skills are clearly above average!. The fact that the order of calculation made a difference made me wonder which of the results was closer to the right answer so I used the Casio calculator here:
https://keisan.casio.com/calculator. This offers up to 130 decimal places. On my orginal post the Haversine2 function is closer than the Haversine function, a comfort to insane programmers everywhere.
Re: Same Calculation, different results
Posted: Fri Apr 24, 2020 4:30 pm
by vwidmer
Re: Same Calculation, different results
Posted: Fri Apr 24, 2020 6:24 pm
by Michael Vogel
Just an idea how errors could be avoided when transforming Haversine1 to Haversine2...
...let's start with some macros...
Code: Select all
Macro _slat_
Sin(Radian(Lat2 - Lat1) / 2)
EndMacro
Macro _slon_
Sin(Radian(Lon2 - Lon1) / 2)
EndMacro
Macro _a_
(_slat_*_slat_ + _slon_* _slon_ * Cos(Radian(Lat1)) * Cos(Radian(Lat2)))
EndMacro
Procedure.d Haversine(PlanetRadius.d, Lat1.d, Lon1.d, Lat2.d, Lon2.d)
ProcedureReturn PlanetRadius * 2 * ATan2(Sqr(1 - _a_), Sqr(_a_))
EndProcedure
Next step is to add an error by adding an asterisk for instance:
Code: Select all
Procedure.d Haversine(PlanetRadius.d, Lat1.d, Lon1.d, Lat2.d, Lon2.d)
ProcedureReturn PlanetRadius * 2 * ATan2(Sqr(1 - _a_), Sqr(_a_)) *
EndProcedure
When you start the compiler now, you'll get an error message with the resulting code which can be easily copied into the Haversine2 routine (without the asterisk)...
Code: Select all
EnableExplicit
#Earth_Radius = 6371.000785
Macro _slat_
Sin(Radian(Lat2 - Lat1) / 2)
EndMacro
Macro _slon_
Sin(Radian(Lon2 - Lon1) / 2)
EndMacro
Macro _a_
(_slat_*_slat_ + _slon_* _slon_ * Cos(Radian(Lat1)) * Cos(Radian(Lat2)))
EndMacro
Procedure.d Haversine(PlanetRadius.d, Lat1.d, Lon1.d, Lat2.d, Lon2.d)
ProcedureReturn PlanetRadius * 2 * ATan2(Sqr(1 - _a_), Sqr(_a_))
EndProcedure
Procedure.d Haversine2(PlanetRadius.d, Lat1.d, Lon1.d, Lat2.d, Lon2.d)
ProcedureReturn PlanetRadius * 2 * ATan2(Sqr(1 - (Sin(Radian(Lat2 - Lat1) / 2)*Sin(Radian(Lat2 - Lat1) / 2) + Sin(Radian(Lon2 - Lon1) / 2)* Sin(Radian(Lon2 - Lon1) / 2) * Cos(Radian(Lat1)) * Cos(Radian(Lat2)))), Sqr((Sin(Radian(Lat2 - Lat1) / 2)*Sin(Radian(Lat2 - Lat1) / 2) + Sin(Radian(Lon2 - Lon1) / 2)* Sin(Radian(Lon2 - Lon1) / 2) * Cos(Radian(Lat1)) * Cos(Radian(Lat2)))))
EndProcedure
Debug Haversine(#Earth_Radius, 51.476852, -0.000500, 55.950558, -3.185556)
Debug Haversine2(#Earth_Radius, 51.476852, -0.000500, 55.950558, -3.185556)
Debug "Difference = " + StrD(Haversine(#Earth_Radius, 51.476852, -0.000500, 55.950558, -3.185556) - Haversine2(#Earth_Radius, 51.476852, -0.000500, 55.950558, -3.185556), 17)