"

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

}

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

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

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

}

}

}

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

}

}

}

#ifndef QD_FMS

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

}

#endif

#ifdef QD_FMS

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

}

#else

}

#endif

#ifdef QD_FMS

}

#else

}

#endif

\

\

	    } else {\

	    }\

	} else {\

	    } else {\

	    }\

	}\

\

    }\

}

\

\

\

		} else {\

		}\

	    } else { \

		} else { \

		} \

	    } \

	} else { \

		} else { \

		} \

	    } else { \

		} else { \

		} \

	    } \

	} \

 \

    } \

}

%}

! !

!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;