# HG changeset patch # User Claus Gittinger # Date 1497824208 -7200 # Node ID 8055a8f0b66ff3a68ba974e0489423b7af9204ce # Parent 18aff5006de03b6cc98ad622e68c774b2e8dbab8 #FEATURE by cg class: QDouble added: #isOne #isZero #ln class: QDouble class added: #infinity #negativeInfinity diff -r 18aff5006de0 -r 8055a8f0b66f QDouble.st --- a/QDouble.st Sun Jun 18 23:20:39 2017 +0200 +++ b/QDouble.st Mon Jun 19 00:16:48 2017 +0200 @@ -1,5 +1,3 @@ -"{ Encoding: utf8 }" - " COPYRIGHT (c) 2017 by eXept Software AG All Rights Reserved @@ -24,590 +22,14 @@ ! !QDouble primitiveDefinitions! -%{ -#include -#include - -#define __USE_ISOC9X 1 -#define __USE_ISOC99 1 -#include - -/* - * on some systems, errno is a macro ... check for it here - */ -#ifndef errno - extern errno; -#endif - -#if !defined (__win32__) -# include -#endif - -#if defined (__aix__) -# include -#endif - -#if defined(__irix__) || defined(__solaris__) || defined(__sunos__) -# include -#endif - -#if defined(__linux__) -# ifndef NAN -# include -# 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; \ - } \ -} - -%} +ry:'Magnitude-Numbers' ! ! !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; - 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); += 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; -} - -inline void -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) - s3 += *c4Ptr; - else - s2 += *c4Ptr; - } else { - s1 = quick_two_sum(s1, *c3Ptr, &s2); - if (s2 != 0.0) - s2 = quick_two_sum(s2, *c4Ptr, &s3); - else - s1 = quick_two_sum(s1, *c4Ptr, &s2); - } - } else { - s0 = quick_two_sum(s0, *c2Ptr, &s1); - if (s1 != 0.0) { - s1 = quick_two_sum(s1, *c3Ptr, &s2); - if (s2 != 0.0) - s2 = quick_two_sum(s2, *c4Ptr, &s3); - else - s1 = quick_two_sum(s1, *c4Ptr, &s2); - } else { - s0 = quick_two_sum(s0, *c3Ptr, &s1); - if (s1 != 0.0) - s1 = quick_two_sum(s1, *c4Ptr, &s2); - else - s0 = quick_two_sum(s0, *c4Ptr, &s1); - } - } - - *c0Ptr = s0; - *c1Ptr = s1; - *c2Ptr = s2; - *c3Ptr = s3; -} - -inline void -three_sum(double *aPtr, double *bPtr, double *cPtr) { - double t1, t2, t3; - t1 = two_sum(*aPtr, *bPtr, &t2); - *aPtr = two_sum(*cPtr, t1, &t3); - *bPtr = two_sum(t2, t3, cPtr); -} - -inline void three_sum2(double *aPtr, double *bPtr, double *cPtr) { - double t1, t2, t3; - t1 = two_sum(*aPtr, *bPtr, &t2); - *aPtr = two_sum(*cPtr, t1, &t3); - *bPtr = t2 + t3; -} -#endif - -#if 0 -/* These are provided to give consistent - interface for double with double-double and quad-double. */ -inline void -sincosh(double t, double &sinh_t, double &cosh_t) { - sinh_t = sinh(t); - cosh_t = cosh(t); -} - -inline double -sqr(double t) { - return t * t; -} - -inline double -to_double(double a) { - return a; -} - -inline int -to_int(double a) { - return static_cast(a); -} -#endif - -%} -! ! + if (s1 ! ! !QDouble class methodsFor:'documentation'! @@ -886,6 +308,12 @@ "Created: / 14-06-2017 / 19:14:49 / cg" ! +infinity + ^ Infinity positive + + "Created: / 18-06-2017 / 23:41:07 / cg" +! + ln10 "return the constant e as quad precision double. (returns approx. 212 bits of precision)" @@ -926,6 +354,12 @@ "Created: / 12-06-2017 / 18:31:34 / cg" ! +negativeInfinity + ^ Infinity negative + + "Created: / 18-06-2017 / 23:40:47 / cg" +! + pi "return the constant pi as quad precision double. (returns approx. 212 bits of precision)" @@ -1816,6 +1250,16 @@ "Modified (comment): / 13-06-2017 / 17:33:19 / cg" ! +ln + self isOne ifTrue:[ ^ self class zero ]. + + " + 1.0 asQDouble ln + " + + "Created: / 18-06-2017 / 23:32:54 / cg" +! + negated ^ self class d0:(self d0) negated @@ -2397,6 +1841,21 @@ ^ self d0 isNaN "Created: / 15-06-2017 / 01:57:35 / cg" +! + +isOne + ^ self d0 = 1.0 + and:[ self d1 = 0.0 + and:[ self d2 = 0.0 + and:[ self d3 = 0.0 ]]] + + "Created: / 18-06-2017 / 23:29:07 / cg" +! + +isZero + ^ self d0 = 0.0 + + "Created: / 18-06-2017 / 23:29:25 / cg" ! ! !QDouble class methodsFor:'documentation'! @@ -2408,3 +1867,4 @@ version_CVS ^ '$Header$' ! ! +