#DOCUMENTATION by cg
class: SharedQueue
comment/format in:
#next
#nextWithTimeout:
"{ 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 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
#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 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; \
} \
}
%}
! !
!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);
} 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<int>(a);
}
#endif
%}
! !
!QDouble class methodsFor:'documentation'!
copyright
"
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.
"
!
documentation
"
ATTENTION: ongoing, unfinished work.
No warranty that this works correctly...
QDoubles represent rational numbers with extended, but still limited precision.
In contrast to Floats (which use the C-compiler's native 64bit 'double' format),
QDoubles give you roughly 200 bit or approx. 60 decimal digits of precision.
Representation:
QDoubles use 4 IEEE doubles, each keeping 53 bits of precision.
A qDouble's value is the sum of those 4 doubles,
and a qDouble keeps this unevaluated sum as its state.
(due to overlap, the final precision is less than 53*4)
Range and Precision of Storage Formats: see LimitedPrecisionReal >> documentation
the number of decmal digits:
QDouble decimalPrecision
LongFloat decimalPrecision
Float decimalPrecision
ShortFloat decimalPrecision
the number of bits:
QDouble precision
LongFloat precision
Float precision
ShortFloat precision
Notice:
when assigning a converted double precision number as in:
qd := 1.0 asQDouble.
you still get only a regular double precision approximation to 0.1
because the error is already inherit in the double.
For a full precision constant, you (currently) need to convert from a string
(because the compilers do not know about then, yet):
qd := QDouble readFrom:'0.1'.
To see the error of the double precision version, compute:
(0.1 asQDouble) - (QDouble readFrom:'0.1')
[author:]
Claus Gittinger
[see also:]
Number
Float ShortFloat LongFloat
Fraction FixedPoint Integer Complex
FloatArray DoubleArray
"
!
examples
"
Floats, LongFloats suffer from loosing bits:
(Float readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
-(Float readFrom:'0.333333333333333333333333333333333333333333333333333333333')
-> 0.0
(Float readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
= (Float readFrom:'0.333333333333333333333333333333333333333333333333333333333')
-> true
(Float readFrom:'0.33333333333333333333333333333333333333333333333333333333333333333333')
= (Float readFrom:'0.3333333333333333333333333333333333333333333333333333333333333333333')
-> true
1000 0110 1000 0101 1000 0101 1000 0101 1000 0101 1000 0101 1101 0101 0011 1111
(Float readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
= (Float readFrom:'0.3333333333333333333333333333333333333333333333333333333333333333333')
(LongFloat readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
-(LongFloat readFrom:'0.333333333333333333333333333333333333333333333333333333333')
-> 0.0
(LongFloat readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
= (LongFloat readFrom:'0.333333333333333333333333333333333333333333333333333333333')
-> 0.0
(QDouble readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
-(QDouble readFrom:'0.333333333333333333333333333333333333333333333333333333333')
(QDouble readFrom:'0.33333333333333333333333333333333333333333333333333333333333')
-(QDouble readFrom:'0.3333333333333333333333333333333333333333333333333333333333')
(QDouble readFrom:'0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333')
-(QDouble readFrom:'0.3333333333333333333333333333333333333333333333333333333333333333333333333333333333333333')
"
! !
!QDouble class methodsFor:'instance creation'!
basicNew
"return a new quad-precision double - here we return 0.0
Notice that numbers are usually NOT created this way ...
It's implemented here to allow things like binary store & load
of floats. (but even this support will go away eventually, it's not
a good idea to store the bits of a float - the reader might have a
totally different representation - so floats should be
binary stored in a device independent format."
%{ /* NOCONTEXT */
#ifdef __SCHTEAM__
ERROR("trying to instantiate a quad double");
#else
OBJ newFloat;
__qNew(newFloat, sizeof(struct __quadDoubleStruct));
__stx_setClass(newFloat, QDouble);
__QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[0] = 0.0;
__QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[1] = 0.0;
__QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[2] = 0.0;
__QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[3] = 0.0;
RETURN (newFloat);
#endif /* not SCHTEAM */
%}
"
self basicNew
"
"Created: / 12-06-2017 / 16:00:38 / cg"
!
d0:d0 d1:d1 d2:d2 d3:d3
"return a new quad-precision double from individual double components"
%{ /* NOCONTEXT */
#ifdef __SCHTEAM__
ERROR("trying to instantiate a quad double");
#else
OBJ newQD;
if (__isFloatLike(d0)
&& __isFloatLike(d1)
&& __isFloatLike(d2)
&& __isFloatLike(d3)) {
__qNew_qdReal(newQD, __floatVal(d0), __floatVal(d1),
__floatVal(d2), __floatVal(d3));
RETURN (newQD);
}
#endif
%}.
self error:'invalid argument'
"
self d0: 3.141592653589793116e+00
d1: 1.224646799147353207e-16
d2: -2.994769809718339666e-33
d3: 1.112454220863365282e-49
"
"Created: / 12-06-2017 / 20:17:14 / cg"
!
fromDoubleArray:aDoubleArray
"return a new quad-precision double from coercing a double array"
%{ /* NOCONTEXT */
#ifdef __SCHTEAM__
ERROR("trying to instantiate a quad double");
#else
OBJ newQD;
if (__isDoubleArray(aDoubleArray)) {
__qNew_qdReal(newQD,
__DoubleArrayInstPtr(aDoubleArray)->d_element[0],
__DoubleArrayInstPtr(aDoubleArray)->d_element[1],
__DoubleArrayInstPtr(aDoubleArray)->d_element[2],
__DoubleArrayInstPtr(aDoubleArray)->d_element[3]);
RETURN (newQD);
}
#endif
%}.
self error:'invalid argument'
"
self fromDoubleArray(DoubleArray
with: 3.141592653589793116e+00
with: 1.224646799147353207e-16
with: -2.994769809718339666e-33
with: 1.112454220863365282e-49)
"
"Created: / 12-06-2017 / 18:25:32 / cg"
!
fromFloat:aFloat
"return a new quad-precision double from coercing aFloat"
%{ /* NOCONTEXT */
#ifdef __SCHTEAM__
ERROR("trying to instantiate a quad double");
#else
double dVal;
OBJ newFloat;
if (__isFloatLike(aFloat)) {
dVal = __floatVal(aFloat);
} else if (__isShortFloat(aFloat)) {
dVal = __shortFloatVal(aFloat);
} else {
goto badArg;
}
__qNew(newFloat, sizeof(struct __quadDoubleStruct));
__stx_setClass(newFloat, QDouble);
__QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[0] = dVal;
__QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[1] = 0.0;
__QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[2] = 0.0;
__QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[3] = 0.0;
RETURN (newFloat);
badArg: ;
#endif
%}.
self error:'invalid argument'
"
self fromFloat:1.0
"
"Created: / 12-06-2017 / 16:06:54 / cg"
!
fromInteger:anInteger
"return a new quad-precision double from coercing anInteger"
%{ /* NOCONTEXT */
#ifndef __SCHTEAM__
OBJ newFloat;
if (__isSmallInteger(anInteger)) {
INT iVal = __smallIntegerVal(anInteger);
double *a;
__qNew(newFloat, sizeof(struct __quadDoubleStruct));
__stx_setClass(newFloat, QDouble);
a = __QuadDoubleInstPtr(newFloat)->d_quadDoubleValue;
a[1] = 0.0;
a[2] = 0.0;
a[3] = 0.0;
// need more than 52bits?
if ((sizeof(INT) > 52)
&& ((iVal > 0xFFFFFFFFFFFFF)
|| (iVal < -0xFFFFFFFFFFFFF))) {
a[0] = (double)(iVal & ~0xFFFFFFFF);
a[1] = (double)(iVal & 0xFFFFFFFF);
// m_renorm4(a[0], a[1], a[2], a[3]);
} else {
a[0] = (double)iVal;
}
RETURN (newFloat);
}
#endif
%}.
^ super fromInteger:anInteger
"
self fromInteger:2
self fromInteger:16rFFFFFFFF -- 32bit 4294967295.0
self fromInteger:16rFFFFFFFFFFFF -- 48bit 281474976710655.0
self fromInteger:16rFFFFFFFFFFFFF -- 52bit 4503599627370495.0
self fromInteger:16rFFFFFFFFFFFFFF -- 56bit 72057594037927935.0
self fromInteger:16rFFFFFFFFFFFFFFF -- 60bit 1152921504606846975.0
self fromInteger:16r1FFFFFFFFFFFFFFF -- 61bit 2305843009213693951.0
self fromInteger:16r3FFFFFFFFFFFFFFF -- 62bit 4611686018427387903.0
self fromInteger:16r7FFFFFFFFFFFFFFF -- 63bit 9223372036854775807.0
self fromInteger:16rFFFFFFFFFFFFFFFF -- 64bit 18446744073709551615.0
"
"Created: / 12-06-2017 / 16:10:10 / cg"
"Modified: / 04-07-2017 / 12:51:52 / cg"
! !
!QDouble class methodsFor:'coercing & converting'!
coerce:aNumber
"convert the argument aNumber into an instance of the receiver's class and return it."
^ aNumber asQDouble
"Created: / 12-06-2017 / 17:13:47 / cg"
"Modified: / 12-06-2017 / 21:09:06 / cg"
! !
!QDouble class methodsFor:'constants'!
NaN
"return a QDouble which represents not-a-Number (i.e. an invalid number)"
NaN isNil ifTrue:[
NaN := self d0:(Float NaN) d1:(Float NaN) d2:(Float NaN) d3:(Float NaN)
].
^ NaN
"Created: / 21-06-2017 / 20:44:57 / cg"
!
e
"return the constant e as quad precision double.
(returns approx. 200 bits of precision)"
E isNil ifTrue:[
E := self fromDoubleArray:(DoubleArray
with: 2.718281828459045091e+00
with: 1.445646891729250158e-16
with: -2.127717108038176765e-33
with: 1.515630159841218954e-49)
].
^ E
"
self e
"
"Created: / 12-06-2017 / 18:29:36 / cg"
!
fmax
"return the constant e as quad precision double.
(returns approx. 200 bits of precision)"
FMax isNil ifTrue:[
FMax := self fromDoubleArray:(DoubleArray
with: 1.797693134862314E+308
with: 9.97920154767359795037e+291
with: 5.53956966280111259858e+275
with: 3.07507889307840487279e+259)
].
^ FMax
"
Float fmax
self fmax
"
"Created: / 14-06-2017 / 19:14:18 / cg"
!
fmin
"return the smallest representable instance of this class"
FMin isNil ifTrue:[
FMin := 1.6259745436952323e-260 asQDouble
].
^ FMin
"
QDouble fmin
Float fmin
"
"Created: / 14-06-2017 / 19:14:49 / cg"
!
infinity
^ Infinity positive
"Created: / 18-06-2017 / 23:41:07 / cg"
!
invFact
"table returning 1/n!!
(for taylor series)"
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
"
1.0 / (3 factorial) 0.166666666666667
1 asQDouble / (3 factorial) - (self invFact at:1)
1 asQDouble / (4 factorial) - (self invFact at:2)
1 asQDouble / (5 factorial) - (self invFact at:3)
1 asQDouble / (6 factorial) - (self invFact at:4)
1 asQDouble / (7 factorial) - (self invFact at:5)
1 asQDouble / (8 factorial) - (self invFact at:6)
1 asQDouble / (9 factorial) - (self invFact at:7)
1 asQDouble / (10 factorial) - (self invFact at:8)
1 asQDouble / (11 factorial) - (self invFact at:9)
1 asQDouble / (12 factorial) - (self invFact at:10)
1 asQDouble / (13 factorial) - (self invFact at:11)
1 asQDouble / (14 factorial) - (self invFact at:12)
1 asQDouble / (15 factorial) - (self invFact at:13)
1 asQDouble / (16 factorial) - (self invFact at:14)
1 asQDouble / (17 factorial) - (self invFact at:15)
"
"Created: / 19-06-2017 / 02:22:23 / cg"
"Modified (comment): / 20-06-2017 / 13:04:54 / cg"
!
ln10
"return the constant e as quad precision double.
(returns approx. 200 bits of precision)"
Ln10 isNil ifTrue:[
Ln10 := self fromDoubleArray:(DoubleArray
with: 2.302585092994045901e+00
with: -2.170756223382249351e-16
with: -9.984262454465776570e-33
with: -4.023357454450206379e-49)
].
^ Ln10
"
self ln10
"
"Created: / 12-06-2017 / 18:32:29 / cg"
!
ln2
"return the constant e as quad precision double.
(returns approx. 200 bits of precision)"
Ln2 isNil ifTrue:[
Ln2 := self fromDoubleArray:(DoubleArray
with: 6.931471805599452862e-01
with: 2.319046813846299558e-17
with: 5.707708438416212066e-34
with: -3.582432210601811423e-50)
].
^ Ln2
"
self ln2
"
"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. 200 bits of precision)"
Pi isNil ifTrue:[
Pi := self fromDoubleArray:(DoubleArray
with: 3.141592653589793116e+00
with: 1.224646799147353207e-16
with: -2.994769809718339666e-33
with: 1.112454220863365282e-49)
].
^ Pi
"
self pi
"
"Created: / 12-06-2017 / 18:27:13 / cg"
!
unity
"return the neutral element for multiplication (1.0) as QDouble"
QDoubleOne isNil ifTrue:[
QDoubleOne := 1.0 asQDouble.
].
^ QDoubleOne
"
self unity
"
"Created: / 15-06-2017 / 11:45:22 / cg"
!
zero
"return the neutral element for addition (0.0) as QDouble"
QDoubleZero isNil ifTrue:[
QDoubleZero := 0.0 asQDouble
].
^ QDoubleZero
"
self zero
"
"Created: / 15-06-2017 / 11:44:13 / cg"
! !
!QDouble class methodsFor:'queries'!
defaultPrintPrecision
"return the number of decimal digits printed by default"
^ 20
"
ShortFloat defaultPrintPrecision
Float defaultPrintPrecision
LongFloat defaultPrintPrecision
QDouble defaultPrintPrecision
"
"Created: / 17-06-2017 / 02:58:51 / cg"
"Modified: / 21-06-2017 / 13:39:08 / cg"
!
epsilon
"return the maximum relative spacing of instances of mySelf
(i.e. the value-delta of the least significant bit)"
^ 1.2154326714572500565324311366323150942261000827598106963711353e-64
"
Float epsilon -> 2.22044604925031E-16
ShortFloat epsilon -> 1.19209289550781E-07
LongFloat epsilon -> 1.0842021724855E-19
QDouble epsilon -> 1.21543267145725E-63
"
"Created: / 12-06-2017 / 18:52:44 / cg"
"Modified: / 22-06-2017 / 15:34:56 / cg"
!
numBitsInExponent
"answer the number of bits in the exponent"
^ Float numBitsInExponent
"
1.0 asQDouble numBitsInExponent
"
"Created: / 12-06-2017 / 11:11:04 / cg"
"Modified (comment): / 28-05-2019 / 08:55:04 / Claus Gittinger"
!
numBitsInMantissa
"answer the number of bits in the mantissa.
Here, a fake number is returned"
^ (Float numBitsInMantissa - 1) * 4
"
1.0 asFloat numBitsInMantissa
1.0 asShortFloat numBitsInMantissa
1.0 asLongFloat numBitsInMantissa
1.0 asQDouble numBitsInMantissa
1.0 asQDouble class numBitsInMantissa
Float numBitsInMantissa
ShortFloat numBitsInMantissa
QDouble numBitsInMantissa
"
"Created: / 12-06-2017 / 11:13:44 / cg"
"Modified (comment): / 20-06-2017 / 11:05:26 / cg"
"Modified (comment): / 28-05-2019 / 09:07:07 / Claus Gittinger"
!
precision
"answer the number of bits in the mantissa"
"/ subtract some due to overlap in the component numbers
^ (Float precision - 2) * 4
"
1.0 class numBitsInMantissa
1.0 asShortFloat class numBitsInMantissa
1.0 asLongFloat class numBitsInMantissa
1.0 asQDouble class numBitsInMantissa
Float numBitsInMantissa
ShortFloat numBitsInMantissa
QDouble numBitsInMantissa
QDouble precision
"
"Created: / 12-06-2017 / 18:49:11 / cg"
"Modified (comment): / 20-06-2017 / 12:59:00 / cg"
!
radix
"answer the radix of a Floats exponent
This is an IEEE float, which is represented as binary"
^ Float radix
"Created: / 12-06-2017 / 18:50:04 / cg"
! !
!QDouble methodsFor:'arithmetic'!
* aNumber
"return the product of the receiver and the argument, aNumber"
^ aNumber productFromQDouble:self
"
(((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) * (QDouble fromFloat:2.0)) asDoubleArray
"
"Created: / 13-06-2017 / 01:00:47 / cg"
"Modified (comment): / 14-06-2017 / 12:08:50 / cg"
!
+ aNumber
"return the sum of the receiver and the argument, aNumber"
^ aNumber sumFromQDouble:self
"
((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) asDoubleArray
"
"Created: / 12-06-2017 / 16:17:46 / cg"
"Modified: / 12-06-2017 / 23:06:22 / cg"
!
- aNumber
"return the sum of the receiver and the argument, aNumber"
^ self + (aNumber negated)
"
(QDouble fromFloat:1e20) asDoubleArray
((QDouble fromFloat:1e20) - (QDouble fromFloat:1.0)) asDoubleArray
(QDouble fromFloat:1e-20) asDoubleArray
((QDouble fromFloat:1e-20) - (QDouble fromFloat:1.0)) asDoubleArray
((QDouble fromFloat:2.0) - (QDouble fromFloat:1.0)) asDoubleArray
((QDouble fromFloat:2.0) - (QDouble fromFloat:1.0) + (QDouble fromFloat:1.0)) asDoubleArray
((QDouble fromFloat:1e-20) - (QDouble fromFloat:1.0) + (QDouble fromFloat:1.0)) asDoubleArray
"
"Created: / 12-06-2017 / 23:41:39 / cg"
"Modified (comment): / 15-06-2017 / 00:34:41 / cg"
!
/ aNumber
"return the quotient of the receiver and the argument, aNumber"
^ aNumber quotientFromQDouble:self
"
((QDouble fromFloat:1e20) / (QDouble fromFloat:2.0)) asDoubleArray
((QDouble fromFloat:1.2345) / (QDouble fromFloat:10.0)) asDoubleArray
((QDouble fromFloat:1.2345) / 10.0) asDoubleArray
"
"Created: / 13-06-2017 / 17:59:09 / cg"
"Modified (comment): / 15-06-2017 / 00:14:26 / cg"
! !
!QDouble methodsFor:'coercing & converting'!
asDoubleArray
^ DoubleArray
with:self d0
with:self d1
with:self d2
with:self d3.
"
(QDouble fromFloat:1.0) asDoubleArray
(QDouble fromFloat:2.0) asDoubleArray
"
"Created: / 12-06-2017 / 18:19:19 / cg"
"Modified (comment): / 13-06-2017 / 17:58:09 / cg"
!
asFloat
^ self d0
"
(QDouble fromFloat:1.0) asFloat
(QDouble fromFloat:2.0) asFloat
"
"Created: / 12-06-2017 / 18:15:27 / cg"
"Modified: / 13-06-2017 / 17:56:50 / cg"
!
asInteger
^ self d0 asInteger
+ self d1 asInteger
+ self d2 asInteger
+ self d3 asInteger
"Created: / 19-06-2017 / 18:07:17 / cg"
!
asQDouble
"return a QDouble with same value as myself."
^ self
"Created: / 15-06-2017 / 12:08:02 / cg"
!
asTrueFraction
self halt.
^ self d0 asTrueFraction
+ self d1 asTrueFraction
+ self d2 asTrueFraction
+ self d3 asTrueFraction
"
1e10 asTrueFraction -> 10000000000
1e20 asTrueFraction -> 100000000000000000000
(1e20 + 1) asTrueFraction -> 100000000000000000000 ouch!!
1e10 asQDouble asTrueFraction -> 10000000000
1e20 asQDouble asTrueFraction -> 100000000000000000000
(1e20 asQDouble + 1) asTrueFraction -> 100000000000000000001
(1e40 asQDouble + 1e20 + 1) asTrueFraction -> 10000000000000000303886028427003666890753
(1e40 asQDouble + 1e20) asTrueFraction
"
"Created: / 20-06-2017 / 11:09:03 / cg"
!
coerce:aNumber
"convert the argument aNumber into an instance of the receiver's class and return it."
^ aNumber asQDouble
"Created: / 12-06-2017 / 17:13:47 / cg"
"Modified: / 12-06-2017 / 21:09:06 / cg"
!
exponent
^ self d0 exponent
"Created: / 20-06-2017 / 11:06:02 / cg"
!
generality
"return the generality value - see ArithmeticValue>>retry:coercing:"
^ 95
"Created: / 12-06-2017 / 17:13:14 / cg"
!
negative
^ self d0 negative
"
(QDouble fromFloat:0.0) negative
(QDouble fromFloat:1.0) negative
(QDouble fromFloat:-1.0) negative
"
"Created: / 13-06-2017 / 01:57:39 / cg"
"Modified: / 13-06-2017 / 17:58:26 / cg"
!
positive
"return true, if the receiver is greater or equal to zero (not negative)"
^ self d0 positive
"
(QDouble fromFloat:1.0) positive
(QDouble fromFloat:-1.0) positive
"
"Created: / 13-06-2017 / 01:56:53 / cg"
"Modified: / 13-06-2017 / 17:58:41 / cg"
"Modified (comment): / 28-05-2019 / 05:55:55 / Claus Gittinger"
! !
!QDouble methodsFor:'comparing'!
< aNumber
"return true, if the argument, aNumber is greater than the receiver"
^ aNumber lessFromQDouble:self
"Created: / 13-06-2017 / 16:58:53 / cg"
!
= aNumber
"return true, if the argument, aNumber has the same value as than the receiver"
^ aNumber equalFromQDouble:self
"Created: / 13-06-2017 / 17:12:09 / cg"
! !
!QDouble methodsFor:'double dispatching'!
differenceFromFloat:aFloat
"aFloat - self"
^ aFloat + (self negated)
"
1.0 - (QDouble fromFloat:1.0)
1e20 - (QDouble fromFloat:1.0)
(1.0 - (QDouble fromFloat:1.0)) asFloat
(1e20 - (QDouble fromFloat:1.0)) asFloat
(1.0 - (QDouble fromFloat:1.0)) asDoubleArray
(1e20 - (QDouble fromFloat:1.0)) asDoubleArray
(1e20 - (QDouble fromFloat:1.0) + 1e-20) asDoubleArray
"
"Created: / 12-06-2017 / 23:38:05 / cg"
!
differenceFromQDouble:aQDouble
"aQDouble - self"
^ aQDouble + (self negated)
"
1.0 - (QDouble fromFloat:1.0)
1e20 - (QDouble fromFloat:1.0)
(1.0 - (QDouble fromFloat:1.0)) asFloat
(1e20 - (QDouble fromFloat:1.0)) asFloat
(1.0 - (QDouble fromFloat:1.0)) asDoubleArray
(1e20 - (QDouble fromFloat:1.0)) asDoubleArray
(1e20 - (QDouble fromFloat:1.0) + 1e-20) asDoubleArray
"
"Created: / 12-06-2017 / 23:38:19 / cg"
!
equalFromQDouble:aQDouble
%{
if (__Class(aQDouble) == QDouble) {
double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
double *b = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
RETURN ((a[0] == b[0]
&& a[1] == b[1]
&& a[2] == b[2]
&& a[3] == b[3]) ? true : false);
}
%}.
^ (aQDouble d0 = self d0)
and:[ (aQDouble d1 = self d1)
and:[ (aQDouble d2 = self d2)
and:[ (aQDouble d3 = self d3) ]]]
"
(QDouble fromFloat:1.0) = (QDouble fromFloat:1.0)
(QDouble fromFloat:1.0) = 1.0
1.0 = (QDouble fromFloat:1.0)
"
"Created: / 13-06-2017 / 03:01:19 / cg"
"Modified: / 13-06-2017 / 18:01:52 / cg"
!
lessFromQDouble:aQDouble
"sent when aQDouble does not know how to compare to the receiver..
Return true if aQDouble < self"
%{
if (__Class(aQDouble) == QDouble) {
double *a = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
double *b = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
// now compare if a < b!
RETURN
((a[0] < b[0] ||
(a[0] == b[0] && (a[1] < b[1] ||
(a[1] == b[1] && (a[2] < b[2] ||
(a[2] == b[2] && a[3] < b[3])))))) ? true : false);
}
%}.
^ super lessFromQDouble:aQDouble
"
(1.0 + 1e-40) > 1.0
((QDouble fromFloat:1.0) + (QDouble fromFloat:1e-40)) > (QDouble fromFloat:1.0)
(QDouble fromFloat:1.0) > (QDouble fromFloat:1.0)
(QDouble fromFloat:1.1) > (QDouble fromFloat:1.0)
(QDouble fromFloat:1.0) > 1.0
(QDouble fromFloat:1.1) > 1.0
1.0 > (QDouble fromFloat:1.0)
"
"Created: / 13-06-2017 / 17:07:47 / cg"
!
productFromFloat:aFloat
%{
if (__isFloatLike(aFloat)) {
double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
double b = __floatVal(aFloat);
double p0, p1, p2, p3;
double q0, q1, q2;
double s0, s1, s2, s3, s4;
OBJ newQD;
int savedCV;
fpu_fix_start(&savedCV);
m_two_prod(p0, a[0], b, q0);
m_two_prod(p1, a[1], b, q1);
m_two_prod(p2, a[2], b, q2);
p3 = a[3] * b;
s0 = p0;
m_two_sum(s1, q0, p1, s2);
m_three_sum(s2, q1, p2);
m_three_sum2(q1, q2, p3);
s3 = q1;
s4 = q2 + p2;
m_renorm5(s0, s1, s2, s3, s4);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, s0, s1, s2, s3);
RETURN( newQD );
}
%}.
^ super productFromFloat:aFloat.
"
loosing bits here:
(1e20+1.0)*2.0
(1e20+1.0)*1e20
(1e40+1.0)*2.0
but not here:
(1.0 asQDouble) * 2.0
((1e20 asQDouble) + (1.0)) * 2.0
((1e20 asQDouble) + (1.0)) * 100.0 10000000000000000000100.0
((1e20 asQDouble) + (1.0)) * 1000.0
((1e40 asQDouble) + (1.0)) * 2.0
2.0 * (QDouble fromFloat:1.0)
2.0 * (QDouble fromFloat:3.0)
(QDouble fromFloat:2.0) * (QDouble fromFloat:3.0)
QDouble ln2 DoubleArray(0.693147180559945 2.3190468138463E-17 5.70770843841621E-34 -3.58243221060181E-50)
2.0 * QDouble ln2 DoubleArray(1.38629436111989 4.6380936276926E-17 1.14154168768324E-33 -7.16486442120362E-50)
QDouble ln2 * 2.0
2.0 * ((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) DoubleArray(2E+20 2.0 0.0 0.0)
((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) * 2.0 DoubleArray(2E+20 4E+20 0.0 0.0)
((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) * (QDouble fromFloat:2.0) DoubleArray(2E+20 4E+20 0.0 0.0)
(QDouble fromFloat:2.0) * ((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) DoubleArray(2E+20 4E+20 0.0 0.0)
(2.0 * ((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0))) - (QDouble fromFloat:1e20) - (QDouble fromFloat:1e20)
(2.0 * (QDouble fromFloat:1.0)) asFloat
(1e20 * (QDouble fromFloat:1.0)) asFloat
(1e20 * (QDouble fromFloat:1.0) * 1e-20) asDoubleArray
"
"Created: / 13-06-2017 / 00:58:56 / cg"
"Modified: / 19-06-2017 / 16:48:18 / cg"
"Modified (comment): / 19-06-2017 / 18:11:43 / cg"
!
productFromQDouble:aQDouble
%{
if (__Class(aQDouble) == QDouble) {
double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
double *b = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
OBJ newQD;
int savedCV;
#define QD_IEEE_MUL
#ifndef QD_IEEE_MUL
// sloppy
double p0, p1, p2, p3, p4, p5;
double q0, q1, q2, q3, q4, q5;
double t0, t1;
double s0, s1, s2;
fpu_fix_start(&savedCV);
m_two_prod(p0, a[0], b[0], q0);
// fprintf(stderr, "%f * %f -> %f, %f\n", a[0], b[0], p0, q0);
m_two_prod(p1, a[0], b[1], q1);
m_two_prod(p2, a[1], b[0], q2);
// fprintf(stderr, "%f * %f -> %f, %f\n", a[0], b[1], p1, q1);
// fprintf(stderr, "%f * %f -> %f, %f\n", a[1], b[0], p2, q2);
m_two_prod(p3, a[0], b[2], q3);
m_two_prod(p4, a[1], b[1], q4);
m_two_prod(p5, a[2], b[0], q5);
// fprintf(stderr, "%f * %f -> %f, %f\n", a[0], b[2], p3, q3);
// fprintf(stderr, "%f * %f -> %f, %f\n", a[1], b[1], p4, q4);
// fprintf(stderr, "%f * %f -> %f, %f\n", a[2], b[0], p5, q5);
/* Start Accumulation */
m_three_sum(p1, p2, q0);
// fprintf(stderr, "mul7: after three_sum: p1:%e p2:%e q0:%e\n", p1, p2, q0);
/* Six-Three Sum of p2, q1, q2, p3, p4, p5. */
m_three_sum(p2, q1, q2);
// fprintf(stderr, "mul8: after three_sum: p2:%e q1:%e q2:%e\n", p2, q1, q2);
m_three_sum(p3, p4, p5);
// fprintf(stderr, "mul9: after three_sum: p3:%e p4:%e p5:%e\n", p3, p4, p5);
/* compute (s0, s1, s2) = (p2, q1, q2) + (p3, p4, p5). */
m_two_sum(s0, p2, p3, t0);
// fprintf(stderr, "mul10: after two_sum: s0:%e p2:%e p3:%e t0:%e\n", s0, p2, p3, t0);
m_two_sum(s1, q1, p4, t1);
// fprintf(stderr, "mul11: after two_sum: s1:%e q1:%e p4:%e t1:%e\n", s1, q1, p4, t1);
s2 = q2 + p5;
m_two_sum(s1, s1, t0, t0);
// fprintf(stderr, "mul12: after two_sum: s1:%e s1:%e t0:%e t0:%e\n", s1, s1, t0, t0);
s2 += (t0 + t1);
/* O(eps^3) order terms */
s1 += a[0]*b[3] + a[1]*b[2] + a[2]*b[1] + a[3]*b[0] + q0 + q3 + q4 + q5;
// fprintf(stderr, "before renorm5: p0:%e p2:%e s0:%e s1:%e s2:%e\n", p0, p1, s0, s1, s2);
m_renorm5(p0, p1, s0, s1, s2);
// fprintf(stderr, "after renorm5: p0:%e p2:%e s0:%e s1:%e s2:%e\n", p0, p1, s0, s1, s2);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, p0, p1, s0, s1);
#else
// accurate
double p0, p1, p2, p3, p4, p5;
double q0, q1, q2, q3, q4, q5;
double p6, p7, p8, p9;
double q6, q7, q8, q9;
double r0, r1;
double t0, t1;
double s0, s1, s2;
fpu_fix_start(&savedCV);
m_two_prod(p0, a[0], b[0], q0);
m_two_prod(p1, a[0], b[1], q1);
m_two_prod(p2, a[1], b[0], q2);
m_two_prod(p3, a[0], b[2], q3);
m_two_prod(p4, a[1], b[1], q4);
m_two_prod(p5, a[2], b[0], q5);
/* Start Accumulation */
m_three_sum(p1, p2, q0);
/* Six-Three Sum of p2, q1, q2, p3, p4, p5. */
m_three_sum(p2, q1, q2);
m_three_sum(p3, p4, p5);
/* compute (s0, s1, s2) = (p2, q1, q2) + (p3, p4, p5). */
m_two_sum(s0, p2, p3, t0);
m_two_sum(s1, q1, p4, t1);
s2 = q2 + p5;
m_two_sum(s1, s1, t0, t0);
s2 += (t0 + t1);
/* O(eps^3) order terms */
m_two_prod(p6, a[0], b[3], q6);
m_two_prod(p7, a[1], b[2], q7);
m_two_prod(p8, a[2], b[1], q8);
m_two_prod(p9, a[3], b[0], q9);
/* Nine-Two-Sum of q0, s1, q3, q4, q5, p6, p7, p8, p9. */
m_two_sum(q0, q0, q3, q3);
m_two_sum(q4, q4, q5, q5);
m_two_sum(p6, p6, p7, p7);
m_two_sum(p8, p8, p9, p9);
/* Compute (t0, t1) = (q0, q3) + (q4, q5). */
m_two_sum(t0, q0, q4, t1);
t1 += (q3 + q5);
/* Compute (r0, r1) = (p6, p7) + (p8, p9). */
m_two_sum(r0, p6, p8, r1);
r1 += (p7 + p9);
/* Compute (q3, q4) = (t0, t1) + (r0, r1). */
m_two_sum(q3, t0, r0, q4);
q4 += (t1 + r1);
/* Compute (t0, t1) = (q3, q4) + s1. */
m_two_sum(t0, q3, s1, t1);
t1 += q4;
/* O(eps^4) terms -- Nine-One-Sum */
t1 += a[1] * b[3] + a[2] * b[2] + a[3] * b[1] + q6 + q7 + q8 + q9 + s2;
// fprintf(stderr, "before accur-renorm5: p0:%e p1:%e s0:%e t0:%e t1:%e\n", p0, p1, s0, t0, t1);
m_renorm5(p0, p1, s0, t0, t1);
// fprintf(stderr, "after accur-renorm5: p0:%e p1:%e s0:%e t0:%e t1:%e\n", p0, p1, s0, t0, t1);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, p0, p1, s0, t0);
#endif
RETURN( newQD );
}
%}.
^ super productFromQDouble:aQDouble.
"
(QDouble fromFloat:1.0) * 2.0
2.0 * (QDouble fromFloat:1.0)
(QDouble fromFloat:1.0) * (QDouble fromFloat:2.0)
1e20 * (QDouble fromFloat:1.0)
(2.0 * (QDouble fromFloat:1.0)) asFloat
(1e20 * (QDouble fromFloat:1.0)) asFloat
(1e20 * (QDouble fromFloat:1.0) * 1e-20) asDoubleArray
( ((QDouble fromFloat:1.0) + (QDouble fromFloat:1e20)) * (QDouble fromFloat:2.0)) asDoubleArray
"
"Created: / 13-06-2017 / 01:06:22 / cg"
"Modified: / 05-07-2017 / 11:07:16 / cg"
!
quotientFromQDouble:aQDouble
"sloppy"
|q0 q1 q2 q3 r|
q0 := aQDouble d0 / self d0.
"/ Stdout showCR:('q0: %1 (a[0]=%2; b[0]=%3)\n' bindWith:q0 with:self d0 with:aQDouble d0).
r := aQDouble - (self * q0).
"/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
q1 := r d0 / self d0.
"/ Stdout showCR:('q1: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q1 with:r d0 with:aQDouble d0).
r := r - (self * q1).
"/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
q2 := r d0 / self d0.
"/ Stdout showCR:('q2: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q2 with:r d0 with:aQDouble d0).
r := r - (self * q2).
"/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
q3 := r d0 / self d0.
"/ Stdout showCR:('q3: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q3 with:r d0 with:aQDouble d0).
r := self class d0:q0 d1:q1 d2:q2 d3:q3.
"/ Stdout showCR:('before renorm: %1\n' bindWith:r asDoubleArray).
r renorm.
"/ Stdout showCR:('after renorm: %1\n' bindWith:r asDoubleArray).
^ r
"
2.0 / (QDouble fromFloat:2.0)
2.0 / (QDouble fromFloat:1.0)
1e20 / (QDouble fromFloat:1.0)
(2.0 / (QDouble fromFloat:1.0)) asFloat
(1e20 / (QDouble fromFloat:1.0)) asFloat
(QDouble fromFloat:2.0) / 2.0
(QDouble fromFloat:1e20) / 2.0
((QDouble fromFloat:1.0) / 2.0) asFloat
((QDouble fromFloat:1e20 / 2.0)) asFloat
((1e20 + (QDouble fromFloat:1.0) + 1e-20) / 2.0) asDoubleArray
((QDouble fromFloat:10.0) quotientFromQDouble: (QDouble fromFloat:1.234)) asDoubleArray
((QDouble fromFloat:1.234) / (QDouble fromFloat:10.0)) asDoubleArray
q0: 1.234000e-01 (a[0]=1.234000e+00; b[0]=1.000000e+01)
a: 1.234000e+00/0.000000e+00/0.000000e+00/0.000000e+00
b: 1.000000e+01/0.000000e+00/0.000000e+00/0.000000e+00
(b * q0): 1.234000e+00/-2.775558e-17/0.000000e+00/0.000000e+00
r: 2.775558e-17/0.000000e+00/0.000000e+00/0.000000e+00
q1: 2.775558e-18 (r[0]=2.775558e-17; b[0]=1.000000e+01)
r: -1.540744e-33/0.000000e+00/0.000000e+00/0.000000e+00
q2: -1.540744e-34 (r[0]=-1.540744e-33; b[0]=1.000000e+01)
r: 8.552847e-50/0.000000e+00/0.000000e+00/0.000000e+00
q3: 8.552847e-51 (r[0]=8.552847e-50; b[0]=1.000000e+01)
before renorm: 1.234000e-01/2.775558e-18/-1.540744e-34/8.552847e-51
after renorm: 1.234000e-01/2.775558e-18/-1.540744e-34/8.552847e-51
1.234/10.0 is: 0.123400 / 0.000000 / -0.000000 / 0.000000
"
"Created: / 13-06-2017 / 17:50:35 / cg"
"Modified (comment): / 15-06-2017 / 01:02:05 / cg"
!
sumFromFloat:aFloat
%{
if (__isFloatLike(aFloat)) {
double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
double b = __floatVal(aFloat);
double c0, c1, c2, c3;
double e;
OBJ newQD;
int savedCV;
fpu_fix_start(&savedCV);
m_two_sum(c0, a[0], b, e);
m_two_sum(c1 ,a[1], e, e);
m_two_sum(c2, a[2], e, e);
m_two_sum(c3, a[3], e, e);
m_renorm5(c0, c1, c2, c3, e);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, c0, c1, c2, c3);
RETURN( newQD );
}
%}.
^ super sumFromFloat:aFloat.
"
1.0 + (QDouble fromFloat:1.0)
1e20 + (QDouble fromFloat:1.0)
(1.0 + (QDouble fromFloat:1.0)) asFloat
(1e20 + (QDouble fromFloat:1.0)) asFloat
(1.0 + (QDouble fromFloat:1.0)) asDoubleArray
(1e20 + (QDouble fromFloat:1.0)) asDoubleArray
(1e20 + (QDouble fromFloat:1.0) + 1e-20) asDoubleArray
"
"Created: / 12-06-2017 / 17:16:41 / cg"
"Modified: / 14-06-2017 / 11:43:47 / cg"
!
sumFromInteger:anInteger
^ self sumFromFloat:(anInteger asFloat)
"
1 + (QDouble fromFloat:1.0)
1e20 asInteger + (QDouble fromFloat:1.0)
(1 + (QDouble fromFloat:1.0)) asFloat
(1e20 asInteger + (QDouble fromFloat:1.0)) asFloat
"
"Created: / 03-07-2017 / 10:35:46 / cg"
!
sumFromQDouble:aQDouble
%{
/* s = quick_three_accum(a, b, c) adds c to the dd-pair (a, b).
* If the result does not fit in two doubles, then the sum is
* output into s and (a,b) contains the remainder. Otherwise
* s is zero and (a,b) contains the sum. */
# define m_quick_three_accum(q3_outRef, q3_a, q3_b, q3_c) \
{ \
double s;\
int za, zb;\
\
m_two_sum(s, q3_b, q3_c, q3_b); \
m_two_sum(s, q3_a, s, q3_a); \
\
za = (q3_a != 0.0);\
zb = (q3_b != 0.0);\
\
if (za && zb) {\
q3_outRef = s;\
} else {\
if (!zb) {\
q3_b = q3_a;\
q3_a = s;\
} else {\
q3_a = s;\
}\
q3_outRef = 0.0;\
}\
}
if (__Class(aQDouble) == QDouble) {
double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
double *b = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
OBJ newQD;
int savedCV;
#define QD_IEEE_ADD
#ifndef xQD_IEEE_ADD
// sloppy_add...
# if 0
double s0, s1, s2, s3;
double t0, t1, t2, t3;
fpu_fix_start(&savedCV);
m_two_sum(s0, a[0], b[0], t0);
m_two_sum(s1, a[1], b[1], t1);
m_two_sum(s2, a[2], b[2], t2);
m_two_sum(s3, a[3], b[3], t3);
m_two_sum(s1, s1, t0, t0);
m_three_sum(s2, t0, t1);
m_three_sum2(s3, t0, t2);
t0 = t0 + t1 + t3;
m_renorm5(s0, s1, s2, s3, t0);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, s0, s1, s2, s3);
# else
/* Same as above, but addition re-organized to minimize
data dependency ... unfortunately some compilers are
not very smart to do this automatically */
double s0, s1, s2, s3;
double t0, t1, t2, t3;
double v0, v1, v2, v3;
double u0, u1, u2, u3;
double w0, w1, w2, w3;
fpu_fix_start(&savedCV);
s0 = a[0] + b[0];
s1 = a[1] + b[1];
s2 = a[2] + b[2];
s3 = a[3] + b[3];
v0 = s0 - a[0];
v1 = s1 - a[1];
v2 = s2 - a[2];
v3 = s3 - a[3];
u0 = s0 - v0;
u1 = s1 - v1;
u2 = s2 - v2;
u3 = s3 - v3;
w0 = a[0] - u0;
w1 = a[1] - u1;
w2 = a[2] - u2;
w3 = a[3] - u3;
u0 = b[0] - v0;
u1 = b[1] - v1;
u2 = b[2] - v2;
u3 = b[3] - v3;
t0 = w0 + u0;
t1 = w1 + u1;
t2 = w2 + u2;
t3 = w3 + u3;
m_two_sum(s1, s1, t0, t0);
m_three_sum(s2, t0, t1);
m_three_sum2(s3, t0, t2);
t0 = t0 + t1 + t3;
/* renormalize */
m_renorm5(s0, s1, s2, s3, t0);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, s0, s1, s2, s3);
# endif
#else
// ieee_add...
int i, j, k;
double s, t;
double u, v; /* double-length accumulator */
double x[4] = {0.0, 0.0, 0.0, 0.0};
fpu_fix_start(&savedCV);
i = j = k = 0;
if (abs(a[i]) > abs(b[j]))
u = a[i++];
else
u = b[j++];
if (abs(a[i]) > abs(b[j]))
v = a[i++];
else
v = b[j++];
m_quick_two_sum(u, u, v, v);
while (k < 4) {
if (i >= 4 && j >= 4) {
x[k] = u;
if (k < 3)
x[++k] = v;
break;
}
if (i >= 4)
t = b[j++];
else if (j >= 4)
t = a[i++];
else if (abs(a[i]) > abs(b[j])) {
t = a[i++];
} else
t = b[j++];
m_quick_three_accum(s, u, v, t);
if (s != 0.0) {
x[k++] = s;
}
}
/* add the rest. */
for (k = i; k < 4; k++)
x[3] += a[k];
for (k = j; k < 4; k++)
x[3] += b[k];
m_renorm4(x[0], x[1], x[2], x[3]);
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, x[0], x[1], x[2], x[3]);
#endif
RETURN(newQD);
}
%}.
^ super sumFromQDouble:aQDouble
"
(QDouble fromFloat:1.0) + (QDouble fromFloat:1.0)
(QDouble fromFloat:1.0) + 1.0
1.0 + (QDouble fromFloat:1.0)
((QDouble fromFloat:1.0) + (QDouble fromFloat:1.0)) asDoubleArray
((QDouble fromFloat:1.0) + 1.0) asDoubleArray
(1.0 + (QDouble fromFloat:1.0)) asDoubleArray
(1e-20 + (QDouble fromFloat:1.0)) asDoubleArray
(1e20 + (QDouble fromFloat:1.0)) asDoubleArray
"
"Created: / 12-06-2017 / 21:15:43 / cg"
"Modified: / 03-07-2017 / 23:09:11 / cg"
! !
!QDouble methodsFor:'inspecting'!
inspectorExtraAttributes
"extra (pseudo instvar) entries to be shown in an inspector."
^ super inspectorExtraAttributes
add:'-{doubles}' -> [ self asDoubleArray printString ];
yourself
"Created: / 12-06-2017 / 23:43:05 / cg"
"Modified (format): / 18-07-2017 / 19:54:48 / cg"
! !
!QDouble methodsFor:'mathematical functions'!
exp
"/ the exp_sloppy code is the algorithm from the original C++ qd package;
"/ however, it is inexact in the 37th digit
"/ Therefore, use the inherited code, which is slow, but more precise.
^ super exp
"
2.0 exp -> 7.38905609893065
2.0 asQDouble exp -> 7.38905609893065022723
2 asQDouble exp printfPrintString:'%70.68f'
-> 7.389056098930650227230427460575007813180315570551847324087127822432669
actual result:
-> 7.3890560989306502272304274605750078131803155705518473240871278225225737960790577633843124850791217947737531612654...
"
"Created: / 19-06-2017 / 01:49:32 / cg"
"Modified (comment): / 22-06-2017 / 14:32:53 / cg"
!
exp_sloppy
"/ this is the algorithm from the qd C++ package;
"/ however, it is inexact in the 37th digit
"/
"/ the inherited exp code is much slower, but more precise.
"/ Strategy: We first reduce the size of x by noting that
"/
"/ exp(kr + m * log(2)) = 2^m * exp(r)^k
"/
"/ where m and k are integers. By choosing m appropriately
"/ we can make |kr| <= log(2) / 2 = 0.347. Then exp(r) is
"/ evaluated using the familiar Taylor series. Reducing the
"/ argument substantially speeds up the convergence. */
|k inv_k d0 m mul_pwr2 mul2 r
s p t thresh eps i|
eps := 1.21543267145725e-63. "/ = 2^-209.
d0 := self d0.
(d0 <= -709.0) ifTrue:[
^ 0.0 asQDouble.
].
(d0 >= 709.0) ifTrue:[
^ Infinity positive
].
(d0 = 0.0) ifTrue:[
^ 1.0 asQDouble.
].
(d0 = 1.0) ifTrue:[
self isOne ifTrue:[
^ self class e
].
].
k := 1.0 ldexp:16.
inv_k := 1.0 / k.
m := (d0 / Float ln2 + 0.5) floor.
mul_pwr2 := [:a :b | QDouble d0:(a d0 * b) d1:(a d1 * b) d2:(a d2 * b) d3:(a d3 * b) ].
mul2 := [:a | QDouble d0:(a d0 * 2.0) d1:(a d1 * 2.0) d2:(a d2 * 2.0) d3:(a d3 * 2.0) ].
r := mul_pwr2 value:(self - (self class ln2 * m)) value:inv_k.
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.
] doWhile:[ (t asFloat abs > thresh) and:[i < 10] ].
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := (mul2 value:s) + s squared.
s := s + 1.0.
^ s ldexp:m asInteger.
"
1.0 exp -> 2.71828182845905
1.0 asQDouble exp -> 2.71828182845905
"
"Created: / 21-06-2017 / 13:48:27 / cg"
!
exp_withAccuracy:epsilon
"compute e^x of the receiver using a taylor series approximation.
This method is only invoked for limitedPrecisionReal classes, which do not compute
exp themself (i.e. QDouble)"
"/ uses taylor series:
"/ x x^2 x^3
"/ e^x = 1 + --- + --- + --- ...
"/ 1!! 2!! 3!!
|x2 facN num den approx delta|
x2 := self asFloat squared.
num := x2.
den := 2.
facN := 2.
approx := self + 1 + (num / den).
[
facN := facN + 1.
den := den * facN.
num := num * self.
delta := num / den.
"/ Transcript showCR:delta.
delta isNaN ifTrue:[self halt:'nan when dividing for delta'. num / den].
delta = 0 ifTrue:[self halt:'zero delta'].
approx := approx + delta.
] doUntil:[delta abs <= epsilon].
^ approx
"
wolfram:
7.389056098930650227230427460575007813180315570551847324087
2 asQDouble exp_withAccuracy:QDouble epsilon 7.38905609893065022723
2 asQDouble exp_withAccuracy:LongFloat epsilon 7.38905609893065022723
2 asQDouble exp_withAccuracy:Float epsilon 7.38905609893065022489
2.718281828459045235360287471352662497757247093699959574966...
1 asQDouble exp_withAccuracy:QDouble epsilon 2.71828182845904523536
1 asQDouble exp_withAccuracy:LongFloat epsilon 2.71828182845904523536
1 asQDouble exp_withAccuracy:Float epsilon 2.71828182845904522671
1.7392749415205010473946813036112352261479840577250084... × 10^18
42 asQDouble exp_withAccuracy:QDouble epsilon NAN
42 asQDouble exp_withAccuracy:LongFloat epsilon 1.73927494152050104739e18
42 asQDouble exp_withAccuracy:Float epsilon 1.73927494152050104739e18
"
"Created: / 04-07-2017 / 11:58:26 / cg"
"Modified: / 10-10-2017 / 16:03:34 / cg"
!
floor
"return the receiver truncated towards negative infinity"
%{
double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
OBJ newQD;
double x0, x1, x2, x3;
x1 = x2 = x3 = 0.0;
x0 =floor(a[0]);
if (x0 == a[0]) {
x1 = floor(a[1]);
if (x1 == a[1]) {
x2 = floor(a[2]);
if (x2 == a[2]) {
x3 = floor(a[3]);
}
}
m_renorm4(x0, x1, x2, x3);
}
__qNew_qdReal(newQD, x0, x1, x2, x3);
RETURN( newQD );
%}.
"
(QDouble fromFloat:4.0) floor
(QDouble fromFloat:4.1) floor
(QDouble fromFloat:0.1) floor
(0.1 + (QDouble fromFloat:1.0)) floor
(1e20 + (QDouble fromFloat:1.0)) floor
(QDouble fromFloat:1.5) floor
(QDouble fromFloat:0.5) floor
(QDouble fromFloat:-0.5) floor
(QDouble fromFloat:-1.5) floor
"
"Created: / 13-06-2017 / 01:52:44 / cg"
"Modified (comment): / 13-06-2017 / 17:33:19 / cg"
!
floor2
"return the receiver truncated towards negative infinity"
|f0 f1 f2 f3|
%{
double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
OBJ newQD;
double x0, x1, x2, x3;
x1 = x2 = x3 = 0.0;
x0 =floor(a[0]);
if (x0 == a[0]) {
x1 = floor(a[1]);
if (x1 == a[1]) {
x2 = floor(a[2]);
if (x2 == a[2]) {
x3 = floor(a[3]);
}
}
m_renorm4(x0, x1, x2, x3);
}
__qMKFLOAT(f0, x0);
__qMKFLOAT(f1, x1);
__qMKFLOAT(f2, x2);
__qMKFLOAT(f3, x3);
%}.
^ (f0 + f1 + f2 + f3) asInteger
"
(QDouble fromFloat:4.0) floor
(QDouble fromFloat:4.1) floor
(QDouble fromFloat:0.1) floor
(0.1 + (QDouble fromFloat:1.0)) floor
(1e20 + (QDouble fromFloat:1.0)) floor
(QDouble fromFloat:1.5) floor
(QDouble fromFloat:0.5) floor
(QDouble fromFloat:-0.5) floor
(QDouble fromFloat:-1.5) floor
"
"Created: / 13-06-2017 / 01:52:44 / cg"
"Modified (comment): / 13-06-2017 / 17:33:19 / cg"
!
ln
"return the natural logarithm of myself.
Raises an exception, if the receiver is less or equal to zero.
Not sure if this is really faster than using a taylor right away:
the three exp-computations at the end are done in qDouble and are tailors themself..."
|d0 x|
^ super ln.
d0 := self d0.
d0 = 1.0 ifTrue:[
self isOne ifTrue:[ ^ self class zero ].
].
d0 > 0.0 ifTrue:[
"/ initial approx.
x := d0 ln asQDouble.
"/ three more iterations of newton...
x := x + (self / (x exp)) - 1.0.
x := x + (self / (x exp)) - 1.0.
x := x + (self / (x exp)) - 1.0.
^ x
].
"/ now done via trapInfinity; was:
"/ d0 = 0.0 ifTrue:[
"/ ^ Infinity negative
"/ ].
"/ if you need -INF for a zero receiver, try Number trapInfinity:[...]
^ self class
raise:(self = 0 ifTrue:[#infiniteResultSignal] ifFalse:[#domainErrorSignal])
receiver:self
selector:#ln
arguments:#()
errorString:'bad receiver in ln (not strictly positive)'
"
-1 ln
-1.0 asQDouble ln
0.0 asQDouble ln
1.0 asQDouble ln
3.0 ln printfPrintString:'%60.58lf'
-> 1.0986122886681097821082175869378261268138885498046875000000'
^
3.0 asQDouble ln printfPrintString:'%60.58f'
-> 1.0986122886681096913952452369225257046474905578227494517347
3.0 asQDouble ln printfPrintString:'%70.68f'
-> 1.09861228866810969139524523692252570464749055782274945173469433364779
(3.0 asQDouble ln_withAccuracy:1e-64) printfPrintString:'%70.68f'
1.09861228866810969139524523692252570464749055782274945173469433364475
(3.0 asQDouble ln_withAccuracy:1e-100) printfPrintString:'%70.68f'
'1.098612288668109691395245236922525704647490557822749451734694333656909'
actual result:
-> 1.0986122886681096913952452369225257046474905578227494517346943336374942932186089668736157548137320887879700290659...
"
"Created: / 18-06-2017 / 23:32:54 / cg"
"Modified: / 04-07-2017 / 11:46:27 / cg"
!
negated
^ self class
d0:(self d0) negated
d1:(self d1) negated
d2:(self d2) negated
d3:(self d3) negated
"
(QDouble fromFloat:1.0) negated
((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) negated asDoubleArray
(((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0))
+ ((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0))) asDoubleArray
"
"Created: / 12-06-2017 / 20:14:55 / cg"
"Modified (comment): / 12-06-2017 / 23:46:57 / cg"
!
sqrt
"Return the square root of the receiver"
"this computes a roughly 65 digits precision result,
using sqrt from the double as an initial guess"
|guess|
guess := self asFloat sqrt asQDouble.
^ self sqrt_withAccuracy:(self epsilon) fromInitialGuess:guess
"FloatD is only correct in roughly 17 digits
2 sqrt printfPrintString:'%50.48f'
-> 1.414213562373095145474621858738828450441360473633
^
QDouble gives you much more:
2 asQDouble sqrt printfPrintString:'%50.48f'
-> 1.414213562373095048801688724209698078569671875377
2 asQDouble sqrt printfPrintString:'%60.58f'
-> 1.4142135623730950488016887242096980785696718753769480731767'
2 asQDouble sqrt printfPrintString:'%70.68f'
-> 1.41421356237309504880168872420969807856967187537694807317667973799602
actual digits:
-> 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309...
"
"Created: / 22-06-2017 / 13:53:47 / cg"
"Modified: / 22-06-2017 / 17:08:06 / mawalch"
"Modified (format): / 03-07-2017 / 12:09:04 / cg"
!
squared
"return receiver * receiver"
%{
double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
double p0, p1, p2, p3, p4, p5;
double q0, q1, q2, q3;
double s0, s1;
double t0, t1;
double t;
OBJ newQD;
// code makes use of some symetries to avoid a few multiplications...
m_two_sqr(p0, a[0], q0);
t = 2.0 * a[0];
m_two_prod(p1, t, a[1], q1);
m_two_prod(p2, t, a[2], q2);
m_two_sqr(p3, a[1], q3);
m_two_sum(p1, q0, p1, q0);
m_two_sum(q0, q0, q1, q1);
m_two_sum(p2, p2, p3, p3);
m_two_sum(s0, q0, p2, t0);
m_two_sum(s1, q1, p3, t1);
m_two_sum(s1, s1, t0, t0);
t0 += t1;
m_quick_two_sum(s1, s1, t0, t0);
m_quick_two_sum(p2, s0, s1, t1);
m_quick_two_sum(p3, t1, t0, q0);
p4 = 2.0 * a[0] * a[3];
p5 = 2.0 * a[1] * a[2];
m_two_sum(p4, p4, p5, p5);
m_two_sum(q2, q2, q3, q3);
m_two_sum(t0, p4, q2, t1);
t1 = t1 + p5 + q3;
m_two_sum(p3, p3, t0, p4);
p4 = p4 + q0 + t1;
m_renorm5(p0, p1, p2, p3, p4);
__qNew_qdReal(newQD, p0, p1, s0, s1);
RETURN( newQD );
%}.
"
(QDouble fromFloat:4.0) squared
(1e20 + (QDouble fromFloat:1.0)) squared
"
"Created: / 13-06-2017 / 01:27:58 / cg"
"Modified: / 22-06-2017 / 14:08:31 / cg"
! !
!QDouble methodsFor:'printing & storing'!
digitsWithPrecision:precision
<resource: #obsolete>
"generate digits and exponent.
if precision is >0, that many digits are generated.
If it is 0 the required number of digits is generated
(but never more than the decimalPrecision, which is 65)"
|numDigits r exp i d out str|
numDigits := precision+1. "/ number of digits
r := self abs.
self d0 = 0.0 ifTrue:[
^ { String new:(precision max:1) withAll:$0 . 0 }
].
out := WriteStream on:(String new:precision+5).
"/ determine approx. exponent
exp := self d0 abs log10 floor.
exp < -300 ifTrue:[
"/ 1e-305 asQDouble
r := r * (10.0 raisedToInteger:300).
r := r / (10.0 raisedToInteger:(exp+300)).
] ifFalse:[
exp > 300 ifTrue:[
"/ 1e305 asQDouble
"/ lexpr(x,exp) = x * 2 ^ exp
self halt.
r := r * (2 raisedTo:-53).
r := r / (10.0 asQDouble raisedTo: exp).
r := r * (2 raisedTo:53).
] ifFalse:[
r := r / (10.0 asQDouble raisedTo:exp).
]
].
"/ Fix exponent if we are off by one
(r >= 10.0) ifTrue:[
r := r / 10.0.
exp := exp + 1.
] ifFalse:[
(r < 1.0) ifTrue:[
r := r * 10.0.
exp := exp - 1.
]
].
((r >= 10.0) or:[ r < 1.0 ]) ifTrue:[
self error:'can''t compute exponent.'.
].
"/
"/ Extract the digits
"/ notice, that the d1,d2 and d3 components might
"/ be negative; therefore characters out of the 0..9 range
"/ might be produced here
"/
i := 1.
[ (precision ~~ 0 and:[ i <= numDigits ])
or:[ (precision == 0 and:[r d0 ~= 0.0]) ]] whileTrue:[
d := r d0 truncated.
r := r - d.
r := r * 10.0.
out nextPut:($0 + d).
i := i + 1.
].
numDigits := i-1.
str := out contents.
"/ Fix out-of-range digits.
numDigits to:2 by:-1 do:[:i |
(str at:i) < $0 ifTrue:[
str at:i-1 put:(str at:i-1) - 1.
str at:i put:(str at:i) + 10.
] ifFalse:[
(str at:i) > $9 ifTrue:[
str at:i-1 put:(str at:i-1) + 1.
str at:i put:(str at:i) - 10.
].
].
].
str first <= $0 ifTrue:[
self error:'non-positive leading digit'
].
"/ Round, handle carry
(str at:numDigits) >= $5 ifTrue:[
str at:numDigits-1 put:(str at:numDigits-1) + 1.
i := numDigits-1.
[i > 1 and:[(str at:i) > $9]] whileTrue:[
str at:i put:(str at:i) - 10.
i := i - 1.
str at:i put:(str at:i) + 1.
]
].
"/ If first digit is 10, shift everything.
str first > $9 ifTrue:[
exp := exp + 1.
str at:1 put:$0.
str := '1',str
].
^ { (str copyTo:numDigits-1) . exp }
"
0 asQDouble digitsWithPrecision:1 -> #('0' 0)
0 asQDouble digitsWithPrecision:0 -> #('0' 0)
1.2345 printfPrintString:'%.4f'
1.2345 asQDouble digitsWithPrecision:5 -> #('12345' 0)
--- but 1.2345 is not really what you think:
1.2345 printfPrintString:'%.20f'
1.2345 asQDouble digitsWithPrecision:20 -> #('12344999999999999307' 0)
12.345 asQDouble digitsWithPrecision:5 -> #('12345' 1)
12345 asQDouble digitsWithPrecision:5 -> #('12345' 4)
12345.1 asQDouble digitsWithPrecision:5 -> #('12345' 4)
12345.9 asQDouble digitsWithPrecision:5 -> #('12346' 4)
1.2345 asQDouble / 10.0 asQDouble
1.2345 asQDouble / 10.0
"
"Created: / 15-06-2017 / 09:10:01 / cg"
"Modified: / 16-06-2017 / 10:01:03 / cg"
!
printOn:aStream
"return a printed representation of the receiver.
Notice:
this code was adapted from an ugly piece of c++ code,
which was obviously hacked.
It does need a rework.
As an alternative, use the printf functions, which should also deal wth QDoubles"
thisContext isRecursive ifTrue:[
aStream nextPutAll:'aQDouble (error while printing)'.
^ self.
].
PrintfScanf printf:'%g' on:aStream argument:self.
"/ self
"/ printOn:aStream precision:40 width:0
"/ fixed:true showPositive:false uppercase:false fillChar:(Character space)
"
(1.2345 asQDouble) printOn:Transcript . Transcript cr.
(1.2345 asFloat) printOn:Transcript . Transcript cr.
(1.2345 asLongFloat) printOn:Transcript . Transcript cr.
(1.2345 asShortFloat) printOn:Transcript . Transcript cr.
((QDouble fromFloat:1.2345) / 10.0) printOn:Transcript
((QDouble fromFloat:1.2345) / 10000.0) printOn:Transcript
((QDouble fromFloat:1.2345) / 1000000000.0) printOn:Transcript
"
"Created: / 15-06-2017 / 01:51:36 / cg"
"Modified (comment): / 21-06-2017 / 09:55:10 / cg"
"Modified: / 05-06-2019 / 20:38:58 / Claus Gittinger"
!
printOn:aStream precision:precisionIn width:width
fixed:fixed showPositive:showPositive uppercase:uppercase fillChar:fillChar
<resource: #obsolete>
"return a printed representation of the receiver.
This is a parametrized entry, which can be used by printf-like functions.
Notice:
this code was adapted from an ugly piece of c++ code,
which was obviously hacked.
It does need a rework.
As an alternative, use the printf functions, which should also deal wth QDoubles
"
"
1.2345 asQDouble printString
12.345 asQDouble printString
12345 asQDouble printString
"
|sgn count delta exp precision|
"/ self d1 = 0.0 ifTrue:[
"/ self d0 printOn:aStream.
"/ ^ self.
"/ ].
count := 0.
sgn := true.
exp := 0.
precision := precisionIn.
self isInfinite ifTrue:[
self < 0 ifTrue:[
aStream nextPut:$-.
count := 1.
] ifFalse:[
showPositive ifTrue:[
aStream nextPut:$+.
count := 1.
] ifFalse:[
sgn := false.
].
].
uppercase ifTrue:[
aStream nextPutAll:'INF'
] ifFalse:[
aStream nextPutAll:'inf'
].
count := count + 3.
] ifFalse:[
self isNaN ifTrue:[
uppercase ifTrue:[
aStream nextPutAll:'NAN'
] ifFalse:[
aStream nextPutAll:'nan'
].
count := count + 3.
] ifFalse:[
self < 0 ifTrue:[
aStream nextPut:$-.
count := count + 1.
] ifFalse:[
showPositive ifTrue:[
aStream nextPut:$+.
count := count + 1.
] ifFalse:[
sgn := false.
].
].
self = 0.0 ifTrue:[
aStream nextPut:$0.
count := count + 1.
precision > 0 ifTrue:[
aStream nextPut:$..
count := count + 1.
precision timesRepeat:[ aStream nextPut:$0 ].
count := count + precision.
].
self halt.
] ifFalse:[
|off d d_width_extra|
"/ non-zero case
off := fixed
ifTrue:[ 1 + self asFloat abs log10 floor asInteger ]
ifFalse:[1].
d := precision + off.
d_width_extra := d.
fixed ifTrue:[
d_width_extra := 40 max:d.
].
"/ highly special case - fixed mode, precision is zero, abs(*this) < 1.0
"/ without this trap a number like 0.9 printed fixed with 0 precision prints as 0
"/ should be rounded to 1.
(fixed and:[ (precision == 0) and:[ (self abs < 1.0) ]]) ifTrue:[
(self abs >= 0.5) ifTrue:[
aStream nextPut:$1
] ifFalse:[
aStream nextPut:$0
].
^ self
].
"/ handle near zero to working precision (but not exactly zero)
(fixed and:[ d <= 0 ]) ifTrue:[
aStream nextPut:$0.
(precision > 0) ifTrue:[
aStream nextPut:$. .
aStream next:precision put:$0.
]
] ifFalse:[
"/ default
|t j|
t := self digitsWithPrecision:(fixed ifTrue:[d_width_extra] ifFalse:[d])+1.
exp := t second.
t := t first.
fixed ifTrue:[
"/ fix the string if it's been computed incorrectly
"/ round here in the decimal string if required
t := self round_string_qd:t at:(d + 1) offset:off.
precision := t at:3.
off := t at:2.
t := t at:1.
(off > 0) ifTrue:[
aStream next:off putAll:t startingAt:1.
(precision > 0) ifTrue:[
aStream nextPut:$. .
aStream next:precision-1 putAll:t startingAt:off+1.
]
] ifFalse:[
aStream nextPutAll:'0.'.
(off < 0) ifTrue:[
aStream next:off negated put:$0.
].
aStream next:d putAll:t startingAt:0.
]
] ifFalse:[
aStream nextPut:(t at:1).
(precision > 0) ifTrue:[
aStream nextPut:$. .
].
aStream next:precision putAll:t startingAt:2.
]
].
].
]
].
"/ trap for improper offset with large values
"/ without this trap, output of values of the for 10^j - 1 fail for j > 28
"/ and are output with the point in the wrong place, leading to a dramatically off value
"/ (fixed and:[ (precision > 0) ]) ifTrue:[
"/ "/ make sure that the value isn't dramatically larger
"/ from_string = atof(s.c_str());
"/
"/ // if this ratio is large, then we've got problems
"/ if( fabs( from_string / this->x[0] ) > 3.0 ){
"/
"/ int point_position;
"/ char temp;
"/
"/ // loop on the string, find the point, move it up one
"/ // don't act on the first character
"/ for(i=1; i < s.length(); i++){
"/ if(s[i] == '.'){
"/ s[i] = s[i-1] ;
"/ s[i-1] = '.' ;
"/ break;
"/ }
"/ }
"/
"/ from_string = atof(s.c_str());
"/ // if this ratio is large, then the string has not been fixed
"/ if( fabs( from_string / this->x[0] ) > 3.0 ){
"/ dd_real::error("Re-rounding unsuccessful in large number fixed point trap.") ;
"/ }
"/ }
"/ }
"/
fixed ifFalse:[
"/ Fill in exponent part
aStream nextPut:(uppercase ifTrue:[$E] ifFalse:[$e]).
aStream print:exp.
].
"/ fill in the blanks
(delta := width-count) > 0 ifTrue:[
self halt.
"/ if (fmt & ios_base::internal) {
"/ if (sgn)
"/ s.insert(static_cast<string::size_type>(1), delta, fill);
"/ else
"/ s.insert(static_cast<string::size_type>(0), delta, fill);
"/ } else if (fmt & ios_base::left) {
"/ s.append(delta, fill);
"/ } else {
"/ s.insert(static_cast<string::size_type>(0), delta, fill);
"/ }
"/ }
].
"Created: / 15-06-2017 / 02:37:31 / cg"
"Modified (comment): / 16-06-2017 / 14:48:30 / cg"
!
round_string_qd:str at:precisionIn offset:offsetIn
<resource: #obsolete>
"returns a triple of: { new-str . new-offset . new-precision }"
"/
"/ Input string must be all digits or errors will occur.
"/
|i numDigits offsetOut precisionOut|
numDigits := precisionIn.
offsetOut := offsetIn.
precisionOut := precisionIn.
"/ Round, handle carry
((str at:numDigits) >= $5) ifTrue:[
str at:numDigits-1 put:(str at:numDigits-1)+1.
i := numDigits-1.
[ i > 1 and:[ (str at:i) > $9] ] whileTrue:[
str at:i put:(str at:i)-10.
i := i - 1.
str at:i put:(str at:i)+1.
]
].
"/ If first digit is 10, shift everything.
(str at:1) > $9 ifTrue:[
"/ e++; // don't modify exponent here
str replaceFrom:2 with:str startingAt:1.
str at:1 put:$1.
str at:2 put:$0.
offsetOut := offsetOut + 1.
precisionOut := precisionOut + 1.
].
^ { (str copyTo:precisionOut) . offsetOut . precisionOut }
"Created: / 16-06-2017 / 10:12:39 / cg"
"Modified (comment): / 16-06-2017 / 11:22:03 / cg"
! !
!QDouble methodsFor:'private accessing'!
d0
"the most significant (and highest valued) 53 bits of precision"
%{
RETURN ( __MKFLOAT(__QuadDoubleInstPtr(self)->d_quadDoubleValue[0]) );
%}
"Created: / 12-06-2017 / 20:15:12 / cg"
"Modified (comment): / 13-06-2017 / 17:59:47 / cg"
!
d1
"the next most significant (and next highest valued) 53 bits of precision"
%{
RETURN ( __MKFLOAT(__QuadDoubleInstPtr(self)->d_quadDoubleValue[1]) );
%}
"Created: / 12-06-2017 / 20:15:12 / cg"
"Modified (comment): / 13-06-2017 / 18:00:00 / cg"
!
d2
%{
RETURN ( __MKFLOAT(__QuadDoubleInstPtr(self)->d_quadDoubleValue[2]) );
%}
"Created: / 12-06-2017 / 20:15:29 / cg"
!
d3
"the least significant (and smallest valued) 53 bits of precision"
%{
RETURN ( __MKFLOAT(__QuadDoubleInstPtr(self)->d_quadDoubleValue[3]) );
%}
"Created: / 12-06-2017 / 20:15:32 / cg"
"Modified (comment): / 13-06-2017 / 18:00:18 / cg"
!
renorm
"destructive renormalization"
%{
double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
double c0, c1, c2, c3;
c0 = a[0];
c1 = a[1];
c2 = a[2];
c3 = a[3];
m_renorm4(c0, c1, c2, c3);
a[0] = c0;
a[1] = c1;
a[2] = c2;
a[3] = c3;
RETURN( self );
%}.
^ self error.
"
(QDouble fromFloat:1.0) renorm
"
"Created: / 13-06-2017 / 18:05:33 / cg"
"Modified: / 15-06-2017 / 00:12:59 / cg"
! !
!QDouble methodsFor:'testing'!
isFinite
^ self d0 isFinite
"Created: / 17-06-2017 / 03:40:30 / cg"
!
isInfinite
^ self d0 isInfinite
"Created: / 15-06-2017 / 01:57:57 / cg"
!
isNaN
^ 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'!
version
^ '$Header$'
!
version_CVS
^ '$Header$'
! !