30位有效数字的浮点数结构解决double数据类型多次累加后的明显的误差

标准的c或者c++的double数据类型只有15位有效数字（好像有这么回事，参看IEEE 754 ），

因此产生了大的数字多次累加后的明显的误差，在财务计算中，这种误差是不能接受的。

参看http://blog.csdn.net/l1t/archive/2004/10/09/129110.aspx

来自Lcc-win32 的src\doubledouble\doubledouble.c

（作者主页：http://members.lycos.co.uk/keithmbriggs/doubledouble.html

上还有一个c++封装的类，有兴趣的人可以去看。）

利用2个double变量构造出一个doubledouble结构，解决了这个问题。

测试办法：把main()函数的内容替换成下面内容。

int main(void)

{

doubledouble a = new_doubledouble(5.0,0.0);

doubledouble b = new_doubledouble(4.0,0.0);

doubledouble c = a*b;

double x1,x2;

char buffer[256];

int i;

a=a-a+99999999.99;

c=c-c+0.0;

for(i=0;i<81920;i++)

{

c=c+a;

}

doubledoubletoa(c,buffer);

printf("sum() =%s\n",buffer);

c=a*81920.0;

doubledoubletoa(c,buffer);

printf("mult()=%s\n",buffer);

x1=c.hi+c.lo;

printf("mult()=%18.2lf\n",x1);

return 0;

}

运行结果：

sum() =8.191999999180799560546875000000000e12

mult()=8.191999999180799560546875000000000e12

mult()= 8191999999180.80

在Lcc-win32 中，可以不加修改地通过编译。

如果在gcc下编译，需要定义宏#define GCC

在gcc和vs.net 2003 中，函数

doubledouble ddexp(const doubledouble& x)

{

...

sum1=y*((((ysq+3960.)*ysq+2162160.)*ysq+302702400.)*ysq+8821612800.);

sum2=(((90.*ysq+110880.)*ysq+30270240.)*ysq+2075673600.)*ysq+17643225600.;

...

}

上述2行不能通过编译，如果不需要这个函数，可以删除。

/* doubledouble.c */

#include <string.h>

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

This software was written by the mathematician Keith Martin Briggs. It implements

an extended precision data type "doubledouble", that is made up of two floating

point numbers. I have rewritten it for lcc-win32 and extensiveley used this code

to test the compiler's operator subclassing module.

The code should compile without problems in any C++ compiler.

---------------------------------------------------------------jacob

Except where otherwise indicated,

This program is free software; you can redistribute it and/or

modify it under the terms of the GNU General Public License

as published by the Free Software Foundation; either version 2

of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,

but WITHOUT ANY WARRANTY; without even the implied warranty of

MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

GNU General Public License for more details.

You should have received a copy of the GNU General Public License

along with This program; if not, write to the Free Software

Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

// 97 Aug 04 KMB added ldexp

// 97 Jul 11 Moved This stuff out of quad.h, created inline.h so it can

// be #included even if we're not inlining, by just "#define inline"

// 97 Jul 12 Added all combos of doubledouble/double/int +-*/. Only a couple actually

// written; most just call the others by swapping arguments. Generates

// equivalent code with a good inlining compiler (tested with g++ 2.7.2).

// - e.g., all subtraction is now done by addition of unary minus.

// - removed if's from doubledouble*int. Zero-branch code is 2.5 faster (tested)

// - generally cleaned up and re-organized the order of the functions,

// now they're grouped nicely by function.

// - tested Newton's division. Works but is terribly slow, because

// it needs to do several doubledouble + and * operations, and doubledouble /

// without it is about the same speed at doubledouble * anyway, so no win.

// - moved recip here as an inline.

// - checked and tested doubledouble/double (BUG?), seems fine.

typedef struct doubledouble {

double hi;

double lo;

} doubledouble;

#ifdef __LCC__

#define x86_FIX unsigned short __old_cw, __new_cw; _asm("\tfnstcw\t%__old_cw"); __new_cw = (__old_cw & ~0x300) | 0x200; _asm ("\tfldcw\t%__new_cw")

#define END_x86_FIX _asm ("\tfldcw\t%__old_cw")

#else

#define x86_FIX unsigned short __old_cw,__new_cw; _asm {fnstcw __old_cw} __new_cw = (__old_cw & ~0x300) | 0x200; _asm { fldcw __new_cw};

#define END_x86_FIX _asm { fldcw __new_cw}

#endif

#ifdef GCC

#undef x86_FIX

#undef END_x86_FIX

#define x86_FIX

#define END_x86_FIX

#endif

static double Split = 134217729.0;

doubledouble new_doubledouble(double x,double y);

doubledouble operator-(const doubledouble &x)

{

doubledouble z;

z.hi = -x.hi;

z.lo = -x.lo;

return z;

}

// Absolute value

doubledouble ddfabs(const doubledouble& x)

{

if (x.hi>=0.0) return x;

else {

doubledouble result;

result.hi = -x.hi;

result.lo = -x.lo;

return result;

}

// Addition

/* (C) W. Kahan 1989

* NOTICE:

* Copyrighted programs may not be translated, used, nor

* reproduced without the author's permission. Normally that

* permission is granted freely for academic and scientific

* purposes subject to the following three requirements:

* 1. This NOTICE and the copyright notices must remain attached

* to the programs and their translations.

* 2. Users of such a program should regard themselves as voluntary

* participants in the author's researches and so are obliged

* to report their experience with the program back to the author.

* 3. Neither the author nor the institution that employs him will

* be held responsible for the consequences of using a program

* for which neither has received payment.

* Would-be commercial users of these programs must correspond

* with the author to arrange terms of payment and warranty.

///////////////////////////////////////////////////////////////

/////////////////// Addition and Subtraction /////////////////

///////////////////////////////////////////////////////////////

// Binary addition

doubledouble operator +(const doubledouble& x,const doubledouble& y)

{

double H, h, T, t, S, s, e, f;

doubledouble z;

x86_FIX;

S = x.hi+y.hi;

T = x.lo+y.lo;

e = S-x.hi;

f = T-x.lo;

// s = S-e;

t = T-f;

s = (y.hi-e)+(x.hi-(S-e));

t = (y.lo-f)+(x.lo-(T-f));

e = s+T;

H = S+e;

h = e+(S-H);

e = t+h;

z.hi = H+e;

z.lo = e+(double)(H-z.hi);

END_x86_FIX;

return z;

}

doubledouble operator +(const double& x,const doubledouble& y)

{

double H, h, S, s, e;

doubledouble z;

x86_FIX;

S = x+y.hi;

e = S-x;

s = S-e;

s = (y.hi-e)+(x-s);

H = S+(s+y.lo);

h = (s+y.lo)+(S-H);

z.hi = H+h;

z.lo = h+(double)(H-z.hi);

END_x86_FIX;

return z;

}

doubledouble operator +(const doubledouble& x,const double& y )

{

return y+x;

}

doubledouble operator +(int x,const doubledouble& y)

{

return (double)(x)+y;

}

doubledouble operator +(doubledouble& x, int y) {

return y+x;

}

// Subtraction

doubledouble operator -(const doubledouble& x,const doubledouble& y)

{

return x+(-y);

}

doubledouble operator -(const double& x, const doubledouble& y)

{

return x+(-y);

}

doubledouble operator -(const int x, const doubledouble& y)

{

return x+(-y);

}

doubledouble operator -(const doubledouble& x, const double& y)

{

return x+(-y);

}

doubledouble operator -(const doubledouble& x,const int y)

{

return x+(-y);

}

////////////////////////// Self-Addition ///////////////////////

doubledouble& operator +=(doubledouble &This,const doubledouble& y)

{

double H, h, T, t, S, s, e, f;

x86_FIX;

S =This.hi+y.hi;

T = This.lo+y.lo;

e = S-This.hi;

f = T-This.lo;

s = S-e;

t = T-f;

s = (y.hi-e)+(This.hi-s);

t = (y.lo-f)+(This.lo-t);

f = s+T;

H = S+f;

h = f+(S-H);

This.hi = H+(t+h);

This.lo = (t+h)+(double)(H-This.hi);

END_x86_FIX;

return This;

}

doubledouble& operator +=(doubledouble &This,const double& y)

{

double H, h, S, s, e, f;

x86_FIX;

S =This.hi+y;

e = S-This.hi;

s = S-e;

s = (y-e)+(This.hi-s);

f = s+This.lo;

H = S+f;

h = f+(S-H);

This.hi = H+h;

This.lo = h+(double)(H-This.hi);

END_x86_FIX;

return This;

}

doubledouble& operator +=(doubledouble &This,const int y)

{

return This += (double)(y);

}

// Self-Subtraction

doubledouble& operator -=(doubledouble &This,const doubledouble& y)

{

return This += -y;

}

doubledouble& operator -=(doubledouble &This,const double& y)

{

return This += -y;

}

doubledouble& operator -=(doubledouble &This,const int y)

{

return This += -y;

}

/////////////////////////////////////////////////////////////

//////////////////// Multiplication ///////////////////////

/////////////////////////////////////////////////////////////

// Binary Multiplication

doubledouble operator *(const doubledouble& x,const doubledouble& y)

{

double hx, tx, hy, ty, C, c;

doubledouble z;

x86_FIX;

C = Split*x.hi;

hx = C-x.hi;

c = Split*y.hi;

hx = C-hx;

tx = x.hi-hx;

hy = c-y.hi;

C = x.hi*y.hi;

hy = c-hy;

ty = y.hi-hy;

c = ((((hx*hy-C)+hx*ty)+tx*hy)+tx*ty)+(x.hi*y.lo+x.lo*y.hi);

z.hi = C+c;

hx=C-z.hi;

z.lo = c+hx;

END_x86_FIX;

return z;

}

// double*doubledouble

doubledouble operator *(const double& x, const doubledouble& y)

{

double hx, tx, hy, ty, C, c;

doubledouble z;

x86_FIX;

C = Split*x;

hx = C-x;

c = Split*y.hi;

hx = C-hx ;

tx = x-hx;

hy = c-y.hi;

C = x*y.hi;

hy = c-hy;

ty = y.hi - hy;

c = ((((hx*hy-C)+hx*ty)+tx*hy)+tx*ty)+x*y.lo;

z.hi = C+c;

z.lo = c+(double)(C-z.hi);

END_x86_FIX;

return z;

}

doubledouble operator *(int x, const doubledouble& y )

{

return ((const double)(x))*y;

}

doubledouble operator *(const doubledouble& x,const double& y )

{

return y*x;

}

doubledouble operator *(const doubledouble& x,const int y )

{

return y*x;

}

// Self-multiplication

doubledouble& operator *=(doubledouble &This,const doubledouble& y )

{

double hx, tx, hy, ty, C, c;

x86_FIX;

C = Split*This.hi;

hx = C-This.hi;

c = Split*y.hi;

hx = C-hx;

tx =This.hi-hx;

hy = c-y.hi;

C =This.hi*y.hi;

hy = c-hy;

ty = y.hi-hy;

c = ((((hx*hy-C)+hx*ty)+tx*hy)+tx*ty)+(This.hi*y.lo+This.lo*y.hi);

hx = C+c;

This.hi = hx;

This.lo = c+(double)(C-hx);

END_x86_FIX;

return This;

}

doubledouble& operator *=(doubledouble &This,const double& y )

{

double hx, tx, hy, ty, C, c;

x86_FIX;

C = Split*This.hi;

hx = C-This.hi;

c = Split*y;

hx = C-hx;

tx =This.hi-hx;

hy = c-y;

C =This.hi*y;

hy = c-hy;

ty = y-hy;

c = ((((hx*hy-C)+hx*ty)+tx*hy)+tx*ty)+(This.lo*y);

hx = C+c;

This.hi = hx;

This.lo = c+(double)(C-hx);

END_x86_FIX;

return This;

}

doubledouble& operator *=(doubledouble &This,const int y )

{

return This *= (const double)(y);

}

////////////////////////////////////////////////////////////////

////////////////////////// Division //////////////////////////

////////////////////////////////////////////////////////////////

// Reciprocal

doubledouble recip(const doubledouble& y) {

double hc, tc, hy, ty, C, c, U, u;

doubledouble z;

x86_FIX;

C = 1.0/y.hi;

c = Split*C;

hc =c-C;

u = Split*y.hi;

hc = c-hc;

tc = C-hc;

hy = u-y.hi;

U = C*y.hi;

hy = u-hy;

ty = y.hi-hy;

u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty;

c = ((((1.0-U)-u))-C*y.lo)/y.hi;

z.hi = C+c;

z.lo = (double)(C-z.hi)+c;

END_x86_FIX;

return z;

}

doubledouble operator /(const doubledouble& x,const doubledouble& y )

{

double hc, tc, hy, ty, C, c, U, u;

doubledouble z;

x86_FIX;

C = x.hi/y.hi;

c =Split*C;

hc =c-C;

u = Split*y.hi;

hc = c-hc;

tc = C-hc;

hy = u-y.hi;

U = C * y.hi;

hy = u-hy;

ty = y.hi-hy;

u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty;

c = ((((x.hi-U)-u)+x.lo)-C*y.lo)/y.hi;

z.hi = C+c;

z.lo = (double)(C-z.hi)+c;

END_x86_FIX;

return z;

}

// double/doubledouble:

doubledouble operator /(const double& x,const doubledouble& y )

{

double hc, tc, hy, ty, C, c, U, u;

doubledouble z;

x86_FIX;

C = x/y.hi;

c = Split*C;

hc =c-C;

u = Split*y.hi;

hc = c-hc;

tc = C-hc;

hy = u-y.hi;

U = C*y.hi;

hy = u-hy;

ty = y.hi-hy;

u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty;

c = ((((x-U)-u))-C*y.lo)/y.hi;

z.hi = C+c;

z.lo = (double)(C-z.hi)+c;

END_x86_FIX;

return z;

}

doubledouble operator /(const doubledouble& x,const double& y )

{

double hc, tc, hy, ty, C, c, U, u;

doubledouble z;

x86_FIX;

C = x.hi/y;

c = Split*C;

hc = c-C;

u = Split*y;

hc = c-hc;

tc = C-hc;

hy = u-y;

U = C*y;

hy = u-hy;

ty = y-hy;

u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty;

c = (((x.hi-U)-u)+x.lo)/y;

z.hi = C+c;

z.lo = (double)(C-z.hi)+c;

END_x86_FIX;

return z;

}

// doubledouble/int

doubledouble operator /(const doubledouble& x,const int y)

{

return x/(const double)(y);

}

doubledouble operator /(const int x,const doubledouble& y)

{

return (const double)(x)/y;

}

// Self-division.

doubledouble& operator /=(doubledouble &This,const doubledouble& y)

{

double hc, tc, hy, ty, C, c, U, u;

x86_FIX;

C =This.hi/y.hi;

c = Split*C;

hc =c-C;

u = Split*y.hi;

hc = c-hc;

tc = C-hc;

hy = u-y.hi;

U = C * y.hi;

hy = u-hy;

ty = y.hi-hy;

u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty;

c = ((((This.hi-U)-u)+This.lo)-C*y.lo)/y.hi;

u = C+c;

This.hi = u;

This.lo = (double)(C-u)+c;

END_x86_FIX;

return This;

}

doubledouble& operator /=(doubledouble &This,const double& y)

{

double hc, tc, hy, ty, C, c, U, u;

x86_FIX;

C =This.hi/y;

c = Split*C;

hc =c-C;

u = Split*y;

hc = c-hc;

tc = C-hc;

hy = u-y;

U = C * y;

hy = u-hy;

ty = y-hy;

u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty;

c = (((This.hi-U)-u)+This.lo)/y;

u = C+c;

This.hi = u;

This.lo = (double)(C-u)+c;

END_x86_FIX;

return This;

}

doubledouble& operator /=(doubledouble &This,const int y)

{

return (This) /=(const double)(y);

}

doubledouble Qcopysign(const doubledouble& x,const double y)

{

if (y>=0.0) return ddfabs(x);

else return -ddfabs(x);

}

doubledouble atodd(const char *s)

{

doubledouble result = {

0.0,0.0 };

int n, sign, ex = 0;

/* eat whitespace */

while (*s == ' ' || *s == '\t' || *s == '\n') s++;

switch (*s) { // get sign of mantissa

case '-':

{

sign = -1;

s++;

break;

}

case '+':

s++; // no break

default:

sign = 1;

}

/* get digits before decimal point */

while (n=(*s++)-'0', n>=0 && n<10) result = 10.0*result+n;

s--;

if (*s == '.') /* get digits after decimal point */ {

s++;

while (n=(*s++)-'0', n>=0 && n<10) {

result = 10.0*result+n;

--ex;

}

s--;

}

if (*s == 'e' || *s == 'E') /* get exponent */ ex+=atoi(++s);

/* exponent adjustment */

while (ex-- > 0) result *= 10.0;

while (++ex < 0) result /= 10.0;

return (sign>=0) ? result : -result;

}

// Square (faster than x*x)

doubledouble sqr(const doubledouble& x)

{

double hx,tx,C,c;

doubledouble r;

x86_FIX;

C=Split*x.hi;

hx=C-x.hi;

hx=C-hx;

tx=x.hi-hx;

C=x.hi*x.hi;

c=((((hx*hx-C)+2.0*hx*tx))+tx*tx)+2.0*x.hi*x.lo;

hx=C+c;

r.hi=hx;

r.lo = c+(C-hx);

END_x86_FIX;

return r;

}

extern double pow(double,double);

// doubledouble^int

doubledouble powint(const doubledouble& u, int c)

{

doubledouble result;

if (c<0) return recip(powint(u,-c));

switch (c) {

case 0:

result.hi = result.lo = 0.0;

if (u.hi == 0.0)

result.lo = pow(0.0,0.0);

else

result.lo = 1.0;

return result;

// let math.h handle 0^0.

case 1:

return u;

case 2:

return sqr(u);

case 3:

return sqr(u)*u;

default:

{ // binary method

int n=c, m;

doubledouble y={

1.0,0.0 }

, z=u;

if (n<0) n=-n;

while (1) {

m=n;

n/=2;

if (n+n!=m) { // if m odd

y*=z;

if (0==n) return y;

}

z=sqr(z);

}

/* not reached */

return u;

}

doubledouble new_doubledouble(double x,double y)

{

doubledouble result;

x86_FIX;

result.hi = x+y;

result.lo = y+(x-result.hi);

END_x86_FIX;

return result;

}

// Floor. See Graham, Knuth and Patashnik `Concrete Mathematics' page 70.

// Greatest integer <= x

// maybe floor_s is better?

doubledouble ddfloor(const doubledouble& x)

{

double fh=floor(x.hi), fl=floor(x.lo);

double t1=x.hi-fh;

double t2=x.lo-fl;

double t3=t1+t2;

int t;

t3 = x.hi-fh+x.lo-fl;

t=(int)(floor(t3));

switch (t) {

case 0:

return new_doubledouble(fh,0.0)+new_doubledouble(fl,0.0);

case 1:

return new_doubledouble(fh,0.0)+new_doubledouble(fl+1,0.0);

// case 2 can only arise when fp({hi})<=1, fp({lo})<=1, but

fp({hi}+{lo})=2

case 2:

return new_doubledouble(fh,0.0)+new_doubledouble(fl+1.0,0.0);

}

// never get here

return new_doubledouble(0.0,0.0);

}

// Signum

int sign(doubledouble& x)

{

if (x.hi>0.0) return 1;

else if (x.hi<0.0) return -1;

else return 0;

}

char *doubledoubletoa(doubledouble &x,char *s)

{

int Digits;

doubledouble ten = new_doubledouble(10.0,0.0),l,y,z;

double q;

int m,n,d,i;

if (x.hi==0.0) {

strcpy(s, "0.0 ");

return s;

}

if (x.hi!=x.hi) {

strcpy(s , "NaN ");

return s;

}

Digits=34;

y=ddfabs(x);

q=log10(y.hi);

n=(int)(floor(q));

if (n<0) n++;

if (n == 0)

l = new_doubledouble(1.0,0.0);

else

l = powint(ten,n);

y = y/l;

if (sign(x)<0) *s++ = '-'; //else s<<" ";

d = Digits>34 ? 34 : Digits;

d = d<3 ? 3 : d;

for (i=1; i<=d; i++) {

if (2==i) *s++ = '.';

m = (int)(ddfloor(y).hi);

if (m<0 || m>9) {

// fprintf(stderr,

// "Internal error in doubledouble.cc: digit=%d",m);

}

sprintf(s,"%d",m);

s += strlen(s);

// z.lo = 0.0;

// z.hi = (double)m;

z = new_doubledouble((double)m,0.0);

y = (y-z)*ten;

if (y.hi<=0.0) break; // x must be an integer

//if (y.hi==0.0) break; // ???????????

}

if (n!=0) {

*s++ = 'e';

sprintf(s,"%d",n);

s += strlen(s);

}

else *s++ = 0;

return s;

}

int operator> (const doubledouble& x, const doubledouble& y)

{

return (x.hi> y.hi) || (x.hi==y.hi && x.lo> y.lo);

}

int operator>=(const doubledouble& x, const doubledouble& y)

{

return (x.hi>y.hi) || (x.hi==y.hi && x.lo>=y.lo);

}

int operator< (const doubledouble& x, const doubledouble& y)

{

return (x.hi< y.hi) || (x.hi==y.hi && x.lo< y.lo);

}

int operator<=(const doubledouble& x, const doubledouble& y)

{

return (x.hi<y.hi) || (x.hi==y.hi && x.lo<=y.lo);

}

int operator==(const doubledouble& x, const doubledouble& y)

{

return x.hi==y.hi && x.lo==y.lo;

}

int operator!=(const doubledouble& x, const doubledouble& y)

{

return x.hi!=y.hi || x.lo!=y.lo;

}

struct ddConstants {

doubledouble Log2;

doubledouble Log10;

doubledouble Pi;

doubledouble TwoPi;

doubledouble Pion2;

doubledouble Pion4;

doubledouble _Pi;

} ddconst;

void InitConstants(void)

{

ddconst.Log2 =atodd("0.6931471805599453094172321214581765680755");

ddconst.Log10=atodd("2.302585092994045684017991454684364207601");

ddconst.Pi =atodd("3.1415926535897932384626433832795028841972");

ddconst.TwoPi=atodd("6.2831853071795864769252867665590057683943");

ddconst.Pion2=atodd("1.5707963267948966192313216916397514420985");

ddconst.Pion4=atodd("0.7853981633974483096156608458198757210493");

ddconst._Pi =atodd("0.3183098861837906715377675267450287240689");

}

doubledouble ddldexp(const doubledouble x, const int exp)

{ // x*2^exp

return new_doubledouble(ldexp(x.hi,exp),ldexp(x.lo,exp));

}

// Another floor. V. Shoup 97 Sep 23...

doubledouble ddfloor_s(const doubledouble& x)

{

double fhi=floor(x.hi);

if (fhi!=x.hi)

return new_doubledouble(fhi,0.0);

else

return new_doubledouble(fhi,floor(x.lo));

}

// rint (round to nearest int)

doubledouble ddrint(const doubledouble& x)

{

const doubledouble &rx = x+new_doubledouble(0.5,0.0);

return ddfloor(rx);

}

doubledouble ddexp(const doubledouble& x)

{

/* Uses some ideas by Alan Miller

Method:

x x.log2(e) nint[x.log2(e)] + frac[x.log2(e)]

e = 2 = 2

iy fy

= 2 . 2

Then

fy y.loge(2)

2 = e

Now y.loge(2) will be less than 0.3466 in absolute value.

This is halved and a Pade' approximant is used to approximate e^x over

the region (-0.1733, +0.1733). This approximation is then squared.

int iy;

doubledouble y,temp,ysq,sum1,sum2;

if (x.hi<-744.4400719213812) return new_doubledouble(0.0,0.0); //

exp(x)<1e-300

y=x/ddconst.Log2;

iy = (int)(ddrint(y)).hi;

temp=new_doubledouble((double)iy,0.0);

y=(y-temp)*ddconst.Log2;

y=ddldexp(y,-1);

ysq=sqr(y);

sum1=y*((((ysq+3960.)*ysq+2162160.)*ysq+302702400.)*ysq+8821612800.);

sum2=(((90.*ysq+110880.)*ysq+30270240.)*ysq+2075673600.)*ysq+17643225600.;

sum2 + sum1 2.sum1

Now approximation = ----------- = 1 + ----------- = 1 + 2.temp

sum2 - sum1 sum2 - sum1

Then (1 + 2.temp)^2 = 4.temp.(1 + temp) + 1

temp=sum1/(sum2-sum1);

y=temp*(temp+1);

y=ddldexp(y,2);

return ddldexp(y+1,iy);

}

// Natural logarithm

doubledouble ddlog(const doubledouble& x)

{ // Newton method

doubledouble s,e;

if (x.hi<0.0) {

fprintf(stderr,"\ndoubledouble: attempt to take log of negative

argument!\n");

return new_doubledouble(0.0,0.0);

}

if (x.hi==0.0) {

fprintf(stderr,"\ndoubledouble: attempt to take log(0)!\n");

return new_doubledouble(0.0,0.0);

}

s=new_doubledouble(log(x.hi),0.0), e=ddexp(s); // s=double approximation to

result

return s+(x-e)/e; // Newton correction, good enough

//doubledouble d=(x-e)/e; return s+d-0.5*sqr(d); // or 2nd order correction

}

void base_and_prec(void) { // Linnainmaa ACM TOMS 7, 272 Thm 3

int p;

x86_FIX;

printf("Base and precision determination by Linnainmaa's method:\n");

{

double U,R,u,r,beta;

U=4.0/3.0;

U-=1;

U*=3;

U-=1;

U=fabs(U);

R=U/2+1;

R-=1;

if (R!=0.0) U=R;

u=2;

u/=3;

u-=0.5;

u*=3;

u-=0.5;

u=fabs(u);

r=u/2+0.5;

r-=0.5;

if (r!=0.0) u=r;

beta=U/u;

p=(int)(-log(u)/log(beta)+0.5);

printf("Type double: base is %g, precision %d\n",beta,p);

}

{

doubledouble ddU,ddR,u,r,beta;

ddU=new_doubledouble(4.0,0.0);

ddU/=new_doubledouble(3.0,0.0);

ddU-=1;

ddU*=3;

ddU-=1;

ddU=ddfabs(ddU);

ddR=ddU/2+1;

ddR-=1;

if (ddR.hi!=0.0) ddU=ddR;

u=new_doubledouble(2.0,0.0);

u/=3.0;

u-=0.5;

u*=3;

u-=0.5;

u=ddfabs(u);

r=u/2+0.5;

r-=0.5;

if (r.hi!=0.0) u=r;

beta=ddU/u;

p=(int)(-ddlog(u)/ddlog(beta)+0.5).hi;

printf("Type doubledouble: base is %g, precision %d\n",beta.hi,p);

}

END_x86_FIX;

}

// Assumption: if it works, don't even bother doing output

int Verify(char msg[], doubledouble correct, doubledouble computed)

{

// default fractional error allowed.

doubledouble delta = new_doubledouble(1.1e-31,0.0);

if(correct == computed)

return 1;

else if(correct.hi == 0.0) // can't do a relative test

{

error:

printf("%s", msg);

printf(" corret.hi %g correct.lo %g", correct.hi,correct.lo);

printf("\ngot hi %g, lo %g\n",computed.hi,computed.lo);

return 0;

}

else if(ddfabs((correct-computed)/correct) <= delta)

return 1;

else goto error;

}

int main(void)

{

doubledouble a = new_doubledouble(5.0,0.0);

doubledouble b = new_doubledouble(4.0,0.0);

doubledouble c = a*b;

double x1,x2;

char buffer[256];

int i;

InitConstants();

base_and_prec();

doubledoubletoa(ddconst.Log2,buffer);

printf("%s\n",buffer);

x1=(pow(2,53)-2)/pow(2,53);

x2=(3.0)/pow(2,55);

Verify("Test -2: should be 0: ",

new_doubledouble(0.0,0.0),

ddfloor(new_doubledouble(x1,0.0)+new_doubledouble(x2,0.0)));

Verify("Test -1: should be 0: ", new_doubledouble(0.0,0.0),

ddfloor(atodd("0.9999999999999997779553950749686919152737")

atodd("0.00000000000000008326672684688674053177237510681152343750")));

for (i=0; i<5;i++) {

c = c*c;

doubledoubletoa(c,buffer);

printf("%s\n",buffer);

a = c - 1.0;

b = new_doubledouble(-1.0,0.0);

a = c + b;

doubledoubletoa(a,buffer);

printf("%s\n",buffer);

}

return 0;

}