#FEATURE by cg
class: QDouble
added:
#isOne
#isZero
#ln
class: QDouble class
added:
#infinity
#negativeInfinity
"
COPYRIGHT (c) 2017 by eXept Software AG
All Rights Reserved
This software is furnished under a license and may be used
only in accordance with the terms of that license and with the
inclusion of the above copyright notice. This software may not
be provided or otherwise made available to, or used by, any
other person. No title to or ownership of the software is
hereby transferred.
"
"{ Package: 'stx:libbasic2' }"
"{ NameSpace: Smalltalk }"
LimitedPrecisionReal variableByteSubclass:#QDouble
instanceVariableNames:''
classVariableNames:'DefaultPrintFormat QDoubleZero QDoubleOne Pi E Ln2 Ln10 Epsilon
FMax FMin'
poolDictionaries:''
category:'Magnitude-Numbers'
!
!QDouble primitiveDefinitions!
ry:'Magnitude-Numbers'
! !
!QDouble primitiveFunctions!
= quick_two_sum(s1, *c3Ptr, &s2);
} else {
s0 = quick_two_sum(s0, *c2Ptr, &s1);
if (s1 ! !
!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.
This does not work correctly at the moment.
QDoubles represent rational numbers with extended, but still limited precision.
In contrast to Floats (which use the C-compilers native 64bit 'double' format),
QDoubles give you 212 bit or approx. 64 decimal digits of precision.
Representation:
QDoubles use 4 IEEE doubles, each keeping 53 bits of precision.
Range and Precision of Storage Formats: see LimitedPrecisionReal >> documentation
Float decimalPrecision
QDouble decimalPrecision
[author:]
Claus Gittinger
[see also:]
Number
Float ShortFloat Fraction FixedPoint Integer Complex
FloatArray DoubleArray
"
! !
!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)) {
__qNew(newFloat, sizeof(struct __quadDoubleStruct));
__stx_setClass(newFloat, QDouble);
__QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[0] = (double)(__smallIntegerVal(anInteger));
__QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[1] = 0.0;
__QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[2] = 0.0;
__QuadDoubleInstPtr(newFloat)->d_quadDoubleValue[3] = 0.0;
RETURN (newFloat);
}
#endif
%}.
^ super fromInteger:anInteger
"
self fromInteger:2
"
"Created: / 12-06-2017 / 16:10:10 / 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'!
e
"return the constant e as quad precision double.
(returns approx. 212 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. 212 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"
!
ln10
"return the constant e as quad precision double.
(returns approx. 212 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. 212 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. 212 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"
^ 10
"
ShortFloat defaultPrintPrecision
Float defaultPrintPrecision
LongFloat defaultPrintPrecision
QDouble defaultPrintPrecision
"
"Created: / 17-06-2017 / 02:58:51 / cg"
!
epsilon
"return the maximum relative spacing of instances of mySelf
(i.e. the value-delta of the least significant bit)"
^ 1.2154326714572500565324311366323150942261000827598106963711353e-63
"
Float epsilon
ShortFloat epsilon
QDouble epsilon
"
"Created: / 12-06-2017 / 18:52:44 / cg"
"Modified: / 14-06-2017 / 19:12:23 / cg"
!
numBitsInExponent
"answer the number of bits in the exponent"
^ Float numBitsInExponent
"
1.0 asQDouble class numBitsInExponent
"
"Created: / 12-06-2017 / 11:11:04 / cg"
!
numBitsInMantissa
"answer the number of bits in the mantissa"
^ Float precision * 4
"
1.0 class numBitsInMantissa
1.0 asShortFloat class numBitsInMantissa
1.0 asLongFloat class numBitsInMantissa
1.0 asDDouble class numBitsInMantissa
1.0 asQDouble class numBitsInMantissa
Float numBitsInMantissa
ShortFloat numBitsInMantissa
QDouble numBitsInMantissa
"
"Created: / 12-06-2017 / 11:13:44 / cg"
"Modified: / 12-06-2017 / 18:48:48 / cg"
!
precision
"answer the number of bits in the mantissa"
^ Float precision * 4
"
1.0 class numBitsInMantissa
1.0 asShortFloat class numBitsInMantissa
1.0 asLongFloat class numBitsInMantissa
1.0 asDDouble class numBitsInMantissa
1.0 asQDouble class numBitsInMantissa
Float numBitsInMantissa
ShortFloat numBitsInMantissa
QDouble numBitsInMantissa
QDouble precision
"
"Created: / 12-06-2017 / 18:49:11 / 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"
!
asQDouble
"return a QDouble with same value as myself."
^ self
"Created: / 15-06-2017 / 12:08:02 / 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"
!
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
^ 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"
! !
!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.
"
(QDouble fromFloat:1.0) productFromFloat:2.0
((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) productFromFloat:2.0
((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) productFromFloat:2e20
2.0 * (QDouble fromFloat:1.0)
2.0 * (QDouble fromFloat:3.0)
2.0 * ((QDouble fromFloat:1e20) + (QDouble fromFloat:1.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: / 14-06-2017 / 11:42:57 / cg"
!
productFromQDouble:aQDouble
%{
if (__Class(aQDouble) == QDouble) {
double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
double *b = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
OBJ newQD;
// sloppy
double p0, p1, p2, p3, p4, p5;
double q0, q1, q2, q3, q4, q5;
double t0, t1;
double s0, s1, s2;
int savedCV;
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);
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: / 14-06-2017 / 16:59:58 / 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"
!
sumFromQDouble:aQDouble
%{
if (__Class(aQDouble) == QDouble) {
double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
double *b = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
OBJ newQD;
#ifndef QD_IEEE_ADD
// sloppy_add...
# if 0
double s0, s1, s2, s3;
double t0, t1, t2, t3;
int savedCV;
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);
RETURN(newQD);
# 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;
int savedCV;
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);
RETURN(newQD);
# 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};
int savedCV;
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++];
u = quick_two_sum(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++];
s = quick_three_accum(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]);
RETURN(newQD);
#endif
}
%}.
^ 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: / 15-06-2017 / 00:49:10 / 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"
! !
!QDouble methodsFor:'mathematical functions'!
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 );
}
__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"
!
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
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"
!
squared
%{
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;
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"
! !
!QDouble methodsFor:'printing & storing'!
digitsWithPrecision:precision
"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"
PrintfScanf printf:'%f' 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
"
"Created: / 15-06-2017 / 01:51:36 / cg"
"Modified: / 17-06-2017 / 03:30:41 / cg"
!
printOn:aStream precision:precisionIn width:width
fixed:fixed showPositive:showPositive uppercase:uppercase fillChar:fillChar
"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"
!
printString
"return a printed representation of the receiver"
^ String streamContents:[:s | self printOn:s]
"Created: / 12-06-2017 / 23:41:04 / cg"
"Modified: / 15-06-2017 / 01:53:46 / cg"
!
round_string_qd:str at:precisionIn offset:offsetIn
"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$'
! !