"

 COPYRIGHT (c) 2017 by eXept Software AG

	      All Rights Reserved

 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

		FMax FMin InvFact'

	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

#ifdef __win32__

# ifndef isinf

#  define isinf(x) \

	((((unsigned int *)(&x))[0] == 0x00000000) && \

	 ((((unsigned int *)(&x))[1] & 0x7FF00000) == 0x7FF00000))

# endif

#endif

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

# ifndef _FPU_EXTENDED

#  define _FPU_EXTENDED 0x0300

# endif

# ifndef _FPU_DOUBLE

#  define _FPU_DOUBLE 0x0200

# endif

# if defined( __win32__ ) && (defined( __BORLANDC__ ) || defined( __VISUALC__ ))

#  define fpu_fix_start(old_cw_ptr)\

    {\

	*old_cw_ptr = _control87(0, 0); \

	_control87(_FPU_DOUBLE, _FPU_EXTENDED);\

    }

#  define fpu_fix_end(old_cw_ptr)\

    {\

	_control87(*old_cw_ptr, _FPU_EXTENDED);\

    }

# 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)\

    {\

	volatile unsigned short cw, new_cw;\

	_FPU_GETCW(cw);\

	new_cw = (cw & ~_FPU_EXTENDED) | _FPU_DOUBLE;\

	_FPU_SETCW(new_cw);\

	*old_cw_ptr = cw;\

    }

#  define fpu_fix_end(old_cw_ptr)\

    {\

	volatile unsigned short cw = *old_cw_ptr;\

	_FPU_SETCW(cw);\

    }

# endif

#endif

struct qd_real {

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

};

struct __quadDoubleStruct {

	STX_OBJ_HEADER

#ifdef __NEED_DOUBLE_ALIGN

	__FILLTYPE_DOUBLE f_filler;

#endif

	double d_quadDoubleValue[4];

};

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

#ifndef __isQuadDouble

# define __isQuadDouble(o) \

	(__Class(o) == @global(QuadDouble))

#endif

#ifndef __qIsQuadDouble

# define __qIsQuadDouble(o) \

	(__qClass(o) == @global(QuadDouble))

#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;   \

}

// sigh: not all compilers (borland) support inline functions;

// therefore we have to use macros...

// sigh2: c-macros are unhygienic - to avoid catching/hiding variable bindings,

// use different names in each macros (i.e. a_xxx)

// 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_1, a_1, b_1, err_1)\

{\

    double s_1 = (a_1) + (b_1);\

    (err_1) = (b_1) - (s_1 - (a_1));\

    (rslt_1) = s_1; \

}

#define m_quick_two_diff(rslt_2, a_2, b_2, err_2)\

{\

    double s_2 = (a_2) - (b_2);\

    (err_2) = ((a_2) - s_2) - (b_2);\

    (rslt_2) = s_2;\

}

#define m_two_sum(rslt_3, a_3, b_3, err_3)\

{\

    double s_3 = a_3 + b_3;\

    double bb_3 = s_3 - a_3;\

    err_3 = (a_3 - (s_3 - bb_3)) + (b_3 - bb_3);\

    rslt_3 = s_3;\

}

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

#define m_two_diff(rslt_4, a_4, b_4, err_4)\

{\

    double s_4 = a_4 - b_4;\

    double bb_4 = s_4 - a_4;\

    err_4 = (a_4 - (s_4 - bb_4)) - (b_4 + bb_4);\

    rslt_4 = s_4;\

}

#define m_three_sum(a_5, b_5, c_5)\

{ \

    double t1_5, t2_5, t3_5; \

    m_two_sum(t1_5, a_5, b_5, t2_5); \

    m_two_sum(a_5, c_5, t1_5, t3_5); \

    m_two_sum(b_5, t2_5, t3_5, c_5); \

}

#define m_three_sum2(a_6, b_6, c_6)\

{\

    double t1_6, t2_6, t3_6;\

    m_two_sum(t1_6, a_6, b_6, t2_6);\

    m_two_sum(a_6, c_6, t1_6, t3_6);\

    b_6 = t2_6 + t3_6;\

}

#ifndef QD_FMS

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

#define m_split(a_7, hi_7, lo_7)\

{\

    double temp_7;\

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

	a_7 *= 3.7252902984619140625e-09;\

	temp_7 = _QD_SPLITTER * a_7;\

	hi_7 = temp_7 - (temp_7 - a_7);\

	lo_7 = a_7 - hi_7;\

	hi_7 *= 268435456.0;\

	lo_7 *= 268435456.0;\

    } else {\

	temp_7 = _QD_SPLITTER * a_7;\

	hi_7 = temp_7 - (temp_7 - a_7);\

	lo_7 = a_7 - hi_7;\

    }\

}

#endif

#ifdef QD_FMS

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

#define m_two_prod(rslt_8, a_8, b_8, err_8)\

{\

    double p_8 = a_8 * b_8;\

    err_8 = QD_FMS(a_8, b_8, p_8);\

    rslt_8 = p_8; \

}

#else

#define m_two_prod(rslt_8, a_8, b_8, err_8)\

{\

    double a_hi_8, a_lo_8, b_hi_8, b_lo_8;\

    double p_8 = a_8 * b_8;\

    m_split(a_8, a_hi_8, a_lo_8);\

    m_split(b_8, b_hi_8, b_lo_8);\

    err_8 = ((a_hi_8 * b_hi_8 - p_8) + a_hi_8 * b_lo_8 + a_lo_8 * b_hi_8) + a_lo_8 * b_lo_8;\

    rslt_8 = p_8; \

}

#endif

#ifdef QD_FMS

#define m_two_sqr(rslt_9, a_9, err_9)\

{\

    double p_9 = a_9 * a_9;\

    err_9 = QD_FMS(a_9, a_9, p_9);\

    rslt_9 = p_9;\

}

#else

#define m_two_sqr(rslt_9, a_9, err_9)\

{\

    double hi_9, lo_9;\

    double q_9 = a_9 * a_9;\

    m_split(a_9, hi_9, lo_9);\

    err_9 = ((hi_9 * hi_9 - q_9) + 2.0 * hi_9 * lo_9) + lo_9 * lo_9;\

    rslt_9 = q_9;\

}

#endif

#define m_renorm4(c0_10, c1_10, c2_10, c3_10)\

{\

    double s0_10, s1_10, s2_10 = 0.0, s3_10 = 0.0;\

\

    if (! isinf(c0_10)) { \

\

	m_quick_two_sum(s0_10, c2_10, c3_10, c3_10);\

	m_quick_two_sum(s0_10, c1_10, s0_10, c2_10);\

	m_quick_two_sum(c0_10, c0_10, s0_10, c1_10);\

\

	s0_10 = c0_10;\

	s1_10 = c1_10;\

	if (s1_10 != 0.0) {\

	     m_quick_two_sum(s1_10, s1_10, c2_10, s2_10);\

	    if (s2_10 != 0.0) {\

		 m_quick_two_sum(s2_10, s2_10, c3_10, s3_10);\

	    } else {\

		 m_quick_two_sum(s1_10, s1_10, c3_10, s2_10);\

	    }\

	} else {\

	    m_quick_two_sum(s0_10, s0_10, c2_10, s1_10);\

	    if (s1_10 != 0.0) {\

		 m_quick_two_sum(s1_10, s1_10, c3_10, s2_10);\

	    } else {\

		 m_quick_two_sum(s0_10, s0_10, c3_10, s1_10);\

	    }\

	}\

\

	c0_10 = s0_10;\

	c1_10 = s1_10;\

	c2_10 = s2_10;\

	c3_10 = s3_10;\

    }\

}

#define m_renorm5(c0_11, c1_11, c2_11, c3_11, c4_11)\

{\

    double s0_11, s1_11, s2_11 = 0.0, s3_11 = 0.0; \

\

    if (! isinf(c0_11)) { \

	m_quick_two_sum(s0_11, c3_11, c4_11, c4_11); \

	m_quick_two_sum(s0_11, c2_11, s0_11, c3_11); \

	m_quick_two_sum(s0_11, c1_11, s0_11, c2_11); \

	m_quick_two_sum(c0_11, c0_11, s0_11, c1_11); \

\

	s0_11 = c0_11; \

	s1_11 = c1_11; \

\

	m_quick_two_sum(s0_11, c0_11, c1_11, s1_11); \

	if (s1_11 != 0.0) { \

	    m_quick_two_sum(s1_11, s1_11, c2_11, s2_11); \

	    if (s2_11 != 0.0) { \

		m_quick_two_sum(s2_11 ,s2_11, c3_11, s3_11); \

		if (s3_11 != 0.0) {\

		    s3_11 += c4_11; \

		} else {\

		    s2_11 += c4_11;\

		}\

	    } else { \

		m_quick_two_sum(s1_11, s1_11, c3_11, s2_11); \

		if (s2_11 != 0.0) {\

		    m_quick_two_sum(s2_11, s2_11, c4_11, s3_11); \

		} else { \

		    m_quick_two_sum(s1_11, s1_11, c4_11, s2_11); \

		} \

	    } \

	} else { \

	    m_quick_two_sum(s0_11,s0_11, c2_11, s1_11); \

	    if (s1_11 != 0.0) { \

		m_quick_two_sum(s1_11,s1_11, c3_11, s2_11); \

		if (s2_11 != 0.0) {\

		    m_quick_two_sum(s2_11,s2_11, c4_11, s3_11); \

		} else { \

		    m_quick_two_sum(s1_11 ,s1_11, c4_11, s2_11); \

		} \

	    } else { \

		m_quick_two_sum(s0_11,s0_11, c3_11, s1_11); \

		if (s1_11 != 0.0) { \

		    m_quick_two_sum(s1_11,s1_11, c4_11, s2_11); \

		} else { \

		    m_quick_two_sum(s0_11,s0_11, c4_11, s1_11); \

		} \

	    } \

	} \

 \

	c0_11 = s0_11; \

	c1_11 = s1_11; \

	c2_11 = s2_11; \

	c3_11 = s3_11; \

    } \

}

%}

! !

!QDouble primitiveFunctions!

%{

#if 0

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

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

inline double

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| */

inline double

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).  */

inline double

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).  */

inline double

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 */

inline void

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). */

inline double

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. */

inline double

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. */

inline double

nint(double d) {

  if (d == floor(d))

    return d;

  return floor(d + 0.5);

}

/* Computes the truncated integer. */

inline double

aint(double d) {

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

}

inline void

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;