Fractalforums

Added timings for perturbed iterations (excluding reference orbit calculations and glitch correction passes) to the table above.

I get an unexpected result using f107. Has someone else had issues with the datatype - or did I make a mistake in the coding below?

I use the following three numbers:

a1: 570949632*2^-32 +          0*2^-64
a2:   1899995*2^-32 + 1073741824*2^-64
a3:  32192101*2^-32 + 3782580368*2^-64

and computed (aa, bb was set to 0 at the start)

aa=aa+a1-a2-a3 = 536857534*2^-32 + 3733612400*2^-64 = 144111597179676279/1152921504606846976
bb=bb-a2-a3+a1 = 536857534*2^-32 + 3733612400*2^-64

which are expectedly identical using WolframAlpha, GNU's __float128 and my fpa datatype using the statement (aa==bb).

But when I use f107 I get "not-identical" (0) as a result.
(APPROX and FMA are both off).

If changing the order of an expression leads to different results, that hints to rounding errors.

But where would those be? a1,a2,a3 are at most 64 bits to the right of the radix, and none to the left. So they can be represented by f107 offering (author's statement) 107 bits mantissa precision. For adding/subtracting them from one another, there's room enough for carry-over.

Can anyone repeat that result?

f107_o aa,a1,a2,a3,bb,two8,two32,tp;

two8=256.0;
two32=1.0;
two32 = two32 * two8;
two32 = two32 * two8;
two32 = two32 * two8;
two32 = two32 * two8;

a1=570949632;
a1 /= two32;

a2=1899995;
a2 /= two32;
tp=1073741824;
tp /= two32;
tp /= two32;
a2 = a2+tp;

a3=32192101;
a3 /= two32;
tp=3782580368;
tp /= two32;
tp /= two32;
a3 = a3+tp;

aa=0.0;
aa=aa+a1-a2-a3;

bb=0.0;
bb=bb-a2-a3+a1;

printf("is aa==bb => %i\n",(aa==bb));


Can't reproduce with "g++ (Debian 8.3.0-6) 8.3.0" on "Intel(R) Core(TM)2 Duo CPU     P7550  @ 2.26GHz".

Test prints 1 in all variations I've tried, with both libqd's dd_real and mrob's f107_o.  Adding -O3 -ffast-math gives incorrect result (smaller double is 0 instead of 6e-18), but the test still prints 1.  I note the smaller double is 0 for all of a1, a2, a3 so they all fit in 53 bits.

Additional code:
#include <stdio.h>
#if 0
#include <qd/dd_real.h>
typedef dd_real R;
#else
#define F107_APPROX 0
#define F107_FMA 0
#include "f107_o.h"
typedef f107_o R;
double to_double(R x)
{
  double r(x);
  return r;
}
#endif

int main(int argc, char **argv)
{
  // declare variables as R
...
printf("a1 = %.18e + %.18e\n", to_double(a1), to_double(a1 - to_double(a1)));
printf("a2 = %.18e + %.18e\n", to_double(a2), to_double(a2 - to_double(a2)));
printf("a3 = %.18e + %.18e\n", to_double(a3), to_double(a3 - to_double(a3)));
printf("aa = %.18e + %.18e\n", to_double(aa), to_double(aa - to_double(aa)));
printf("bb = %.18e + %.18e\n", to_double(bb), to_double(bb - to_double(bb)));
...
return 0;
}

@claude: Thanks very much for your coding and testing effort!

After reading you couldn't reproduce my results, I did a new clean install of the TDM compiler, and now it works fine.
Not sure what happened there, next time I'll do that as a first step.