Almond/libmandel/qd-2.3.22/src/qd_real.cpp

/*
 * src/qd_real.cc
 *
 * This work was supported by the Director, Office of Science, Division
 * of Mathematical, Information, and Computational Sciences of the
 * U.S. Department of Energy under contract number DE-AC03-76SF00098.
 *
 * Copyright (c) 2000-2007
 *
 * Contains implementation of non-inlined functions of quad-double
 * package.  Inlined functions are found in qd_inline.h (in include directory).
 */
#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <iomanip>
#include <string>

#include "config.h"
#include <qd/qd_real.h>
#include "util.h"

#include <qd/bits.h>

#ifndef QD_INLINE
#include <qd/qd_inline.h>
#endif

using std::cout;
using std::cerr;
using std::endl;
using std::istream;
using std::ostream;
using std::ios_base;
using std::string;
using std::setw;

using namespace qd;

void qd_real::error(const char *msg) {
  if (msg) { cerr << "ERROR " << msg << endl; }
}

/********** Multiplications **********/

qd_real nint(const qd_real &a) {
  double x0, x1, x2, x3;

  x0 = nint(a[0]);
  x1 = x2 = x3 = 0.0;

  if (x0 == a[0]) {
    /* First double is already an integer. */
    x1 = nint(a[1]);

    if (x1 == a[1]) {
      /* Second double is already an integer. */
      x2 = nint(a[2]);
      
      if (x2 == a[2]) {
        /* Third double is already an integer. */
        x3 = nint(a[3]);
      } else {
        if (std::abs(x2 - a[2]) == 0.5 && a[3] < 0.0) {
          x2 -= 1.0;
        }
      }

    } else {
      if (std::abs(x1 - a[1]) == 0.5 && a[2] < 0.0) {
          x1 -= 1.0;
      }
    }

  } else {
    /* First double is not an integer. */
      if (std::abs(x0 - a[0]) == 0.5 && a[1] < 0.0) {
          x0 -= 1.0;
      }
  }
  
  renorm(x0, x1, x2, x3);
  return qd_real(x0, x1, x2, x3);
}

qd_real floor(const qd_real &a) {
  double x0, x1, x2, x3;
  x1 = x2 = x3 = 0.0;
  x0 = std::floor(a[0]);

  if (x0 == a[0]) {
    x1 = std::floor(a[1]);
    
    if (x1 == a[1]) {
      x2 = std::floor(a[2]);

      if (x2 == a[2]) {
        x3 = std::floor(a[3]);
      }
    }

    renorm(x0, x1, x2, x3);
    return qd_real(x0, x1, x2, x3);
  }

  return qd_real(x0, x1, x2, x3);
}

qd_real ceil(const qd_real &a) {
  double x0, x1, x2, x3;
  x1 = x2 = x3 = 0.0;
  x0 = std::ceil(a[0]);

  if (x0 == a[0]) {
    x1 = std::ceil(a[1]);
    
    if (x1 == a[1]) {
      x2 = std::ceil(a[2]);

      if (x2 == a[2]) {
        x3 = std::ceil(a[3]);
      }
    }

    renorm(x0, x1, x2, x3);
    return qd_real(x0, x1, x2, x3);
  }

  return qd_real(x0, x1, x2, x3);
}



/********** Divisions **********/
/* quad-double / double */
qd_real operator/(const qd_real &a, double b) {
  /* Strategy:  compute approximate quotient using high order
     doubles, and then correct it 3 times using the remainder.
     (Analogous to long division.)                             */
  double t0, t1;
  double q0, q1, q2, q3;
  qd_real r;

  q0 = a[0] / b;  /* approximate quotient */

  /* Compute the remainder  a - q0 * b */
  t0 = two_prod(q0, b, t1);
  r = a - dd_real(t0, t1);

  /* Compute the first correction */
  q1 = r[0] / b;
  t0 = two_prod(q1, b, t1);
  r -= dd_real(t0, t1);

  /* Second correction to the quotient. */
  q2 = r[0] / b;
  t0 = two_prod(q2, b, t1);
  r -= dd_real(t0, t1);

  /* Final correction to the quotient. */
  q3 = r[0] / b;

  renorm(q0, q1, q2, q3);
  return qd_real(q0, q1, q2, q3);
}

qd_real::qd_real(const char *s) {
  if (qd_real::read(s, *this)) {
    qd_real::error("(qd_real::qd_real): INPUT ERROR.");
    *this = qd_real::_nan;
  }
}

qd_real &qd_real::operator=(const char *s) {
  if (qd_real::read(s, *this)) {
    qd_real::error("(qd_real::operator=): INPUT ERROR.");
    *this = qd_real::_nan;
  }
  return *this;
}

istream &operator>>(istream &s, qd_real &qd) {
  char str[255];
  s >> str;
  qd = qd_real(str);
  return s;
}

ostream &operator<<(ostream &os, const qd_real &qd) {
  bool showpos = (os.flags() & ios_base::showpos) != 0;
  bool uppercase = (os.flags() & ios_base::uppercase) != 0;
  return os << qd.to_string(os.precision(), os.width(), os.flags(), 
      showpos, uppercase, os.fill());
}

/* Read a quad-double from s. */
int qd_real::read(const char *s, qd_real &qd) {
  const char *p = s;
  char ch;
  int sign = 0;
  int point = -1;  /* location of decimal point */
  int nd = 0;      /* number of digits read */
  int e = 0;       /* exponent. */
  bool done = false;
  qd_real r = 0.0;  /* number being read */

  /* Skip any leading spaces */
  while (*p == ' ') p++;

  while (!done && (ch = *p) != '\0') {
    if (ch >= '0' && ch <= '9') {
      /* It's a digit */
      int d = ch - '0';
      r *= 10.0;
      r += static_cast<double>(d);
      nd++;
    } else {
      /* Non-digit */
      switch (ch) {
      case '.':
        if (point >= 0)
          return -1;   /* we've already encountered a decimal point. */
        point = nd;
        break;
      case '-':
      case '+':
        if (sign != 0 || nd > 0)
          return -1;  /* we've already encountered a sign, or if its
                            not at first position. */
        sign = (ch == '-') ? -1 : 1;
        break;
      case 'E':
      case 'e':
        int nread;
        nread = std::sscanf(p+1, "%d", &e);
        done = true;
        if (nread != 1)
          return -1;  /* read of exponent failed. */
        break;
      case ' ':
        done = true;
        break;
      default:
        return -1;
      
      }
    }

    p++;
  }



  /* Adjust exponent to account for decimal point */
  if (point >= 0) {
    e -= (nd - point);
  }

  /* Multiply the the exponent */
  if (e != 0) {
    r *= (qd_real(10.0) ^ e);
  }

  qd = (sign < 0) ? -r : r;
  return 0;
}

void qd_real::to_digits(char *s, int &expn, int precision) const {
  int D = precision + 1;  /* number of digits to compute */

  qd_real r = abs(*this);
  int e;  /* exponent */
  int i, d;

  if (x[0] == 0.0) {
    /* this == 0.0 */
    expn = 0;
    for (i = 0; i < precision; i++) s[i] = '0';
    return;
  }

  /* First determine the (approximate) exponent. */
  e = static_cast<int>(std::floor(std::log10(std::abs(x[0]))));

  if (e < -300) {
    r *= qd_real(10.0) ^ 300;
    r /= qd_real(10.0) ^ (e + 300);
  } else if (e > 300) {
    r = ldexp(r, -53);
    r /= qd_real(10.0) ^ e;
    r = ldexp(r, 53);
  } else {
    r /= qd_real(10.0) ^ e;
  }

  /* Fix exponent if we are off by one */
  if (r >= 10.0) {
    r /= 10.0;
    e++;
  } else if (r < 1.0) {
    r *= 10.0;
    e--;
  }

  if (r >= 10.0 || r < 1.0) {
    qd_real::error("(qd_real::to_digits): can't compute exponent.");
    return;
  }

  /* Extract the digits */
  for (i = 0; i < D; i++) {
    d = static_cast<int>(r[0]);
    r -= d;
    r *= 10.0;

    s[i] = static_cast<char>(d + '0');
  }

  /* Fix out of range digits. */
  for (i = D-1; i > 0; i--) {
    if (s[i] < '0') {
      s[i-1]--;
      s[i] += 10;
    } else if (s[i] > '9') {
      s[i-1]++;
      s[i] -= 10;
    }
  }

  if (s[0] <= '0') {
    qd_real::error("(qd_real::to_digits): non-positive leading digit.");
    return;
  }

  /* Round, handle carry */
  if (s[D-1] >= '5') {
    s[D-2]++;

    i = D-2;
    while (i > 0 && s[i] > '9') {
      s[i] -= 10;
      s[--i]++;
    }
  }

  /* If first digit is 10, shift everything. */
  if (s[0] > '9') { 
    e++; 
    for (i = precision; i >= 2; i--) s[i] = s[i-1]; 
    s[0] = '1';
    s[1] = '0';
  }

  s[precision] = 0;
  expn = e;
}

/* Writes the quad-double number into the character array s of length len.
   The integer d specifies how many significant digits to write.
   The string s must be able to hold at least (d+8) characters.  
   showpos indicates whether to use the + sign, and uppercase indicates
   whether the E or e is to be used for the exponent. */
void qd_real::write(char *s, int len, int precision, 
    bool showpos, bool uppercase) const {
  string str = to_string(precision, 0, ios_base::scientific, showpos, uppercase);
  strncpy(s, str.c_str(), len-1);
  s[len-1] = 0;
}

void round_string_qd(char *s, int precision, int *offset){
	/*
	 Input string must be all digits or errors will occur.
	 */

	int i;
	int D = precision ;

	/* Round, handle carry */
	  if (D>0 && s[D] >= '5') {
	    s[D-1]++;

	    i = D-1;
	    while (i > 0 && s[i] > '9') {
	      s[i] -= 10;
	      s[--i]++;
	    }
	  }

	  /* If first digit is 10, shift everything. */
	  if (s[0] > '9') {
	    // e++; // don't modify exponent here
	    for (i = precision; i >= 1; i--) s[i+1] = s[i];
	    s[0] = '1';
	    s[1] = '0';

	    (*offset)++ ; // now offset needs to be increased by one
	    precision++ ;
	  }

	  s[precision] = 0; // add terminator for array
}


string qd_real::to_string(int precision, int width, ios_base::fmtflags fmt, 
    bool showpos, bool uppercase, char fill) const {
  string s;
  bool fixed = (fmt & ios_base::fixed) != 0;
  bool sgn = true;
  int i, e = 0;

  if (isinf()) {
    if (*this < 0.0)
      s += '-';
    else if (showpos)
      s += '+';
    else
      sgn = false;
    s += uppercase ? "INF" : "inf";
  } else if (isnan()) {
    s = uppercase ? "NAN" : "nan";
    sgn = false;
  } else {
    if (*this < 0.0)
      s += '-';
    else if (showpos)
      s += '+';
    else
      sgn = false;

    if (*this == 0.0) {
      /* Zero case */
      s += '0';
      if (precision > 0) {
        s += '.';
        s.append(precision, '0');
      }
    } else {
      /* Non-zero case */
      int off = (fixed ? (1 + to_int(floor(log10(abs(*this))))) : 1);
      int d = precision + off;

      int d_with_extra = d;
      if(fixed)
    	  d_with_extra = std::max(120, d); // longer than the max accuracy for DD

      // highly special case - fixed mode, precision is zero, abs(*this) < 1.0
      // without this trap a number like 0.9 printed fixed with 0 precision prints as 0
      // should be rounded to 1.
      if(fixed && (precision == 0) && (abs(*this) < 1.0)){
    	  if(abs(*this) >= 0.5)
    		  s += '1';
    	  else
    		  s += '0';

    	  return s;
      }

      // handle near zero to working precision (but not exactly zero)
      if (fixed && d <= 0) {
        s += '0';
        if (precision > 0) {
          s += '.';
          s.append(precision, '0');
        }
      } else {  // default

        char *t ; // = new char[d+1];
        int j;

        if(fixed){
        	t = new char[d_with_extra+1];
        	to_digits(t, e, d_with_extra);
        }
        else{
        	t = new char[d+1];
        	to_digits(t, e, d);
        }

        off = e + 1;

        if (fixed) {
          // fix the string if it's been computed incorrectly
          // round here in the decimal string if required
          round_string_qd(t, d, &off);

          if (off > 0) {
            for (i = 0; i < off; i++) s += t[i];
            if (precision > 0) {
              s += '.';
              for (j = 0; j < precision; j++, i++) s += t[i];
            }
          } else {
            s += "0.";
            if (off < 0) s.append(-off, '0');
            for (i = 0; i < d; i++) s += t[i];
          }
        } else {
          s += t[0];
          if (precision > 0) s += '.';

          for (i = 1; i <= precision; i++)
            s += t[i];

        }
		delete [] t;
      }
    }

    // trap for improper offset with large values
    // without this trap, output of values of the for 10^j - 1 fail for j > 28
    // and are output with the point in the wrong place, leading to a dramatically off value
    if(fixed && (precision > 0)){
    	// make sure that the value isn't dramatically larger
    	double from_string = atof(s.c_str());

    	// if this ratio is large, then we've got problems
    	if( fabs( from_string / this->x[0] ) > 3.0 ){

    		int point_position;
    		char temp;

    		// loop on the string, find the point, move it up one
    		// don't act on the first character
    		for(i=1; i < s.length(); i++){
    			if(s[i] == '.'){
    				s[i] = s[i-1] ;
    				s[i-1] = '.' ;
    				break;
    			}
    		}

        	from_string = atof(s.c_str());
        	// if this ratio is large, then the string has not been fixed
        	if( fabs( from_string / this->x[0] ) > 3.0 ){
        		dd_real::error("Re-rounding unsuccessful in large number fixed point trap.") ;
        	}
    	}
    }

    if (!fixed) {
      /* Fill in exponent part */
      s += uppercase ? 'E' : 'e';
      append_expn(s, e);
    }
  }

  /* Fill in the blanks */
  int len = s.length();
  if (len < width) {
    int delta = width - len;
    if (fmt & ios_base::internal) {
      if (sgn)
        s.insert(static_cast<string::size_type>(1), delta, fill);
      else
        s.insert(static_cast<string::size_type>(0), delta, fill);
    } else if (fmt & ios_base::left) {
      s.append(delta, fill);
    } else {
      s.insert(static_cast<string::size_type>(0), delta, fill);
    }
  }

  return s;
}

/* Computes  qd^n, where n is an integer. */
qd_real pow(const qd_real &a, int n) {
  if (n == 0)
    return 1.0;

  qd_real r = a;   /* odd-case multiplier */
  qd_real s = 1.0;  /* current answer */
  int N = std::abs(n);

  if (N > 1) {

    /* Use binary exponentiation. */
    while (N > 0) {
      if (N % 2 == 1) {
        /* If odd, multiply by r. Note eventually N = 1, so this
         eventually executes. */
        s *= r;
      }
      N /= 2;
      if (N > 0)
        r = sqr(r);
    }

  } else {
    s = r;
  }

  if (n < 0)
    return (1.0 / s);

  return s;
}

qd_real pow(const qd_real &a, const qd_real &b) {
  return exp(b * log(a));
}

qd_real npwr(const qd_real &a, int n) {
  return pow(a, n);
}

/* Debugging routines */
void qd_real::dump_bits(const string &name, std::ostream &os) const {
  string::size_type len = name.length();
  if (len > 0) {
    os << name << " = ";
    len += 3;
  }
  os << "[ ";
  len += 2;
  for (int j = 0; j < 4; j++) {
    if (j > 0) for (string::size_type i = 0; i < len; i++) os << ' ';
    print_double_info(os, x[j]);
    if (j < 3)
      os << endl;
    else
      os << " ]" << endl;
  }
}

void qd_real::dump(const string &name, std::ostream &os) const {
  std::ios_base::fmtflags old_flags = os.flags();
  std::streamsize old_prec = os.precision(19);
  os << std::scientific;

  string::size_type len = name.length();
  if (len > 0) {
    os << name << " = ";
    len += 3;
  }
  os << "[ ";
  len += 2;
  os << setw(27) << x[0] << ", " << setw(26) << x[1] << "," << endl;
  for (string::size_type i = 0; i < len; i++) os << ' ';
  os << setw(27) << x[2] << ", " << setw(26) << x[3] << "  ]" << endl;

  os.precision(old_prec);
  os.flags(old_flags);
}

/* Divisions */
/* quad-double / double-double */
qd_real qd_real::sloppy_div(const qd_real &a, const dd_real &b) {
  double q0, q1, q2, q3;
  qd_real r;
  qd_real qd_b(b);

  q0 = a[0] / b._hi();
  r = a - q0 * qd_b;

  q1 = r[0] / b._hi();
  r -= (q1 * qd_b);

  q2 = r[0] / b._hi();
  r -= (q2 * qd_b);

  q3 = r[0] / b._hi();

  ::renorm(q0, q1, q2, q3);
  return qd_real(q0, q1, q2, q3);
}

qd_real qd_real::accurate_div(const qd_real &a, const dd_real &b) {
  double q0, q1, q2, q3, q4;
  qd_real r;
  qd_real qd_b(b);

  q0 = a[0] / b._hi();
  r = a - q0 * qd_b;

  q1 = r[0] / b._hi();
  r -= (q1 * qd_b);

  q2 = r[0] / b._hi();
  r -= (q2 * qd_b);

  q3 = r[0] / b._hi();
  r -= (q3 * qd_b);

  q4 = r[0] / b._hi();

  ::renorm(q0, q1, q2, q3, q4);
  return qd_real(q0, q1, q2, q3);
}

/* quad-double / quad-double */
qd_real qd_real::sloppy_div(const qd_real &a, const qd_real &b) {
  double q0, q1, q2, q3;

  qd_real r;

  q0 = a[0] / b[0];
  r = a - (b * q0);

  q1 = r[0] / b[0];
  r -= (b * q1);

  q2 = r[0] / b[0];
  r -= (b * q2);

  q3 = r[0] / b[0];

  ::renorm(q0, q1, q2, q3);

  return qd_real(q0, q1, q2, q3);
}

qd_real qd_real::accurate_div(const qd_real &a, const qd_real &b) {
  double q0, q1, q2, q3;

  qd_real r;

  q0 = a[0] / b[0];
  r = a - (b * q0);

  q1 = r[0] / b[0];
  r -= (b * q1);

  q2 = r[0] / b[0];
  r -= (b * q2);

  q3 = r[0] / b[0];

  r -= (b * q3);
  double q4 = r[0] / b[0];

  ::renorm(q0, q1, q2, q3, q4);

  return qd_real(q0, q1, q2, q3);
}

QD_API qd_real sqrt(const qd_real &a) {
  /* Strategy:  

     Perform the following Newton iteration:

       x' = x + (1 - a * x^2) * x / 2;
       
     which converges to 1/sqrt(a), starting with the
     double precision approximation to 1/sqrt(a).
     Since Newton's iteration more or less doubles the
     number of correct digits, we only need to perform it 
     twice.
  */

  if (a.is_zero())
    return 0.0;

  if (a.is_negative()) {
    qd_real::error("(qd_real::sqrt): Negative argument.");
    return qd_real::_nan;
  }

  qd_real r = (1.0 / std::sqrt(a[0]));
  qd_real h = mul_pwr2(a, 0.5);

  r += ((0.5 - h * sqr(r)) * r);
  r += ((0.5 - h * sqr(r)) * r);
  r += ((0.5 - h * sqr(r)) * r);

  r *= a;
  return r;
}


/* Computes the n-th root of a */
qd_real nroot(const qd_real &a, int n) {
  /* Strategy:  Use Newton's iteration to solve
     
        1/(x^n) - a = 0

     Newton iteration becomes

        x' = x + x * (1 - a * x^n) / n

     Since Newton's iteration converges quadratically, 
     we only need to perform it twice.

   */
	if (n <= 0) {
		qd_real::error("(qd_real::nroot): N must be positive.");
		return qd_real::_nan;
	}

	if (n % 2 == 0 && a.is_negative()) {
		qd_real::error("(qd_real::nroot): Negative argument.");
		return qd_real::_nan;
	}

	if (n == 1) {
		return a;
	}
	if (n == 2) {
		return sqrt(a);
	}
	if (a.is_zero()) {
		return qd_real(0.0);
	}


	/* Note  a^{-1/n} = exp(-log(a)/n) */
	qd_real r = abs(a);
	qd_real x = std::exp(-std::log(r.x[0]) / n);

	/* Perform Newton's iteration. */
	double dbl_n = static_cast<double>(n);
	x += x * (1.0 - r * npwr(x, n)) / dbl_n;
	x += x * (1.0 - r * npwr(x, n)) / dbl_n;
	x += x * (1.0 - r * npwr(x, n)) / dbl_n;
	if (a[0] < 0.0){
		x = -x;
	}
	return 1.0 / x;
}

static const int n_inv_fact = 15;
static const qd_real inv_fact[n_inv_fact] = {
  qd_real( 1.66666666666666657e-01,  9.25185853854297066e-18,
           5.13581318503262866e-34,  2.85094902409834186e-50),
  qd_real( 4.16666666666666644e-02,  2.31296463463574266e-18,
           1.28395329625815716e-34,  7.12737256024585466e-51),
  qd_real( 8.33333333333333322e-03,  1.15648231731787138e-19,
           1.60494162032269652e-36,  2.22730392507682967e-53),
  qd_real( 1.38888888888888894e-03, -5.30054395437357706e-20,
          -1.73868675534958776e-36, -1.63335621172300840e-52),
  qd_real( 1.98412698412698413e-04,  1.72095582934207053e-22,
           1.49269123913941271e-40,  1.29470326746002471e-58),
  qd_real( 2.48015873015873016e-05,  2.15119478667758816e-23,
           1.86586404892426588e-41,  1.61837908432503088e-59),
  qd_real( 2.75573192239858925e-06, -1.85839327404647208e-22,
           8.49175460488199287e-39, -5.72661640789429621e-55),
  qd_real( 2.75573192239858883e-07,  2.37677146222502973e-23,
          -3.26318890334088294e-40,  1.61435111860404415e-56),
  qd_real( 2.50521083854417202e-08, -1.44881407093591197e-24,
           2.04267351467144546e-41, -8.49632672007163175e-58),
  qd_real( 2.08767569878681002e-09, -1.20734505911325997e-25,
           1.70222792889287100e-42,  1.41609532150396700e-58),
  qd_real( 1.60590438368216133e-10,  1.25852945887520981e-26,
          -5.31334602762985031e-43,  3.54021472597605528e-59),
  qd_real( 1.14707455977297245e-11,  2.06555127528307454e-28,
           6.88907923246664603e-45,  5.72920002655109095e-61),
  qd_real( 7.64716373181981641e-13,  7.03872877733453001e-30,
          -7.82753927716258345e-48,  1.92138649443790242e-64),
  qd_real( 4.77947733238738525e-14,  4.39920548583408126e-31,
          -4.89221204822661465e-49,  1.20086655902368901e-65),
  qd_real( 2.81145725434552060e-15,  1.65088427308614326e-31,
          -2.87777179307447918e-50,  4.27110689256293549e-67)
};

qd_real exp(const qd_real &a) {
  /* Strategy:  We first reduce the size of x by noting that
     
          exp(kr + m * log(2)) = 2^m * exp(r)^k

     where m and k are integers.  By choosing m appropriately
     we can make |kr| <= log(2) / 2 = 0.347.  Then exp(r) is 
     evaluated using the familiar Taylor series.  Reducing the 
     argument substantially speeds up the convergence.       */  

  const double k = ldexp(1.0, 16);
  const double inv_k = 1.0 / k;

  if (a[0] <= -709.0)
    return 0.0;

  if (a[0] >=  709.0)
    return qd_real::_inf;

  if (a.is_zero())
    return 1.0;

  if (a.is_one())
    return qd_real::_e;

  double m = std::floor(a.x[0] / qd_real::_log2.x[0] + 0.5);
  qd_real r = mul_pwr2(a - qd_real::_log2 * m, inv_k);
  qd_real s, p, t;
  double thresh = inv_k * qd_real::_eps;

  p = sqr(r);
  s = r + mul_pwr2(p, 0.5);
  int i = 0;
  do {
    p *= r;
    t = p * inv_fact[i++];
    s += t;
  } while (std::abs(to_double(t)) > thresh && i < 9);

  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s = mul_pwr2(s, 2.0) + sqr(s);
  s += 1.0;
  return ldexp(s, static_cast<int>(m));
}

/* Logarithm.  Computes log(x) in quad-double precision.
   This is a natural logarithm (i.e., base e).            */
qd_real log(const qd_real &a) {
  /* Strategy.  The Taylor series for log converges much more
     slowly than that of exp, due to the lack of the factorial
     term in the denominator.  Hence this routine instead tries
     to determine the root of the function

         f(x) = exp(x) - a

     using Newton iteration.  The iteration is given by

         x' = x - f(x)/f'(x) 
            = x - (1 - a * exp(-x))
            = x + a * exp(-x) - 1.
           
     Two iteration is needed, since Newton's iteration 
     approximately doubles the number of digits per iteration. */

  if (a.is_one()) {
    return 0.0;
  }

  if (a[0] <= 0.0) {
    qd_real::error("(qd_real::log): Non-positive argument.");
    return qd_real::_nan;
  }

  if (a[0] == 0.0) {
    return -qd_real::_inf;
  }

  qd_real x = std::log(a[0]);   /* Initial approximation */

  x = x + a * exp(-x) - 1.0;
  x = x + a * exp(-x) - 1.0;
  x = x + a * exp(-x) - 1.0;

  return x;
}

qd_real log10(const qd_real &a) {
  return log(a) / qd_real::_log10;
}

static const qd_real _pi1024 = qd_real(
    3.067961575771282340e-03, 1.195944139792337116e-19,
   -2.924579892303066080e-36, 1.086381075061880158e-52);

/* Table of sin(k * pi/1024) and cos(k * pi/1024). */
static const qd_real sin_table [] = {
  qd_real( 3.0679567629659761e-03, 1.2690279085455925e-19,
       5.2879464245328389e-36, -1.7820334081955298e-52),
  qd_real( 6.1358846491544753e-03, 9.0545257482474933e-20,
       1.6260113133745320e-37, -9.7492001208767410e-55),
  qd_real( 9.2037547820598194e-03, -1.2136591693535934e-19,
       5.5696903949425567e-36, 1.2505635791936951e-52),
  qd_real( 1.2271538285719925e-02, 6.9197907640283170e-19,
       -4.0203726713435555e-36, -2.0688703606952816e-52),
  qd_real( 1.5339206284988102e-02, -8.4462578865401696e-19,
       4.6535897505058629e-35, -1.3923682978570467e-51),
  qd_real( 1.8406729905804820e-02, 7.4195533812833160e-19,
       3.9068476486787607e-35, 3.6393321292898614e-52),
  qd_real( 2.1474080275469508e-02, -4.5407960207688566e-19,
       -2.2031770119723005e-35, 1.2709814654833741e-51),
  qd_real( 2.4541228522912288e-02, -9.1868490125778782e-20,
       4.8706148704467061e-36, -2.8153947855469224e-52),
  qd_real( 2.7608145778965743e-02, -1.5932358831389269e-18,
       -7.0475416242776030e-35, -2.7518494176602744e-51),
  qd_real( 3.0674803176636626e-02, -1.6936054844107918e-20,
       -2.0039543064442544e-36, -1.6267505108658196e-52),
  qd_real( 3.3741171851377587e-02, -2.0096074292368340e-18,
       -1.3548237016537134e-34, 6.5554881875899973e-51),
  qd_real( 3.6807222941358832e-02, 6.1060088803529842e-19,
       -4.0448721259852727e-35, -2.1111056765671495e-51),
  qd_real( 3.9872927587739811e-02, 4.6657453481183289e-19,
       3.4119333562288684e-35, 2.4007534726187511e-51),
  qd_real( 4.2938256934940820e-02, 2.8351940588660907e-18,
       1.6991309601186475e-34, 6.8026536098672629e-51),
  qd_real( 4.6003182130914630e-02, -1.1182813940157788e-18,
       7.5235020270378946e-35, 4.1187304955493722e-52),
  qd_real( 4.9067674327418015e-02, -6.7961037205182801e-19,
       -4.4318868124718325e-35, -9.9376628132525316e-52),
  qd_real( 5.2131704680283324e-02, -2.4243695291953779e-18,
       -1.3675405320092298e-34, -8.3938137621145070e-51),
  qd_real( 5.5195244349689941e-02, -1.3340299860891103e-18,
       -3.4359574125665608e-35, 1.1911462755409369e-51),
  qd_real( 5.8258264500435759e-02, 2.3299905496077492e-19,
       1.9376108990628660e-36, -5.1273775710095301e-53),
  qd_real( 6.1320736302208578e-02, -5.1181134064638108e-19,
       -4.2726335866706313e-35, 2.6368495557440691e-51),
  qd_real( 6.4382630929857465e-02, -4.2325997000052705e-18,
       3.3260117711855937e-35, 1.4736267706718352e-51),
  qd_real( 6.7443919563664065e-02, -6.9221796556983636e-18,
       1.5909286358911040e-34, -7.8828946891835218e-51),
  qd_real( 7.0504573389613870e-02, -6.8552791107342883e-18,
       -1.9961177630841580e-34, 2.0127129580485300e-50),
  qd_real( 7.3564563599667426e-02, -2.7784941506273593e-18,
       -9.1240375489852821e-35, -1.9589752023546795e-51),
  qd_real( 7.6623861392031492e-02, 2.3253700287958801e-19,
       -1.3186083921213440e-36, -4.9927872608099673e-53),
  qd_real( 7.9682437971430126e-02, -4.4867664311373041e-18,
       2.8540789143650264e-34, 2.8491348583262741e-51),
  qd_real( 8.2740264549375692e-02, 1.4735983530877760e-18,
       3.7284093452233713e-35, 2.9024430036724088e-52),
  qd_real( 8.5797312344439894e-02, -3.3881893830684029e-18,
       -1.6135529531508258e-34, 7.7294651620588049e-51),
  qd_real( 8.8853552582524600e-02, -3.7501775830290691e-18,
       3.7543606373911573e-34, 2.2233701854451859e-50),
  qd_real( 9.1908956497132724e-02, 4.7631594854274564e-18,
       1.5722874642939344e-34, -4.8464145447831456e-51),
  qd_real( 9.4963495329639006e-02, -6.5885886400417564e-18,
       -2.1371116991641965e-34, 1.3819370559249300e-50),
  qd_real( 9.8017140329560604e-02, -1.6345823622442560e-18,
       -1.3209238810006454e-35, -3.5691060049117942e-52),
  qd_real( 1.0106986275482782e-01, 3.3164325719308656e-18,
       -1.2004224885132282e-34, 7.2028828495418631e-51),
  qd_real( 1.0412163387205457e-01, 6.5760254085385100e-18,
       1.7066246171219214e-34, -4.9499340996893514e-51),
  qd_real( 1.0717242495680884e-01, 6.4424044279026198e-18,
       -8.3956976499698139e-35, -4.0667730213318321e-51),
  qd_real( 1.1022220729388306e-01, -5.6789503537823233e-19,
       1.0380274792383233e-35, 1.5213997918456695e-52),
  qd_real( 1.1327095217756435e-01, 2.7100481012132900e-18,
       1.5323292999491619e-35, 4.9564432810360879e-52),
  qd_real( 1.1631863091190477e-01, 1.0294914877509705e-18,
       -9.3975734948993038e-35, 1.3534827323719708e-52),
  qd_real( 1.1936521481099137e-01, -3.9500089391898506e-18,
       3.5317349978227311e-34, 1.8856046807012275e-51),
  qd_real( 1.2241067519921620e-01, 2.8354501489965335e-18,
       1.8151655751493305e-34, -2.8716592177915192e-51),
  qd_real( 1.2545498341154623e-01, 4.8686751763148235e-18,
       5.9878105258097936e-35, -3.3534629098722107e-51),
  qd_real( 1.2849811079379317e-01, 3.8198603954988802e-18,
       -1.8627501455947798e-34, -2.4308161133527791e-51),
  qd_real( 1.3154002870288312e-01, -5.0039708262213813e-18,
       -1.2983004159245552e-34, -4.6872034915794122e-51),
  qd_real( 1.3458070850712620e-01, -9.1670359171480699e-18,
       1.5916493007073973e-34, 4.0237002484366833e-51),
  qd_real( 1.3762012158648604e-01, 6.6253255866774482e-18,
       -2.3746583031401459e-34, -9.3703876173093250e-52),
  qd_real( 1.4065823933284924e-01, -7.9193932965524741e-18,
       6.0972464202108397e-34, 2.4566623241035797e-50),
  qd_real( 1.4369503315029444e-01, 1.1472723016618666e-17,
       -5.1884954557576435e-35, -4.2220684832186607e-51),
  qd_real( 1.4673047445536175e-01, 3.7269471470465677e-18,
       3.7352398151250827e-34, -4.0881822289508634e-51),
  qd_real( 1.4976453467732151e-01, 8.0812114131285151e-18,
       1.2979142554917325e-34, 9.9380667487736254e-51),
  qd_real( 1.5279718525844344e-01, -7.6313573938416838e-18,
       5.7714690450284125e-34, -3.7731132582986687e-50),
  qd_real( 1.5582839765426523e-01, 3.0351307187678221e-18,
       -1.0976942315176184e-34, 7.8734647685257867e-51),
  qd_real( 1.5885814333386145e-01, -4.0163200573859079e-18,
       -9.2840580257628812e-35, -2.8567420029274875e-51),
  qd_real( 1.6188639378011183e-01, 1.1850519643573528e-17,
       -5.0440990519162957e-34, 3.0510028707928009e-50),
  qd_real( 1.6491312048996992e-01, -7.0405288319166738e-19,
       3.3211107491245527e-35, 8.6663299254686031e-52),
  qd_real( 1.6793829497473117e-01, 5.4284533721558139e-18,
       -3.3263339336181369e-34, -1.8536367335123848e-50),
  qd_real( 1.7096188876030122e-01, 9.1919980181759094e-18,
       -6.7688743940982606e-34, -1.0377711384318389e-50),
  qd_real( 1.7398387338746382e-01, 5.8151994618107928e-18,
       -1.6751014298301606e-34, -6.6982259797164963e-51),
  qd_real( 1.7700422041214875e-01, 6.7329300635408167e-18,
       2.8042736644246623e-34, 3.6786888232793599e-51),
  qd_real( 1.8002290140569951e-01, 7.9701826047392143e-18,
       -7.0765920110524977e-34, 1.9622512608461784e-50),
  qd_real( 1.8303988795514095e-01, 7.7349918688637383e-18,
       -4.4803769968145083e-34, 1.1201148793328890e-50),
  qd_real( 1.8605515166344666e-01, -1.2564893007679552e-17,
       7.5953844248530810e-34, -3.8471695132415039e-51),
  qd_real( 1.8906866414980622e-01, -7.6208955803527778e-18,
       -4.4792298656662981e-34, -4.4136824096645007e-50),
  qd_real( 1.9208039704989244e-01, 4.3348343941174903e-18,
       -2.3404121848139937e-34, 1.5789970962611856e-50),
  qd_real( 1.9509032201612828e-01, -7.9910790684617313e-18,
       6.1846270024220713e-34, -3.5840270918032937e-50),
  qd_real( 1.9809841071795359e-01, -1.8434411800689445e-18,
       1.4139031318237285e-34, 1.0542811125343809e-50),
  qd_real( 2.0110463484209190e-01, 1.1010032669300739e-17,
       -3.9123576757413791e-34, 2.4084852500063531e-51),
  qd_real( 2.0410896609281687e-01, 6.0941297773957752e-18,
       -2.8275409970449641e-34, 4.6101008563532989e-51),
  qd_real( 2.0711137619221856e-01, -1.0613362528971356e-17,
       2.2456805112690884e-34, 1.3483736125280904e-50),
  qd_real( 2.1011183688046961e-01, 1.1561548476512844e-17,
       6.0355905610401254e-34, 3.3329909618405675e-50),
  qd_real( 2.1311031991609136e-01, 1.2031873821063860e-17,
       -3.4142699719695635e-34, -1.2436262780241778e-50),
  qd_real( 2.1610679707621952e-01, -1.0111196082609117e-17,
       7.2789545335189643e-34, -2.9347540365258610e-50),
  qd_real( 2.1910124015686980e-01, -3.6513812299150776e-19,
       -2.3359499418606442e-35, 3.1785298198458653e-52),
  qd_real( 2.2209362097320354e-01, -3.0337210995812162e-18,
       6.6654668033632998e-35, 2.0110862322656942e-51),
  qd_real( 2.2508391135979283e-01, 3.9507040822556510e-18,
       2.4287993958305375e-35, 5.6662797513020322e-52),
  qd_real( 2.2807208317088573e-01, 8.2361837339258012e-18,
       6.9786781316397937e-34, -6.4122962482639504e-51),
  qd_real( 2.3105810828067111e-01, 1.0129787149761869e-17,
       -6.9359234615816044e-34, -2.8877355604883782e-50),
  qd_real( 2.3404195858354343e-01, -6.9922402696101173e-18,
       -5.7323031922750280e-34, 5.3092579966872727e-51),
  qd_real( 2.3702360599436720e-01, 8.8544852285039918e-18,
       1.3588480826354134e-34, 1.0381022520213867e-50),
  qd_real( 2.4000302244874150e-01, -1.2137758975632164e-17,
       -2.6448807731703891e-34, -1.9929733800670473e-51),
  qd_real( 2.4298017990326390e-01, -8.7514315297196632e-18,
       -6.5723260373079431e-34, -1.0333158083172177e-50),
  qd_real( 2.4595505033579462e-01, -1.1129044052741832e-17,
       4.3805998202883397e-34, 1.2219399554686291e-50),
  qd_real( 2.4892760574572018e-01, -8.1783436100020990e-18,
       5.5666875261111840e-34, 3.8080473058748167e-50),
  qd_real( 2.5189781815421697e-01, -1.7591436032517039e-17,
       -1.0959681232525285e-33, 5.6209426020232456e-50),
  qd_real( 2.5486565960451457e-01, -1.3602299806901461e-19,
       -6.0073844642762535e-36, -3.0072751311893878e-52),
  qd_real( 2.5783110216215899e-01, 1.8480038630879957e-17,
       3.3201664714047599e-34, -5.5547819290576764e-51),
  qd_real( 2.6079411791527551e-01, 4.2721420983550075e-18,
       5.6782126934777920e-35, 3.1428338084365397e-51),
  qd_real( 2.6375467897483140e-01, -1.8837947680038700e-17,
       1.3720129045754794e-33, -8.2763406665966033e-50),
  qd_real( 2.6671275747489837e-01, 2.0941222578826688e-17,
       -1.1303466524727989e-33, 1.9954224050508963e-50),
  qd_real( 2.6966832557291509e-01, 1.5765657618133259e-17,
       -6.9696142173370086e-34, -4.0455346879146776e-50),
  qd_real( 2.7262135544994898e-01, 7.8697166076387850e-18,
       6.6179388602933372e-35, -2.7642903696386267e-51),
  qd_real( 2.7557181931095814e-01, 1.9320328962556582e-17,
       1.3932094180100280e-33, 1.3617253920018116e-50),
  qd_real( 2.7851968938505312e-01, -1.0030273719543544e-17,
       7.2592115325689254e-34, -1.0068516296655851e-50),
  qd_real( 2.8146493792575800e-01, -1.2322299641274009e-17,
       -1.0564788706386435e-34, 7.5137424251265885e-51),
  qd_real( 2.8440753721127182e-01, 2.2209268510661475e-17,
       -9.1823095629523708e-34, -5.2192875308892218e-50),
  qd_real( 2.8734745954472951e-01, 1.5461117367645717e-17,
       -6.3263973663444076e-34, -2.2982538416476214e-50),
  qd_real( 2.9028467725446239e-01, -1.8927978707774251e-17,
       1.1522953157142315e-33, 7.4738655654716596e-50),
  qd_real( 2.9321916269425863e-01, 2.2385430811901833e-17,
       1.3662484646539680e-33, -4.2451325253996938e-50),
  qd_real( 2.9615088824362384e-01, -2.0220736360876938e-17,
       -7.9252212533920413e-35, -2.8990577729572470e-51),
  qd_real( 2.9907982630804048e-01, 1.6701181609219447e-18,
       8.6091151117316292e-35, 3.9931286230012102e-52),
  qd_real( 3.0200594931922808e-01, -1.7167666235262474e-17,
       2.3336182149008069e-34, 8.3025334555220004e-51),
  qd_real( 3.0492922973540243e-01, -2.2989033898191262e-17,
       -1.4598901099661133e-34, 3.7760487693121827e-51),
  qd_real( 3.0784964004153487e-01, 2.7074088527245185e-17,
       1.2568858206899284e-33, 7.2931815105901645e-50),
  qd_real( 3.1076715274961147e-01, 2.0887076364048513e-17,
       -3.0130590791065942e-34, 1.3876739009935179e-51),
  qd_real( 3.1368174039889146e-01, 1.4560447299968912e-17,
       3.6564186898011595e-34, 1.1654264734999375e-50),
  qd_real( 3.1659337555616585e-01, 2.1435292512726283e-17,
       1.2338169231377316e-33, 3.3963542100989293e-50),
  qd_real( 3.1950203081601569e-01, -1.3981562491096626e-17,
       8.1730000697411350e-34, -7.7671096270210952e-50),
  qd_real( 3.2240767880106985e-01, -4.0519039937959398e-18,
       3.7438302780296796e-34, 8.7936731046639195e-51),
  qd_real( 3.2531029216226293e-01, 7.9171249463765892e-18,
       -6.7576622068146391e-35, 2.3021655066929538e-51),
  qd_real( 3.2820984357909255e-01, -2.6693140719641896e-17,
       7.8928851447534788e-34, 2.5525163821987809e-51),
  qd_real( 3.3110630575987643e-01, -2.7469465474778694e-17,
       -1.3401245916610206e-33, 6.5531762489976163e-50),
  qd_real( 3.3399965144200938e-01, 2.2598986806288142e-17,
       7.8063057192586115e-34, 2.0427600895486683e-50),
  qd_real( 3.3688985339222005e-01, -4.2000940033475092e-19,
       -2.9178652969985438e-36, -1.1597376437036749e-52),
  qd_real( 3.3977688440682685e-01, 6.6028679499418282e-18,
       1.2575009988669683e-34, 2.5569067699008304e-51),
  qd_real( 3.4266071731199438e-01, 1.9261518449306319e-17,
       -9.2754189135990867e-34, 8.5439996687390166e-50),
  qd_real( 3.4554132496398904e-01, 2.7251143672916123e-17,
       7.0138163601941737e-34, -1.4176292197454015e-50),
  qd_real( 3.4841868024943456e-01, 3.6974420514204918e-18,
       3.5532146878499996e-34, 1.9565462544501322e-50),
  qd_real( 3.5129275608556715e-01, -2.2670712098795844e-17,
       -1.6994216673139631e-34, -1.2271556077284517e-50),
  qd_real( 3.5416352542049040e-01, -1.6951763305764860e-17,
       1.2772331777814617e-33, -3.3703785435843310e-50),
  qd_real( 3.5703096123343003e-01, -4.8218191137919166e-19,
       -4.1672436994492361e-35, -7.1531167149364352e-52),
  qd_real( 3.5989503653498817e-01, -1.7601687123839282e-17,
       1.3375125473046791e-33, 7.9467815593584340e-50),
  qd_real( 3.6275572436739723e-01, -9.1668352663749849e-18,
       -7.4317843956936735e-34, -2.0199582511804564e-50),
  qd_real( 3.6561299780477385e-01, 1.6217898770457546e-17,
       1.1286970151961055e-33, -7.1825287318139010e-50),
  qd_real( 3.6846682995337232e-01, 1.0463640796159268e-17,
       2.0554984738517304e-35, 1.0441861305618769e-51),
  qd_real( 3.7131719395183754e-01, 3.4749239648238266e-19,
       -7.5151053042866671e-37, -2.8153468438650851e-53),
  qd_real( 3.7416406297145799e-01, 8.0114103761962118e-18,
       5.3429599813406052e-34, 1.0351378796539210e-50),
  qd_real( 3.7700741021641826e-01, -2.7255302041956930e-18,
       6.3646586445018137e-35, 8.3048657176503559e-52),
  qd_real( 3.7984720892405116e-01, 9.9151305855172370e-18,
       4.8761409697224886e-34, 1.4025084000776705e-50),
  qd_real( 3.8268343236508978e-01, -1.0050772696461588e-17,
       -2.0605316302806695e-34, -1.2717724698085205e-50),
  qd_real( 3.8551605384391885e-01, 1.5177665396472313e-17,
       1.4198230518016535e-33, 5.8955167159904235e-50),
  qd_real( 3.8834504669882630e-01, -1.0053770598398717e-17,
       7.5942999255057131e-34, -3.1967974046654219e-50),
  qd_real( 3.9117038430225387e-01, 1.7997787858243995e-17,
       -1.0613482402609856e-33, -5.4582148817791032e-50),
  qd_real( 3.9399204006104810e-01, 9.7649241641239336e-18,
       -2.1233599441284617e-34, -5.5529836795340819e-51),
  qd_real( 3.9680998741671031e-01, 2.0545063670840126e-17,
       6.1347058801922842e-34, 1.0733788150636430e-50),
  qd_real( 3.9962419984564684e-01, -1.5065497476189372e-17,
       -9.9653258881867298e-34, -5.7524323712725355e-50),
  qd_real( 4.0243465085941843e-01, 1.0902619339328270e-17,
       7.3998528125989765e-34, 2.2745784806823499e-50),
  qd_real( 4.0524131400498986e-01, 9.9111401942899884e-18,
       -2.5169070895434648e-34, 9.2772984818436573e-53),
  qd_real( 4.0804416286497869e-01, -7.0006015137351311e-18,
       -1.4108207334268228e-34, 1.5175546997577136e-52),
  qd_real( 4.1084317105790397e-01, -2.4219835190355499e-17,
       -1.1418902925313314e-33, -2.0996843165093468e-50),
  qd_real( 4.1363831223843456e-01, -1.0393984940597871e-17,
       -1.1481681174503880e-34, -2.0281052851028680e-51),
  qd_real( 4.1642956009763721e-01, -2.5475580413131732e-17,
       -3.4482678506112824e-34, 7.1788619351865480e-51),
  qd_real( 4.1921688836322396e-01, -4.2232463750110590e-18,
       -3.6053023045255790e-34, -2.2209673210025631e-50),
  qd_real( 4.2200027079979968e-01, 4.3543266994128527e-18,
       3.1734310272251190e-34, -1.3573247980738668e-50),
  qd_real( 4.2477968120910881e-01, 2.7462312204277281e-17,
       -4.6552847802111948e-34, 6.5961781099193122e-51),
  qd_real( 4.2755509343028208e-01, 9.4111898162954726e-18,
       -1.7446682426598801e-34, -2.2054492626480169e-51),
  qd_real( 4.3032648134008261e-01, 2.2259686974092690e-17,
       8.5972591314085075e-34, -2.9420897889003020e-50),
  qd_real( 4.3309381885315196e-01, 1.1224283329847517e-17,
       5.3223748041075651e-35, 5.3926192627014212e-51),
  qd_real( 4.3585707992225547e-01, 1.6230515450644527e-17,
       -6.4371449063579431e-35, -6.9102436481386757e-51),
  qd_real( 4.3861623853852766e-01, -2.0883315831075090e-17,
       -1.4259583540891877e-34, 6.3864763590657077e-52),
  qd_real( 4.4137126873171667e-01, 2.2360783886964969e-17,
       1.1864769603515770e-34, -3.8087003266189232e-51),
  qd_real( 4.4412214457042926e-01, -2.4218874422178315e-17,
       2.2205230838703907e-34, 9.2133035911356258e-51),
  qd_real( 4.4686884016237421e-01, -1.9222136142309382e-17,
       -4.4425678589732049e-35, -1.3673609292149535e-51),
  qd_real( 4.4961132965460660e-01, 4.8831924232035243e-18,
       2.7151084498191381e-34, -1.5653993171613154e-50),
  qd_real( 4.5234958723377089e-01, -1.4827977472196122e-17,
       -7.6947501088972324e-34, 1.7656856882031319e-50),
  qd_real( 4.5508358712634384e-01, -1.2379906758116472e-17,
       5.5289688955542643e-34, -8.5382312840209386e-51),
  qd_real( 4.5781330359887723e-01, -8.4554254922295949e-18,
       -6.3770394246764263e-34, 3.1778253575564249e-50),
  qd_real( 4.6053871095824001e-01, 1.8488777492177872e-17,
       -1.0527732154209725e-33, 3.3235593490947102e-50),
  qd_real( 4.6325978355186020e-01, -7.3514924533231707e-18,
       6.7175396881707035e-34, 3.9594127612123379e-50),
  qd_real( 4.6597649576796618e-01, -3.3023547778235135e-18,
       3.4904677050476886e-35, 3.4483855263874246e-51),
  qd_real( 4.6868882203582796e-01, -2.2949251681845054e-17,
       -1.1364757641823658e-33, 6.8840522501918612e-50),
  qd_real( 4.7139673682599764e-01, 6.5166781360690130e-18,
       2.9457546966235984e-34, -6.2159717738836630e-51),
  qd_real( 4.7410021465055002e-01, -8.1451601548978075e-18,
       -3.4789448555614422e-34, -1.1681943974658508e-50),
  qd_real( 4.7679923006332214e-01, -1.0293515338305794e-17,
       -3.6582045008369952e-34, 1.7424131479176475e-50),
  qd_real( 4.7949375766015301e-01, 1.8419999662684771e-17,
       -1.3040838621273312e-33, 1.0977131822246471e-50),
  qd_real( 4.8218377207912277e-01, -2.5861500925520442e-17,
       -6.2913197606500007e-36, 4.0802359808684726e-52),
  qd_real( 4.8486924800079112e-01, -1.8034004203262245e-17,
       -3.5244276906958044e-34, -1.7138318654749246e-50),
  qd_real( 4.8755016014843594e-01, 1.4231090931273653e-17,
       -1.8277733073262697e-34, -1.5208291790429557e-51),
  qd_real( 4.9022648328829116e-01, -5.1496145643440404e-18,
       -3.6903027405284104e-34, 1.5172940095151304e-50),
  qd_real( 4.9289819222978404e-01, -1.0257831676562186e-18,
       6.9520817760885069e-35, -2.4260961214090389e-51),
  qd_real( 4.9556526182577254e-01, -9.4323241942365362e-18,
       3.1212918657699143e-35, 4.2009072375242736e-52),
  qd_real( 4.9822766697278187e-01, -1.6126383830540798e-17,
       -1.5092897319298871e-33, 1.1049298890895917e-50),
  qd_real( 5.0088538261124083e-01, -3.9604015147074639e-17,
       -2.2208395201898007e-33, 1.3648202735839417e-49),
  qd_real( 5.0353838372571758e-01, -1.6731308204967497e-17,
       -1.0140233644074786e-33, 4.0953071937671477e-50),
  qd_real( 5.0618664534515534e-01, -4.8321592986493711e-17,
       9.2858107226642252e-34, 4.2699802401037005e-50),
  qd_real( 5.0883014254310699e-01, 4.7836968268014130e-17,
       -1.0727022928806035e-33, 2.7309374513672757e-50),
  qd_real( 5.1146885043797041e-01, -1.3088001221007579e-17,
       4.0929033363366899e-34, -3.7952190153477926e-50),
  qd_real( 5.1410274419322177e-01, -4.5712707523615624e-17,
       1.5488279442238283e-33, -2.5853959305521130e-50),
  qd_real( 5.1673179901764987e-01, 8.3018617233836515e-18,
       5.8251027467695202e-34, -2.2812397190535076e-50),
  qd_real( 5.1935599016558964e-01, -5.5331248144171145e-17,
       -3.1628375609769026e-35, -2.4091972051188571e-51),
  qd_real( 5.2197529293715439e-01, -4.6555795692088883e-17,
       4.6378980936850430e-34, -3.3470542934689532e-51),
  qd_real( 5.2458968267846895e-01, -4.3068869040082345e-17,
       -4.2013155291932055e-34, -1.5096069926700274e-50),
  qd_real( 5.2719913478190139e-01, -4.2202983480560619e-17,
       8.5585916184867295e-34, 7.9974339336732307e-50),
  qd_real( 5.2980362468629472e-01, -4.8067841706482342e-17,
       5.8309721046630296e-34, -8.9740761521756660e-51),
  qd_real( 5.3240312787719801e-01, -4.1020306135800895e-17,
       -1.9239996374230821e-33, -1.5326987913812184e-49),
  qd_real( 5.3499761988709726e-01, -5.3683132708358134e-17,
       -1.3900569918838112e-33, 2.7154084726474092e-50),
  qd_real( 5.3758707629564551e-01, -2.2617365388403054e-17,
       -5.9787279033447075e-34, 3.1204419729043625e-51),
  qd_real( 5.4017147272989285e-01, 2.7072447965935839e-17,
       1.1698799709213829e-33, -5.9094668515881500e-50),
  qd_real( 5.4275078486451589e-01, 1.7148261004757101e-17,
       -1.3525905925200870e-33, 4.9724411290727323e-50),
  qd_real( 5.4532498842204646e-01, -4.1517817538384258e-17,
       -1.5318930219385941e-33, 6.3629921101413974e-50),
  qd_real( 5.4789405917310019e-01, -2.4065878297113363e-17,
       -3.5639213669362606e-36, -2.6013270854271645e-52),
  qd_real( 5.5045797293660481e-01, -8.3319903015807663e-18,
       -2.3058454035767633e-34, -2.1611290432369010e-50),
  qd_real( 5.5301670558002758e-01, -4.7061536623798204e-17,
       -1.0617111545918056e-33, -1.6196316144407379e-50),
  qd_real( 5.5557023301960218e-01, 4.7094109405616768e-17,
       -2.0640520383682921e-33, 1.2290163188567138e-49),
  qd_real( 5.5811853122055610e-01, 1.3481176324765226e-17,
       -5.5016743873011438e-34, -2.3484822739335416e-50),
  qd_real( 5.6066157619733603e-01, -7.3956418153476152e-18,
       3.9680620611731193e-34, 3.1995952200836223e-50),
  qd_real( 5.6319934401383409e-01, 2.3835775146854829e-17,
       1.3511793173769814e-34, 9.3201311581248143e-51),
  qd_real( 5.6573181078361323e-01, -3.4096079596590466e-17,
       -1.7073289744303546e-33, 8.9147089975404507e-50),
  qd_real( 5.6825895267013160e-01, -5.0935673642769248e-17,
       -1.6274356351028249e-33, 9.8183151561702966e-51),
  qd_real( 5.7078074588696726e-01, 2.4568151455566208e-17,
       -1.2844481247560350e-33, -1.8037634376936261e-50),
  qd_real( 5.7329716669804220e-01, 8.5176611669306400e-18,
       -6.4443208788026766e-34, 2.2546105543273003e-50),
  qd_real( 5.7580819141784534e-01, -3.7909495458942734e-17,
       -2.7433738046854309e-33, 1.1130841524216795e-49),
  qd_real( 5.7831379641165559e-01, -2.6237691512372831e-17,
       1.3679051680738167e-33, -3.1409808935335900e-50),
  qd_real( 5.8081395809576453e-01, 1.8585338586613408e-17,
       2.7673843114549181e-34, 1.9605349619836937e-50),
  qd_real( 5.8330865293769829e-01, 3.4516601079044858e-18,
       1.8065977478946306e-34, -6.3953958038544646e-51),
  qd_real( 5.8579785745643886e-01, -3.7485501964311294e-18,
       2.7965403775536614e-34, -7.1816936024157202e-51),
  qd_real( 5.8828154822264533e-01, -2.9292166725006846e-17,
       -2.3744954603693934e-33, -1.1571631191512480e-50),
  qd_real( 5.9075970185887428e-01, -4.7013584170659542e-17,
       2.4808417611768356e-33, 1.2598907673643198e-50),
  qd_real( 5.9323229503979980e-01, 1.2892320944189053e-17,
       5.3058364776359583e-34, 4.1141674699390052e-50),
  qd_real( 5.9569930449243336e-01, -1.3438641936579467e-17,
       -6.7877687907721049e-35, -5.6046937531684890e-51),
  qd_real( 5.9816070699634227e-01, 3.8801885783000657e-17,
       -1.2084165858094663e-33, -4.0456610843430061e-50),
  qd_real( 6.0061647938386897e-01, -4.6398198229461932e-17,
       -1.6673493003710801e-33, 5.1982824378491445e-50),
  qd_real( 6.0306659854034816e-01, 3.7323357680559650e-17,
       2.7771920866974305e-33, -1.6194229649742458e-49),
  qd_real( 6.0551104140432555e-01, -3.1202672493305677e-17,
       1.2761267338680916e-33, -4.0859368598379647e-50),
  qd_real( 6.0794978496777363e-01, 3.5160832362096660e-17,
       -2.5546242776778394e-34, -1.4085313551220694e-50),
  qd_real( 6.1038280627630948e-01, -2.2563265648229169e-17,
       1.3185575011226730e-33, 8.2316691420063460e-50),
  qd_real( 6.1281008242940971e-01, -4.2693476568409685e-18,
       2.5839965886650320e-34, 1.6884412005622537e-50),
  qd_real( 6.1523159058062682e-01, 2.6231417767266950e-17,
       -1.4095366621106716e-33, 7.2058690491304558e-50),
  qd_real( 6.1764730793780398e-01, -4.7478594510902452e-17,
       -7.2986558263123996e-34, -3.0152327517439154e-50),
  qd_real( 6.2005721176328921e-01, -2.7983410837681118e-17,
       1.1649951056138923e-33, -5.4539089117135207e-50),
  qd_real( 6.2246127937414997e-01, 5.2940728606573002e-18,
       -4.8486411215945827e-35, 1.2696527641980109e-52),
  qd_real( 6.2485948814238634e-01, 3.3671846037243900e-17,
       -2.7846053391012096e-33, 5.6102718120012104e-50),
  qd_real( 6.2725181549514408e-01, 3.0763585181253225e-17,
       2.7068930273498138e-34, -1.1172240309286484e-50),
  qd_real( 6.2963823891492698e-01, 4.1115334049626806e-17,
       -1.9167473580230747e-33, 1.1118424028161730e-49),
  qd_real( 6.3201873593980906e-01, -4.0164942296463612e-17,
       -7.2208643641736723e-34, 3.7828920470544344e-50),
  qd_real( 6.3439328416364549e-01, 1.0420901929280035e-17,
       4.1174558929280492e-34, -1.4464152986630705e-51),
  qd_real( 6.3676186123628420e-01, 3.1419048711901611e-17,
       -2.2693738415126449e-33, -1.6023584204297388e-49),
  qd_real( 6.3912444486377573e-01, 1.2416796312271043e-17,
       -6.2095419626356605e-34, 2.7762065999506603e-50),
  qd_real( 6.4148101280858316e-01, -9.9883430115943310e-18,
       4.1969230376730128e-34, 5.6980543799257597e-51),
  qd_real( 6.4383154288979150e-01, -3.2084798795046886e-17,
       -1.2595311907053305e-33, -4.0205885230841536e-50),
  qd_real( 6.4617601298331639e-01, -2.9756137382280815e-17,
       -1.0275370077518259e-33, 8.0852478665893014e-51),
  qd_real( 6.4851440102211244e-01, 3.9870270313386831e-18,
       1.9408388509540788e-34, -5.1798420636193190e-51),
  qd_real( 6.5084668499638088e-01, 3.9714670710500257e-17,
       2.9178546787002963e-34, 3.8140635508293278e-51),
  qd_real( 6.5317284295377676e-01, 8.5695642060026238e-18,
       -6.9165322305070633e-34, 2.3873751224185395e-50),
  qd_real( 6.5549285299961535e-01, 3.5638734426385005e-17,
       1.2695365790889811e-33, 4.3984952865412050e-50),
  qd_real( 6.5780669329707864e-01, 1.9580943058468545e-17,
       -1.1944272256627192e-33, 2.8556402616436858e-50),
  qd_real( 6.6011434206742048e-01, -1.3960054386823638e-19,
       6.1515777931494047e-36, 5.3510498875622660e-52),
  qd_real( 6.6241577759017178e-01, -2.2615508885764591e-17,
       5.0177050318126862e-34, 2.9162532399530762e-50),
  qd_real( 6.6471097820334490e-01, -3.6227793598034367e-17,
       -9.0607934765540427e-34, 3.0917036342380213e-50),
  qd_real( 6.6699992230363747e-01, 3.5284364997428166e-17,
       -1.0382057232458238e-33, 7.3812756550167626e-50),
  qd_real( 6.6928258834663612e-01, -5.4592652417447913e-17,
       -2.5181014709695152e-33, -1.6867875999437174e-49),
  qd_real( 6.7155895484701844e-01, -4.0489037749296692e-17,
       3.1995835625355681e-34, -1.4044414655670960e-50),
  qd_real( 6.7382900037875604e-01, 2.3091901236161086e-17,
       5.7428037192881319e-34, 1.1240668354625977e-50),
  qd_real( 6.7609270357531592e-01, 3.7256902248049466e-17,
       1.7059417895764375e-33, 9.7326347795300652e-50),
  qd_real( 6.7835004312986147e-01, 1.8302093041863122e-17,
       9.5241675746813072e-34, 5.0328101116133503e-50),
  qd_real( 6.8060099779545302e-01, 2.8473293354522047e-17,
       4.1331805977270903e-34, 4.2579030510748576e-50),
  qd_real( 6.8284554638524808e-01, -1.2958058061524531e-17,
       1.8292386959330698e-34, 3.4536209116044487e-51),
  qd_real( 6.8508366777270036e-01, 2.5948135194645137e-17,
       -8.5030743129500702e-34, -6.9572086141009930e-50),
  qd_real( 6.8731534089175916e-01, -5.5156158714917168e-17,
       1.1896489854266829e-33, -7.8505896218220662e-51),
  qd_real( 6.8954054473706694e-01, -1.5889323294806790e-17,
       9.1242356240205712e-34, 3.8315454152267638e-50),
  qd_real( 6.9175925836415775e-01, 2.7406078472410668e-17,
       1.3286508943202092e-33, 1.0651869129580079e-51),
  qd_real( 6.9397146088965400e-01, 7.4345076956280137e-18,
       7.5061528388197460e-34, -1.5928000240686583e-50),
  qd_real( 6.9617713149146299e-01, -4.1224081213582889e-17,
       -3.1838716762083291e-35, -3.9625587412119131e-51),
  qd_real( 6.9837624940897280e-01, 4.8988282435667768e-17,
       1.9134010413244152e-33, 2.6161153243793989e-50),
  qd_real( 7.0056879394324834e-01, 3.1027960192992922e-17,
       9.5638250509179997e-34, 4.5896916138107048e-51),
  qd_real( 7.0275474445722530e-01, 2.5278294383629822e-18,
       -8.6985561210674942e-35, -5.6899862307812990e-51),
  qd_real( 7.0493408037590488e-01, 2.7608725585748502e-17,
       2.9816599471629137e-34, 1.1533044185111206e-50),
  qd_real( 7.0710678118654757e-01, -4.8336466567264567e-17,
       2.0693376543497068e-33, 2.4677734957341755e-50)
};

static const qd_real cos_table [] = {
  qd_real( 9.9999529380957619e-01, -1.9668064285322189e-17,
       -6.3053955095883481e-34, 5.3266110855726731e-52),
  qd_real( 9.9998117528260111e-01, 3.3568103522895585e-17,
       -1.4740132559368063e-35, 9.8603097594755596e-52),
  qd_real( 9.9995764455196390e-01, -3.1527836866647287e-17,
       2.6363251186638437e-33, 1.0007504815488399e-49),
  qd_real( 9.9992470183914450e-01, 3.7931082512668012e-17,
       -8.5099918660501484e-35, -4.9956973223295153e-51),
  qd_real( 9.9988234745421256e-01, -3.5477814872408538e-17,
       1.7102001035303974e-33, -1.0725388519026542e-49),
  qd_real( 9.9983058179582340e-01, 1.8825140517551119e-17,
       -5.1383513457616937e-34, -3.8378827995403787e-50),
  qd_real( 9.9976940535121528e-01, 4.2681177032289012e-17,
       1.9062302359737099e-33, -6.0221153262881160e-50),
  qd_real( 9.9969881869620425e-01, -2.9851486403799753e-17,
       -1.9084787370733737e-33, 5.5980260344029202e-51),
  qd_real( 9.9961882249517864e-01, -4.1181965521424734e-17,
       2.0915365593699916e-33, 8.1403390920903734e-50),
  qd_real( 9.9952941750109314e-01, 2.0517917823755591e-17,
       -4.7673802585706520e-34, -2.9443604198656772e-50),
  qd_real( 9.9943060455546173e-01, 3.9644497752257798e-17,
       -2.3757223716722428e-34, -1.2856759011361726e-51),
  qd_real( 9.9932238458834954e-01, -4.2858538440845682e-17,
       3.3235101605146565e-34, -8.3554272377057543e-51),
  qd_real( 9.9920475861836389e-01, 9.1796317110385693e-18,
       5.5416208185868570e-34, 8.0267046717615311e-52),
  qd_real( 9.9907772775264536e-01, 2.1419007653587032e-17,
       -7.9048203318529618e-34, -5.3166296181112712e-50),
  qd_real( 9.9894129318685687e-01, -2.0610641910058638e-17,
       -1.2546525485913485e-33, -7.5175888806157064e-50),
  qd_real( 9.9879545620517241e-01, -1.2291693337075465e-17,
       2.4468446786491271e-34, 1.0723891085210268e-50),
  qd_real( 9.9864021818026527e-01, -4.8690254312923302e-17,
       -2.9470881967909147e-34, -1.3000650761346907e-50),
  qd_real( 9.9847558057329477e-01, -2.2002931182778795e-17,
       -1.2371509454944992e-33, -2.4911225131232065e-50),
  qd_real( 9.9830154493389289e-01, -5.1869402702792278e-17,
       1.0480195493633452e-33, -2.8995649143155511e-50),
  qd_real( 9.9811811290014918e-01, 2.7935487558113833e-17,
       2.4423341255830345e-33, -6.7646699175334417e-50),
  qd_real( 9.9792528619859600e-01, 1.7143659778886362e-17,
       5.7885840902887460e-34, -9.2601432603894597e-51),
  qd_real( 9.9772306664419164e-01, -2.6394475274898721e-17,
       -1.6176223087661783e-34, -9.9924942889362281e-51),
  qd_real( 9.9751145614030345e-01, 5.6007205919806937e-18,
       -5.9477673514685690e-35, -1.4166807162743627e-54),
  qd_real( 9.9729045667869021e-01, 9.1647695371101735e-18,
       6.7824134309739296e-34, -8.6191392795543357e-52),
  qd_real( 9.9706007033948296e-01, 1.6734093546241963e-17,
       -1.3169951440780028e-33, 1.0311048767952477e-50),
  qd_real( 9.9682029929116567e-01, 4.7062820708615655e-17,
       2.8412041076474937e-33, -8.0006155670263622e-50),
  qd_real( 9.9657114579055484e-01, 1.1707179088390986e-17,
       -7.5934413263024663e-34, 2.8474848436926008e-50),
  qd_real( 9.9631261218277800e-01, 1.1336497891624735e-17,
       3.4002458674414360e-34, 7.7419075921544901e-52),
  qd_real( 9.9604470090125197e-01, 2.2870031707670695e-17,
       -3.9184839405013148e-34, -3.7081260416246375e-50),
  qd_real( 9.9576741446765982e-01, -2.3151908323094359e-17,
       -1.6306512931944591e-34, -1.5925420783863192e-51),
  qd_real( 9.9548075549192694e-01, 3.2084621412226554e-18,
       -4.9501292146013023e-36, -2.7811428850878516e-52),
  qd_real( 9.9518472667219693e-01, -4.2486913678304410e-17,
       1.3315510772504614e-33, 6.7927987417051888e-50),
  qd_real( 9.9487933079480562e-01, 4.2130813284943662e-18,
       -4.2062597488288452e-35, 2.5157064556087620e-51),
  qd_real( 9.9456457073425542e-01, 3.6745069641528058e-17,
       -3.0603304105471010e-33, 1.0397872280487526e-49),
  qd_real( 9.9424044945318790e-01, 4.4129423472462673e-17,
       -3.0107231708238066e-33, 7.4201582906861892e-50),
  qd_real( 9.9390697000235606e-01, -1.8964849471123746e-17,
       -1.5980853777937752e-35, -8.5374807150597082e-52),
  qd_real( 9.9356413552059530e-01, 2.9752309927797428e-17,
       -4.5066707331134233e-34, -3.3548191633805036e-50),
  qd_real( 9.9321194923479450e-01, 3.3096906261272262e-17,
       1.5592811973249567e-33, 1.4373977733253592e-50),
  qd_real( 9.9285041445986510e-01, -1.4094517733693302e-17,
       -1.1954558131616916e-33, 1.8761873742867983e-50),
  qd_real( 9.9247953459870997e-01, 3.1093055095428906e-17,
       -1.8379594757818019e-33, -3.9885758559381314e-51),
  qd_real( 9.9209931314219180e-01, -3.9431926149588778e-17,
       -6.2758062911047230e-34, -1.2960929559212390e-50),
  qd_real( 9.9170975366909953e-01, -2.3372891311883661e-18,
       2.7073298824968591e-35, -1.2569459441802872e-51),
  qd_real( 9.9131085984611544e-01, -2.5192111583372105e-17,
       -1.2852471567380887e-33, 5.2385212584310483e-50),
  qd_real( 9.9090263542778001e-01, 1.5394565094566704e-17,
       -1.0799984133184567e-33, 2.7451115960133595e-51),
  qd_real( 9.9048508425645709e-01, -5.5411437553780867e-17,
       -1.4614017210753585e-33, -3.8339374397387620e-50),
  qd_real( 9.9005821026229712e-01, -1.7055485906233963e-17,
       1.3454939685758777e-33, 7.3117589137300036e-50),
  qd_real( 9.8962201746320089e-01, -5.2398217968132530e-17,
       1.3463189211456219e-33, 5.8021640554894872e-50),
  qd_real( 9.8917650996478101e-01, -4.0987309937047111e-17,
       -4.4857560552048437e-34, -3.9414504502871125e-50),
  qd_real( 9.8872169196032378e-01, -1.0976227206656125e-17,
       3.2311342577653764e-34, 9.6367946583575041e-51),
  qd_real( 9.8825756773074946e-01, 2.7030607784372632e-17,
       7.7514866488601377e-35, 2.1019644956864938e-51),
  qd_real( 9.8778414164457218e-01, -2.3600693397159021e-17,
       -1.2323283769707861e-33, 3.0130900716803339e-50),
  qd_real( 9.8730141815785843e-01, -5.2332261255715652e-17,
       -2.7937644333152473e-33, 1.2074160567958408e-49),
  qd_real( 9.8680940181418553e-01, -5.0287214351061075e-17,
       -2.2681526238144461e-33, 4.4003694320169133e-50),
  qd_real( 9.8630809724459867e-01, -2.1520877103013341e-17,
       1.1866528054187716e-33, -7.8532199199813836e-50),
  qd_real( 9.8579750916756748e-01, -5.1439452979953012e-17,
       2.6276439309996725e-33, 7.5423552783286347e-50),
  qd_real( 9.8527764238894122e-01, 2.3155637027900207e-17,
       -7.5275971545764833e-34, 1.0582231660456094e-50),
  qd_real( 9.8474850180190421e-01, 1.0548144061829957e-17,
       2.8786145266267306e-34, -3.6782210081466112e-51),
  qd_real( 9.8421009238692903e-01, 4.7983922627050691e-17,
       2.2597419645070588e-34, 1.7573875814863400e-50),
  qd_real( 9.8366241921173025e-01, 1.9864948201635255e-17,
       -1.0743046281211033e-35, 1.7975662796558100e-52),
  qd_real( 9.8310548743121629e-01, 4.2170007522888628e-17,
       8.2396265656440904e-34, -8.0803700139096561e-50),
  qd_real( 9.8253930228744124e-01, 1.5149580813777224e-17,
       -4.1802771422186237e-34, -2.2150174326226160e-50),
  qd_real( 9.8196386910955524e-01, 2.1108443711513084e-17,
       -1.5253013442896054e-33, -6.8388082079337969e-50),
  qd_real( 9.8137919331375456e-01, 1.3428163260355633e-17,
       -6.5294290469962986e-34, 2.7965412287456268e-51),
  qd_real( 9.8078528040323043e-01, 1.8546939997825006e-17,
       -1.0696564445530757e-33, 6.6668174475264961e-50),
  qd_real( 9.8018213596811743e-01, -3.6801786963856159e-17,
       6.3245171387992842e-34, 1.8600176137175971e-50),
  qd_real( 9.7956976568544052e-01, 1.5573991584990420e-17,
       -1.3401066029782990e-33, -1.7263702199862149e-50),
  qd_real( 9.7894817531906220e-01, -2.3817727961148053e-18,
       -1.0694750370381661e-34, -8.2293047196087462e-51),
  qd_real( 9.7831737071962765e-01, -2.1623082233344895e-17,
       1.0970403012028032e-33, 7.7091923099369339e-50),
  qd_real( 9.7767735782450993e-01, 5.0514136167059628e-17,
       -1.3254751701428788e-33, 7.0161254312124538e-50),
  qd_real( 9.7702814265775439e-01, -4.3353875751555997e-17,
       5.4948839831535478e-34, -9.2755263105377306e-51),
  qd_real( 9.7636973133002114e-01, 9.3093931526213780e-18,
       -4.1184949155685665e-34, -3.1913926031393690e-50),
  qd_real( 9.7570213003852857e-01, -2.5572556081259686e-17,
       -9.3174244508942223e-34, -8.3675863211646863e-51),
  qd_real( 9.7502534506699412e-01, 2.6642660651899135e-17,
       1.7819392739353853e-34, -3.3159625385648947e-51),
  qd_real( 9.7433938278557586e-01, 2.3041221476151512e-18,
       1.0758686005031430e-34, 5.1074116432809478e-51),
  qd_real( 9.7364424965081198e-01, -5.1729808691005871e-17,
       -1.5508473005989887e-33, -1.6505125917675401e-49),
  qd_real( 9.7293995220556018e-01, -3.1311211122281800e-17,
       -2.6874087789006141e-33, -2.1652434818822145e-51),
  qd_real( 9.7222649707893627e-01, 3.6461169785938221e-17,
       3.0309636883883133e-33, -1.2702716907967306e-51),
  qd_real( 9.7150389098625178e-01, -7.9865421122289046e-18,
       -4.3628417211263380e-34, 3.4307517798759352e-51),
  qd_real( 9.7077214072895035e-01, -4.7992163325114922e-17,
       3.0347528910975783e-33, 8.5989199506479701e-50),
  qd_real( 9.7003125319454397e-01, 1.8365300348428844e-17,
       -1.4311097571944918e-33, 8.5846781998740697e-51),
  qd_real( 9.6928123535654853e-01, -4.5663660261927896e-17,
       9.6147526917239387e-34, 8.1267605207871330e-51),
  qd_real( 9.6852209427441727e-01, 4.9475074918244771e-17,
       2.8558738351911241e-33, 6.2948422316507461e-50),
  qd_real( 9.6775383709347551e-01, -4.5512132825515820e-17,
       -1.4127617988719093e-33, -8.4620609089704578e-50),
  qd_real( 9.6697647104485207e-01, 3.8496228837337864e-17,
       -5.3881631542745647e-34, -3.5221863171458959e-50),
  qd_real( 9.6619000344541250e-01, 5.1298840401665493e-17,
       1.4564075904769808e-34, 1.0095973971377432e-50),
  qd_real( 9.6539444169768940e-01, -2.3745389918392156e-17,
       5.9221515590053862e-34, -3.8811192556231094e-50),
  qd_real( 9.6458979328981276e-01, -3.4189470735959786e-17,
       2.2982074155463522e-33, -4.5128791045607634e-50),
  qd_real( 9.6377606579543984e-01, 2.6463950561220029e-17,
       -2.9073234590199323e-36, -1.2938328629395601e-52),
  qd_real( 9.6295326687368388e-01, 8.9341960404313634e-18,
       -3.9071244661020126e-34, 1.6212091116847394e-50),
  qd_real( 9.6212140426904158e-01, 1.5236770453846305e-17,
       -1.3050173525597142e-33, 7.9016122394092666e-50),
  qd_real( 9.6128048581132064e-01, 2.0933955216674039e-18,
       1.0768607469015692e-34, -5.9453639304361774e-51),
  qd_real( 9.6043051941556579e-01, 2.4653904815317185e-17,
       -1.3792169410906322e-33, -4.7726598378506903e-51),
  qd_real( 9.5957151308198452e-01, 1.1000640085000957e-17,
       -4.2036030828223975e-34, 4.0023704842606573e-51),
  qd_real( 9.5870347489587160e-01, -4.3685014392372053e-17,
       2.2001800662729131e-33, -1.0553721324358075e-49),
  qd_real( 9.5782641302753291e-01, -1.7696710075371263e-17,
       1.9164034110382190e-34, 8.1489235071754813e-51),
  qd_real( 9.5694033573220882e-01, 4.0553869861875701e-17,
       -1.7147013364302149e-33, 2.5736745295329455e-50),
  qd_real( 9.5604525134999641e-01, 3.7705045279589067e-17,
       1.9678699997347571e-33, 8.5093177731230180e-50),
  qd_real( 9.5514116830577067e-01, 5.0088652955014668e-17,
       -2.6983181838059211e-33, 1.0102323575596493e-49),
  qd_real( 9.5422809510910567e-01, -3.7545901690626874e-17,
       1.4951619241257764e-33, -8.2717333151394973e-50),
  qd_real( 9.5330604035419386e-01, -2.5190738779919934e-17,
       -1.4272239821134379e-33, -4.6717286809283155e-50),
  qd_real( 9.5237501271976588e-01, -2.0269300462299272e-17,
       -1.0635956887246246e-33, -3.5514537666487619e-50),
  qd_real( 9.5143502096900834e-01, 3.1350584123266695e-17,
       -2.4824833452737813e-33, 9.5450335525380613e-51),
  qd_real( 9.5048607394948170e-01, 1.9410097562630436e-17,
       -8.1559393949816789e-34, -1.0501209720164562e-50),
  qd_real( 9.4952818059303667e-01, -7.5544151928043298e-18,
       -5.1260245024046686e-34, 1.8093643389040406e-50),
  qd_real( 9.4856134991573027e-01, 2.0668262262333232e-17,
       -5.9440730243667306e-34, 1.4268853111554300e-50),
  qd_real( 9.4758559101774109e-01, 4.3417993852125991e-17,
       -2.7728667889840373e-34, 5.5709160196519968e-51),
  qd_real( 9.4660091308328353e-01, 3.5056800210680730e-17,
       9.8578536940318117e-34, 6.6035911064585197e-50),
  qd_real( 9.4560732538052128e-01, 4.6019102478523738e-17,
       -6.2534384769452059e-34, 1.5758941215779961e-50),
  qd_real( 9.4460483726148026e-01, 8.8100545476641165e-18,
       5.2291695602757842e-34, -3.3487256018407123e-50),
  qd_real( 9.4359345816196039e-01, -2.4093127844404214e-17,
       1.0283279856803939e-34, -2.3398232614531355e-51),
  qd_real( 9.4257319760144687e-01, 1.3235564806436886e-17,
       -5.7048262885386911e-35, 3.9947050442753744e-51),
  qd_real( 9.4154406518302081e-01, -2.7896379547698341e-17,
       1.6273236356733898e-33, -5.3075944708471203e-51),
  qd_real( 9.4050607059326830e-01, 2.8610421567116268e-17,
       2.9261501147538827e-33, -2.6849867690896925e-50),
  qd_real( 9.3945922360218992e-01, -7.0152867943098655e-18,
       -5.6395693818011210e-34, 3.5568142678987651e-50),
  qd_real( 9.3840353406310806e-01, 5.4242545044795490e-17,
       -1.9039966607859759e-33, -1.5627792988341215e-49),
  qd_real( 9.3733901191257496e-01, -3.6570926284362776e-17,
       -1.1902940071273247e-33, -1.1215082331583223e-50),
  qd_real( 9.3626566717027826e-01, -1.3013766145497654e-17,
       5.2229870061990595e-34, -3.3972777075634108e-51),
  qd_real( 9.3518350993894761e-01, -3.2609395302485065e-17,
       -8.1813015218875245e-34, 5.5642140024928139e-50),
  qd_real( 9.3409255040425887e-01, 4.4662824360767511e-17,
       -2.5903243047396916e-33, 8.1505209004343043e-50),
  qd_real( 9.3299279883473885e-01, 4.2041415555384355e-17,
       9.0285896495521276e-34, 5.3019984977661259e-50),
  qd_real( 9.3188426558166815e-01, -4.0785944377318095e-17,
       1.7631450298754169e-33, 2.5776403305507453e-50),
  qd_real( 9.3076696107898371e-01, 1.9703775102838329e-17,
       6.5657908718278205e-34, -1.9480347966259524e-51),
  qd_real( 9.2964089584318121e-01, 5.1282530016864107e-17,
       2.3719739891916261e-34, -1.7230065426917127e-50),
  qd_real( 9.2850608047321559e-01, -2.3306639848485943e-17,
       -7.7799084333208503e-34, -5.8597558009300305e-50),
  qd_real( 9.2736252565040111e-01, -2.7677111692155437e-17,
       2.2110293450199576e-34, 2.0349190819680613e-50),
  qd_real( 9.2621024213831138e-01, -3.7303754586099054e-17,
       2.0464457809993405e-33, 1.3831799631231817e-49),
  qd_real( 9.2504924078267758e-01, 6.0529447412576159e-18,
       -8.8256517760278541e-35, 1.8285462122388328e-51),
  qd_real( 9.2387953251128674e-01, 1.7645047084336677e-17,
       -5.0442537321586818e-34, -4.0478677716823890e-50),
  qd_real( 9.2270112833387852e-01, 5.2963798918539814e-17,
       -5.7135699628876685e-34, 3.0163671797219087e-50),
  qd_real( 9.2151403934204190e-01, 4.1639843390684644e-17,
       1.1891485604702356e-33, 2.0862437594380324e-50),
  qd_real( 9.2031827670911059e-01, -2.7806888779036837e-17,
       2.7011013677071274e-33, 1.1998578792455499e-49),
  qd_real( 9.1911385169005777e-01, -2.6496484622344718e-17,
       6.5403604763461920e-34, -2.8997180201186078e-50),
  qd_real( 9.1790077562139050e-01, -3.9074579680849515e-17,
       2.3004636541490264e-33, 3.9851762744443107e-50),
  qd_real( 9.1667905992104270e-01, -4.1733978698287568e-17,
       1.2094444804381172e-33, 4.9356916826097816e-50),
  qd_real( 9.1544871608826783e-01, -1.3591056692900894e-17,
       5.9923027475594735e-34, 2.1403295925962879e-50),
  qd_real( 9.1420975570353069e-01, -3.6316182527814423e-17,
       -1.9438819777122554e-33, 2.8340679287728316e-50),
  qd_real( 9.1296219042839821e-01, -4.7932505228039469e-17,
       -1.7753551889428638e-33, 4.0607782903868160e-51),
  qd_real( 9.1170603200542988e-01, -2.6913273175034130e-17,
       -5.1928101916162528e-35, 1.1338175936090630e-51),
  qd_real( 9.1044129225806725e-01, -5.0433041673313820e-17,
       1.0938746257404305e-33, 9.5378272084170731e-51),
  qd_real( 9.0916798309052238e-01, -3.6878564091359894e-18,
       2.9951330310507693e-34, -1.2225666136919926e-50),
  qd_real( 9.0788611648766626e-01, -4.9459964301225840e-17,
       -1.6599682707075313e-33, -5.1925202712634716e-50),
  qd_real( 9.0659570451491533e-01, 3.0506718955442023e-17,
       -1.4478836557141204e-33, 1.8906373784448725e-50),
  qd_real( 9.0529675931811882e-01, -4.1153099826889901e-17,
       2.9859368705184223e-33, 5.1145293917439211e-50),
  qd_real( 9.0398929312344334e-01, -6.6097544687484308e-18,
       1.2728013034680357e-34, -4.3026097234014823e-51),
  qd_real( 9.0267331823725883e-01, -1.9250787033961483e-17,
       1.3242128993244527e-33, -5.2971030688703665e-50),
  qd_real( 9.0134884704602203e-01, -1.3524789367698682e-17,
       6.3605353115880091e-34, 3.6227400654573828e-50),
  qd_real( 9.0001589201616028e-01, -5.0639618050802273e-17,
       1.0783525384031576e-33, 2.8130016326515111e-50),
  qd_real( 8.9867446569395382e-01, 2.6316906461033013e-17,
       3.7003137047796840e-35, -2.3447719900465938e-51),
  qd_real( 8.9732458070541832e-01, -3.6396283314867290e-17,
       -2.3611649895474815e-33, 1.1837247047900082e-49),
  qd_real( 8.9596624975618511e-01, 4.9025099114811813e-17,
       -1.9440489814795326e-33, -1.7070486667767033e-49),
  qd_real( 8.9459948563138270e-01, -1.7516226396814919e-17,
       -1.3200670047246923e-33, -1.5953009884324695e-50),
  qd_real( 8.9322430119551532e-01, -4.1161239151908913e-18,
       2.5380253805715999e-34, 4.2849455510516192e-51),
  qd_real( 8.9184070939234272e-01, 4.6690228137124547e-18,
       1.6150254286841982e-34, -3.9617448820725012e-51),
  qd_real( 8.9044872324475788e-01, 1.1781931459051803e-17,
       -1.3346142209571930e-34, -9.4982373530733431e-51),
  qd_real( 8.8904835585466457e-01, -1.1164514966766675e-17,
       -3.4797636107798736e-34, -1.5605079997040631e-50),
  qd_real( 8.8763962040285393e-01, 1.2805091918587960e-17,
       3.9948742059584459e-35, 3.8940716325338136e-51),
  qd_real( 8.8622253014888064e-01, -6.7307369600274315e-18,
       1.2385593432917413e-34, 2.0364014759133320e-51),
  qd_real( 8.8479709843093779e-01, -9.4331469628972690e-18,
       -5.7106541478701439e-34, 1.8260134111907397e-50),
  qd_real( 8.8336333866573158e-01, 1.5822643380255127e-17,
       -7.8921320007588250e-34, -1.4782321016179836e-50),
  qd_real( 8.8192126434835505e-01, -1.9843248405890562e-17,
       -7.0412114007673834e-34, -1.0636770169389104e-50),
  qd_real( 8.8047088905216075e-01, 1.6311096602996350e-17,
       -5.7541360594724172e-34, -4.0128611862170021e-50),
  qd_real( 8.7901222642863353e-01, -4.7356837291118011e-17,
       1.4388771297975192e-33, -2.9085554304479134e-50),
  qd_real( 8.7754529020726124e-01, 5.0113311846499550e-17,
       2.8382769008739543e-34, 1.5550640393164140e-50),
  qd_real( 8.7607009419540660e-01, 5.8729024235147677e-18,
       2.7941144391738458e-34, -1.8536073846509828e-50),
  qd_real( 8.7458665227817611e-01, -5.7216617730397065e-19,
       -2.9705811503689596e-35, 8.7389593969796752e-52),
  qd_real( 8.7309497841829009e-01, 7.8424672990129903e-18,
       -4.8685015839797165e-34, -2.2815570587477527e-50),
  qd_real( 8.7159508665595109e-01, -5.5272998038551050e-17,
       -2.2104090204984907e-33, -9.7749763187643172e-50),
  qd_real( 8.7008699110871146e-01, -4.1888510868549968e-17,
       7.0900185861878415e-34, 3.7600251115157260e-50),
  qd_real( 8.6857070597134090e-01, 2.7192781689782903e-19,
       -1.6710140396932428e-35, -1.2625514734637969e-51),
  qd_real( 8.6704624551569265e-01, 3.0267859550930567e-18,
       -1.1559438782171572e-34, -5.3580556397808012e-52),
  qd_real( 8.6551362409056909e-01, -6.3723113549628899e-18,
       2.3725520321746832e-34, 1.5911880348395175e-50),
  qd_real( 8.6397285612158670e-01, 4.1486355957361607e-17,
       2.2709976932210266e-33, -8.1228385659479984e-50),
  qd_real( 8.6242395611104050e-01, 3.7008992527383130e-17,
       5.2128411542701573e-34, 2.6945600081026861e-50),
  qd_real( 8.6086693863776731e-01, -3.0050048898573656e-17,
       -8.8706183090892111e-34, 1.5005320558097301e-50),
  qd_real( 8.5930181835700836e-01, 4.2435655816850687e-17,
       7.6181814059912025e-34, -3.9592127850658708e-50),
  qd_real( 8.5772861000027212e-01, -4.8183447936336620e-17,
       -1.1044130517687532e-33, -8.7400233444645562e-50),
  qd_real( 8.5614732837519447e-01, 9.1806925616606261e-18,
       5.6328649785951470e-34, 2.3326646113217378e-51),
  qd_real( 8.5455798836540053e-01, -1.2991124236396092e-17,
       1.2893407722948080e-34, -3.6506925747583053e-52),
  qd_real( 8.5296060493036363e-01, 2.7152984251981370e-17,
       7.4336483283120719e-34, 4.2162417622350668e-50),
  qd_real( 8.5135519310526520e-01, -5.3279874446016209e-17,
       2.2281156380919942e-33, -4.0281886404138477e-50),
  qd_real( 8.4974176800085244e-01, 5.1812347659974015e-17,
       3.0810626087331275e-33, -2.5931308201994965e-50),
  qd_real( 8.4812034480329723e-01, 1.8762563415239981e-17,
       1.4048773307919617e-33, -2.4915221509958691e-50),
  qd_real( 8.4649093877405213e-01, -4.7969419958569345e-17,
       -2.7518267097886703e-33, -7.3518959727313350e-50),
  qd_real( 8.4485356524970712e-01, -4.3631360296879637e-17,
       -2.0307726853367547e-33, 4.3097229819851761e-50),
  qd_real( 8.4320823964184544e-01, 9.6536707005959077e-19,
       2.8995142431556364e-36, 9.6715076811480284e-53),
  qd_real( 8.4155497743689844e-01, -3.4095465391321557e-17,
       -8.4130208607579595e-34, -4.9447283960568686e-50),
  qd_real( 8.3989379419599952e-01, -1.6673694881511411e-17,
       -1.4759184141750289e-33, -7.5795098161914058e-50),
  qd_real( 8.3822470555483808e-01, -3.5560085052855026e-17,
       1.1689791577022643e-33, -5.8627347359723411e-50),
  qd_real( 8.3654772722351201e-01, -2.0899059027066533e-17,
       -9.8104097821002585e-35, -3.1609177868229853e-51),
  qd_real( 8.3486287498638001e-01, 4.6048430609159657e-17,
       -5.1827423265239912e-34, -7.0505343435504109e-51),
  qd_real( 8.3317016470191319e-01, 1.3275129507229764e-18,
       4.8589164115370863e-35, 4.5422281300506859e-51),
  qd_real( 8.3146961230254524e-01, 1.4073856984728024e-18,
       4.6951315383980830e-35, 5.1431906049905658e-51),
  qd_real( 8.2976123379452305e-01, -2.9349109376485597e-18,
       1.1496917934149818e-34, 3.5186665544980233e-51),
  qd_real( 8.2804504525775580e-01, -4.4196593225871532e-17,
       2.7967864855211251e-33, 1.0030777287393502e-49),
  qd_real( 8.2632106284566353e-01, -5.3957485453612902e-17,
       6.8976896130138550e-34, 3.8106164274199196e-50),
  qd_real( 8.2458930278502529e-01, -2.6512360488868275e-17,
       1.6916964350914386e-34, 6.7693974813562649e-51),
  qd_real( 8.2284978137582632e-01, 1.5193019034505495e-17,
       9.6890547246521685e-34, 5.6994562923653264e-50),
  qd_real( 8.2110251499110465e-01, 3.0715131609697682e-17,
       -1.7037168325855879e-33, -1.1149862443283853e-49),
  qd_real( 8.1934752007679701e-01, -4.8200736995191133e-17,
       -1.5574489646672781e-35, -9.5647853614522216e-53),
  qd_real( 8.1758481315158371e-01, -1.4883149812426772e-17,
       -7.8273262771298917e-34, 4.1332149161031594e-50),
  qd_real( 8.1581441080673378e-01, 8.2652693782130871e-18,
       -2.3028778135179471e-34, 1.5102071387249843e-50),
  qd_real( 8.1403632970594841e-01, -5.2127351877042624e-17,
       -1.9047670611316360e-33, -1.6937269585941507e-49),
  qd_real( 8.1225058658520388e-01, 3.1054545609214803e-17,
       2.2649541922707251e-34, -7.4221684154649405e-51),
  qd_real( 8.1045719825259477e-01, 2.3520367349840499e-17,
       -7.7530070904846341e-34, -7.2792616357197140e-50),
  qd_real( 8.0865618158817498e-01, 9.3251597879721674e-18,
       -7.1823301933068394e-34, 2.3925440846132106e-50),
  qd_real( 8.0684755354379922e-01, 4.9220603766095546e-17,
       2.9796016899903487e-33, 1.5220754223615788e-49),
  qd_real( 8.0503133114296355e-01, 5.1368289568212149e-17,
       6.3082807402256524e-34, 7.3277646085129827e-51),
  qd_real( 8.0320753148064494e-01, -3.3060609804814910e-17,
       -1.2242726252420433e-33, 2.8413673268630117e-50),
  qd_real( 8.0137617172314024e-01, -2.0958013413495834e-17,
       -4.3798162198006931e-34, 2.0235690497752515e-50),
  qd_real( 7.9953726910790501e-01, 2.0356723822005431e-17,
       -9.7448513696896360e-34, 5.3608109599696008e-52),
  qd_real( 7.9769084094339116e-01, -4.6730759884788944e-17,
       2.3075897077191757e-33, 3.1605567774640253e-51),
  qd_real( 7.9583690460888357e-01, -3.0062724851910721e-17,
       -2.2496210832042235e-33, -6.5881774117183040e-50),
  qd_real( 7.9397547755433717e-01, -7.4194631759921416e-18,
       2.4124341304631069e-34, -4.9956808616244972e-51),
  qd_real( 7.9210657730021239e-01, -3.7087850202326467e-17,
       -1.4874457267228264e-33, 2.9323097289153505e-50),
  qd_real( 7.9023022143731003e-01, 2.3056905954954492e-17,
       1.4481080533260193e-33, -7.6725237057203488e-50),
  qd_real( 7.8834642762660623e-01, 3.4396993154059708e-17,
       1.7710623746737170e-33, 1.7084159098417402e-49),
  qd_real( 7.8645521359908577e-01, -9.7841429939305265e-18,
       3.3906063272445472e-34, 5.7269505320382577e-51),
  qd_real( 7.8455659715557524e-01, -8.5627965423173476e-18,
       -2.1106834459001849e-34, -1.6890322182469603e-50),
  qd_real( 7.8265059616657573e-01, 9.0745866975808825e-18,
       6.7623847404278666e-34, -1.7173237731987271e-50),
  qd_real( 7.8073722857209449e-01, -9.9198782066678806e-18,
       -2.1265794012162715e-36, 3.0772165598957647e-54),
  qd_real( 7.7881651238147598e-01, -2.4891385579973807e-17,
       6.7665497024807980e-35, -6.5218594281701332e-52),
  qd_real( 7.7688846567323244e-01, 7.7418602570672864e-18,
       -5.9986517872157897e-34, 3.0566548232958972e-50),
  qd_real( 7.7495310659487393e-01, -5.2209083189826433e-17,
       -9.6653593393686612e-34, 3.7027750076562569e-50),
  qd_real( 7.7301045336273699e-01, -3.2565907033649772e-17,
       1.3860807251523929e-33, -3.9971329917586022e-50),
  qd_real( 7.7106052426181382e-01, -4.4558442347769265e-17,
       -2.9863565614083783e-33, -6.8795262083596236e-50),
  qd_real( 7.6910333764557959e-01, 5.1546455184564817e-17,
       2.6142829553524292e-33, -1.6199023632773298e-49),
  qd_real( 7.6713891193582040e-01, -1.8885903683750782e-17,
       -1.3659359331495433e-33, -2.2538834962921934e-50),
  qd_real( 7.6516726562245896e-01, -3.2707225612534598e-17,
       1.1177117747079528e-33, -3.7005182280175715e-50),
  qd_real( 7.6318841726338127e-01, 2.6314748416750748e-18,
       1.4048039063095910e-34, 8.9601886626630321e-52),
  qd_real( 7.6120238548426178e-01, 3.5315510881690551e-17,
       1.2833566381864357e-33, 8.6221435180890613e-50),
  qd_real( 7.5920918897838807e-01, -3.8558842175523123e-17,
       2.9720241208332759e-34, -1.2521388928220163e-50),
  qd_real( 7.5720884650648457e-01, -1.9909098777335502e-17,
       3.9409283266158482e-34, 2.0744254207802976e-50),
  qd_real( 7.5520137689653655e-01, -1.9402238001823017e-17,
       -3.7756206444727573e-34, -2.1212242308178287e-50),
  qd_real( 7.5318679904361252e-01, -3.7937789838736540e-17,
       -6.7009539920231559e-34, -6.7128562115050214e-51),
  qd_real( 7.5116513190968637e-01, 4.3499761158645868e-17,
       2.5227718971102212e-33, -6.5969709212757102e-50),
  qd_real( 7.4913639452345937e-01, -4.4729078447011889e-17,
       -2.4206025249983768e-33, 1.1336681351116422e-49),
  qd_real( 7.4710060598018013e-01, 1.1874824875965430e-17,
       2.1992523849833518e-34, 1.1025018564644483e-50),
  qd_real( 7.4505778544146595e-01, 1.5078686911877863e-17,
       8.0898987212942471e-34, 8.2677958765323532e-50),
  qd_real( 7.4300795213512172e-01, -2.5144629669719265e-17,
       7.1128989512526157e-34, 3.0181629077821220e-50),
  qd_real( 7.4095112535495911e-01, -1.4708616952297345e-17,
       -4.9550433827142032e-34, 3.1434132533735671e-50),
  qd_real( 7.3888732446061511e-01, 3.4324874808225091e-17,
       -1.3706639444717610e-33, -3.3520827530718938e-51),
  qd_real( 7.3681656887736990e-01, -2.8932468101656295e-17,
       -3.4649887126202378e-34, -1.8484474476291476e-50),
  qd_real( 7.3473887809596350e-01, -3.4507595976263941e-17,
       -2.3718000676666409e-33, -3.9696090387165402e-50),
  qd_real( 7.3265427167241282e-01, 1.8918673481573520e-17,
       -1.5123719544119886e-33, -9.7922152011625728e-51),
  qd_real( 7.3056276922782759e-01, -2.9689959904476928e-17,
       -1.1276871244239744e-33, -3.0531520961539007e-50),
  qd_real( 7.2846439044822520e-01, 1.1924642323370718e-19,
       5.9001892316611011e-36, 1.2178089069502704e-52),
  qd_real( 7.2635915508434601e-01, -3.1917502443460542e-17,
       7.7047912412039396e-34, 4.1455880160182123e-50),
  qd_real( 7.2424708295146689e-01, 2.9198471334403004e-17,
       2.3027324968739464e-33, -1.2928820533892183e-51),
  qd_real( 7.2212819392921535e-01, -2.3871262053452047e-17,
       1.0636125432862273e-33, -4.4598638837802517e-50),
  qd_real( 7.2000250796138165e-01, -2.5689658854462333e-17,
       -9.1492566948567925e-34, 4.4403780801267786e-50),
  qd_real( 7.1787004505573171e-01, 2.7006476062511453e-17,
       -2.2854956580215348e-34, 9.1726903890287867e-51),
  qd_real( 7.1573082528381871e-01, -5.1581018476410262e-17,
       -1.3736271349300259e-34, -1.2734611344111297e-50),
  qd_real( 7.1358486878079364e-01, -4.2342504403133584e-17,
       -4.2690366101617268e-34, -2.6352370883066522e-50),
  qd_real( 7.1143219574521643e-01, 7.9643298613856813e-18,
       2.9488239510721469e-34, 1.6985236437666356e-50),
  qd_real( 7.0927282643886569e-01, -3.7597359110245730e-17,
       1.0613125954645119e-34, 8.9465480185486032e-51),
  qd_real( 7.0710678118654757e-01, -4.8336466567264567e-17,
       2.0693376543497068e-33, 2.4677734957341755e-50)
};

/* Computes sin(a) and cos(a) using Taylor series.
   Assumes |a| <= pi/2048.                           */
static void sincos_taylor(const qd_real &a, 
                          qd_real &sin_a, qd_real &cos_a) {
  const double thresh = 0.5 * qd_real::_eps * std::abs(to_double(a));
  qd_real p, s, t, x;

  if (a.is_zero()) {
    sin_a = 0.0;
    cos_a = 1.0;
    return;
  }

  x = -sqr(a);
  s = a;
  p = a;
  int i = 0;
  do {
    p *= x;
    t = p * inv_fact[i];
    s += t;
    i += 2;
  } while (i < n_inv_fact && std::abs(to_double(t)) > thresh);

  sin_a = s;
  cos_a = sqrt(1.0 - sqr(s));
}

static qd_real sin_taylor(const qd_real &a) {
  const double thresh = 0.5 * qd_real::_eps * std::abs(to_double(a));
  qd_real p, s, t, x;

  if (a.is_zero()) {
    return 0.0;
  }

  x = -sqr(a);
  s = a;
  p = a;
  int i = 0;
  do {
    p *= x;
    t = p * inv_fact[i];
    s += t;
    i += 2;
  } while (i < n_inv_fact && std::abs(to_double(t)) > thresh);

  return s;
}

static qd_real cos_taylor(const qd_real &a) {
  const double thresh = 0.5 * qd_real::_eps;
  qd_real p, s, t, x;

  if (a.is_zero()) {
    return 1.0;
  }

  x = -sqr(a);
  s = 1.0 + mul_pwr2(x, 0.5);
  p = x;
  int i = 1;
  do {
    p *= x;
    t = p * inv_fact[i];
    s += t;
    i += 2;
  } while (i < n_inv_fact && std::abs(to_double(t)) > thresh);

  return s;
}

qd_real sin(const qd_real &a) {

  /* Strategy.  To compute sin(x), we choose integers a, b so that

       x = s + a * (pi/2) + b * (pi/1024)

     and |s| <= pi/2048.  Using a precomputed table of
     sin(k pi / 1024) and cos(k pi / 1024), we can compute
     sin(x) from sin(s) and cos(s).  This greatly increases the
     convergence of the sine Taylor series.                          */

  if (a.is_zero()) {
    return 0.0;
  }

  // approximately reduce modulo 2*pi
  qd_real z = nint(a / qd_real::_2pi);
  qd_real r = a - qd_real::_2pi * z;

  // approximately reduce modulo pi/2 and then modulo pi/1024
  double q = std::floor(r.x[0] / qd_real::_pi2[0] + 0.5);
  qd_real t = r - qd_real::_pi2 * q;
  int j = static_cast<int>(q);
  q = std::floor(t.x[0] / _pi1024[0] + 0.5);
  t -= _pi1024 * q;
  int k = static_cast<int>(q);
  int abs_k = std::abs(k);

  if (j < -2 || j > 2) {
    qd_real::error("(qd_real::sin): Cannot reduce modulo pi/2.");
    return qd_real::_nan;
  }

  if (abs_k > 256) {
    qd_real::error("(qd_real::sin): Cannot reduce modulo pi/1024.");
    return qd_real::_nan;
  }

  if (k == 0) {
    switch (j) {
      case 0:
        return sin_taylor(t);
      case 1:
        return cos_taylor(t);
      case -1:
        return -cos_taylor(t);
      default:
        return -sin_taylor(t);
    }
  }

  qd_real sin_t, cos_t;
  qd_real u = cos_table[abs_k-1];
  qd_real v = sin_table[abs_k-1];
  sincos_taylor(t, sin_t, cos_t);

  if (j == 0) {
    if (k > 0) {
      r = u * sin_t + v * cos_t;
    } else {
      r = u * sin_t - v * cos_t;
    }
  } else if (j == 1) {
    if (k > 0) {
      r = u * cos_t - v * sin_t;
    } else {
      r = u * cos_t + v * sin_t;
    }
  } else if (j == -1) {
    if (k > 0) {
      r = v * sin_t - u * cos_t;
    } else {
      r = - u * cos_t - v * sin_t;
    }
  } else {
    if (k > 0) {
      r = - u * sin_t - v * cos_t;
    } else {
      r = v * cos_t - u * sin_t;
    }
  }

  return r;
}

qd_real cos(const qd_real &a) {

  if (a.is_zero()) {
    return 1.0;
  }

  // approximately reduce modulo 2*pi
  qd_real z = nint(a / qd_real::_2pi);
  qd_real r = a - qd_real::_2pi * z;

  // approximately reduce modulo pi/2 and then modulo pi/1024
  double q = std::floor(r.x[0] / qd_real::_pi2.x[0] + 0.5);
  qd_real t = r - qd_real::_pi2 * q;
  int j = static_cast<int>(q);
  q = std::floor(t.x[0] / _pi1024.x[0] + 0.5);
  t -= _pi1024 * q;
  int k = static_cast<int>(q);
  int abs_k = std::abs(k);

  if (j < -2 || j > 2) {
    qd_real::error("(qd_real::cos): Cannot reduce modulo pi/2.");
    return qd_real::_nan;
  }

  if (abs_k > 256) {
    qd_real::error("(qd_real::cos): Cannot reduce modulo pi/1024.");
    return qd_real::_nan;
  }

  if (k == 0) {
    switch (j) {
      case 0:
        return cos_taylor(t);
      case 1:
        return -sin_taylor(t);
      case -1:
        return sin_taylor(t);
      default:
        return -cos_taylor(t);
    }
  }

  qd_real sin_t, cos_t;
  sincos_taylor(t, sin_t, cos_t);

  qd_real u = cos_table[abs_k-1];
  qd_real v = sin_table[abs_k-1];

  if (j == 0) {
    if (k > 0) {
      r = u * cos_t - v * sin_t;
    } else {
      r = u * cos_t + v * sin_t;
    }
  } else if (j == 1) {
    if (k > 0) {
      r = - u * sin_t - v * cos_t;
    } else {
      r = v * cos_t - u * sin_t;
    }
  } else if (j == -1) {
    if (k > 0) {
      r = u * sin_t + v * cos_t;
    } else {
      r = u * sin_t - v * cos_t;
    }
  } else {
    if (k > 0) {
      r = v * sin_t - u * cos_t;
    } else {
      r = - u * cos_t - v * sin_t;
    }
  }

  return r;
}

void sincos(const qd_real &a, qd_real &sin_a, qd_real &cos_a) {

  if (a.is_zero()) {
    sin_a = 0.0;
    cos_a = 1.0;
    return;
  }

  // approximately reduce by 2*pi
  qd_real z = nint(a / qd_real::_2pi);
  qd_real t = a - qd_real::_2pi * z;

  // approximately reduce by pi/2 and then by pi/1024.
  double q = std::floor(t.x[0] / qd_real::_pi2.x[0] + 0.5);
  t -= qd_real::_pi2 * q;
  int j = static_cast<int>(q);
  q = std::floor(t.x[0] / _pi1024.x[0] + 0.5);
  t -= _pi1024 * q;
  int k = static_cast<int>(q);
  int abs_k = std::abs(k);

  if (j < -2 || j > 2) {
    qd_real::error("(qd_real::sincos): Cannot reduce modulo pi/2.");
    cos_a = sin_a = qd_real::_nan;
    return;
  }

  if (abs_k > 256) {
    qd_real::error("(qd_real::sincos): Cannot reduce modulo pi/1024.");
    cos_a = sin_a = qd_real::_nan;
    return;
  }

  qd_real sin_t, cos_t;
  sincos_taylor(t, sin_t, cos_t);

  if (k == 0) {
    if (j == 0) {
      sin_a = sin_t;
      cos_a = cos_t;
    } else if (j == 1) {
      sin_a = cos_t;
      cos_a = -sin_t;
    } else if (j == -1) {
      sin_a = -cos_t;
      cos_a = sin_t;
    } else {
      sin_a = -sin_t;
      cos_a = -cos_t;
    }
    return;
  }

  qd_real u = cos_table[abs_k-1];
  qd_real v = sin_table[abs_k-1];

  if (j == 0) {
    if (k > 0) {
      sin_a = u * sin_t + v * cos_t;
      cos_a = u * cos_t - v * sin_t;
    } else {
      sin_a = u * sin_t - v * cos_t;
      cos_a = u * cos_t + v * sin_t;
    }
  } else if (j == 1) {
    if (k > 0) {
      cos_a = - u * sin_t - v * cos_t;
      sin_a = u * cos_t - v * sin_t;
    } else {
      cos_a = v * cos_t - u * sin_t;
      sin_a = u * cos_t + v * sin_t;
    }
  } else if (j == -1) {
    if (k > 0) {
      cos_a = u * sin_t + v * cos_t;
      sin_a =  v * sin_t - u * cos_t;
    } else {
      cos_a = u * sin_t - v * cos_t;
      sin_a = - u * cos_t - v * sin_t;
    }
  } else {
    if (k > 0) {
      sin_a = - u * sin_t - v * cos_t;
      cos_a = v * sin_t - u * cos_t;
    } else {
      sin_a = v * cos_t - u * sin_t;
      cos_a = - u * cos_t - v * sin_t;
    }
  }
}

qd_real atan(const qd_real &a) {
  return atan2(a, qd_real(1.0));
}

qd_real atan2(const qd_real &y, const qd_real &x) {
  /* Strategy: Instead of using Taylor series to compute 
     arctan, we instead use Newton's iteration to solve
     the equation

        sin(z) = y/r    or    cos(z) = x/r

     where r = sqrt(x^2 + y^2).
     The iteration is given by

        z' = z + (y - sin(z)) / cos(z)          (for equation 1)
        z' = z - (x - cos(z)) / sin(z)          (for equation 2)

     Here, x and y are normalized so that x^2 + y^2 = 1.
     If |x| > |y|, then first iteration is used since the 
     denominator is larger.  Otherwise, the second is used.
  */

  if (x.is_zero()) {
    
    if (y.is_zero()) {
      /* Both x and y is zero. */
      qd_real::error("(qd_real::atan2): Both arguments zero.");
      return qd_real::_nan;
    }

    return (y.is_positive()) ? qd_real::_pi2 : -qd_real::_pi2;
  } else if (y.is_zero()) {
    return (x.is_positive()) ? qd_real(0.0) : qd_real::_pi;
  }

  if (x == y) {
    return (y.is_positive()) ? qd_real::_pi4 : -qd_real::_3pi4;
  }

  if (x == -y) {
    return (y.is_positive()) ? qd_real::_3pi4 : -qd_real::_pi4;
  }

  qd_real r = sqrt(sqr(x) + sqr(y));
  qd_real xx = x / r;
  qd_real yy = y / r;

  /* Compute double precision approximation to atan. */
  qd_real z = std::atan2(to_double(y), to_double(x));
  qd_real sin_z, cos_z;

  if (std::abs(xx.x[0]) > std::abs(yy.x[0])) {
    /* Use Newton iteration 1.  z' = z + (y - sin(z)) / cos(z)  */
    sincos(z, sin_z, cos_z);
    z += (yy - sin_z) / cos_z;
    sincos(z, sin_z, cos_z);
    z += (yy - sin_z) / cos_z;
    sincos(z, sin_z, cos_z);
    z += (yy - sin_z) / cos_z;
  } else {
    /* Use Newton iteration 2.  z' = z - (x - cos(z)) / sin(z)  */
    sincos(z, sin_z, cos_z);
    z -= (xx - cos_z) / sin_z;
    sincos(z, sin_z, cos_z);
    z -= (xx - cos_z) / sin_z;
    sincos(z, sin_z, cos_z);
    z -= (xx - cos_z) / sin_z;
  }

  return z;
}


qd_real drem(const qd_real &a, const qd_real &b) {
  qd_real n = nint(a/b);
  return (a - n * b);
}

qd_real divrem(const qd_real &a, const qd_real &b, qd_real &r) {
  qd_real n = nint(a/b);
  r = a - n * b;
  return n;
}

qd_real tan(const qd_real &a) {
  qd_real s, c;
  sincos(a, s, c);
  return s/c;
}

qd_real asin(const qd_real &a) {
  qd_real abs_a = abs(a);

  if (abs_a > 1.0) {
    qd_real::error("(qd_real::asin): Argument out of domain.");
    return qd_real::_nan;
  }

  if (abs_a.is_one()) {
    return (a.is_positive()) ? qd_real::_pi2 : -qd_real::_pi2;
  }

  return atan2(a, sqrt(1.0 - sqr(a)));
}

qd_real acos(const qd_real &a) {
  qd_real abs_a = abs(a);

  if (abs_a > 1.0) {
    qd_real::error("(qd_real::acos): Argument out of domain.");
    return qd_real::_nan;
  }

  if (abs_a.is_one()) {
    return (a.is_positive()) ? qd_real(0.0) : qd_real::_pi;
  }

  return atan2(sqrt(1.0 - sqr(a)), a);
}
 
qd_real sinh(const qd_real &a) {
  if (a.is_zero()) {
    return 0.0;
  }

  if (abs(a) > 0.05) {
    qd_real ea = exp(a);
    return mul_pwr2(ea - inv(ea), 0.5);
  }

  /* Since a is small, using the above formula gives
     a lot of cancellation.   So use Taylor series. */
  qd_real s = a;
  qd_real t = a;
  qd_real r = sqr(t);
  double m = 1.0;
  double thresh = std::abs(to_double(a) * qd_real::_eps);

  do {
    m += 2.0;
    t *= r;
    t /= (m-1) * m;

    s += t;
  } while (abs(t) > thresh);

  return s;
}

qd_real cosh(const qd_real &a) {
  if (a.is_zero()) {
    return 1.0;
  }

  qd_real ea = exp(a);
  return mul_pwr2(ea + inv(ea), 0.5);
}

qd_real tanh(const qd_real &a) {
  if (a.is_zero()) {
    return 0.0;
  }

  if (std::abs(to_double(a)) > 0.05) {
    qd_real ea = exp(a);
    qd_real inv_ea = inv(ea);
    return (ea - inv_ea) / (ea + inv_ea);
  } else {
    qd_real s, c;
    s = sinh(a);
    c = sqrt(1.0 + sqr(s));
    return s / c;
  }
}

void sincosh(const qd_real &a, qd_real &s, qd_real &c) {
  if (std::abs(to_double(a)) <= 0.05) {
    s = sinh(a);
    c = sqrt(1.0 + sqr(s));
  } else {
    qd_real ea = exp(a);
    qd_real inv_ea = inv(ea);
    s = mul_pwr2(ea - inv_ea, 0.5);
    c = mul_pwr2(ea + inv_ea, 0.5);
  }
}

qd_real asinh(const qd_real &a) {
  return log(a + sqrt(sqr(a) + 1.0));
}

qd_real acosh(const qd_real &a) {
  if (a < 1.0) {
    qd_real::error("(qd_real::acosh): Argument out of domain.");
    return qd_real::_nan;
  }

  return log(a + sqrt(sqr(a) - 1.0));
}

qd_real atanh(const qd_real &a) {
  if (abs(a) >= 1.0) {
    qd_real::error("(qd_real::atanh): Argument out of domain.");
    return qd_real::_nan;
  }

  return mul_pwr2(log((1.0 + a) / (1.0 - a)), 0.5);
}

QD_API qd_real fmod(const qd_real &a, const qd_real &b) {
  qd_real n = aint(a / b);
  return (a - b * n);
}

QD_API qd_real qdrand() {
  static const double m_const = 4.6566128730773926e-10;  /* = 2^{-31} */
  double m = m_const;
  qd_real r = 0.0;
  double d;

  /* Strategy:  Generate 31 bits at a time, using lrand48 
     random number generator.  Shift the bits, and repeat
     7 times. */

  for (int i = 0; i < 7; i++, m *= m_const) {
    d = std::rand() * m;
    r += d;
  }

  return r;
}


/* polyeval(c, n, x)
   Evaluates the given n-th degree polynomial at x.
   The polynomial is given by the array of (n+1) coefficients. */
qd_real polyeval(const qd_real *c, int n, const qd_real &x) {
  /* Just use Horner's method of polynomial evaluation. */
  qd_real r = c[n];
  
  for (int i = n-1; i >= 0; i--) {
    r *= x;
    r += c[i];
  }

  return r;
}

/* polyroot(c, n, x0)
   Given an n-th degree polynomial, finds a root close to 
   the given guess x0.  Note that this uses simple Newton
   iteration scheme, and does not work for multiple roots.  */
QD_API qd_real polyroot(const qd_real *c, int n, 
    const qd_real &x0, int max_iter, double thresh) {
  qd_real x = x0;
  qd_real f;
  qd_real *d = new qd_real[n];
  bool conv = false;
  int i;
  double max_c = std::abs(to_double(c[0]));
  double v;

  if (thresh == 0.0) thresh = qd_real::_eps;

  /* Compute the coefficients of the derivatives. */
  for (i = 1; i <= n; i++) {
    v = std::abs(to_double(c[i]));
    if (v > max_c) max_c = v;
    d[i-1] = c[i] * static_cast<double>(i);
  }
  thresh *= max_c;

  /* Newton iteration. */
  for (i = 0; i < max_iter; i++) {
    f = polyeval(c, n, x);

    if (abs(f) < thresh) {
      conv = true;
      break;
    }
    x -= (f / polyeval(d, n-1, x));
  }
  delete [] d;

  if (!conv) {
    qd_real::error("(qd_real::polyroot): Failed to converge.");
    return qd_real::_nan;
  }

  return x;
}

qd_real qd_real::debug_rand() {
  if (std::rand() % 2 == 0)
    return qdrand();

  int expn = 0;
  qd_real a = 0.0;
  double d;
  for (int i = 0; i < 4; i++) {
    d = std::ldexp(std::rand() / static_cast<double>(RAND_MAX), -expn);
    a += d;
    expn = expn + 54 + std::rand() % 200;
  }
  return a;
}