Another take on decimal data types

[kdevel]

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 at 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]