1 // Written in the D programming language.
2
3 /**
4 This is a submodule of $(MREF std, math).
5
6 It contains classical algebraic functions like `abs`, `sqrt`, and `poly`.
7
8 Copyright: Copyright The D Language Foundation 2000 - 2011.
9 License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
10 Authors: $(HTTP digitalmars.com, Walter Bright), Don Clugston,
11 Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger
12 Source: $(PHOBOSSRC std/math/algebraic.d)
13
14 Macros:
15 TABLE_SV = <table border="1" cellpadding="4" cellspacing="0">
16 <caption>Special Values</caption>
17 $0</table>
18 NAN = $(RED NAN)
19 POWER = $1<sup>$2</sup>
20 SUB = $1<sub>$2</sub>
21 PLUSMN = ±
22 INFIN = ∞
23 PLUSMNINF = ±∞
24 LT = <
25
26 */
27
28 module std.math.algebraic;
29
30 static import core.math;
31 static import core.stdc.math;
32 import std.traits : CommonType, isFloatingPoint, isIntegral, isSigned, Unqual;
33
34 /***********************************
35 * Calculates the absolute value of a number.
36 *
37 * Params:
38 * Num = (template parameter) type of number
39 * x = real number value
40 *
41 * Returns:
42 * The absolute value of the number. If floating-point or integral,
43 * the return type will be the same as the input.
44 *
45 * Limitations:
46 * Does not work correctly for signed intergal types and value `Num`.min.
47 */
48 auto abs(Num)(Num x) @nogc pure nothrow
49 if ((is(immutable Num == immutable short) || is(immutable Num == immutable byte)) ||
50 (is(typeof(Num.init >= 0)) && is(typeof(-Num.init))))
51 {
52 static if (isFloatingPoint!(Num))
53 return fabs(x);
54 else
55 {
56 static if (is(immutable Num == immutable short) || is(immutable Num == immutable byte))
57 return x >= 0 ? x : cast(Num) -int(x);
58 else
59 return x >= 0 ? x : -x;
60 }
61 }
62
63 /// ditto
64 @safe pure nothrow @nogc unittest
65 {
66 import std.math.traits : isIdentical, isNaN;
67
68 assert(isIdentical(abs(-0.0L), 0.0L));
69 assert(isNaN(abs(real.nan)));
70 assert(abs(-real.infinity) == real.infinity);
71 assert(abs(-56) == 56);
72 assert(abs(2321312L) == 2321312L);
73 }
74
75 @safe pure nothrow @nogc unittest
76 {
77 short s = -8;
78 byte b = -8;
79 assert(abs(s) == 8);
80 assert(abs(b) == 8);
81 immutable(byte) c = -8;
82 assert(abs(c) == 8);
83 }
84
85 @safe pure nothrow @nogc unittest
86 {
87 import std.meta : AliasSeq;
88 static foreach (T; AliasSeq!(float, double, real))
89 {{
90 T f = 3;
91 assert(abs(f) == f);
92 assert(abs(-f) == f);
93 }}
94 }
95
96 // see https://issues.dlang.org/show_bug.cgi?id=20205
97 // to avoid falling into the trap again
98 @safe pure nothrow @nogc unittest
99 {
100 assert(50 - abs(-100) == -50);
101 }
102
103 // https://issues.dlang.org/show_bug.cgi?id=19162
104 @safe unittest
105 {
Vector(T,int size)106 struct Vector(T, int size)
107 {
108 T x, y, z;
109 }
110
abs(T,int size)111 static auto abs(T, int size)(auto ref const Vector!(T, size) v)
112 {
113 return v;
114 }
115 Vector!(int, 3) v;
116 assert(abs(v) == v);
117 }
118
119 /*******************************
120 * Returns |x|
121 *
122 * $(TABLE_SV
123 * $(TR $(TH x) $(TH fabs(x)))
124 * $(TR $(TD $(PLUSMN)0.0) $(TD +0.0) )
125 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) )
126 * )
127 */
pragma(inline,true)128 pragma(inline, true)
129 real fabs(real x) @safe pure nothrow @nogc { return core.math.fabs(x); }
130
131 ///ditto
pragma(inline,true)132 pragma(inline, true)
133 double fabs(double x) @safe pure nothrow @nogc { return core.math.fabs(x); }
134
135 ///ditto
pragma(inline,true)136 pragma(inline, true)
137 float fabs(float x) @safe pure nothrow @nogc { return core.math.fabs(x); }
138
139 ///
140 @safe unittest
141 {
142 import std.math.traits : isIdentical;
143
144 assert(isIdentical(fabs(0.0f), 0.0f));
145 assert(isIdentical(fabs(-0.0f), 0.0f));
146 assert(fabs(-10.0f) == 10.0f);
147
148 assert(isIdentical(fabs(0.0), 0.0));
149 assert(isIdentical(fabs(-0.0), 0.0));
150 assert(fabs(-10.0) == 10.0);
151
152 assert(isIdentical(fabs(0.0L), 0.0L));
153 assert(isIdentical(fabs(-0.0L), 0.0L));
154 assert(fabs(-10.0L) == 10.0L);
155 }
156
157 @safe unittest
158 {
159 real function(real) pfabs = &fabs;
160 assert(pfabs != null);
161 }
162
163 @safe pure nothrow @nogc unittest
164 {
165 float f = fabs(-2.0f);
166 assert(f == 2);
167
168 double d = fabs(-2.0);
169 assert(d == 2);
170
171 real r = fabs(-2.0L);
172 assert(r == 2);
173 }
174
175 /***************************************
176 * Compute square root of x.
177 *
178 * $(TABLE_SV
179 * $(TR $(TH x) $(TH sqrt(x)) $(TH invalid?))
180 * $(TR $(TD -0.0) $(TD -0.0) $(TD no))
181 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD yes))
182 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no))
183 * )
184 */
pragma(inline,true)185 pragma(inline, true)
186 float sqrt(float x) @nogc @safe pure nothrow { return core.math.sqrt(x); }
187
188 /// ditto
pragma(inline,true)189 pragma(inline, true)
190 double sqrt(double x) @nogc @safe pure nothrow { return core.math.sqrt(x); }
191
192 /// ditto
pragma(inline,true)193 pragma(inline, true)
194 real sqrt(real x) @nogc @safe pure nothrow { return core.math.sqrt(x); }
195
196 ///
197 @safe pure nothrow @nogc unittest
198 {
199 import std.math.operations : feqrel;
200 import std.math.traits : isNaN;
201
202 assert(sqrt(2.0).feqrel(1.4142) > 16);
203 assert(sqrt(9.0).feqrel(3.0) > 16);
204
205 assert(isNaN(sqrt(-1.0f)));
206 assert(isNaN(sqrt(-1.0)));
207 assert(isNaN(sqrt(-1.0L)));
208 }
209
210 @safe unittest
211 {
212 // https://issues.dlang.org/show_bug.cgi?id=5305
213 float function(float) psqrtf = &sqrt;
214 assert(psqrtf != null);
215 double function(double) psqrtd = &sqrt;
216 assert(psqrtd != null);
217 real function(real) psqrtr = &sqrt;
218 assert(psqrtr != null);
219
220 //ctfe
221 enum ZX80 = sqrt(7.0f);
222 enum ZX81 = sqrt(7.0);
223 enum ZX82 = sqrt(7.0L);
224 }
225
226 @safe pure nothrow @nogc unittest
227 {
228 float f = sqrt(2.0f);
229 assert(fabs(f * f - 2.0f) < .00001);
230
231 double d = sqrt(2.0);
232 assert(fabs(d * d - 2.0) < .00001);
233
234 real r = sqrt(2.0L);
235 assert(fabs(r * r - 2.0) < .00001);
236 }
237
238 /***************
239 * Calculates the cube root of x.
240 *
241 * $(TABLE_SV
242 * $(TR $(TH $(I x)) $(TH cbrt(x)) $(TH invalid?))
243 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) )
244 * $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes) )
245 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no) )
246 * )
247 */
cbrt(real x)248 real cbrt(real x) @trusted nothrow @nogc
249 {
250 version (CRuntime_Microsoft)
251 {
252 import std.math.traits : copysign;
253 import std.math.exponential : exp2;
254
255 version (INLINE_YL2X)
256 return copysign(exp2(core.math.yl2x(fabs(x), 1.0L/3.0L)), x);
257 else
258 return core.stdc.math.cbrtl(x);
259 }
260 else
261 return core.stdc.math.cbrtl(x);
262 }
263
264 ///
265 @safe unittest
266 {
267 import std.math.operations : feqrel;
268
269 assert(cbrt(1.0).feqrel(1.0) > 16);
270 assert(cbrt(27.0).feqrel(3.0) > 16);
271 assert(cbrt(15.625).feqrel(2.5) > 16);
272 }
273
274 /***********************************************************************
275 * Calculates the length of the
276 * hypotenuse of a right-angled triangle with sides of length x and y.
277 * The hypotenuse is the value of the square root of
278 * the sums of the squares of x and y:
279 *
280 * sqrt($(POWER x, 2) + $(POWER y, 2))
281 *
282 * Note that hypot(x, y), hypot(y, x) and
283 * hypot(x, -y) are equivalent.
284 *
285 * $(TABLE_SV
286 * $(TR $(TH x) $(TH y) $(TH hypot(x, y)) $(TH invalid?))
287 * $(TR $(TD x) $(TD $(PLUSMN)0.0) $(TD |x|) $(TD no))
288 * $(TR $(TD $(PLUSMNINF)) $(TD y) $(TD +$(INFIN)) $(TD no))
289 * $(TR $(TD $(PLUSMNINF)) $(TD $(NAN)) $(TD +$(INFIN)) $(TD no))
290 * )
291 */
292 T hypot(T)(const T x, const T y) @safe pure nothrow @nogc
293 if (isFloatingPoint!T)
294 {
295 // Scale x and y to avoid underflow and overflow.
296 // If one is huge and the other tiny, return the larger.
297 // If both are huge, avoid overflow by scaling by 2^^-N.
298 // If both are tiny, avoid underflow by scaling by 2^^N.
299 import core.math : fabs, sqrt;
300 import std.math : floatTraits, RealFormat;
301
302 alias F = floatTraits!T;
303
304 T u = fabs(x);
305 T v = fabs(y);
306 if (!(u >= v)) // check for NaN as well.
307 {
308 v = u;
309 u = fabs(y);
310 if (u == T.infinity) return u; // hypot(inf, nan) == inf
311 if (v == T.infinity) return v; // hypot(nan, inf) == inf
312 }
313
314 static if (F.realFormat == RealFormat.ieeeSingle)
315 {
316 enum SQRTMIN = 0x1p-60f;
317 enum SQRTMAX = 0x1p+60f;
318 enum SCALE_UNDERFLOW = 0x1p+90f;
319 enum SCALE_OVERFLOW = 0x1p-90f;
320 }
321 else static if (F.realFormat == RealFormat.ieeeDouble ||
322 F.realFormat == RealFormat.ieeeExtended53 ||
323 F.realFormat == RealFormat.ibmExtended)
324 {
325 enum SQRTMIN = 0x1p-450L;
326 enum SQRTMAX = 0x1p+500L;
327 enum SCALE_UNDERFLOW = 0x1p+600L;
328 enum SCALE_OVERFLOW = 0x1p-600L;
329 }
330 else static if (F.realFormat == RealFormat.ieeeExtended ||
331 F.realFormat == RealFormat.ieeeQuadruple)
332 {
333 enum SQRTMIN = 0x1p-8000L;
334 enum SQRTMAX = 0x1p+8000L;
335 enum SCALE_UNDERFLOW = 0x1p+10000L;
336 enum SCALE_OVERFLOW = 0x1p-10000L;
337 }
338 else
339 assert(0, "hypot not implemented");
340
341 // Now u >= v, or else one is NaN.
342 T ratio = 1.0;
343 if (v >= SQRTMAX)
344 {
345 // hypot(huge, huge) -- avoid overflow
346 ratio = SCALE_UNDERFLOW;
347 u *= SCALE_OVERFLOW;
348 v *= SCALE_OVERFLOW;
349 }
350 else if (u <= SQRTMIN)
351 {
352 // hypot (tiny, tiny) -- avoid underflow
353 // This is only necessary to avoid setting the underflow
354 // flag.
355 ratio = SCALE_OVERFLOW;
356 u *= SCALE_UNDERFLOW;
357 v *= SCALE_UNDERFLOW;
358 }
359
360 if (u * T.epsilon > v)
361 {
362 // hypot (huge, tiny) = huge
363 return u;
364 }
365
366 // both are in the normal range
367 return ratio * sqrt(u*u + v*v);
368 }
369
370 ///
371 @safe unittest
372 {
373 import std.math.operations : feqrel;
374
375 assert(hypot(1.0, 1.0).feqrel(1.4142) > 16);
376 assert(hypot(3.0, 4.0).feqrel(5.0) > 16);
377 assert(hypot(real.infinity, 1.0L) == real.infinity);
378 assert(hypot(real.infinity, real.nan) == real.infinity);
379 }
380
381 @safe unittest
382 {
383 import std.math.operations : feqrel;
384
385 assert(hypot(1.0f, 1.0f).feqrel(1.4142f) > 16);
386 assert(hypot(3.0f, 4.0f).feqrel(5.0f) > 16);
387 assert(hypot(float.infinity, 1.0f) == float.infinity);
388 assert(hypot(float.infinity, float.nan) == float.infinity);
389
390 assert(hypot(1.0L, 1.0L).feqrel(1.4142L) > 16);
391 assert(hypot(3.0L, 4.0L).feqrel(5.0L) > 16);
392 assert(hypot(double.infinity, 1.0) == double.infinity);
393 assert(hypot(double.infinity, double.nan) == double.infinity);
394 }
395
396 @safe unittest
397 {
398 import std.math.operations : feqrel;
399 import std.math.traits : isIdentical;
400 import std.meta : AliasSeq;
401
402 static foreach (T; AliasSeq!(float, double, real))
403 {{
404 static T[3][] vals = // x,y,hypot
405 [
406 [ 0.0, 0.0, 0.0],
407 [ 0.0, -0.0, 0.0],
408 [ -0.0, -0.0, 0.0],
409 [ 3.0, 4.0, 5.0],
410 [ -300, -400, 500],
411 [0.0, 7.0, 7.0],
412 [9.0, 9*T.epsilon, 9.0],
413 [88/(64*sqrt(T.min_normal)), 105/(64*sqrt(T.min_normal)), 137/(64*sqrt(T.min_normal))],
414 [88/(128*sqrt(T.min_normal)), 105/(128*sqrt(T.min_normal)), 137/(128*sqrt(T.min_normal))],
415 [3*T.min_normal*T.epsilon, 4*T.min_normal*T.epsilon, 5*T.min_normal*T.epsilon],
416 [ T.min_normal, T.min_normal, sqrt(2.0L)*T.min_normal],
417 [ T.max/sqrt(2.0L), T.max/sqrt(2.0L), T.max],
418 [ T.infinity, T.nan, T.infinity],
419 [ T.nan, T.infinity, T.infinity],
420 [ T.nan, T.nan, T.nan],
421 [ T.nan, T.max, T.nan],
422 [ T.max, T.nan, T.nan],
423 ];
424 for (int i = 0; i < vals.length; i++)
425 {
426 T x = vals[i][0];
427 T y = vals[i][1];
428 T z = vals[i][2];
429 T h = hypot(x, y);
430 assert(isIdentical(z,h) || feqrel(z, h) >= T.mant_dig - 1);
431 }
432 }}
433 }
434
435 /***********************************************************************
436 * Calculates the distance of the point (x, y, z) from the origin (0, 0, 0)
437 * in three-dimensional space.
438 * The distance is the value of the square root of the sums of the squares
439 * of x, y, and z:
440 *
441 * sqrt($(POWER x, 2) + $(POWER y, 2) + $(POWER z, 2))
442 *
443 * Note that the distance between two points (x1, y1, z1) and (x2, y2, z2)
444 * in three-dimensional space can be calculated as hypot(x2-x1, y2-y1, z2-z1).
445 *
446 * Params:
447 * x = floating point value
448 * y = floating point value
449 * z = floating point value
450 *
451 * Returns:
452 * The square root of the sum of the squares of the given arguments.
453 */
454 T hypot(T)(const T x, const T y, const T z) @safe pure nothrow @nogc
455 if (isFloatingPoint!T)
456 {
457 import core.math : fabs, sqrt;
458 import std.math.operations : fmax;
459 const absx = fabs(x);
460 const absy = fabs(y);
461 const absz = fabs(z);
462
463 // Scale all parameters to avoid overflow.
464 const ratio = fmax(absx, fmax(absy, absz));
465 if (ratio == 0.0)
466 return ratio;
467
468 return ratio * sqrt((absx / ratio) * (absx / ratio)
469 + (absy / ratio) * (absy / ratio)
470 + (absz / ratio) * (absz / ratio));
471 }
472
473 ///
474 @safe unittest
475 {
476 import std.math.operations : isClose;
477
478 assert(isClose(hypot(1.0, 2.0, 2.0), 3.0));
479 assert(isClose(hypot(2.0, 3.0, 6.0), 7.0));
480 assert(isClose(hypot(1.0, 4.0, 8.0), 9.0));
481 }
482
483 @safe unittest
484 {
485 import std.meta : AliasSeq;
486 import std.math.traits : isIdentical;
487 import std.math.operations : isClose;
488 static foreach (T; AliasSeq!(float, double, real))
489 {{
490 static T[4][] vals = [
491 [ 0.0L, 0.0L, 0.0L, 0.0L ],
492 [ 0.0L, 1.0L, 1.0L, sqrt(2.0L) ],
493 [ 1.0L, 1.0L, 1.0L, sqrt(3.0L) ],
494 [ 1.0L, 2.0L, 2.0L, 3.0L ],
495 [ 2.0L, 3.0L, 6.0L, 7.0L ],
496 [ 1.0L, 4.0L, 8.0L, 9.0L ],
497 [ 4.0L, 4.0L, 7.0L, 9.0L ],
498 [ 12.0L, 16.0L, 21.0L, 29.0L ],
499 [ 1e+8L, 1.0L, 1e-8L, 1e+8L ],
500 [ 1.0L, 1e+8L, 1e-8L, 1e+8L ],
501 [ 1e-8L, 1.0L, 1e+8L, 1e+8L ],
502 [ 1e-2L, 1e-4L, 1e-4L, 0.010000999950004999375L ],
503 [ 2147483647.0L, 2147483647.0L, 2147483647.0L, 3719550785.027307813987L ]
504 ];
505 for (int i = 0; i < vals.length; i++)
506 {
507 T x = vals[i][0];
508 T y = vals[i][1];
509 T z = vals[i][2];
510 T r = vals[i][3];
511 T a = hypot(x, y, z);
512 assert(isIdentical(r, a) || isClose(r, a));
513 }
514 }}
515 }
516
517 /***********************************
518 * Evaluate polynomial A(x) = $(SUB a, 0) + $(SUB a, 1)x + $(SUB a, 2)$(POWER x,2) +
519 * $(SUB a,3)$(POWER x,3); ...
520 *
521 * Uses Horner's rule A(x) = $(SUB a, 0) + x($(SUB a, 1) + x($(SUB a, 2) +
522 * x($(SUB a, 3) + ...)))
523 * Params:
524 * x = the value to evaluate.
525 * A = array of coefficients $(SUB a, 0), $(SUB a, 1), etc.
526 */
527 Unqual!(CommonType!(T1, T2)) poly(T1, T2)(T1 x, in T2[] A) @trusted pure nothrow @nogc
528 if (isFloatingPoint!T1 && isFloatingPoint!T2)
529 in
530 {
531 assert(A.length > 0);
532 }
533 do
534 {
535 static if (is(immutable T2 == immutable real))
536 {
537 return polyImpl(x, A);
538 }
539 else
540 {
541 return polyImplBase(x, A);
542 }
543 }
544
545 /// ditto
546 Unqual!(CommonType!(T1, T2)) poly(T1, T2, int N)(T1 x, ref const T2[N] A) @safe pure nothrow @nogc
547 if (isFloatingPoint!T1 && isFloatingPoint!T2 && N > 0 && N <= 10)
548 {
549 // statically unrolled version for up to 10 coefficients
550 typeof(return) r = A[N - 1];
551 static foreach (i; 1 .. N)
552 {
553 r *= x;
554 r += A[N - 1 - i];
555 }
556 return r;
557 }
558
559 ///
560 @safe nothrow @nogc unittest
561 {
562 real x = 3.1L;
563 static real[] pp = [56.1L, 32.7L, 6];
564
565 assert(poly(x, pp) == (56.1L + (32.7L + 6.0L * x) * x));
566 }
567
568 @safe nothrow @nogc unittest
569 {
570 double x = 3.1;
571 static double[] pp = [56.1, 32.7, 6];
572 double y = x;
573 y *= 6.0;
574 y += 32.7;
575 y *= x;
576 y += 56.1;
577 assert(poly(x, pp) == y);
578 }
579
580 @safe unittest
581 {
582 static assert(poly(3.0, [1.0, 2.0, 3.0]) == 34);
583 }
584
585 private Unqual!(CommonType!(T1, T2)) polyImplBase(T1, T2)(T1 x, in T2[] A) @trusted pure nothrow @nogc
586 if (isFloatingPoint!T1 && isFloatingPoint!T2)
587 {
588 ptrdiff_t i = A.length - 1;
589 typeof(return) r = A[i];
590 while (--i >= 0)
591 {
592 r *= x;
593 r += A[i];
594 }
595 return r;
596 }
597
598 version (linux) version = GenericPosixVersion;
599 else version (FreeBSD) version = GenericPosixVersion;
600 else version (OpenBSD) version = GenericPosixVersion;
601 else version (Solaris) version = GenericPosixVersion;
602 else version (DragonFlyBSD) version = GenericPosixVersion;
603
polyImpl(real x,in real[]A)604 private real polyImpl(real x, in real[] A) @trusted pure nothrow @nogc
605 {
606 version (D_InlineAsm_X86)
607 {
608 if (__ctfe)
609 {
610 return polyImplBase(x, A);
611 }
612 version (Windows)
613 {
614 // BUG: This code assumes a frame pointer in EBP.
615 asm pure nothrow @nogc // assembler by W. Bright
616 {
617 // EDX = (A.length - 1) * real.sizeof
618 mov ECX,A[EBP] ; // ECX = A.length
619 dec ECX ;
620 lea EDX,[ECX][ECX*8] ;
621 add EDX,ECX ;
622 add EDX,A+4[EBP] ;
623 fld real ptr [EDX] ; // ST0 = coeff[ECX]
624 jecxz return_ST ;
625 fld x[EBP] ; // ST0 = x
626 fxch ST(1) ; // ST1 = x, ST0 = r
627 align 4 ;
628 L2: fmul ST,ST(1) ; // r *= x
629 fld real ptr -10[EDX] ;
630 sub EDX,10 ; // deg--
631 faddp ST(1),ST ;
632 dec ECX ;
633 jne L2 ;
634 fxch ST(1) ; // ST1 = r, ST0 = x
635 fstp ST(0) ; // dump x
636 align 4 ;
637 return_ST: ;
638 }
639 }
640 else version (GenericPosixVersion)
641 {
642 asm pure nothrow @nogc // assembler by W. Bright
643 {
644 // EDX = (A.length - 1) * real.sizeof
645 mov ECX,A[EBP] ; // ECX = A.length
646 dec ECX ;
647 lea EDX,[ECX*8] ;
648 lea EDX,[EDX][ECX*4] ;
649 add EDX,A+4[EBP] ;
650 fld real ptr [EDX] ; // ST0 = coeff[ECX]
651 jecxz return_ST ;
652 fld x[EBP] ; // ST0 = x
653 fxch ST(1) ; // ST1 = x, ST0 = r
654 align 4 ;
655 L2: fmul ST,ST(1) ; // r *= x
656 fld real ptr -12[EDX] ;
657 sub EDX,12 ; // deg--
658 faddp ST(1),ST ;
659 dec ECX ;
660 jne L2 ;
661 fxch ST(1) ; // ST1 = r, ST0 = x
662 fstp ST(0) ; // dump x
663 align 4 ;
664 return_ST: ;
665 }
666 }
667 else version (OSX)
668 {
669 asm pure nothrow @nogc // assembler by W. Bright
670 {
671 // EDX = (A.length - 1) * real.sizeof
672 mov ECX,A[EBP] ; // ECX = A.length
673 dec ECX ;
674 lea EDX,[ECX*8] ;
675 add EDX,EDX ;
676 add EDX,A+4[EBP] ;
677 fld real ptr [EDX] ; // ST0 = coeff[ECX]
678 jecxz return_ST ;
679 fld x[EBP] ; // ST0 = x
680 fxch ST(1) ; // ST1 = x, ST0 = r
681 align 4 ;
682 L2: fmul ST,ST(1) ; // r *= x
683 fld real ptr -16[EDX] ;
684 sub EDX,16 ; // deg--
685 faddp ST(1),ST ;
686 dec ECX ;
687 jne L2 ;
688 fxch ST(1) ; // ST1 = r, ST0 = x
689 fstp ST(0) ; // dump x
690 align 4 ;
691 return_ST: ;
692 }
693 }
694 else
695 {
696 static assert(0);
697 }
698 }
699 else
700 {
701 return polyImplBase(x, A);
702 }
703 }
704
705 /**
706 * Gives the next power of two after `val`. `T` can be any built-in
707 * numerical type.
708 *
709 * If the operation would lead to an over/underflow, this function will
710 * return `0`.
711 *
712 * Params:
713 * val = any number
714 *
715 * Returns:
716 * the next power of two after `val`
717 */
718 T nextPow2(T)(const T val)
719 if (isIntegral!T)
720 {
721 return powIntegralImpl!(PowType.ceil)(val);
722 }
723
724 /// ditto
725 T nextPow2(T)(const T val)
726 if (isFloatingPoint!T)
727 {
728 return powFloatingPointImpl!(PowType.ceil)(val);
729 }
730
731 ///
732 @safe @nogc pure nothrow unittest
733 {
734 assert(nextPow2(2) == 4);
735 assert(nextPow2(10) == 16);
736 assert(nextPow2(4000) == 4096);
737
738 assert(nextPow2(-2) == -4);
739 assert(nextPow2(-10) == -16);
740
741 assert(nextPow2(uint.max) == 0);
742 assert(nextPow2(uint.min) == 0);
743 assert(nextPow2(size_t.max) == 0);
744 assert(nextPow2(size_t.min) == 0);
745
746 assert(nextPow2(int.max) == 0);
747 assert(nextPow2(int.min) == 0);
748 assert(nextPow2(long.max) == 0);
749 assert(nextPow2(long.min) == 0);
750 }
751
752 ///
753 @safe @nogc pure nothrow unittest
754 {
755 assert(nextPow2(2.1) == 4.0);
756 assert(nextPow2(-2.0) == -4.0);
757 assert(nextPow2(0.25) == 0.5);
758 assert(nextPow2(-4.0) == -8.0);
759
760 assert(nextPow2(double.max) == 0.0);
761 assert(nextPow2(double.infinity) == double.infinity);
762 }
763
764 @safe @nogc pure nothrow unittest
765 {
766 assert(nextPow2(ubyte(2)) == 4);
767 assert(nextPow2(ubyte(10)) == 16);
768
769 assert(nextPow2(byte(2)) == 4);
770 assert(nextPow2(byte(10)) == 16);
771
772 assert(nextPow2(short(2)) == 4);
773 assert(nextPow2(short(10)) == 16);
774 assert(nextPow2(short(4000)) == 4096);
775
776 assert(nextPow2(ushort(2)) == 4);
777 assert(nextPow2(ushort(10)) == 16);
778 assert(nextPow2(ushort(4000)) == 4096);
779 }
780
781 @safe @nogc pure nothrow unittest
782 {
783 foreach (ulong i; 1 .. 62)
784 {
785 assert(nextPow2(1UL << i) == 2UL << i);
786 assert(nextPow2((1UL << i) - 1) == 1UL << i);
787 assert(nextPow2((1UL << i) + 1) == 2UL << i);
788 assert(nextPow2((1UL << i) + (1UL<<(i-1))) == 2UL << i);
789 }
790 }
791
792 @safe @nogc pure nothrow unittest
793 {
794 import std.math.traits : isNaN;
795 import std.meta : AliasSeq;
796
797 static foreach (T; AliasSeq!(float, double, real))
798 {{
799 enum T subNormal = T.min_normal / 2;
800
801 static if (subNormal) assert(nextPow2(subNormal) == T.min_normal);
802
803 assert(nextPow2(T(0.0)) == 0.0);
804
805 assert(nextPow2(T(2.0)) == 4.0);
806 assert(nextPow2(T(2.1)) == 4.0);
807 assert(nextPow2(T(3.1)) == 4.0);
808 assert(nextPow2(T(4.0)) == 8.0);
809 assert(nextPow2(T(0.25)) == 0.5);
810
811 assert(nextPow2(T(-2.0)) == -4.0);
812 assert(nextPow2(T(-2.1)) == -4.0);
813 assert(nextPow2(T(-3.1)) == -4.0);
814 assert(nextPow2(T(-4.0)) == -8.0);
815 assert(nextPow2(T(-0.25)) == -0.5);
816
817 assert(nextPow2(T.max) == 0);
818 assert(nextPow2(-T.max) == 0);
819
820 assert(nextPow2(T.infinity) == T.infinity);
821 assert(nextPow2(T.init).isNaN);
822 }}
823 }
824
825 // https://issues.dlang.org/show_bug.cgi?id=15973
826 @safe @nogc pure nothrow unittest
827 {
828 assert(nextPow2(uint.max / 2) == uint.max / 2 + 1);
829 assert(nextPow2(uint.max / 2 + 2) == 0);
830 assert(nextPow2(int.max / 2) == int.max / 2 + 1);
831 assert(nextPow2(int.max / 2 + 2) == 0);
832 assert(nextPow2(int.min + 1) == int.min);
833 }
834
835 /**
836 * Gives the last power of two before `val`. $(T) can be any built-in
837 * numerical type.
838 *
839 * Params:
840 * val = any number
841 *
842 * Returns:
843 * the last power of two before `val`
844 */
845 T truncPow2(T)(const T val)
846 if (isIntegral!T)
847 {
848 return powIntegralImpl!(PowType.floor)(val);
849 }
850
851 /// ditto
852 T truncPow2(T)(const T val)
853 if (isFloatingPoint!T)
854 {
855 return powFloatingPointImpl!(PowType.floor)(val);
856 }
857
858 ///
859 @safe @nogc pure nothrow unittest
860 {
861 assert(truncPow2(3) == 2);
862 assert(truncPow2(4) == 4);
863 assert(truncPow2(10) == 8);
864 assert(truncPow2(4000) == 2048);
865
866 assert(truncPow2(-5) == -4);
867 assert(truncPow2(-20) == -16);
868
869 assert(truncPow2(uint.max) == int.max + 1);
870 assert(truncPow2(uint.min) == 0);
871 assert(truncPow2(ulong.max) == long.max + 1);
872 assert(truncPow2(ulong.min) == 0);
873
874 assert(truncPow2(int.max) == (int.max / 2) + 1);
875 assert(truncPow2(int.min) == int.min);
876 assert(truncPow2(long.max) == (long.max / 2) + 1);
877 assert(truncPow2(long.min) == long.min);
878 }
879
880 ///
881 @safe @nogc pure nothrow unittest
882 {
883 assert(truncPow2(2.1) == 2.0);
884 assert(truncPow2(7.0) == 4.0);
885 assert(truncPow2(-1.9) == -1.0);
886 assert(truncPow2(0.24) == 0.125);
887 assert(truncPow2(-7.0) == -4.0);
888
889 assert(truncPow2(double.infinity) == double.infinity);
890 }
891
892 @safe @nogc pure nothrow unittest
893 {
894 assert(truncPow2(ubyte(3)) == 2);
895 assert(truncPow2(ubyte(4)) == 4);
896 assert(truncPow2(ubyte(10)) == 8);
897
898 assert(truncPow2(byte(3)) == 2);
899 assert(truncPow2(byte(4)) == 4);
900 assert(truncPow2(byte(10)) == 8);
901
902 assert(truncPow2(ushort(3)) == 2);
903 assert(truncPow2(ushort(4)) == 4);
904 assert(truncPow2(ushort(10)) == 8);
905 assert(truncPow2(ushort(4000)) == 2048);
906
907 assert(truncPow2(short(3)) == 2);
908 assert(truncPow2(short(4)) == 4);
909 assert(truncPow2(short(10)) == 8);
910 assert(truncPow2(short(4000)) == 2048);
911 }
912
913 @safe @nogc pure nothrow unittest
914 {
915 foreach (ulong i; 1 .. 62)
916 {
917 assert(truncPow2(2UL << i) == 2UL << i);
918 assert(truncPow2((2UL << i) + 1) == 2UL << i);
919 assert(truncPow2((2UL << i) - 1) == 1UL << i);
920 assert(truncPow2((2UL << i) - (2UL<<(i-1))) == 1UL << i);
921 }
922 }
923
924 @safe @nogc pure nothrow unittest
925 {
926 import std.math.traits : isNaN;
927 import std.meta : AliasSeq;
928
929 static foreach (T; AliasSeq!(float, double, real))
930 {
931 assert(truncPow2(T(0.0)) == 0.0);
932
933 assert(truncPow2(T(4.0)) == 4.0);
934 assert(truncPow2(T(2.1)) == 2.0);
935 assert(truncPow2(T(3.5)) == 2.0);
936 assert(truncPow2(T(7.0)) == 4.0);
937 assert(truncPow2(T(0.24)) == 0.125);
938
939 assert(truncPow2(T(-2.0)) == -2.0);
940 assert(truncPow2(T(-2.1)) == -2.0);
941 assert(truncPow2(T(-3.1)) == -2.0);
942 assert(truncPow2(T(-7.0)) == -4.0);
943 assert(truncPow2(T(-0.24)) == -0.125);
944
945 assert(truncPow2(T.infinity) == T.infinity);
946 assert(truncPow2(T.init).isNaN);
947 }
948 }
949
950 private enum PowType
951 {
952 floor,
953 ceil
954 }
955
pragma(inline,true)956 pragma(inline, true)
957 private T powIntegralImpl(PowType type, T)(T val)
958 {
959 import core.bitop : bsr;
960
961 if (val == 0 || (type == PowType.ceil && (val > T.max / 2 || val == T.min)))
962 return 0;
963 else
964 {
965 static if (isSigned!T)
966 return cast(Unqual!T) (val < 0 ? -(T(1) << bsr(0 - val) + type) : T(1) << bsr(val) + type);
967 else
968 return cast(Unqual!T) (T(1) << bsr(val) + type);
969 }
970 }
971
powFloatingPointImpl(PowType type,T)972 private T powFloatingPointImpl(PowType type, T)(T x)
973 {
974 import std.math.traits : copysign, isFinite;
975 import std.math.exponential : frexp;
976
977 if (!x.isFinite)
978 return x;
979
980 if (!x)
981 return x;
982
983 int exp;
984 auto y = frexp(x, exp);
985
986 static if (type == PowType.ceil)
987 y = core.math.ldexp(cast(T) 0.5, exp + 1);
988 else
989 y = core.math.ldexp(cast(T) 0.5, exp);
990
991 if (!y.isFinite)
992 return cast(T) 0.0;
993
994 y = copysign(y, x);
995
996 return y;
997 }
998