add sum!"decimal", wip: dynamic struct (#339)

* wip: dynamic struct * wip * wip * add decimal sum

Author: 9il

meson.build

@@ -69,6 +69,7 @@ sources_list = [
  'mir/series',
  'mir/small_array',
  'mir/small_string',
+ 'mir/string_map',
  'mir/type_info',
 ]
 

source/mir/bignum/decimal.d

@@ -752,6 +752,59 @@ struct Decimal(size_t maxSize64)
  {
  return exponent == exponent.max && !coefficient.length;
  }
+
+ ///
+ ref opOpAssign(string op, size_t rhsMaxSize64)(ref const Decimal!rhsMaxSize64 rhs) @safe pure return
+ if (op == "+" || op == "-")
+ {
+ BigInt!rhsMaxSize64 rhsCopy;
+ BigIntView!(const size_t) rhsView;
+ auto expDiff = exponent - rhs.exponent;
+ if (expDiff >= 0)
+ {
+ exponent = rhs.exponent;
+ coefficient.mulPow5(expDiff);
+ coefficient.opOpAssign!"<<"(expDiff);
+ rhsView = rhs.coefficient.view;
+ }
+ else
+ {
+ rhsCopy.copyFrom(rhs.coefficient.view);
+ rhsCopy.mulPow5(-expDiff);
+ rhsCopy.opOpAssign!"<<"(-expDiff);
+ rhsView = rhsCopy.view;
+ }
+ coefficient.opOpAssign!op(rhsView);
+ return this;
+ }
+
+ static if (maxSize64 == 3)
+ ///
+ version(mir_bignum_test) @safe pure @nogc unittest
+ {
+ import std.stdio;
+ auto a = Decimal!1("777.7");
+ auto b = Decimal!1("777");
+ import mir.format;
+ assert(stringBuf() << cast(double)a - cast(double)b << getData == "0.7000000000000455");
+ a -= b;
+ assert(stringBuf() << a << getData == "0.7");
+
+ a = Decimal!1("-777.7");
+ b = Decimal!1("777");
+ a += b;
+ assert(stringBuf() << a << getData == "-0.7");
+
+ a = Decimal!1("777.7");
+ b = Decimal!1("-777");
+ a += b;
+ assert(stringBuf() << a << getData == "0.7");
+
+ a = Decimal!1("777");
+ b = Decimal!1("777.7");
+ a -= b;
+ assert(stringBuf() << a << getData == "-0.7");
+ }
 }
 
 ///

source/mir/bignum/integer.d

@@ -324,22 +324,30 @@ struct BigInt(size_t maxSize64)
  Returns:
  overflow
  +/
- bool opOpAssign(string op)(ref const BigInt rhs)
+ bool opOpAssign(string op, size_t rhsMaxSize64)(ref const BigInt!rhsMaxSize64 rhs)
  @safe pure nothrow @nogc
  if (op == "+" || op == "-")
  {
- int diff = length - rhs.length;
+ return opOpAssign!op(rhs.view);
+ }
+
+ /// ditto
+ bool opOpAssign(string op)(BigIntView!(const size_t) rhs)
+ @safe pure nothrow @nogc
+ if (op == "+" || op == "-")
+ {
+ sizediff_t diff = length - rhs.coefficients.length;
  if (diff < 0)
  {
  auto oldLength = length;
- length = rhs.length;
+ length = cast(int)rhs.coefficients.length;
  view.unsigned.leastSignificantFirst[oldLength .. $] = 0;
  }
  else
- if (rhs.length == 0)
+ if (rhs.coefficients.length == 0)
  return false;
  auto thisView = view;
- auto overflow = thisView.opOpAssign!op(rhs.view);
+ auto overflow = thisView.opOpAssign!op(rhs);
  this.sign = thisView.sign;
  if (overflow)
  {
@@ -519,7 +527,7 @@ struct BigInt(size_t maxSize64)
  +/
  bool copyFrom(W, WordEndian endian)(BigIntView!(const W, endian) view)
  {
- static if (W.sizeof > size_t.sizeof)
+ static if (W.sizeof > size_t.sizeof && endian == TargetEndian)
  {
  return this.copyFrom(cast(BigIntView!(const size_t))view);
  }

source/mir/math/sum.d

@@ -20,6 +20,23 @@ unittest
  assert(r == ar.sum!"kbn");
  assert(r == ar.sum!"kb2");
  assert(r == ar.sum!"precise");
+ assert(r == ar.sum!"decimal");
+}
+
+/// Decimal precise summation
+version(mir_test)
+unittest
+{
+ auto ar = [777.7, -777];
+ assert(ar.sum!"decimal" == 0.7);
+
+ // The exact binary reuslt is 0.7000000000000455
+ assert(ar[0] + ar[1] == 0.7000000000000455);
+ assert(ar.sum!"fast" == 0.7000000000000455);
+ assert(ar.sum!"kahan" == 0.7000000000000455);
+ assert(ar.sum!"kbn" == 0.7000000000000455);
+ assert(ar.sum!"kb2" == 0.7000000000000455);
+ assert(ar.sum!"precise" == 0.7000000000000455);
 }
 
 ///
@@ -111,6 +128,7 @@ unittest
  assert(r == ar.sum!"kb2");
  }
  assert(r == ar.sum!"precise");
+ assert(r == ar.sum!"decimal");
 }
 
 ///
@@ -317,7 +335,7 @@ enum Summation
  Precise summation algorithm.
  The value of the sum is rounded to the nearest representable
  floating-point number using the $(LUCKY round-half-to-even rule).
- The result can differ from the exact value on `X86`, `nextDown(proir) <= result && result <= nextUp(proir)`.
+ The result can differ from the exact value on 32bit `x86`, `nextDown(proir) <= result && result <= nextUp(proir)`.
  The current implementation re-establish special value semantics across iterations (i.e. handling ±inf).
 
  References: $(LINK2 http://www.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps,
@@ -335,6 +353,17 @@ enum Summation
  +/
  precise,
 
+ /++
+ Precise decimal summation algorithm.
+
+ The elements of the sum are converted to a shortest decimal representation that being converted back would result the same floating-point number.
+ The resulting decimal elements are summed without rounding.
+ The decimal sum is converted back to a binary floating point representation using round-half-to-even rule.
+
+ See_also: The $(HTTPS https://github.com/ulfjack/ryu, Ryu algorithm)
+ +/
+ decimal,
+
  /++
  $(WEB en.wikipedia.org/wiki/Kahan_summation, Kahan summation) algorithm.
  +/
@@ -611,6 +640,13 @@ struct Summator(T, Summation summation)
  F s = summationInitValue!F;
  }
  else
+ static if (summation == Summation.decimal)
+ {
+ import mir.bignum.decimal;
+ Decimal!256 s;
+ T ss = 0;
+ }
+ else
  static assert(0, "Unsupported summation type for std.numeric.Summator.");
 
 
@@ -675,6 +711,17 @@ public:
  s = n;
  }
  else
+ static if (summation == Summation.decimal)
+ {
+ ss = 0;
+ if (!(-n.infinity < n && n < n.infinity))
+ {
+ ss = n;
+ n = 0;
+ }
+ s = n;
+ }
+ else
  static assert(0);
  }
 
@@ -895,6 +942,20 @@ public:
  s += n;
  }
  else
+ static if (summation == summation.decimal)
+ {
+ import mir.bignum.internal.ryu.generic_128: genericBinaryToDecimal;
+ if (-n.infinity < n && n < n.infinity)
+ {
+ auto decimal = genericBinaryToDecimal(n);
+ s += decimal;
+ }
+ else
+ {
+ ss += n;
+ }
+ }
+ else
  static assert(0);
  }
 
@@ -1134,7 +1195,12 @@ public:
  return s;
  }
  else
- static assert(0);
+ static if (summation == Summation.decimal)
+ {
+ return cast(T) s + ss;
+ }
+ else
+ static assert(0);
  }
 
  version(none)
@@ -1311,6 +1377,11 @@ public:
  s = rhs;
  }
  else
+ static if (summation == summation.decimal)
+ {
+ __ctor(rhs);
+ }
+ else
  static assert(0);
  }
 
@@ -1671,13 +1742,21 @@ template sum(F, Summation summation = Summation.appropriate)
  ///
  F sum(Range r)
  {
- static if (isComplex!F && summation == Summation.precise)
+ static if (isComplex!F && (summation == Summation.precise || summation == Summation.decimal))
  {
  return sum(r, summationInitValue!F);
  }
  else
  {
- Summator!(F, ResolveSummationType!(summation, Range, sumType!Range)) sum;
+ static if (summation == Summation.decimal)
+ {
+ Summator!(F, summation) sum = void;
+ sum = 0;
+ }
+ else
+ {
+ Summator!(F, ResolveSummationType!(summation, Range, sumType!Range)) sum;
+ }
  sum.put(r.move);
  return sum.sum;
  }
@@ -1686,11 +1765,22 @@ template sum(F, Summation summation = Summation.appropriate)
  ///
  F sum(Range r, F seed)
  {
- static if (isComplex!F && summation == Summation.precise)
+ static if (isComplex!F && (summation == Summation.precise || summation == Summation.decimal))
  {
  alias T = typeof(F.init.re);
- auto sumRe = Summator!(T, Summation.precise)(seed.re);
- auto sumIm = Summator!(T, Summation.precise)(seed.im);
+ static if (summation == Summation.decimal)
+ {
+ Summator!(T, summation) sumRe = void;
+ sumRe = seed.re;
+
+ Summator!(T, summation) sumIm = void;
+ sumIm = seed.im;
+ }
+ else
+ {
+ auto sumRe = Summator!(T, Summation.precise)(seed.re);
+ auto sumIm = Summator!(T, Summation.precise)(seed.im);
+ }
  import mir.ndslice.slice: isSlice;
  static if (isSlice!Range)
  {
@@ -1713,7 +1803,15 @@ template sum(F, Summation summation = Summation.appropriate)
  }
  else
  {
- auto sum = Summator!(F, ResolveSummationType!(summation, Range, F))(seed);
+ static if (summation == Summation.decimal)
+ {
+ Summator!(F, summation) sum = void;
+ sum = seed;
+ }
+ else
+ {
+ auto sum = Summator!(F, ResolveSummationType!(summation, Range, F))(seed);
+ }
  sum.put(r.move);
  return sum.sum;
  }
@@ -1723,7 +1821,7 @@ template sum(F, Summation summation = Summation.appropriate)
  ///
  F sum(scope const F[] r...)
  {
- static if (isComplex!F && summation == Summation.precise)
+ static if (isComplex!F && (summation == Summation.precise || summation == Summation.decimal))
  {
  return sum(r, summationInitValue!F);
  }

source/mir/ndslice/sorting.d

@@ -660,8 +660,8 @@ version(mir_test)
  assert(j == 3);
  auto upperBound = a[j .. $];
 
- assert(a.transitionIndex(a[$-1]) == a.length - 1);
- assert(a.transitionIndex!"a <= b"(a[$-1]) == a.length);
+ assert(a.transitionIndex(a[$ - 1]) == a.length - 1);
+ assert(a.transitionIndex!"a <= b"(a[$ - 1]) == a.length);
 }
 
 /++