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