"

 COPYRIGHT (c) 2017 by eXept Software AG

 This software is furnished under a license and may be used

 only in accordance with the terms of that license and with the

 inclusion of the above copyright notice.   This software may not

 be provided or otherwise made available to, or used by, any

 other person.  No title to or ownership of the software is

 hereby transferred.

"

"{ Package: 'stx:libbasic2' }"

"{ NameSpace: Smalltalk }"

LimitedPrecisionReal variableByteSubclass:#QDouble

	instanceVariableNames:''

	classVariableNames:'DefaultPrintFormat QDoubleZero QDoubleOne Pi E Ln2 Ln10 Epsilon'

	poolDictionaries:''

	category:'Magnitude-Numbers'

!

!QDouble primitiveDefinitions!

%{

#include <stdio.h>

#include <errno.h>

#define __USE_ISOC9X 1

#define __USE_ISOC99 1

#include <math.h>

/*

 * on some systems, errno is a macro ... check for it here

 */

#ifndef errno

 extern errno;

#endif

#if !defined (__win32__)

# include <locale.h>

#endif

#if defined (__aix__)

# include <float.h>

#endif

#if defined(__irix__) || defined(__solaris__) || defined(__sunos__)

# include <nan.h>

#endif

#if defined(__linux__)

# ifndef NAN

#  include <bits/nan.h>

# endif

#endif

#if defined(__x86__) || defined(__x86_64__)

# ifndef _FPU_EXTENDED

#  define _FPU_EXTENDED 0x0300

# endif

# ifndef _FPU_DOUBLE

#  define _FPU_DOUBLE 0x0200

# endif

#  define fpu_fix_start(old_cw_ptr)\

    {\

    }

#  define fpu_fix_end(old_cw_ptr)\

    {\

    }

# else // assume MINGW, GCC or CLANG

#  ifndef _FPU_GETCW

#   define _FPU_GETCW(x) asm volatile ("fnstcw %0":"=m" (x));

#  endif

#  ifndef _FPU_SETCW

#   define _FPU_SETCW(x) asm volatile ("fldcw %0": :"m" (x));

#  endif

#  define fpu_fix_start(old_cw_ptr)\

    {\

    }

#  define fpu_fix_end(old_cw_ptr)\

    {\

    }

# endif

#endif

struct qd_real {

    double x[4];    /* The Components. */

};

struct __quadDoubleStruct {

#ifdef __NEED_DOUBLE_ALIGN

#endif

};

#define __QuadDoubleInstPtr(obj)      ((struct __quadDoubleStruct *)(__objPtr(obj)))

#ifndef __isQuadDouble

# define __isQuadDouble(o) \

#endif

#ifndef __qIsQuadDouble

# define __qIsQuadDouble(o) \

#endif

#define __qNew_qdReal(newQD, d0,d1,d2,d3) { \

    __qNew(newQD, sizeof(struct __quadDoubleStruct)); \

    __stx_setClass(newQD, QDouble);                \

    __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[0] = d0;   \

    __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[1] = d1;   \

    __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[2] = d2;   \

    __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[3] = d3;   \

}

// qd_real(c0, c1, c2, c3);

#define _QD_SPLITTER 134217729.0               // = 2^27 + 1

#define _QD_SPLIT_THRESH 6.69692879491417e+299 // = 2^996

#define m_quick_two_sum(rslt, a, b, err) {\

  double s = a + b;\

  err = b - (s - a);\

  rslt = s; \

}

#define m_quick_two_diff(rslt, a, b, err) {\

  double s = a - b;\

  err = (a - s) - b;\

  rslt = s;\

}

#define m_two_sum(rslt, a, b, err) {\

  double s = a + b;\

  double bb = s - a;\

  err = (a - (s - bb)) + (b - bb);\

  rslt = s;\

}

/* Computes fl(a-b) and err(a-b).  */

#define m_two_diff(rslt, a, b, err) {\

  double s = a - b;\

  double bb = s - a;\

  err = (a - (s - bb)) - (b + bb);\

  rslt = s;\

}

#define m_three_sum(a, b, c) { \

  double t1, t2, t3; \

  m_two_sum(t1, a, b, t2); \

  m_two_sum(a, c, t1, t3); \

  m_two_sum(b, t2, t3, c); \

}

#define m_three_sum2(a, b, c) {\

  double t1, t2, t3;\

  m_two_sum(t1, a, b, t2);\

  m_two_sum(a, c, t1, t3);\

  b = t2 + t3;\

}

#ifndef QD_FMS

/* Computes high word and lo word of a */

#define m_split(a, hi, lo) {\

  double temp;\

  if (a > _QD_SPLIT_THRESH || a < -_QD_SPLIT_THRESH) {\

    a *= 3.7252902984619140625e-09;\

    temp = _QD_SPLITTER * a;\

    hi = temp - (temp - a);\

    lo = a - hi;\

    hi *= 268435456.0;\

    lo *= 268435456.0;\

  } else {\

    temp = _QD_SPLITTER * a;\

    hi = temp - (temp - a);\

    lo = a - hi;\

  }\

}

#endif

#ifdef QD_FMS

/* Computes fl(a*b) and err(a*b). */

#define m_two_prod(rslt, a, b, err) {\

  double p = a * b;\

  err = QD_FMS(a, b, p);\

  rslt = p; \

}

#else

#define m_two_prod(rslt, a, b, err) {\

  double a_hi, a_lo, b_hi, b_lo;\

  double p = a * b;\

  m_split(a, a_hi, a_lo);\

  m_split(b, b_hi, b_lo);\

  err = ((a_hi * b_hi - p) + a_hi * b_lo + a_lo * b_hi) + a_lo * b_lo;\

  rslt = p; \

}

#endif

#ifdef QD_FMS

#define m_two_sqr(rslt, a, err) {\

  double p = a * a;\

  err = QD_FMS(a, a, p);\

  rslt = p;\

}

#else

#define m_two_sqr(rslt, a, err) {\

  double hi, lo;\

  double q = a * a;\

  m_split(a, hi, lo);\

  err = ((hi * hi - q) + 2.0 * hi * lo) + lo * lo;\

  rslt = q;\

}

#endif

#define m_renorm4(c0, c1, c2, c3) {\

  double s0, s1, s2 = 0.0, s3 = 0.0;\

\

    if (! isinf(c0)) { \

\

\

\

    }\

}

#define m_renorm5(c0, c1, c2, c3, c4) { \

    double s0, s1, s2 = 0.0, s3 = 0.0; \

\

    if (! isinf(c0)) { \

\

\

 \

    } \

}

%}

! !

!QDouble primitiveFunctions!

%{

/*********** Basic Functions ************/

/* Computes fl(a+b) and err(a+b).  Assumes |a| >= |b|. */

quick_two_sum(double a, double b, double *errPtr) {

  double s = a + b;

  *errPtr = b - (s - a);

  return s;

}

/* Computes fl(a-b) and err(a-b).  Assumes |a| >= |b| */

quick_two_diff(double a, double b, double *errPtr) {

  double s = a - b;

  *errPtr = (a - s) - b;

  return s;

}

/* Computes fl(a+b) and err(a+b).  */

two_sum(double a, double b, double *errPtr) {

  double s = a + b;

  double bb = s - a;

  *errPtr = (a - (s - bb)) + (b - bb);

  return s;

}

/* Computes fl(a-b) and err(a-b).  */

two_diff(double a, double b, double *errPtr) {

  double s = a - b;

  double bb = s - a;

  *errPtr = (a - (s - bb)) - (b + bb);

  return s;

}

#ifndef QD_FMS

/* Computes high word and lo word of a */

split(double a, double *hiPtr, double *loPtr) {

  double temp;

  if (a > _QD_SPLIT_THRESH || a < -_QD_SPLIT_THRESH) {

    a *= 3.7252902984619140625e-09;  // 2^-28

    temp = _QD_SPLITTER * a;

    *hiPtr = temp - (temp - a);

    *loPtr = a - *hiPtr;

    *hiPtr *= 268435456.0;          // 2^28

    *loPtr *= 268435456.0;          // 2^28

  } else {

    temp = _QD_SPLITTER * a;

    *hiPtr = temp - (temp - a);

    *loPtr = a - *hiPtr;

  }

}

#endif

/* Computes fl(a*b) and err(a*b). */

two_prod(double a, double b, double *errPtr) {

#ifdef QD_FMS

  double p = a * b;

  *errPtr = QD_FMS(a, b, p);

  return p;

#else

  double a_hi, a_lo, b_hi, b_lo;

  double p = a * b;

  split(a, &a_hi, &a_lo);

  split(b, &b_hi, &b_lo);

  *errPtr = ((a_hi * b_hi - p) + a_hi * b_lo + a_lo * b_hi) + a_lo * b_lo;

  return p;

#endif

}

/* Computes fl(a*a) and err(a*a).  Faster than the above method. */

two_sqr(double a, double *errPtr) {

#ifdef QD_FMS

  double p = a * a;

  *errPtr = QD_FMS(a, a, p);

  return p;

#else

  double hi, lo;

  double q = a * a;

  split(a, &hi, &lo);

  *errPtr = ((hi * hi - q) + 2.0 * hi * lo) + lo * lo;

  return q;

#endif

}

/* Computes the nearest integer to d. */

nint(double d) {

  if (d == floor(d))

    return d;

  return floor(d + 0.5);

}

/* Computes the truncated integer. */

aint(double d) {

  return (d >= 0.0) ? floor(d) : ceil(d);

}

renorm4(double *c0Ptr, double *c1Ptr, double *c2Ptr, double *c3Ptr) {

  double s0, s1, s2 = 0.0, s3 = 0.0;

  double c0 = *c0Ptr;

  if (isinf(c0)) return;

  s0 = quick_two_sum(*c2Ptr, *c3Ptr, c3Ptr);

  s0 = quick_two_sum(*c1Ptr, s0, c2Ptr);

  c0 = quick_two_sum(c0, s0, c1Ptr);

  s0 = c0;

  s1 = *c1Ptr;

  if (s1 != 0.0) {

    s1 = quick_two_sum(s1, *c2Ptr, &s2);

    if (s2 != 0.0)

      s2 = quick_two_sum(s2, *c3Ptr, &s3);

    else

      s1 = quick_two_sum(s1, *c3Ptr, &s2);

  } else {

    s0 = quick_two_sum(s0, *c2Ptr, &s1);

    if (s1 != 0.0)

      s1 = quick_two_sum(s1, *c3Ptr, &s2);

    else

      s0 = quick_two_sum(s0, *c3Ptr, &s1);

  }

  *c0Ptr = s0;

  *c1Ptr = s1;

  *c2Ptr = s2;

  *c3Ptr = s3;

}

renorm5(double *c0Ptr, double *c1Ptr, double *c2Ptr, double *c3Ptr, double *c4Ptr) {

  double s0, s1, s2 = 0.0, s3 = 0.0;

  if (isinf(*c0Ptr)) return;

  s0 = quick_two_sum(*c3Ptr, *c4Ptr, c4Ptr);

  s0 = quick_two_sum(*c2Ptr, s0, c3Ptr);

  s0 = quick_two_sum(*c1Ptr, s0, c2Ptr);

  *c0Ptr = quick_two_sum(*c0Ptr, s0, c1Ptr);

  s0 = *c0Ptr;

  s1 = *c1Ptr;

  s0 = quick_two_sum(*c0Ptr, *c1Ptr, &s1);

  if (s1 != 0.0) {

    s1 = quick_two_sum(s1, *c2Ptr, &s2);

    if (s2 != 0.0) {

      s2 =quick_two_sum(s2, *c3Ptr, &s3);

      if (s3 != 0.0)

      else

    } else {

      s1 = quick_two_sum(s1, *c3Ptr, &s2);

      if (s2 != 0.0)

      else