Two_Sum(_l, _i, _m, _2); \  Two_Sum(_1, _k, _i, x4); \  Two_Sum(_2, _i, _k, x5); \  Two_Sum(_m, _k, x7, x6)/* An expansion of length two can be squared more quickly than finding the   *//*   product of two different expansions of length two, and the result is    *//*   guaranteed to have no more than six (rather than eight) components.     */#define Two_Square(a1, a0, x5, x4, x3, x2, x1, x0) \  Square(a0, _j, x0); \  _0 = a0 + a0; \  Two_Product(a1, _0, _k, _1); \  Two_One_Sum(_k, _1, _j, _l, _2, x1); \  Square(a1, _j, _1); \  Two_Two_Sum(_j, _1, _l, _2, x5, x4, x3, x2)/* splitter = 2^ceiling(p / 2) + 1.  Used to split floats in half.           */static REAL splitter;static REAL epsilon;         /* = 2^(-p).  Used to estimate roundoff errors. *//* A set of coefficients used to calculate maximum roundoff errors.          */static REAL resulterrbound;static REAL ccwerrboundA, ccwerrboundB, ccwerrboundC;static REAL o3derrboundA, o3derrboundB, o3derrboundC;static REAL iccerrboundA, iccerrboundB, iccerrboundC;static REAL isperrboundA, isperrboundB, isperrboundC;/*****************************************************************************//*                                                                           *//*  doubleprint()   Print the bit representation of a double.                *//*                                                                           *//*  Useful for debugging exact arithmetic routines.                          *//*                                                                           *//*****************************************************************************//*void doubleprint(number)double number;{  unsigned long long no;  unsigned long long sign, expo;  int exponent;  int i, bottomi;  no = *(unsigned long long *) &number;  sign = no & 0x8000000000000000ll;  expo = (no >> 52) & 0x7ffll;  exponent = (int) expo;  exponent = exponent - 1023;  if (sign) {    printf("-");  } else {    printf(" ");  }  if (exponent == -1023) {    printf(      "0.0000000000000000000000000000000000000000000000000000_     (   )");  } else {    printf("1.");    bottomi = -1;    for (i = 0; i < 52; i++) {      if (no & 0x0008000000000000ll) {        printf("1");        bottomi = i;      } else {        printf("0");      }      no <<= 1;    }    printf("_%d  (%d)", exponent, exponent - 1 - bottomi);  }}*//*****************************************************************************//*                                                                           *//*  floatprint()   Print the bit representation of a float.                  *//*                                                                           *//*  Useful for debugging exact arithmetic routines.                          *//*                                                                           *//*****************************************************************************//*void floatprint(number)float number;{  unsigned no;  unsigned sign, expo;  int exponent;  int i, bottomi;  no = *(unsigned *) &number;  sign = no & 0x80000000;  expo = (no >> 23) & 0xff;  exponent = (int) expo;  exponent = exponent - 127;  if (sign) {    printf("-");  } else {    printf(" ");  }  if (exponent == -127) {    printf("0.00000000000000000000000_     (   )");  } else {    printf("1.");    bottomi = -1;    for (i = 0; i < 23; i++) {      if (no & 0x00400000) {        printf("1");        bottomi = i;      } else {        printf("0");      }      no <<= 1;    }    printf("_%3d  (%3d)", exponent, exponent - 1 - bottomi);  }}*//*****************************************************************************//*                                                                           *//*  expansion_print()   Print the bit representation of an expansion.        *//*                                                                           *//*  Useful for debugging exact arithmetic routines.                          *//*                                                                           *//*****************************************************************************//*void expansion_print(elen, e)int elen;REAL *e;{  int i;  for (i = elen - 1; i >= 0; i--) {    REALPRINT(e[i]);    if (i > 0) {      printf(" +\n");    } else {      printf("\n");    }  }}*//*****************************************************************************//*                                                                           *//*  doublerand()   Generate a double with random 53-bit significand and a    *//*                 random exponent in [0, 511].                              *//*                                                                           *//*****************************************************************************//*double doublerand(){  double result;  double expo;  long a, b, c;  long i;  a = random();  b = random();  c = random();  result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8);  for (i = 512, expo = 2; i <= 131072; i *= 2, expo = expo * expo) {    if (c & i) {      result *= expo;    }  }  return result;}*//*****************************************************************************//*                                                                           *//*  narrowdoublerand()   Generate a double with random 53-bit significand    *//*                       and a random exponent in [0, 7].                    *//*                                                                           *//*****************************************************************************//*double narrowdoublerand(){  double result;  double expo;  long a, b, c;  long i;  a = random();  b = random();  c = random();  result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8);  for (i = 512, expo = 2; i <= 2048; i *= 2, expo = expo * expo) {    if (c & i) {      result *= expo;    }  }  return result;}*//*****************************************************************************//*                                                                           *//*  uniformdoublerand()   Generate a double with random 53-bit significand.  *//*                                                                           *//*****************************************************************************//*double uniformdoublerand(){  double result;  long a, b;  a = random();  b = random();  result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8);  return result;}*//*****************************************************************************//*                                                                           *//*  floatrand()   Generate a float with random 24-bit significand and a      *//*                random exponent in [0, 63].                                *//*                                                                           *//*****************************************************************************//*float floatrand(){  float result;  float expo;  long a, c;  long i;  a = random();  c = random();  result = (float) ((a - 1073741824) >> 6);  for (i = 512, expo = 2; i <= 16384; i *= 2, expo = expo * expo) {    if (c & i) {      result *= expo;    }  }  return result;}*//*****************************************************************************//*                                                                           *//*  narrowfloatrand()   Generate a float with random 24-bit significand and  *//*                      a random exponent in [0, 7].                         *//*                                                                           *//*****************************************************************************//*float narrowfloatrand(){  float result;  float expo;  long a, c;  long i;  a = random();  c = random();  result = (float) ((a - 1073741824) >> 6);  for (i = 512, expo = 2; i <= 2048; i *= 2, expo = expo * expo) {    if (c & i) {      result *= expo;    }  }  return result;}*//*****************************************************************************//*                                                                           *//*  uniformfloatrand()   Generate a float with random 24-bit significand.    *//*                                                                           *//*****************************************************************************//*float uniformfloatrand(){  float result;  long a;  a = random();  result = (float) ((a - 1073741824) >> 6);  return result;}*//*****************************************************************************//*                                                                           *//*  exactinit()   Initialize the variables used for exact arithmetic.        *//*                                                                           *//*  `epsilon' is the largest power of two such that 1.0 + epsilon = 1.0 in   *//*  floating-point arithmetic.  `epsilon' bounds the relative roundoff       *//*  error.  It is used for floating-point error analysis.                    *//*                                                                           *//*  `splitter' is used to split floating-point numbers into two half-        *//*  length significands for exact multiplication.                            *//*                                                                           *//*  I imagine that a highly optimizing compiler might be too smart for its   *//*  own good, and somehow cause this routine to fail, if it pretends that    *//*  floating-point arithmetic is too much like real arithmetic.              *//*                                                                           *//*  Don't change this routine unless you fully understand it.                *//*                                                                           *//*****************************************************************************/REAL exactinit(){  REAL half;  REAL check, lastcheck;  int every_other;#ifdef LINUX  int cword;#endif /* LINUX */#ifdef CPU86#ifdef SINGLE  _control87(_PC_24, _MCW_PC); /* Set FPU control word for single precision. */#else /* not SINGLE */  _control87(_PC_53, _MCW_PC); /* Set FPU control word for double precision. */#endif /* not SINGLE */#endif /* CPU86 */#ifdef LINUX#ifdef SINGLE  /*  cword = 4223; */  cword = 4210;                 /* set FPU control word for single precision */#else /* not SINGLE */  /*  cword = 4735; */  cword = 4722;                 /* set FPU control word for double precision */#endif /* not SINGLE */  _FPU_SETCW(cword);#endif /* LINUX */  every_other = 1;  half = 0.5;  epsilon = 1.0;  splitter = 1.0;  check = 1.0;  /* Repeatedly divide `epsilon' by two until it is too small to add to    */  /*   one without causing roundoff.  (Also check if the sum is equal to   */  /*   the previous sum, for machines that round up instead of using exact */  /*   rounding.  Not that this library will work on such machines anyway. */  do {    lastcheck = check;