Same Calculation, different results

[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)

[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

[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.

[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.

[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.

[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?

[Little John]

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).

[Helle]

PureBasic don´t use FPU-instructions, it uses the msvcrt.dll (routines without FPU).

[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.

[vwidmer]

not sure if you seen this post already

http://www.purebasic.fr/english/viewtop ... 13&t=47853

[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)