"{ Encoding: utf8 }"

"

 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 DefaultPrintPrecision E Epsilon FMax FMin

		InvFact Ln10 Ln2 NaN Pi QDoubleOne QDoubleZero'

	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

/*

 * fpu_fix_start:

 *     Set the round-to-double flag, and save the old control word in old_cw.

 * fpu_fix_end:

 *     Restore the control word.

 */

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

#  if defined( __BORLANDC__ )

#   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, 0xFFFF);\

    }

#  else

#   define fpu_fix_start(old_cw_ptr) {  \

        *old_cw_ptr = _control87(0, 0); \

        _control87(0x00010000, 0x00030000);\

    }

#   define fpu_fix_end(old_cw_ptr) { \

        _control87(*old_cw_ptr, 0xFFFFFFFF);\

    }

#  endif

# 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 // defined(__x86__) || defined(__x86_64__)

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

    double* __d__;  \

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

    __stx_setClass(newQD, QDouble);                  \

    __d__ = __QDoubleInstPtr(newQD)->d_qDoubleValue; \

    __d__[0] = d0;   \

    __d__[1] = d1;   \

    __d__[2] = d2;   \

    __d__[3] = d3;   \

}

#if 0

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

#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 v_3 = s_3 - (a_3);\

    (err_3) = ((a_3) - (s_3 - v_3)) + ((b_3) - v_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;\

    double thi_7, tlo_7;\

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

        (a_7) *= 3.7252902984619140625e-09;\

        temp_7 = _QD_SPLITTER * (a_7);\

        thi_7 = temp_7 - (temp_7 - (a_7));\

        tlo_7 = (a_7) - thi_7;\

        thi_7 *= 268435456.0;\

        tlo_7 *= 268435456.0;\

    } else {\

        temp_7 = _QD_SPLITTER * (a_7);\

        thi_7 = temp_7 - (temp_7 - (a_7));\

        tlo_7 = (a_7) - thi_7;\

    }\

    (hi_7) = thi_7; \

    (lo_7) = tlo_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; \

    } \

}

#endif // 0

%}

! !

!QDouble primitiveFunctions!

%{

#ifdef __BORLANDC__

# define INLINE /* */

#else

# define INLINE inline

#endif

#define QD_IEEE_ADD

// routines

// fast_two_sum : s + e = a + b, s = fl(a + b), e = err(a + b), assumption |a|>|b|, Dekker('71)

// two_sum      : s + e = a + b, s = fl(a + b), e = err(a + b), Knuth('69)

// three_sum    : d + e + f = a + b + c, def are nonoverlapping, Bailey

// three_sum2   : d + e = a + b + c, de are nonoverlapping, Bailey

// two_prod     : s + e = a * b, s = fl(a * b), e = err(a * b), Verkamp and Dekker

// sqr          : s + e = a^2, s = fl(a * a), e = err(a * a)

// renorm       : renormalization algorithm

// qd_add_s  : qd + double

// qd_add_qd : qd + qd (sloppy add)

// s_sub_qd  : double - qd

// qd_sub_qd : qd - qd

// s_mul_qd  : double * qd

// qd_mul_qd : qd * qd

// qd_div_qd : qd / qd

// qd_sqr    : qd ^ 2

// qd_sqrt   : square root (scalar)

// qd_pow    : qd ^ n (n integer)

static INLINE void

fast_two_sum(double *s, double *e, double a, double b)

{

    double t;

    s[0] = t = a + b;

    e[0] = b - (t - a);

}

static INLINE void

two_sum(double *s, double *e, double a, double b)

{

    double t, bb;

    s[0] = t = a + b;

    bb = t - a;

    e[0] = (a - (t - bb)) + (b - bb);

}

static INLINE void

three_sum(double *d, double *e, double *f, double a, double b, double c)

{

    double t1,t2,t3,v;

    d[0] = 0.0; e[0] = 0.0; f[0] = 0.0;

    t1= a + b;

    v = t1 - a;

    t2= (a - (t1 - v))+(b - v);

    d[0] = t1 + c;

    v = d[0] - t1;

    t3= (t1 - (d[0] - v))+(c - v);

    e[0] = t2 + t3;

    v = e[0] - t2;

    f[0] = (t2 - (e[0] - v))+(t3 - v);

}

static INLINE void

three_sum2(double *d, double *e, double a, double b, double c) {

    double t1,t2,t3,v;

    d[0] = 0.0; e[0] = 0.0;

    t1= a + b;

    v = t1 - a;

    t2= (a - (t1 - v))+(b - v);

    d[0] = t1 + c;

    v = d[0] - t1;

    t3= (t1 - (d[0] - v))+(c - v);

    e[0] = t2 + t3;

}

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

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

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

static INLINE void

split(double a, double *hi, double *lo) {

    double temp;

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

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

        temp = _QD_SPLITTER * a;

        *hi = temp - (temp - a);

        *lo = a - *hi;

        *hi *= 268435456.0;          // 2^28

        *lo *= 268435456.0;          // 2^28

    } else {

        temp = _QD_SPLITTER * a;

        *hi = temp - (temp - a);

        *lo = a - *hi;

    }

}

#if 0

static INLINE void

two_prod(double *p, double *e, double a, double b)

{

    double t,ah,al,bh,bl;

    p[0] = a * b;

    t = 134217729 * a;       // splitter: 2^27 + 1

    ah = t -(t - a);

    al = a - ah;

    t = 134217729 * b;       // splitter: 2^27 + 1

    bh = t -(t - b);