I am debugging with gdb and I find that the quadruple precision numbers are not displayed properly. The code I am testing with is:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
extern "C" {
#include "quadmath.h"
}
int main ()
{
char buf[128];
int n;
double x = 4.0;
double y = x + 1e-5;
double z = 0.0;
z = (exp(x) - exp(y))/(x-y);
printf("In double precision : %.32e\n", z);
double zz = exp(x/2)*exp(y/2)*2*sinh((x-y)/2)/(x-y);
printf("In double precision using sinh: %.32e\n", zz);
__float128 qx = x;
__float128 qy = y;
__float128 qz = 0.0q;
qz = (expq(qx) - expq(qy))/(qx-qy);
n = quadmath_snprintf(buf, sizeof buf, "%.32Qe", qz);
if ((size_t) n < sizeof buf)
printf ("In quadruple precision : %s\n", buf);
return 0;
}
(Note the extern to work around a libquadmath bug.)
I compile with: g++ a.C -g -lquadmath
Then I get:
Breakpoint 1, main () at a.C:19
19 printf("In double precision using sinh: %.32e\n", zz);
(gdb) n
In double precision using sinh: 5.45984230248043687083736585918814e+01
21 __float128 qx = x;
(gdb)
22 __float128 qy = y;
(gdb) p qx
$1 = 0
which appears to be a gdb error, because the final result appears to be ok. My question simply is whether this is supposed to work? Is there gdb support for quads? Do I need to do something special to get gdb work under these circumstances?
The point of the program is to demonstrate the catastrophic cancellation and that using sinh provides the required (double precision) accuracy without quad arithmetic:
In double precision : 5.45984230249338509111112216487527e+01
In double precision using sinh: 5.45984230248043687083736585918814e+01
In quadruple precision : 5.45984230248043659065767385111870e+01
Thanks!
P.S.: Could not tag with quadruple-precision, I do not have enough reputation to create new tags.