--- a/QDouble.st Mon Jun 19 09:15:41 2017 +0200
+++ b/QDouble.st Mon Jun 19 09:57:38 2017 +0200
@@ -22,17 +22,577 @@
!
!QDouble primitiveDefinitions!
-ry:'Magnitude-Numbers'
+
+%{
+#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
+
+#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)) { \
+\
+ m_quick_two_sum(s0, c2, c3, c3);\
+ m_quick_two_sum(s0, c1, s0, c2);\
+ m_quick_two_sum(c0, c0, s0, c1);\
+\
+ s0 = c0;\
+ s1 = c1;\
+ if (s1 != 0.0) {\
+ m_quick_two_sum(s1, s1, c2, s2);\
+ if (s2 != 0.0) {\
+ m_quick_two_sum(s2, s2, c3, s3);\
+ } else {\
+ m_quick_two_sum(s1, s1, c3, s2);\
+ }\
+ } else {\
+ m_quick_two_sum(s0, s0, c2, s1);\
+ if (s1 != 0.0) {\
+ m_quick_two_sum(s1, s1, c3, s2);\
+ } else {\
+ m_quick_two_sum(s0, s0, c3, s1);\
+ }\
+ }\
+\
+ c0 = s0;\
+ c1 = s1;\
+ c2 = s2;\
+ c3 = s3;\
+ }\
+}
+
+#define m_renorm5(c0, c1, c2, c3, c4) { \
+ double s0, s1, s2 = 0.0, s3 = 0.0; \
+\
+ if (! isinf(c0)) { \
+ m_two_sum(s0, c3, c4, c4); \
+ m_quick_two_sum(s0, c2, s0, c3); \
+ m_quick_two_sum(s0, c1, s0, c2); \
+ m_quick_two_sum(c0, c0, s0, c1); \
+\
+ s0 = c0; \
+ s1 = c1; \
+\
+ m_two_sum(s0, c0, c1, s1); \
+ if (s1 != 0.0) { \
+ m_quick_two_sum(s1, s1, c2, s2); \
+ if (s2 != 0.0) { \
+ m_quick_two_sum(s2 ,s2, c3, s3); \
+ if (s3 != 0.0) {\
+ s3 += c4; \
+ } else {\
+ s2 += c4;\
+ }\
+ } else { \
+ m_quick_two_sum(s1, s1, c3, s2); \
+ if (s2 != 0.0) {\
+ m_quick_two_sum(s2, s2, c4, s3); \
+ } else { \
+ m_quick_two_sum(s1, s1, c4, s2); \
+ } \
+ } \
+ } else { \
+ m_quick_two_sum(s0,s0, c2, s1); \
+ if (s1 != 0.0) { \
+ m_quick_two_sum(s1,s1, c3, s2); \
+ if (s2 != 0.0) {\
+ m_quick_two_sum(s2,s2, c4, s3); \
+ } else { \
+ m_quick_two_sum(s1 ,s1, c4, s2); \
+ } \
+ } else { \
+ m_quick_two_sum(s0,s0, c3, s1); \
+ if (s1 != 0.0) { \
+ m_quick_two_sum(s1,s1, c4, s2); \
+ } else { \
+ m_quick_two_sum(s0,s0, c4, s1); \
+ } \
+ } \
+ } \
+ \
+ c0 = s0; \
+ c1 = s1; \
+ c2 = s2; \
+ c3 = s3; \
+ } \
+}
+
+%}
! !
!QDouble primitiveFunctions!
-= quick_two_sum(s1, *c3Ptr, &s2);
+%{
+
+#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);
} else {
s0 = quick_two_sum(s0, *c2Ptr, &s1);
- if (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<int>(a);
+}
+#endif
+
+%}
+! !
!QDouble class methodsFor:'documentation'!
+
copyright
"
COPYRIGHT (c) 2017 by eXept Software AG
@@ -316,22 +876,22 @@
invFact
InvFact isNil ifTrue:[
- InvFact := Array new:15.
- InvFact at:1 put:(self d0:1.66666666666666657e-01 d1:9.25185853854297066e-18 d2:5.13581318503262866e-34 d3:2.85094902409834186e-50).
- InvFact at:2 put:(self d0:4.16666666666666644e-02 d1:2.31296463463574266e-18 d2:1.28395329625815716e-34 d3:7.12737256024585466e-51).
- InvFact at:3 put:(self d0:8.33333333333333322e-03 d1:1.15648231731787138e-19 d2:1.60494162032269652e-36 d3:2.22730392507682967e-53).
- InvFact at:4 put:(self d0:1.38888888888888894e-03 d1:-5.30054395437357706e-20 d2:-1.73868675534958776e-36 d3:-1.63335621172300840e-52).
- InvFact at:5 put:(self d0:1.98412698412698413e-04 d1:1.72095582934207053e-22 d2:1.49269123913941271e-40 d3:1.29470326746002471e-58).
- InvFact at:6 put:(self d0:2.48015873015873016e-05 d1:2.15119478667758816e-23 d2:1.86586404892426588e-41 d3:1.61837908432503088e-59).
- InvFact at:7 put:(self d0:2.75573192239858925e-06 d1:-1.85839327404647208e-22 d2:8.49175460488199287e-39 d3:-5.72661640789429621e-55).
- InvFact at:8 put:(self d0:2.75573192239858883e-07 d1:2.37677146222502973e-23 d2:-3.26318890334088294e-40 d3:1.61435111860404415e-56).
- InvFact at:9 put:(self d0:2.50521083854417202e-08 d1:-1.44881407093591197e-24 d2:2.04267351467144546e-41 d3:-8.49632672007163175e-58).
- InvFact at:10 put:(self d0:2.08767569878681002e-09 d1:-1.20734505911325997e-25 d2:1.70222792889287100e-42 d3:1.41609532150396700e-58).
- InvFact at:11 put:(self d0:1.60590438368216133e-10 d1:1.25852945887520981e-26 d2:-5.31334602762985031e-43 d3:3.54021472597605528e-59).
- InvFact at:12 put:(self d0:1.14707455977297245e-11 d1:2.06555127528307454e-28 d2:6.88907923246664603e-45 d3:5.72920002655109095e-61).
- InvFact at:13 put:(self d0:7.64716373181981641e-13 d1:7.03872877733453001e-30 d2:-7.82753927716258345e-48 d3:1.92138649443790242e-64).
- InvFact at:14 put:(self d0:4.77947733238738525e-14 d1:4.39920548583408126e-31 d2:-4.89221204822661465e-49 d3:1.20086655902368901e-65).
- InvFact at:15 put:(self d0:2.81145725434552060e-15 d1:1.65088427308614326e-31 d2:-2.87777179307447918e-50 d3:4.27110689256293549e-67).
+ InvFact := Array new:15.
+ InvFact at:1 put:(self d0:1.66666666666666657e-01 d1:9.25185853854297066e-18 d2:5.13581318503262866e-34 d3:2.85094902409834186e-50).
+ InvFact at:2 put:(self d0:4.16666666666666644e-02 d1:2.31296463463574266e-18 d2:1.28395329625815716e-34 d3:7.12737256024585466e-51).
+ InvFact at:3 put:(self d0:8.33333333333333322e-03 d1:1.15648231731787138e-19 d2:1.60494162032269652e-36 d3:2.22730392507682967e-53).
+ InvFact at:4 put:(self d0:1.38888888888888894e-03 d1:-5.30054395437357706e-20 d2:-1.73868675534958776e-36 d3:-1.63335621172300840e-52).
+ InvFact at:5 put:(self d0:1.98412698412698413e-04 d1:1.72095582934207053e-22 d2:1.49269123913941271e-40 d3:1.29470326746002471e-58).
+ InvFact at:6 put:(self d0:2.48015873015873016e-05 d1:2.15119478667758816e-23 d2:1.86586404892426588e-41 d3:1.61837908432503088e-59).
+ InvFact at:7 put:(self d0:2.75573192239858925e-06 d1:-1.85839327404647208e-22 d2:8.49175460488199287e-39 d3:-5.72661640789429621e-55).
+ InvFact at:8 put:(self d0:2.75573192239858883e-07 d1:2.37677146222502973e-23 d2:-3.26318890334088294e-40 d3:1.61435111860404415e-56).
+ InvFact at:9 put:(self d0:2.50521083854417202e-08 d1:-1.44881407093591197e-24 d2:2.04267351467144546e-41 d3:-8.49632672007163175e-58).
+ InvFact at:10 put:(self d0:2.08767569878681002e-09 d1:-1.20734505911325997e-25 d2:1.70222792889287100e-42 d3:1.41609532150396700e-58).
+ InvFact at:11 put:(self d0:1.60590438368216133e-10 d1:1.25852945887520981e-26 d2:-5.31334602762985031e-43 d3:3.54021472597605528e-59).
+ InvFact at:12 put:(self d0:1.14707455977297245e-11 d1:2.06555127528307454e-28 d2:6.88907923246664603e-45 d3:5.72920002655109095e-61).
+ InvFact at:13 put:(self d0:7.64716373181981641e-13 d1:7.03872877733453001e-30 d2:-7.82753927716258345e-48 d3:1.92138649443790242e-64).
+ InvFact at:14 put:(self d0:4.77947733238738525e-14 d1:4.39920548583408126e-31 d2:-4.89221204822661465e-49 d3:1.20086655902368901e-65).
+ InvFact at:15 put:(self d0:2.81145725434552060e-15 d1:1.65088427308614326e-31 d2:-2.87777179307447918e-50 d3:4.27110689256293549e-67).
].
^ InvFact
@@ -1237,27 +1797,27 @@
"/ argument substantially speeds up the convergence. */
|k inv_k d0 m mul_pwr2 r
- s p t thresh _eps i|
-
- _eps := 1.21543267145725e-63. "/ = 2^-209.
-
+ s p t thresh eps i|
+
+ eps := 1.21543267145725e-63. "/ = 2^-209.
+
k := 1.0 ldexp:16.
inv_k := 1.0 / k.
d0 := self d0.
(d0 <= -709.0) ifTrue:[
- ^ 0.0 asQDouble.
- ].
- (d0 >= 709.0) ifTrue:[
- ^ Infinity positive
- ].
+ ^ 0.0 asQDouble.
+ ].
+ (d0 >= 709.0) ifTrue:[
+ ^ Infinity positive
+ ].
self isZero ifTrue:[
- ^ 1.0 asQDouble.
- ].
+ ^ 1.0 asQDouble.
+ ].
self isOne ifTrue:[
- ^ self class e
+ ^ self class e
].
"/ double m = std::floor(a.x[0] / qd_real::_log2.x[0] + 0.5);
@@ -1267,16 +1827,16 @@
"/ qd_real r = mul_pwr2(a - qd_real::_log2 * m, inv_k);
r := mul_pwr2 value:(self - (self class ln2 * m)) value:inv_k.
- thresh := inv_k * _eps.
+ thresh := inv_k * eps.
p := r squared.
s := r + (mul_pwr2 value:p value:0.5).
i := 1.
[
- p := p * r.
- t := p * (self class invFact at:i).
- i := i+1.
- s := s + t.
+ p := p * r.
+ t := p * (self class invFact at:i).
+ i := i+1.
+ s := s + t.
] doWhile:[ (t asFloat abs > thresh) and:[i < 10] ]. "/ (std::abs(to_double(t)) > thresh && i < 9);
s := (mul_pwr2 value:s value:2.0) + s squared.
@@ -1359,18 +1919,18 @@
d0 := self d0.
d0 = 1.0 ifTrue:[
- self isOne ifTrue:[ ^ self class zero ].
+ self isOne ifTrue:[ ^ self class zero ].
].
d0 < 0.0 ifTrue:[
- ^ self class
- raise:#domainErrorSignal
- receiver:self
- selector:#ln
- arguments:#()
- errorString:'bad receiver in ln'
+ ^ self class
+ raise:#domainErrorSignal
+ receiver:self
+ selector:#ln
+ arguments:#()
+ errorString:'bad receiver in ln'
].
d0 = 0.0 ifTrue:[
- ^ Infinity negative
+ ^ Infinity negative
].
"/ initial approx.
@@ -1997,4 +2557,3 @@
version_CVS
^ '$Header$'
! !
-