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

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
Copyright (C) 1997 Keith Martin Briggs

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. Us