Getting odd floating point results.

• Please include the declarations of lat1, lat2, long1, long2, and delta.
• How is BETA defined?
• Please provide example values for the variables.
• Given those values, what results do you expect?
• Given those values, what results do you get?

#define BETA  6356752.0

float procedure( float lat1, float long1, float lat2, float long2 ) {
        float delta;
        delta = ALPHA * cos(lat1);
        return ((lat2-lat1)*BETA*(lat2-lat1)*BETA + (long2-long1)*delta*(long2-long1)*delta);
}

lat1 = 0.727414575
long1 = -1.24152067
lat2 = 0.727336559
long2 = -1.24135643

I expect this to return approximately 737991.335
and instead it returns 6356752.0 (same as BETA).

What is ALPHA?

Sorry, #define ALPHA 6378137.0

Well, you can track the process of computation. But it involves single-stepping in disassembler, watching intermediate results in registers and converting them from binary to human-readable representation.

"... instead it returns 6356752.0 (same as BETA)."

For me, it does not return that value or the value you expect, but rather 857998.3.

I suggest you break the long return expression into small intermediate expressions with values assigned to separate float variables so you can see the intermediate results.

I'd have to look into it further, but I wonder if the large difference in magnitude between intermediate values combined with the limited resolution of a 32-bit float is causing problems.

I just double checked my hand calculations... the result you got is correct, please disregard the one I initially supplied.

It still does not work for me, I can only assume some other condition is impacting the proper execution or compiling of this code.

FWIW, here's the code I used to test:

#include <reg52.h>
#include <math.h>
#include <intrins.h>

#define ALPHA 6378137.0
#define BETA  6356752.0

float procedure( float lat1, float long1, float lat2, float long2 ) {
    float delta;
    delta = ALPHA * cos(lat1);
    return ((lat2-lat1)*BETA*(lat2-lat1)*BETA + (long2-long1)*delta*(long2-long1)*delta);
}

void main(void)
{
    float f;

    f = procedure(0.727414575, -1.24152067, 0.727336559, -1.24135643);

    for (;;) {
        _nop_();
    }
}

I would be very careful about cancellation errors with just single-point precision.

But most important - what are you trying to compute? You are squaring the difference between two coordinates and multiplying with a constant. Can your full formula be rewritten to reduce the problem with cancellations?

And what are your constants? Are they of the correct magnitude? I think you have an error in a unit somewhere. Do you have any web reference to the formula you are trying to compute?

I would be very careful about cancellation errors with just single-point precision.

Seconded. We're looking at an attempt to turn 10^{-4} relative differences into results, using a number format that's only capable of handling 10^{-7} on a good day. That mean's you're down to three credible decimals before the main computation even starts.

The fact that the same formula, evaluated in double-precision on a PC, yields a result of 858106, shows just how deep into cancellation this calculation is.

But most important - what are you trying to compute?

That part is pretty obvious once you look at the names and numbers. The "lat"s and "long"s are geographical latitudes and longitudes in radians. The planet these are used as coordinates on is Earth (ALPHA and BETA are radii of the Earth in meters, in the directions to the equator and the poles respectively). delta is the radius of a circle of latitude = lat1.

The whole thing is an attempt at computing the squared surface distance between two positions on an (idealized, ellipsoidal) Earth, by doing a Pythagoras on the arc lengths corresponding to the differences in longitude and latitude, respectively.

The attempt is flawed because spherical trigonometry is nowhere near as simple as this. This shows in several defects, the most obvious of which is a bias: lat1 is preferred over lat2 to compute delta, for no good reason. It gets a lot worse if you consider what will happen if you use it across the international date line. Or over any sizeable distance at all, for that matter.

In other words, this formula is doomed by cancellation error on short distances, and by conceptual fallacy on large ones.

Hans-Bernhard is in fact correct...

If you take the square root of the result of the provided code and you should get an approximate distance between two geographical coordinates.

The intended application for this code is small distances, as a more accurate computation will not work with single precision floating point numbers.

The issue of cancellation error strikes me as somewhat odd. I am porting code from SDCC over to keil and was getting good results while running SDCC. I'm not an expert on compilers by any means, but perhaps the C51 compiler is take some assumptions in the way that it computes results that make it more vulnerable to this type of error. Without knowing the exact differences in methods between the two compilers I can't say for sure.

The real question then is how can I force the compiler to execute each operation in an order that minimizes the relative differences?

But the big thing here is of course that earth is not flat, so pythagoras is no fun to use except for very small areas - unless the goal is to bild a straight tunnel (which would feel like a downhill slope followed by an uphill slope if really long).

The wikipedia entry for evaluating the great-circle distance is:
en.wikipedia.org/.../Great-circle_distance

But the simplification of computing using an averaged earth radius instead of caring about earth being oblongated will still give _way_ better results than the naive pythagoras.

What quality that is needed comes back to the original question: what is the evaluation intended to be used for?

Note that a PC using single-precision arithmetic will not compute in the same way that the Keil code does. When Keil multiply two float values they do it in single precision. A PC will have extended precision for intermediate operations. Most PC compilers then have compilation flags to control if the generated code should follow the IEEE regulations for storing back intermediate results into limited-sized variables or if the compiler should keep the intermediate values in the floating-point registers (which normally are 80-bit large).

But in your case, you will get very large numbers because you do not divide with 2*pi before you square.

earth is not flat, so pythagoras is no fun

It's not quite that bad. As I said, the formula uses Pythagoras on arc lengths instead of straight lines (secants). I.e. roads, not tunnels.

But yes, a less incorrect function would have to go from arcs to secants, use Pythagoras, then convert the resulting secant back into a surface arc.

Without knowing the exact differences in mehods between the two compilers I can't say for sure.

You can't really know "the methods" as such. What you can know is what code the compilers actually generated for this task. It will take quite some guessing of what vendor library function does what, and possibly you'll have to decode intermediate results into human-readable numbers. In the end you'll probably know more than you ever really wanted to, but it can be done.

The real question then is how can I force the compiler to execute each operation in an order that minimizes the relative differences?

Break up the computation into individual steps. Simplify it (computing stuff twice is wasteful). Turn down optimization, optionally by using "volatile" qualifiers. Print out intermediate results for inspection.