@@ -10,19 +10,23 @@ module mir.math.numeric;
1010import mir.math.common;
1111import mir.primitives;
1212import mir.primitives: isInputRange;
13- import std.traits : CommonType, Unqual, isIterable, ForeachType, isPointer, isNumeric ;
13+ import std.traits : CommonType, Unqual, isIterable, ForeachType, isPointer;
1414import mir.internal.utility: isFloatingPoint;
1515
16- deprecated ( " Use mir.bignum.fp.Fp!128 instead. " )
16+ // /
1717struct ProdAccumulator (T)
1818 if (isFloatingPoint! T)
1919{
2020 alias F = Unqual! T;
2121
22+ // /
2223 long exp = 1L ;
24+ // /
2325 F x = cast (F) 0.5 ;
26+ // /
2427 alias mantissa = x;
2528
29+ // /
2630 @safe pure @nogc nothrow
2731 this (F value)
2832 {
@@ -33,13 +37,15 @@ struct ProdAccumulator(T)
3337 this .exp = lexp;
3438 }
3539
40+ // /
3641 @safe pure @nogc nothrow
3742 this (long exp, F x)
3843 {
3944 this .exp = exp;
4045 this .x = x;
4146 }
4247
48+ // /
4349 @safe pure @nogc nothrow
4450 void put (U)(U e)
4551 if (is (U : T))
@@ -60,6 +66,7 @@ struct ProdAccumulator(T)
6066 }
6167 }
6268
69+ // /
6370 @safe pure @nogc nothrow
6471 void put (ProdAccumulator! T value)
6572 {
@@ -72,6 +79,7 @@ struct ProdAccumulator(T)
7279 }
7380 }
7481
82+ // /
7583 void put (Range )(Range r)
7684 if (isIterable! Range )
7785 {
@@ -101,6 +109,7 @@ struct ProdAccumulator(T)
101109 }
102110 }
103111
112+ // /
104113 @safe pure @nogc nothrow
105114 T prod () const scope @property
106115 {
@@ -112,27 +121,74 @@ struct ProdAccumulator(T)
112121 return ldexp (mantissa, e);
113122 }
114123
124+ // /
115125 @safe pure @nogc nothrow
116126 ProdAccumulator! T ldexp (long exp) const
117127 {
118128 return typeof (return )(this .exp + exp, mantissa);
119129 }
120130
131+ // ///
121132 alias opOpAssign (string op : " *" ) = put;
122133
134+ // /
123135 @safe pure @nogc nothrow
124136 ProdAccumulator! T opUnary (string op : " -" )() const
125137 {
126138 return typeof (return )(exp, - mantissa);
127139 }
128140
141+ // /
129142 @safe pure @nogc nothrow
130143 ProdAccumulator! T opUnary (string op : " +" )() const
131144 {
132145 return typeof (return )(exp, + mantissa);
133146 }
134147}
135148
149+ // /
150+ version (mir_test)
151+ @safe pure nothrow
152+ unittest
153+ {
154+ import mir.ndslice.slice: sliced;
155+
156+ ProdAccumulator! float x;
157+ x.put([1 , 2 , 3 ].sliced);
158+ assert (x.prod == 6f );
159+ x.put(4 );
160+ assert (x.prod == 24f );
161+ }
162+
163+ version (mir_test)
164+ @safe pure @nogc nothrow
165+ unittest
166+ {
167+ import mir.ndslice.slice: sliced;
168+
169+ static immutable a = [1 , 2 , 3 ];
170+ ProdAccumulator! float x;
171+ x.put(a);
172+ assert (x.prod == 6f );
173+ x.put(4 );
174+ assert (x.prod == 24f );
175+ static assert (is (typeof (x.prod) == float ));
176+ }
177+
178+ version (mir_test)
179+ @safe pure nothrow
180+ unittest
181+ {
182+ import mir.ndslice.slice: sliced;
183+
184+ ProdAccumulator! double x;
185+ x.put([1.0 , 2.0 , 3.0 ]);
186+ assert (x.prod == 6.0 );
187+ x.put(4.0 );
188+ assert (x.prod == 24.0 );
189+ static assert (is (typeof (x.prod) == double ));
190+ }
191+
136192package template prodType(T)
137193{
138194 import mir.math.sum: elementType;
@@ -171,48 +227,30 @@ See_also:
171227$(MREF mir, algorithm, iteration, reduce)
172228$(MREF mir, algorithm, iteration, fold)
173229+/
174- template prod (uint size)
230+ F prod (F, Range )(Range r)
231+ if (isFloatingPoint! F && isIterable! Range )
175232{
176- import mir.bignum.fp: Fp;
177-
178- Fp! size prod (Range )(Range r)
179- if (isIterable! Range )
180- {
181- import core.lifetime : move;
182- import mir.algorithm.iteration: each;
183- Fp! size prod = 1LU;
184- each! ((e) {
185- prod *= Fp! 128 (e);
186- })(move(r));
187- return prod;
188- }
189- }
233+ import core.lifetime : move;
190234
191- // / ditto
192- template prod (F)
193- if (isFloatingPoint! F)
194- {
195- F prod (Range )(Range r)
196- if (isFloatingPoint! F && isIterable! Range )
197- {
198- import core.lifetime : move;
199- return cast (F) .prod! 128 (r);
200- }
235+ ProdAccumulator! F prod;
236+ prod.put(r.move);
237+ return prod.prod;
201238}
202239
203-
204- deprecated (" Use prod!128 instead." )
240+ /+ +
241+ Params:
242+ r = finite iterable range
243+ exp = value of exponent
244+ Returns:
245+ The mantissa, such that the product equals the mantissa times 2^^exp
246+ +/
205247F prod (F, Range )(Range r, ref long exp)
206248 if (isFloatingPoint! F && isIterable! Range )
207249{
208250 import core.lifetime : move;
209251
210- import core.lifetime : move;
211-
212- auto p = .prod! 128 (r.move);
213-
214252 ProdAccumulator! F prod;
215- prod.put(cast (F) p );
253+ prod.put(r.move );
216254 exp = prod.exp;
217255 return prod.x;
218256}
@@ -228,10 +266,17 @@ prodType!Range prod(Range)(Range r)
228266{
229267 import core.lifetime : move;
230268
231- return cast (typeof (return )) .prod! 128 (r.move);
269+ alias F = typeof (return );
270+ return .prod! (F, Range )(r.move);
232271}
233272
234- deprecated (" Use prod!128 instead." )
273+ /+ +
274+ Params:
275+ r = finite iterable range
276+ exp = value of exponent
277+ Returns:
278+ The mantissa, such that the product equals the mantissa times 2^^exp
279+ +/
235280prodType! Range prod (Range )(Range r, ref long exp)
236281 if (isIterable! Range )
237282{
@@ -249,7 +294,10 @@ Returns:
249294+/
250295prodType! T prod (T)(scope const T[] ar... )
251296{
252- return cast (typeof (return )) .prod! 128 (ar);
297+ alias F = typeof (return );
298+ ProdAccumulator! F prod;
299+ prod.put(ar);
300+ return prod.prod;
253301}
254302
255303// / Product of arbitrary inputs
@@ -273,6 +321,11 @@ unittest
273321 auto r = [l, l, l, s, s, s, 0.8 * 2.0 ^^ 10 ];
274322
275323 assert (r.prod == 0.8 * 2.0 ^^ 10 );
324+
325+ // Can get the mantissa and exponent
326+ long e;
327+ assert (r.prod(e).approxEqual(0.8 ));
328+ assert (e == 10 );
276329}
277330
278331// / Product of vector
@@ -291,6 +344,10 @@ unittest
291344 auto r = [l, l, l, s, s, s, u, u, u].sliced;
292345
293346 assert (r.prod == reduce! " a * b" (1.0 , [u, u, u]));
347+
348+ long e;
349+ assert (r.prod(e).approxEqual(reduce! " a * b" (1.0 , [c, c, c])));
350+ assert (e == 30 );
294351}
295352
296353// / Product of matrix
@@ -312,6 +369,10 @@ unittest
312369 ].fuse;
313370
314371 assert (r.prod == reduce! " a * b" (1.0 , [u, u, u]));
372+
373+ long e;
374+ assert (r.prod(e) == reduce! " a * b" (1.0 , [c, c, c]));
375+ assert (e == 30 );
315376}
316377
317378// / Column prod of matrix
@@ -391,6 +452,10 @@ unittest
391452 static immutable result2 = [c, c, c];
392453
393454 assert (r.sliced.prod.approxEqual(reduce! " a * b" (1.0 , result1)));
455+
456+ long e;
457+ assert (r.sliced.prod(e).approxEqual(reduce! " a * b" (1.0 , result2)));
458+ assert (e == 30 );
394459}
395460
396461version (mir_test)
@@ -410,27 +475,29 @@ unittest
410475 static immutable result2 = [c, c, c];
411476
412477 assert (r.sliced.prod! double .approxEqual(reduce! " a * b" (1.0 , result1)));
478+
479+ long e;
480+ assert (r.sliced.prod! double (e).approxEqual(reduce! " a * b" (1.0 , result2)));
481+ assert (e == 30 );
413482}
414483
415484/+ +
416485Compute the sum of binary logarithms of the input range $(D r).
417486The error of this method is much smaller than with a naive sum of log2.
418487+/
419488Unqual! (DeepElementType! Range ) sumOfLog2s (Range )(Range r)
420- if (isNumeric ! (DeepElementType! Range ))
489+ if (isFloatingPoint ! (DeepElementType! Range ))
421490{
422- import mir.bignum.fp: fp_log2 ;
423- import core.lifetime : move ;
424- return prod ! 128 (move(r)).fp_log2 ! ( typeof ( return ) );
491+ long exp = 0 ;
492+ auto x = .prod(r, exp) ;
493+ return exp + log2(x );
425494}
426495
427496// /
428497version (mir_test)
429498@safe pure
430499unittest
431500{
432- import mir.format: text;
433- import mir.bignum.fp;
434501 alias isNaN = x => x != x;
435502
436503 assert (sumOfLog2s(new double [0 ]) == 0 );
@@ -450,7 +517,7 @@ Compute the sum of binary logarithms of the input range $(D r).
450517The error of this method is much smaller than with a naive sum of log.
451518+/
452519Unqual! (DeepElementType! Range ) sumOfLogs (Range )(Range r)
453- if (isNumeric ! (DeepElementType! Range ))
520+ if (isFloatingPoint ! (DeepElementType! Range ))
454521{
455522 import core.lifetime : move;
456523 import mir.math.constant: LN2 ;
0 commit comments