On Saturday, 13 January 2018 at 19:28:40 UTC, rumbu wrote:
On Saturday, 13 January 2018 at 18:37:10 UTC, kdevel wrote:
I get large numerical dicrepancies and an exception:
That's because you are mixing floating point and decimal.
Just to take one example: double 1.1 cannot be represented
exactly as floating point and it's in fact
1.10000002384185791015625.
Sure. double has 53 mantissa bits but float only 24. But that is
not my point.
I am *not* talking about the 1e-18 but about these lines:
1.10 0.891207360061435339970703
0.878166666666666700439167 -1e-02
1.20 0.932039085967226349689098
0.932735999999999982229600 1e-04
1.30 0.963558185417192964729825
0.964774416666666678037840
1e-03
and
5.70 -0.550685542597637763537807
-25.165500000000002545915000 -
At 5.7 the value is obviously out of range.
sin is calculated using a Taylor series: sin(x) = x - x^3/3! +
x^5/5! - x^7/7! ... and so on. Raising to power all that junk
after 1.1 will lead finally to error. If you really want to
find out sin(1.1) using decimal, try sin(decimal128("1.1")) or
sin(decimal128(1)/10);
Chapter 9.2 of IEEE-754-2008 says the domain of sin(x) is (-inf,
inf). So if the argument x is outside the radius of convergence
of the Taylor series x must be reduced modulo 2pi such that it
fits.
For exact values like sin(1.0), I let you decide which one is
more exact:
Wolfram Alpha: 0.8414709848078965066525023216302989
real: 0.8414709848078965066645910000000000
double: 0.8414709848078965048700000000000000
float: 0.8414709568023681640600000000000000
decimal128: 0.8414709848078965066329679978908351
decimal64: 0.8414709848078965000000000000000000
decimal32: 0.8414710000000000000000000000000000
Anyway, I wouldn't call a difference at the 18th digit a "large
discrepancy" when we are talking about irrational numbers.
See above.
Regarding the exception, I cannot reproduce it, but I'll look
into it.
Nice.
Thank you for your exhaustive testing :)
There's more to come (sorry for the greek symbols, could not
resist to check this utf-8 feature):
unity.d
```
import std.stdio;
import std.typecons;
import std.range;
import std.math;
import decimal;
immutable size_t N = 100;
void unity (T) ()
{
writeln ("\n=== ", T.stringof, " ===\n");
immutable one = T (1);
immutable two = T (2);
immutable π = atan (one) * 4;
writefln!"π = <%30.24f>" (π);
foreach (i; iota(N + 1)) {
auto φ = two * π * i / N;
auto sinφ = sin (φ);
auto cosφ = cos (φ);
auto unity = sinφ * sinφ + cosφ * cosφ;
auto δ = one - unity;
writeln ("φ = <", φ, ">, δ = <", δ, ">");
}
}
void main ()
{
// unity!float;
// unity!double;
// unity!real;
// unity!decimal32;
unity!decimal64;
unity!decimal128;
}
```
Selected lines from the output produced here (Linux, DMD64 D
Compiler v2.077.1):
=== Decimal!64 ===
π = < 3.141592653589793000000000>
φ = <0>, δ = <0>
φ = <0.0628319>, δ = <0>
φ = <0.125664>, δ = <1e-16> possibly okay
:
φ = <0.942478>, δ = <2.41721e-06> too large
:
φ = <2.45044>, δ = <-6.57811e-07> ditto
φ = <2.51327>, δ = <-3.83012e-10> ditto
φ = <2.57611>, δ = <0.479476> ditto
:
:
:
=== Decimal!128 ===
π = < 3.141592653589793238612055>
φ = <0>, δ = <0>
φ = <1e-02>, δ = <-1.73780e-22>
:
Something went wrong with printing φ. And also this program
crashes:
core.exception.RangeError@decimal/package.d(6652): Range violation
----------------
??:? _d_arrayboundsp [0x80bc65f]
??:? void decimal.sinkFloat!(char,
decimal.integrals.unsigned!(128).unsigned).sinkFloat(ref
const(std.format.FormatSpec!(char).FormatSpec), void
delegate(const(char)[]),
const(decimal.integrals.unsigned!(128).unsigned), const(int),
const(bool), const(decimal.RoundingMode), const(bool)) [0x80b2976]
??:? void decimal.sinkGeneral!(char,
decimal.integrals.unsigned!(128).unsigned).sinkGeneral(ref
const(std.format.FormatSpec!(char).FormatSpec), void
delegate(const(char)[]),
const(decimal.integrals.unsigned!(128).unsigned), const(int),
const(bool), const(decimal.RoundingMode)) [0x80b33a4]
??:? void decimal.sinkDecimal!(char,
decimal.Decimal!(128).Decimal).sinkDecimal(ref
const(std.format.FormatSpec!(char).FormatSpec), void
delegate(const(char)[]), ref
const(decimal.Decimal!(128).Decimal),
const(decimal.RoundingMode)) [0x80b241a]
??:? const void
decimal.Decimal!(128).Decimal.toString!(char).toString(scope void
delegate(const(char)[]), std.format.FormatSpec!(char).FormatSpec)
[0x80b22bb]
??:? void
std.format.formatObject!(std.stdio.File.LockingTextWriter,
decimal.Decimal!(128).Decimal, char).formatObject(ref
std.stdio.File.LockingTextWriter, ref
decimal.Decimal!(128).Decimal, ref
const(std.format.FormatSpec!(char).FormatSpec)) [0x80b73cf]
??:? void
std.format.formatValue!(std.stdio.File.LockingTextWriter,
decimal.Decimal!(128).Decimal, char).formatValue(ref
std.stdio.File.LockingTextWriter, ref
decimal.Decimal!(128).Decimal, ref
const(std.format.FormatSpec!(char).FormatSpec)) [0x80b7393]
??:? uint
std.format.formattedWrite!(std.stdio.File.LockingTextWriter,
char, decimal.Decimal!(128).Decimal).formattedWrite(ref
std.stdio.File.LockingTextWriter, const(char[]),
decimal.Decimal!(128).Decimal) [0x80b6efd]
??:? void std.stdio.File.write!(immutable(char)[],
decimal.Decimal!(128).Decimal, immutable(char)[],
decimal.Decimal!(128).Decimal, immutable(char)[],
char).write(immutable(char)[], decimal.Decimal!(128).Decimal,
immutable(char)[], decimal.Decimal!(128).Decimal,
immutable(char)[], char) [0x8093176]
??:? void std.stdio.writeln!(immutable(char)[],
decimal.Decimal!(128).Decimal, immutable(char)[],
decimal.Decimal!(128).Decimal,
immutable(char)[]).writeln(immutable(char)[],
decimal.Decimal!(128).Decimal, immutable(char)[],
decimal.Decimal!(128).Decimal, immutable(char)[]) [0x80930c8]
??:? void unity.unity!(decimal.Decimal!(128).Decimal).unity()
[0x8092f62]
??:? _Dmain [0x80905ec]
