'Nuther Challenge to the Code Monkeys

[Lester Thorpe]

Here is another stumper for the dumb fuck code monkeys.
Watch as they rattle their cage -- and scream for Python --
just like the previous time.

The following C program computes the value of an expression
using three different parameters. The expression is:

33375/100*y^6+x^2*(11*x^2*y^2-y^6-121*y^4-2)+55/10*y^8+x/(2*y)

The exact answer is:

-54767/66192 = -0.8273960599468214

The C program gives these results:

# ./cancellation1
Float: inf
Double: 5.764608e+17
Long Double: 5.927442e+39

# ./cancellation2
Float: inf
Double: 1.172604e+00
Long Double: 5.927442e+39

# ./cancellation3
Float: inf
Double: 1.172604e+00
Long Double: 5.927442e+39

Holy double fucking moley! What the fuck is going on here!

Code monkey want banana? Code monkey solve problem.

Otherwise, code monkey will be gassed.

Ahahahahahahahahahahahahaha!


=============================
Begin C Program
=============================

// Compile three ways with:
// gcc -O0 -march=native -mfpmath=387 -fexcess-precision=standard -o cancellation1 cancellation.c
// gcc -O0 -march=native -mfpmath=sse -o cancellation2 cancellation.c
// gcc -O0 -march=native -mfpmath=sse -ffloat-store -o cancellation3 cancellation.c
//
// The expression to compute:
// 333.75 * y^6 + x^2 * (11 * x^2 * y^2 − y^6 − 121 * y^4 − 2) + 5.5 * y^8 + x/2y
//
// Let x=77617, y=3096

#include <stdio.h>

int main()
{
float xf, yf, fxyf;
double x,y, fxy;
long double xl, yl, fxyl;
x=77617.0; xf=(float)x; xl=(long double)x;
y=33096.0; yf=(float)x; yl=(long double)x;

fxyf=333.75*yf*yf*yf*yf*yf*yf+5.5*yf*yf*yf*yf*yf*yf*yf*yf+xf*xf*(11.0*xf*xf*yf*yf-yf*yf*yf*yf*yf*yf-121.0*yf*yf*yf*yf-2.0)+xf/(2.0*yf);
printf("%-14s %e\n", "Float:", fxyf);

fxy=333.75*y*y*y*y*y*y+5.5*y*y*y*y*y*y*y*y+x*x*(11.0*x*x*y*y-y*y*y*y*y*y-121.0*y*y*y*y-2.0)+x/(2.0*y);
printf("%-14s %e\n", "Double:", fxy);

fxyl=333.75*yl*yl*yl*yl*yl*yl+5.5*yl*yl*yl*yl*yl*yl*yl*yl+xl*xl*(11.0*xl*xl*yl*yl-yl*yl*yl*yl*yl*yl-121.0*yl*yl*yl*yl-2.0)+xl/(2.0*yl);
printf("%-14s %Le\n", "Long Double:", fxyl);

return(0);
}

=============================
End C Program
=============================

[Physfitfreak]

Only the case of double float works for this problem. In the long double
and float you've taken y and x equal. I don't know why.

Correct it and repost, then we can comment.



--
This email has been checked for viruses by Avast antivirus software.
www.avast.com

[Lester Thorpe]

On Sun, 18 Feb 2024 16:07:35 -0600, Physfitfreak wrote:

>
> Only the case of double float works for this problem. In the long double
> and float you've taken y and x equal. I don't know why.
>
> Correct it and repost, then we can comment.
>

I do apologize for my typographical error. The corrected section is:

x=77617.0; xf=(float)x; xl=(long double)x;
y=33096.0; yf=(float)y; yl=(long double)y;

The corrected results are:

[~]# ./cancellation1
Float: 5.764608e+17
Double: 5.764608e+17
Long Double: 5.764608e+17
[~]# ./cancellation2
Float: 7.089312e+29
Double: 1.172604e+00
Long Double: 5.764608e+17
[~]# ./cancellation3
Float: 7.089312e+29
Double: 1.172604e+00
Long Double: 5.764608e+17

The correct answer is still -0.8273960599468214.

Try this with gcc. You will get the same results.

Note: I chose to do exponentiation, e.g. y^6, by using
repeated multiplication because it is possibly better
than pow(y, 6).

[Physfitfreak]

What do you mean by "cancellation1" and "cancellation2" and
"cancellation3"?

[Farley Flud]

On Mon, 19 Feb 2024 21:57:09 -0600, Physfitfreak wrote:

>
>
> What do you mean by "cancellation1" and "cancellation2" and
> "cancellation3"?
>

It's a hint.

In this example the old nemesis of "catastrophic cancellation"
bites again.

Using the bc calculator (Linux only), we get 55/10*y^8 =

7917111340668961361101134701524942848

The other part, 33375/100*y^6+x^2*(11*x^2*y^2-y^6-121*y^4-2) =

-7917111340668961361101134701524942850

Adding gives:

-2

But with all floating point formats adding gives:

0

The equation now becomes:

-2 + x/(2*y)

But with all FP formats we have:

0 + x/(2*y)


The moral of the story is that the programmer has to be
damned careful or else his results might be totally FUBAR.

Unless one wants to become a skilled numerical analyst,
which takes years, then one should use proven math packages.

[DFS]

On 2/20/2024 4:39 AM, Lameass Larry wrote:

> Using the bc calculator (Linux only)

https://github.com/gavinhoward/bc
https://gnuwin32.sourceforge.net/packages/bc.htm

[rbowman]

On Tue, 20 Feb 2024 09:39:42 +0000, Farley Flud wrote:

> Unless one wants to become a skilled numerical analyst, which takes
> years, then one should use proven math packages.

That's something we can agree on. When I'm doing a reprojection from SPCS
to WGS84 going down the rabbit hole isn't that productive. Even when
getting the azimuth and distance between two points I'm happy the use -lm.


int
get_azimuth_and_distance(double lat1, double lon1, double lat2, double
lon2,

double* azimuth, double* distance)
{
double rlat1, rlon1;
double rlat2, rlon2;
double pi;
double a;
double d;

pi = 4.0 * atan(1.0);
rlat1 = (lat1*pi)/180.0;
rlon1 = (lon1*pi)/180.0;
rlat2 = (lat2*pi)/180.0;
rlon2 = (lon2*pi)/180.0;

d = acos(sin(rlat1)*sin(rlat2)+cos(rlat1)*cos(rlat2)*cos(rlon1-
rlon2));
if (isnan(d)) {
return -1;
}
if (sin(rlon2-rlon1) < 0) {
a =
acos((sin(rlat2)-sin(rlat1)*cos(d))/
(sin(d)*cos(rlat1)));
}
else {
a =
2*pi-acos((sin(rlat2)-sin(rlat1)*cos(d))/
(sin(d)*cos(rlat1)));
}
if (isnan(a)) {
return -1;
}

*azimuth = (a*(180.0/pi));
*distance = ((d*180.0*60)/pi) * (6076.11549/5280.0);
return 0;
}

[Physfitfreak]

Try this same problem with hoc, using mathcw :) I talked about it with
Bowman last night. Both programming software and the library are free to
download from Nelson Beebe's online site. hoc is very much like c so
there should be no learning curve.

[chrisv]

I don't know if that needed to be tested, or if a clean compile was
good enough to know that it will work correctly.

According to some shameless liar, some dumb fscking snit.

--
'[chrisv] literally said it was "stupid" to test the code vs relying
on compiler warnings.' - DumFSck, lying shamelessly

[%]

[rbowman]

On Tue, 20 Feb 2024 15:03:07 -0600, chrisv wrote:

> I don't know if that needed to be tested, or if a clean compile was good
> enough to know that it will work correctly.

I did some testing but it wasn't a rocket science grade problem. In the
best Physfitfreak mode,

You are a dispatcher in a 911 call center and you know the coordinates of
your desk are 38.246915, -104.501669. A caller reports a car wreck and the
ali feed from his cellphone shows the location at 38.169221, -103.062460.
How far away is he and in what direction?

It's a quick and dirty test of whether it's your problem. That's an
ongoing problem. About 10 years ago there was a big hullabaloo when a 911
call from the west side of the US Virgin Islands wound up in San Juan.
Leaving the language difficulties aside obviously there was nothing SJ
could do but frantically try to contact someone in the VI.

That's a dramatic case. Usually it's only the next county. There are mode
precise algorithms but it goes back to Engineering 101 -- how good is good
enough? Few consumer grade GPS receivers are more accurate than 9 meters.
Given that you're fooling yourself with 6 decimal points.

[Physfitfreak]

On 2/20/2024 12:03 PM, rbowman wrote:
> d = acos(sin(rlat1)*sin(rlat2)+cos(rlat1)*cos(rlat2)*cos(rlon1-
> rlon2));


Two questions.

1. Did you derive the above relation by yourself?

2. So you just threw the above expression at the computer so the
compiler would figure it out what to do with it?

I'm asking the second question because if you just took the relation
from some calculus book or somewhere on the web that had used calculus
to come to that relation, that doesn't mean you can then just use it
like that in your program :)

Here is a famous quote by two of the experts (I just came across it
yesterday in Beebe's book):

"MANY PITFALLS AWAIT A SYSTEMS PROGRAMMER WHO ATTEMPTS TO IMPLEMENT
BASIC FUNCTION ROUTINES USING INFORMATION GLEANED FROM A CALCULUS TEXT.

— W. J. CODY, JR. AND W. WAITE (1980)"


Hehe :) Beebe begins his "Introduction" section of his book with that
quote.

And I asked the first question because it can be a nice challenge
problem here to solve and derive the answer, and then some more
challenge to code it well for computation.

I did a somewhat similar problem, way back, for a friend who wanted to
know toward which direction one should perform the daily prayers, when
one is at a particular spot on Earth. The direction of prayer must
always be along the shortest distance to Mecca.

I solved it for any spot on earth and wrote a FORTRAN code for it so he
could use from then on, no matter in which city he was. Then I had to
convert the code to BASIC cause he didn't have a FORTRAN compiler on his
pc. That was the first time I learned BASIC. All he needed to do was to
give it the longitude and latitude of where he was. It would give the
angle, clockwise beginning from strict north, along which prayer would
be "accepted" by god. Hehe :)

It was "fun".

In Iran, the direction of prayer is already included in the design of
the house one is living. The houses there are always facing the shortest
distance to Mecca. So no calculations for those mortals is needed. But
in the USA it can get pretty tricky.

[Physfitfreak]

On 2/20/2024 4:02 PM, rbowman wrote:
> You are a dispatcher in a 911 call center and you know the coordinates of
> your desk are 38.246915, -104.501669. A caller reports a car wreck and the
> ali feed from his cellphone shows the location at 38.169221, -103.062460.
> How far away is he and in what direction?

If Earth is taken as a sphere, it is easy and fast calculation in polar
coordinates. If it's taken as an ellipsoid, then the radial distance to
the point on the path, itself, becomes a variable. So the easy solution
for the sphere case and the method to calculate it fast won't apply.

For Ellipsoid, I'd just use the equation of ellipsoid in polar
coordinates as a given, then would do the line integral from one
coordinate on it to the other coordinate. If not in college algebra
text, then definitely in the calculus text, material to know how to do
that line integral can be found.

It's a very fun problem :)

[Physfitfreak]

And before anything, I'd check the equation of ellipsoid in Cartesian
coordinates as well to see whether it yields easier integrations over a
path on its surface. Then will choose whether to use Cartesian or polar
coordinates for this problem.

I give it a good chance that the integrations necessary to perform will
not be easy enough to be done analytically, so appropriate numerical
methods for computing the integrations should be employed to compute
that distance.

But the formula you wrote in your earlier post for the distance looked
like somebody had performed the integrations analytically, so this makes
me suspicious that they may have used a sphere to model the Earth.

[Physfitfreak]

On 2/20/2024 3:39 AM, Farley Flud wrote:

Also before I begin working on this you need to use more brackets () to
clarify what expression you have in mind:

33375/100*y^6+x^2*(11*x^2*y^2-y^6-121*y^4-2)+55/10*y^8+x/(2*y)

In the above expression it is not clear you mean to start it as
(33375/100)*y^6 + ... or you mean 33375/(100*y^6) + ...

Same ambiguity about the term 55/10*y^8.

[rbowman]

On Tue, 20 Feb 2024 18:19:53 -0600, Physfitfreak wrote:

> If Earth is taken as a sphere, it is easy and fast calculation in polar
> coordinates. If it's taken as an ellipsoid, then the radial distance to
> the point on the path, itself, becomes a variable. So the easy solution
> for the sphere case and the method to calculate it fast won't apply.

https://community.esri.com/t5/coordinate-reference-systems-blog/distance-
on-an-ellipsoid-vincenty-s-formulae/ba-p/902053

The haversine formula is in one of the links and falls into the category
of 'close enough' in my example. There's also a brief explanation of what
a haversine is, back to the wonderful days of six place tables.

Reprojection from one coordinate system to another gets a bit chewier.

https://geodesy.noaa.gov/PC_PROD/SPCS83/

You can download the source in Fortran 66, complete with some comments by
Vincenty. Note that there is also a NAD 27 version. There is a short
discussion of datums here:

https://geodesy.noaa.gov/datums/faq.shtml

It's non-trivial. Garmin GPS receivers will show NAD 27 coordinates and a
guy used those to post a cache on geocaching.com. The default is 83. The
difference on the ground at N 47 is about 200' which made finding the
cache interesting.

https://www.e-education.psu.edu/geog862/node/1794

has several pages of discussion. For even more fun

https://www.esri.com/about/newsroom/arcuser/moving-from-static-spatial-
reference-systems-in-2022/

afaik they didn't make the 2022 target date and I haven't heard much
lately. "The hard part will be dealing with the many exabytes of data
found in legacy datasets." no shit...

When you consider the distance from Timbuktu to Kalamazoo will add a time
dependency based on the motion of tectonic plates some formulae will need
revamping.

I'm enough of a nerd to find the big picture interesting rather than
delving into the intricacies of how cos() works.

Out of curiosity,

https://www.rome2rio.com/s/Timbuktu/Kalamazoo

"We are unable to find a way to get between Timbuktu and Kalamazoo".
Damn, and I was planning my summer vacation...

[rbowman]

On Tue, 20 Feb 2024 17:52:23 -0600, Physfitfreak wrote:


> 1. Did you derive the above relation by yourself?

Hell, no.

> 2. So you just threw the above expression at the computer so the
> compiler would figure it out what to do with it?

You betcha

> I'm asking the second question because if you just took the relation
> from some calculus book or somewhere on the web that had used calculus
> to come to that relation, that doesn't mean you can then just use it
> like that in your program

Ever hear about standing on the shoulders of giants? There are well known
formulae that I don't feel a need to reinvent. I never felt a burning
desire to redo Djikstra's shortest path algorithm either. You have to
realize there are people in the world concerned with larger problems.

fwiw, I also believe the pytorch library can do a better job of
calculating gradient descent that I can. The same can be said for most
third party libraries. I *could* write an XML parser but I'll link in the
expat library.

https://en.wikipedia.org/wiki/Expat_(software)

Note: because it has a MIT license I CAN use it. If it was GPL I wouldn't
touch it except for hobby projects.

[Lester Thorpe]

On Tue, 20 Feb 2024 20:06:21 -0600, Physfitfreak wrote:

>
> In the above expression it is not clear you mean to start it as

> + ... or you mean 33375/(100*y^6) + ...
>
> Same ambiguity about the term 55/10*y^8.
>

There is no ambiguity. The standard rules of mathematical operator
precedence apply. Therefore it is (33375/100)*y^6 and (55/10)*y^8.
That is how C will process the expressions.

[%]

sometimes assworms just hit send without adding anything

[DFS]

On 2/21/2024 4:45 AM, % wrote:
> % wrote:
>> chrisv wrote:

<snip rbowman code>

>>> I don't know if that needed to be tested, or if a clean compile was
>>> good enough to know that it will work correctly.
>>>
>>> According to some shameless liar, some dumb fscking snit.
>>>
>>
> sometimes assworms just hit send without adding anything

shitv has been extremely butthurt for nearly 10 years because he made
the ridiculous, ignorant claim that it's "stupid" to test code vs
relying on compiler warnings.

<d6plp958n3o545aqg...@4ax.com>

He can't get over it.

[Chris Ahlstrom]

% wrote this copyrighted missive and expects royalties:

> % wrote:

>> chrisv wrote:
>>>
>>> I don't know if that needed to be tested, or if a clean compile was
>>> good enough to know that it will work correctly.
>>>
>>> According to some shameless liar, some dumb fscking snit.

He was referring to a troll he quoted in his message signature, lied
about something chrisv said.

> sometimes assworms just hit send without adding anything

--
You are dishonest, but never to the point of hurting a friend.

[%]

Chris Ahlstrom wrote:
> % wrote this copyrighted missive and expects royalties:
>
>> % wrote:
>>> chrisv wrote:
>>>>
>>>> I don't know if that needed to be tested, or if a clean compile was
>>>> good enough to know that it will work correctly.
>>>>
>>>> According to some shameless liar, some dumb fscking snit.
>
> He was referring to a troll he quoted in his message signature, lied
> about something chrisv said.
>
>> sometimes assworms just hit send without adding anything
>

but i move the post to another place

[Relf]

Your code doesn't consider altitude ?

[rbowman]

On Wed, 21 Feb 2024 10:11:09 -0800 (Seattle), Relf wrote:

> Your code doesn't consider altitude ?

No, nor does it consider topographic features. It's 13 miles from here to
the intersection of Graves Creek Rd and US 12 -- if you're a crow.

[Physfitfreak]

I know that. But before spending time working on a possibly messy
problem, as a matter of principle, one should make sure the initiator of
the problem is also aware of that :)

[Physfitfreak]

On 2/20/2024 6:19 PM, Physfitfreak wrote:

> If it's taken as an ellipsoid, then the radial distance to the point on
> the path, itself, becomes a variable.

More than just that comes into the calculations, of course. The shortest
path on an ellipsoid between two points on it, does _not_, generally,
fall on any circle. So the solution path is weird because you're not
strolling down your path on a spherical surface.

The path doesn't even fall on a concentric ellipse, and not even on an
ellipse that's not concentric with the ellipsoid. So one is left to hit
the shit and do a line integral in the most general case.

But how can you find the shortest path's equation so you could integrate
and find its length? You have to find a general equation of any path
between the two points first, then calculate the line integral of that
general path to find the general distance, then minimize that general
distance by taking partial derivatives of that diabolic expression with
respect to the three coordinates, and equate the result to zero.

So you're left with a very nasty PDE to solve. Surprise! You're now in
physics, and are going to taste it too.

So you'd begin reviewing the available numerical methods to deal with
that particular form of PDE. If you're lucky, you will find one, and
will use it to turn the diabolic equation into a difference equation,
then and only then, you begin coding it carefully.

I give it 999999 out of a million chance that you'll find yourself
getting pure trash results if you use the public gcc to run your
program. You will then RUN, with LOVE and DEVOTION, towards Nelson
Beebe, and implore to use his library and especially his methods
explained in his book for handling basic math functions, and use them
for _every_ math function you find in those pages of expressions that
constitute your PDE.

And only then, you get results that will be in most cases good enough to
use. But hopefully not after waiting a week for each run of the
computation to end, but perhaps a few hours. Here, you've taken care of
the accuracy well enough, but the difference between asking for _one_
extra bit of accuracy, and not asking it, can be the difference between
two hours and two weeks of waiting for the run to end. So you became
more confident that you were correct in throwing your gcc away and
replaced them with mathcw codes.

So let's say you spent three months and did all that.

Now. Are you done? Of course not. Hehe :) Now, you only have the
equation of the shortest line between two general points on the
ellipsoid; i.e. a polynomial in three variables and 6 parameters for the
two points' coordinates. You give the coordinates of the two points to
it, and it gives you the equation of the shortest path between those
particular two points.

Now you're ready for the line integral. You'll find the appropriate
numerical method for it, and use it to compute the distance.

Are you done? Of course not. Even if from point 1 you can see the point
2 in your horizon, the shortest path will not be used if you begin
walking straight towards the destination that you see. The correct path
will be different, sometimes going a bit to the right of the line of
sight of the destination, and sometimes going a bit to the left of that
line of sight. But magically always ending up right at the destination
quicker than the guy who took the direct path stretching between him and
the destination he was seeing, assuming both you and him walked with
same speed!...

Even then, you're left with some cases that results aren't good enough.
There will always remain some gremlins in weird cases of the two points.

So after all that, you've tasted a case of "physics". Well, if others
had not done that work before you, you could get a master's degree in
physics with your work too :-)

I know of one "student" (he was a wealthy young guy who already owned
his own company, and his own one-engine airplane, etc) in graduate
school who solved just one nasty PDE resulting from a situation in space
science, and got his masters with it and left. If he'd stay and do two
more of such computational works, he'd get his PhD with them without
even passing the qualifying exams. But he didn't have much time. He was
always fighting fatigue. Sometimes when he was in my office asking a
question, I'd notice his eyes would begin the involuntary motion of a
micro-sleep, as he was listening to my answer. He was quite a guy
though. An enthusiast.

[Henry George]

On 2/20/2024 5:47 PM, Physfitfreak wrote:
>
> And before anything, I'd check the equation of ellipsoid in Cartesian
> coordinates as well to see whether it yields easier integrations over a
> path on its surface. Then will choose whether to use Cartesian or polar
> coordinates for this problem.

For a SPCS projection Cartesian coordinates are good enough. That's the
whole idea of a projection for a limited area. There is distortion like
the old flatten the orange peel analogy so you select the center wisely.
Missouri, for example, has 3 SPCS zones centered on St. Louis, Jefferson
City, and Kansas City MO.

It doesn't work well for larger areas. Most internet sites like google
use 'web mercator' which is essentially WGS84. That means when you take
a photo with a camera or phone with a GPS the EXIF will contain the
latitude and longitude, or in other words it's geotagged. Since google,
osm, esri, and other base maps are also lat/lon the photo can be
accurately positioned.

Sde note: if you're paranoid scrub the metadata before publicly
publishing a photo.