cexp error
Bruce Evans
brde at optusnet.com.au
Sat Sep 15 19:57:46 UTC 2012
On Sat, 15 Sep 2012, Stephen Montgomery-Smith wrote:
> On 09/15/2012 07:00 AM, Bruce Evans wrote:
>> ...
>> I have Kahan's program which does searches using the continued fraction
>> method (das sent it to me in 2008), but don't really understand it, and
>> in particular don't know how to hack it to give the bound for the 2**45
>> range. Saved results show that I didn't finish checking with it in 2008
>> either.
>
> I am just making a guess here. But I am thinking that Kahan's program might
> be based around the following type of argument. First, find integers p, q
> and r so that
> |pi - p/q| < 1/r
> The idea is that r should be much bigger than q. There is a theorem that
> says you can find a p, q and r so that r=q^2. In any case, I believe
> continued fractions is the way to find these kinds of numbers.
> ...
> Probably the argument used by Kahan or das@ is sharper than this.
google for nearpi.c (2001 version in ~wkahan).
It is poorly engineered, with no support for long doubles, no makefile,
hard-coded configuration in the source file, and no (?) usage info in
the source file. It doesn't quite compile, and after fixing the syntax
error it doesn't work for me.
The following patch (by das I think) makes it mostly work. It still doesn't
work on i386, but it seems to work on amd64. Its hard-coded values are
still inadequate for ld128, and it dumps core on sparc64.
% --- nearpi.c~ Tue Aug 7 13:56:21 2001
% +++ nearpi.c Fri Feb 29 10:18:51 2008
% @@ -1,2 +1,5 @@
% +#include <stdio.h>
% +#include <math.h>
% +#include <tgmath.h>
tgmath.h is critical (I think it is mainly a hack to avoid editing all
the hard-coded double function calls), but not available on my old i386
FreeBSD installations.
% /*
% * Nearpi, a C program to exhibit large floating-point numbers
% @@ -427,11 +430,12 @@
% * Global macro definitions.
% */
% +#include <float.h>
%
% -# define hex( double ) *(1 + ((long *) &double)), *((long *) &double)
% +# define hex( double ) *(1 + ((int *) &double)), *((int *) &double)
long wasn't so long in 1983, and *((int *)&foo is still 1983'ish.
% # define sgn(a) (a >= 0 ? 1 : -1)
% -# define MAX_k 2500
% -# define D 56
% -# define MAX_EXP 127
That was only for floats? 56 for vaxes?
% -# define THRESHOLD 2.22e-16
% +# define MAX_k 260000
% +# define D 64
Probably needs to be 113 for sparc64.
% +# define MAX_EXP 16383
% +# define THRESHOLD (LDBL_EPSILON/512)//2.22e-16
%
% /*
% @@ -441,5 +445,5 @@
% int CFlength, /* length of CF including terminator */
% binade;
% -double e,
% +long double e,
% f; /* [e,f] range of D-bit unsigned int of f;
% form 1X...X */
It had little chance of working for long doubles with none declared.
% @@ -453,5 +457,5 @@
% {
% int k; /* subscript variable */
% - double i[MAX_k],
% + long double i[MAX_k],
% j[MAX_k]; /* i and j are continued fractions
% (coeffs) */
% @@ -576,12 +580,12 @@
%
% dbleCF (i, j)
% -double i[],
% +long double i[],
% j[];
% {
% - register double k,
% + long double k,
% l,
% l0,
% j0;
% - register int n,
% + int n,
% m;
% n = 1;
% @@ -596,5 +600,5 @@
% break;
% };
% - modf (l / 2, &l0);
% + modfl (l / 2, &l0);
% l = l - l0 - l0;
% k = i[n + 1];
Removing the register declarations fixes &l0.
% @@ -679,8 +683,8 @@
%
% input (i)
% -double i[];
% +long double i[];
% {
% int k;
% - double j[MAX_k];
% + long double j[MAX_k];
%
% /*
% @@ -693,5 +697,5 @@
% {
% k = k + 1;
% - scanf ("%E", &i[k]);
% + scanf ("%LE", &i[k]);
% } while (i[k] >= 0);
%
% @@ -741,5 +745,5 @@
%
% nearPiOver2 (i)
% -double i[];
% +long double i[];
% {
% int k, /* subscript for recurrences (see
% @@ -747,5 +751,5 @@
% K; /* like k , but used during cancel. elim.
% */
% - double p[MAX_k], /* product of the q's (see
% + long double p[MAX_k], /* product of the q's (see
% handout) */
% q[MAX_k], /* successive tail evals of CF (see
% @@ -769,14 +773,6 @@
% m0, /* the mantissa of z as a D-bit integer
% */
% - x, /* the reduced argument (see
% + x; /* the reduced argument (see
% handout) */
% - ldexp (), /* sys routine to multiply by a power of
% - two */
% - fabs (), /* sys routine to compute FP absolute
% - value */
% - floor (), /* sys routine to compute greatest int <=
% - value */
% - ceil (); /* sys routine to compute least int >=
% - value */
%
% /*
% @@ -787,5 +783,5 @@
% */
%
% - q[CFlength] = 1.0E + 30;
% + q[CFlength] = 1.0E+30;
%
% for (k = CFlength - 1; k >= 0; k = k - 1)
The syntax error.
BTW, running clog() sources through indent(1) creates large diffs, all
wrong, mainly because indent(1) doesn't understand C99 (?) hex constants
or even C90 F suffixes. It mangles 0x1p30 to something like the above.
% @@ -907,12 +903,16 @@
% * Set z = m0 * 2 ^ (binade+1-D) .
% */
% - z = ldexp (m0, binade + 1 - D);
% + z = ldexpl (m0, binade + 1 - D);
%
% /*
% * Print z (hex), z (dec), m0 (dec), binade+1-D, x (hex), x (dec).
% */
% +#if 1
% printf ("\
% -%08x %08x Z=%22.16E M=%17.17G L+1-%d=%3d %08x %08x x=%23.16E\n",
% +%08x %08x Z=%22.16LE M=%17.17LG L+1-%d=%3d %08x %08x x=%23.16LE\n",
The 16s and 17s may be inadeqate even for double precision.
% hex (z), z, m0, D, binade + 1 - D, hex (x), x);
% +#endif
% + printf("%22.22LgL\tevalf(frac(%Lf*2^%d/(2*Pi)));\n", z, m0, binade +1-D);
22 is adequate for ld80. Perhaps 2 too many (DECIMAL_DIG is 21, but the format
gives 1 extra before the point). ld128 needs 36, 37 or 38.
% + fflush(stdout);
%
% }
So these changes are routine.
You also need an input file with the continued fraction expansion of Pi
(3\n7\n15\n1\n292\n...). das sent me one with 25542 terms.
I generated the following test vectors from the output from this, and
recently cleaned it up:
% #include <float.h>
% #include <math.h>
% #include <stdio.h>
%
% static void
% test(long double x)
% {
% printf("% .*Le % .*Le\n", DECIMAL_DIG, cosl(x), DECIMAL_DIG, sinl(x));
% }
%
% int
% main(void)
% {
% test(14529431823429108538.L * 0x1p98L);
% test(10823840566309211864.L * 0x1p230L);
% test(15670828281847873939.L * 0x1p271L);
% test(15670828281847873939.L * 0x1p272L);
% test(15670828281847873939.L * 0x1p273L);
% test(10396564734341006689.L * 0x1p1217L);
% test(16875986960654748944.L * 0x1p1443L);
% test(16279065367393374496.L * 0x1p1775L);
% test(15294540618075682663.L * 0x1p1842L);
% test(11494079389130875706.L * 0x1p1913L);
% test(15091839851549065173.L * 0x1p1957L);
% test(10192914972581512681.L * 0x1p2006L);
% test(16218102405256695008.L * 0x1p3694L);
% test(15775002102157321889.L * 0x1p4389L);
% test(14742866506737243196.L * 0x1p4896L);
% test(13638149838171524804.L * 0x1p5401L);
% test(15805865952080730952.L * 0x1p5750L);
% test(15449633091134885233.L * 0x1p6509L);
% test(16314936504709115283.L * 0x1p6547L);
% test(18129316451163514697.L * 0x1p6868L);
% test(18129316451163514697.L * 0x1p6869L);
% test(17033015872180140473.L * 0x1p7070L);
% test(17033015872180140473.L * 0x1p7071L);
% test(17033015872180140473.L * 0x1p7072L);
% test(18164632514172025418.L * 0x1p7858L);
% test(16386729801286219183.L * 0x1p8190L);
% test(16386729801286219183.L * 0x1p8191L);
% test(16386729801286219183.L * 0x1p8192L);
% test(15435742767874430341.L * 0x1p8426L);
% test(15435742767874430341.L * 0x1p8427L);
% test(15435742767874430341.L * 0x1p8428L);
% test(15435742767874430341.L * 0x1p8429L);
% test(10290136764324338438.L * 0x1p8445L);
% test(13831617603129791506.L * 0x1p9240L);
% test(13800459120262740616.L * 0x1p9610L);
% test(16069825796649748274.L * 0x1p10154L);
% test(11122677606925071542.L * 0x1p10431L);
% test(16362320716166409536.L * 0x1p10513L);
% test(17476981849448541921.L * 0x1p10531L);
% test(17476981849448541921.L * 0x1p10532L);
% test(17476981849448541921.L * 0x1p10533L);
% test(17476981849448541921.L * 0x1p10534L);
% test(17476981849448541921.L * 0x1p10535L);
% test(17476981849448541921.L * 0x1p10536L);
I wonder what or who generated sequences like this.
% test(12570838853826129195.L * 0x1p10955L);
% test(12570838853826129195.L * 0x1p10956L);
% test(12452308034733593216.L * 0x1p10981L);
% test(17180372579887641792.L * 0x1p11409L);
% test(12672796772319852249.L * 0x1p11779L);
% test(12672796772319852249.L * 0x1p11780L);
% test(12672796772319852249.L * 0x1p11781L);
% test(12371399102696562706.L * 0x1p11922L);
% test(12827592373563150124.L * 0x1p12574L);
% test(10411496087288248933.L * 0x1p12587L);
% test(10411496087288248933.L * 0x1p12588L);
% test(10411496087288248933.L * 0x1p12589L);
% test(10411496087288248933.L * 0x1p12590L);
% test(12374443443369762001.L * 0x1p13142L);
% test(11840786828664647755.L * 0x1p13996L);
% test(9868908899506664085.L * 0x1p14270L);
% test(9868908899506664085.L * 0x1p14271L);
% test(11711452751349741921.L * 0x1p14307L);
% test(11711452751349741921.L * 0x1p14308L);
% test(17278682452196525512.L * 0x1p14444L);
% test(15177699244736331235.L * 0x1p15776L);
% test(14227603679013389146.L * 0x1p15865L);
% return (0);
% }
The nearness can be read off from the value of cosl() or sinl(). All values
in the table are near to within a small multiple of 1e-23, with the worst
case:
test(17476981849448541921. * 2^10531);
producing:
-1.823402780633777070226e-23 -1.000000000000000000000e+00
pari agrees with amd64 and sparc64 libm for cosl() and sinl() on these values.
The above program is too naive to work on i386 (gcc truncates all the long
constants to 53 bits).
Bruce
More information about the freebsd-numerics
mailing list