?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fortran: Fix PR 47485.
[gcc.git] / libphobos / src / std / math / exponential.d
blob7a72f428cd72a61285f1daec3f4e0a34f1272814
1 // Written in the D programming language.
2
3 /**
4 This is a submodule of $(MREF std, math).
5
6 It contains several exponential and logarithm functions.
7
8 Copyright: Copyright The D Language Foundation 2000 - 2011.
9 D implementations of exp, expm1, exp2, log, log10, log1p, and log2
10 functions are based on the CEPHES math library, which is Copyright
11 (C) 2001 Stephen L. Moshier $(LT)steve@moshier.net$(GT) and are
12 incorporated herein by permission of the author. The author reserves
13 the right to distribute this material elsewhere under different
14 copying permissions. These modifications are distributed here under
15 the following terms:
16 License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
17 Authors: $(HTTP digitalmars.com, Walter Bright), Don Clugston,
18 Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger
19 Source: $(PHOBOSSRC std/math/exponential.d)
20
21 Macros:
22 TABLE_SV = <table border="1" cellpadding="4" cellspacing="0">
23 <caption>Special Values</caption>
24 $0</table>
25 NAN = $(RED NAN)
26 PLUSMN = &plusmn;
27 INFIN = &infin;
28 PLUSMNINF = &plusmn;&infin;
29 LT = &lt;
30 GT = &gt;
31 */
32
33 module std.math.exponential;
34
35 import std.traits : isFloatingPoint, isIntegral, isSigned, isUnsigned, Largest, Unqual;
36
37 static import core.math;
38 static import core.stdc.math;
39
40 version (DigitalMars)
41 {
42 version (OSX) { } // macOS 13 (M1) has issues emulating instruction
43 else version = INLINE_YL2X; // x87 has opcodes for these
44 }
45
46 version (D_InlineAsm_X86) version = InlineAsm_X86_Any;
47 version (D_InlineAsm_X86_64) version = InlineAsm_X86_Any;
48
49 version (InlineAsm_X86_Any) version = InlineAsm_X87;
50 version (InlineAsm_X87)
51 {
52 static assert(real.mant_dig == 64);
53 version (CRuntime_Microsoft) version = InlineAsm_X87_MSVC;
54 }
55
56 version (D_HardFloat)
57 {
58 // FloatingPointControl.clearExceptions() depends on version IeeeFlagsSupport
59 version (IeeeFlagsSupport) version = FloatingPointControlSupport;
60 }
61
62 /**
63 * Compute the value of x $(SUPERSCRIPT n), where n is an integer
64 */
65 Unqual!F pow(F, G)(F x, G n) @nogc @trusted pure nothrow
66 if (isFloatingPoint!(F) && isIntegral!(G))
67 {
68 import std.traits : Unsigned;
69
70 real p = 1.0, v = void;
71 Unsigned!(Unqual!G) m = n;
72
73 if (n < 0)
74 {
75 if (n == -1) return 1 / x;
76
77 m = cast(typeof(m))(0 - n);
78 v = p / x;
79 }
80 else
81 {
82 switch (n)
83 {
84 case 0:
85 return 1.0;
86 case 1:
87 return x;
88 case 2:
89 return x * x;
90 default:
91 }
92
93 v = x;
94 }
95
96 while (1)
97 {
98 if (m & 1)
99 p *= v;
100 m >>= 1;
101 if (!m)
102 break;
103 v *= v;
104 }
105 return p;
106 }
107
108 ///
109 @safe pure nothrow @nogc unittest
110 {
111 import std.math.operations : feqrel;
112
113 assert(pow(2.0, 5) == 32.0);
114 assert(pow(1.5, 9).feqrel(38.4433) > 16);
115 assert(pow(real.nan, 2) is real.nan);
116 assert(pow(real.infinity, 2) == real.infinity);
117 }
118
119 @safe pure nothrow @nogc unittest
120 {
121 import std.math.operations : isClose, feqrel;
122
123 // Make sure it instantiates and works properly on immutable values and
124 // with various integer and float types.
125 immutable real x = 46;
126 immutable float xf = x;
127 immutable double xd = x;
128 immutable uint one = 1;
129 immutable ushort two = 2;
130 immutable ubyte three = 3;
131 immutable ulong eight = 8;
132
133 immutable int neg1 = -1;
134 immutable short neg2 = -2;
135 immutable byte neg3 = -3;
136 immutable long neg8 = -8;
137
138
139 assert(pow(x,0) == 1.0);
140 assert(pow(xd,one) == x);
141 assert(pow(xf,two) == x * x);
142 assert(pow(x,three) == x * x * x);
143 assert(pow(x,eight) == (x * x) * (x * x) * (x * x) * (x * x));
144
145 assert(pow(x, neg1) == 1 / x);
146
147 assert(isClose(pow(xd, neg2), cast(double) (1 / (x * x)), 1e-15));
148 assert(isClose(pow(xf, neg8), cast(float) (1 / ((x * x) * (x * x) * (x * x) * (x * x))), 1e-15));
149
150 assert(feqrel(pow(x, neg3), 1 / (x * x * x)) >= real.mant_dig - 1);
151 }
152
153 @safe @nogc nothrow unittest
154 {
155 import std.math.operations : isClose;
156
157 assert(isClose(pow(2.0L, 10L), 1024, 1e-18));
158 }
159
160 // https://issues.dlang.org/show_bug.cgi?id=21601
161 @safe @nogc nothrow pure unittest
162 {
163 // When reals are large enough the results of pow(b, e) can be
164 // calculated correctly, if b is of type float or double and e is
165 // not too large.
166 static if (real.mant_dig >= 64)
167 {
168 // expected result: 3.790e-42
169 assert(pow(-513645318757045764096.0f, -2) > 0.0);
170
171 // expected result: 3.763915357831797e-309
172 assert(pow(-1.6299717435255677e+154, -2) > 0.0);
173 }
174 }
175
176 @safe @nogc nothrow unittest
177 {
178 import std.math.operations : isClose;
179 import std.math.traits : isInfinity;
180
181 static float f1 = 19100.0f;
182 static float f2 = 0.000012f;
183
184 assert(isClose(pow(f1,9), 3.3829868e+38f));
185 assert(isInfinity(pow(f1,10)));
186 assert(pow(f2,9) > 0.0f);
187 assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal));
188
189 static double d1 = 21800.0;
190 static double d2 = 0.000012;
191
192 assert(isClose(pow(d1,71), 1.0725339442974e+308));
193 assert(isInfinity(pow(d1,72)));
194 assert(pow(d2,65) > 0.0f);
195 assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal));
196
197 static if (real.mant_dig == 64) // x87
198 {
199 static real r1 = 21950.0L;
200 static real r2 = 0.000011L;
201
202 assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L));
203 assert(isInfinity(pow(r1,1137)));
204 assert(pow(r2,998) > 0.0L);
205 assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal));
206 }
207 }
208
209 @safe @nogc nothrow pure unittest
210 {
211 import std.math.operations : isClose;
212
213 enum f1 = 19100.0f;
214 enum f2 = 0.000012f;
215
216 static assert(isClose(pow(f1,9), 3.3829868e+38f));
217 static assert(pow(f1,10) > float.max);
218 static assert(pow(f2,9) > 0.0f);
219 static assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal));
220
221 enum d1 = 21800.0;
222 enum d2 = 0.000012;
223
224 static assert(isClose(pow(d1,71), 1.0725339442974e+308));
225 static assert(pow(d1,72) > double.max);
226 static assert(pow(d2,65) > 0.0f);
227 static assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal));
228
229 static if (real.mant_dig == 64) // x87
230 {
231 enum r1 = 21950.0L;
232 enum r2 = 0.000011L;
233
234 static assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L));
235 static assert(pow(r1,1137) > real.max);
236 static assert(pow(r2,998) > 0.0L);
237 static assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal));
238 }
239 }
240
241 /**
242 * Compute the power of two integral numbers.
243 *
244 * Params:
245 * x = base
246 * n = exponent
247 *
248 * Returns:
249 * x raised to the power of n. If n is negative the result is 1 / pow(x, -n),
250 * which is calculated as integer division with remainder. This may result in
251 * a division by zero error.
252 *
253 * If both x and n are 0, the result is 1.
254 *
255 * Throws:
256 * If x is 0 and n is negative, the result is the same as the result of a
257 * division by zero.
258 */
259 typeof(Unqual!(F).init * Unqual!(G).init) pow(F, G)(F x, G n) @nogc @safe pure nothrow
260 if (isIntegral!(F) && isIntegral!(G))
261 {
262 import std.traits : isSigned;
263
264 typeof(return) p, v = void;
265 Unqual!G m = n;
266
267 static if (isSigned!(F))
268 {
269 if (x == -1) return cast(typeof(return)) (m & 1 ? -1 : 1);
270 }
271 static if (isSigned!(G))
272 {
273 if (x == 0 && m <= -1) return x / 0;
274 }
275 if (x == 1) return 1;
276 static if (isSigned!(G))
277 {
278 if (m < 0) return 0;
279 }
280
281 switch (m)
282 {
283 case 0:
284 p = 1;
285 break;
286
287 case 1:
288 p = x;
289 break;
290
291 case 2:
292 p = x * x;
293 break;
294
295 default:
296 v = x;
297 p = 1;
298 while (1)
299 {
300 if (m & 1)
301 p *= v;
302 m >>= 1;
303 if (!m)
304 break;
305 v *= v;
306 }
307 break;
308 }
309 return p;
310 }
311
312 ///
313 @safe pure nothrow @nogc unittest
314 {
315 assert(pow(2, 3) == 8);
316 assert(pow(3, 2) == 9);
317
318 assert(pow(2, 10) == 1_024);
319 assert(pow(2, 20) == 1_048_576);
320 assert(pow(2, 30) == 1_073_741_824);
321
322 assert(pow(0, 0) == 1);
323
324 assert(pow(1, -5) == 1);
325 assert(pow(1, -6) == 1);
326 assert(pow(-1, -5) == -1);
327 assert(pow(-1, -6) == 1);
328
329 assert(pow(-2, 5) == -32);
330 assert(pow(-2, -5) == 0);
331 assert(pow(cast(double) -2, -5) == -0.03125);
332 }
333
334 @safe pure nothrow @nogc unittest
335 {
336 immutable int one = 1;
337 immutable byte two = 2;
338 immutable ubyte three = 3;
339 immutable short four = 4;
340 immutable long ten = 10;
341
342 assert(pow(two, three) == 8);
343 assert(pow(two, ten) == 1024);
344 assert(pow(one, ten) == 1);
345 assert(pow(ten, four) == 10_000);
346 assert(pow(four, 10) == 1_048_576);
347 assert(pow(three, four) == 81);
348 }
349
350 // https://issues.dlang.org/show_bug.cgi?id=7006
351 @safe pure nothrow @nogc unittest
352 {
353 assert(pow(5, -1) == 0);
354 assert(pow(-5, -1) == 0);
355 assert(pow(5, -2) == 0);
356 assert(pow(-5, -2) == 0);
357 assert(pow(-1, int.min) == 1);
358 assert(pow(-2, int.min) == 0);
359
360 assert(pow(4294967290UL,2) == 18446744022169944100UL);
361 assert(pow(0,uint.max) == 0);
362 }
363
364 /**Computes integer to floating point powers.*/
365 real pow(I, F)(I x, F y) @nogc @trusted pure nothrow
366 if (isIntegral!I && isFloatingPoint!F)
367 {
368 return pow(cast(real) x, cast(Unqual!F) y);
369 }
370
371 ///
372 @safe pure nothrow @nogc unittest
373 {
374 assert(pow(2, 5.0) == 32.0);
375 assert(pow(7, 3.0) == 343.0);
376 assert(pow(2, real.nan) is real.nan);
377 assert(pow(2, real.infinity) == real.infinity);
378 }
379
380 /**
381 * Calculates x$(SUPERSCRIPT y).
382 *
383 * $(TABLE_SV
384 * $(TR $(TH x) $(TH y) $(TH pow(x, y))
385 * $(TH div 0) $(TH invalid?))
386 * $(TR $(TD anything) $(TD $(PLUSMN)0.0) $(TD 1.0)
387 * $(TD no) $(TD no) )
388 * $(TR $(TD |x| $(GT) 1) $(TD +$(INFIN)) $(TD +$(INFIN))
389 * $(TD no) $(TD no) )
390 * $(TR $(TD |x| $(LT) 1) $(TD +$(INFIN)) $(TD +0.0)
391 * $(TD no) $(TD no) )
392 * $(TR $(TD |x| $(GT) 1) $(TD -$(INFIN)) $(TD +0.0)
393 * $(TD no) $(TD no) )
394 * $(TR $(TD |x| $(LT) 1) $(TD -$(INFIN)) $(TD +$(INFIN))
395 * $(TD no) $(TD no) )
396 * $(TR $(TD +$(INFIN)) $(TD $(GT) 0.0) $(TD +$(INFIN))
397 * $(TD no) $(TD no) )
398 * $(TR $(TD +$(INFIN)) $(TD $(LT) 0.0) $(TD +0.0)
399 * $(TD no) $(TD no) )
400 * $(TR $(TD -$(INFIN)) $(TD odd integer $(GT) 0.0) $(TD -$(INFIN))
401 * $(TD no) $(TD no) )
402 * $(TR $(TD -$(INFIN)) $(TD $(GT) 0.0, not odd integer) $(TD +$(INFIN))
403 * $(TD no) $(TD no))
404 * $(TR $(TD -$(INFIN)) $(TD odd integer $(LT) 0.0) $(TD -0.0)
405 * $(TD no) $(TD no) )
406 * $(TR $(TD -$(INFIN)) $(TD $(LT) 0.0, not odd integer) $(TD +0.0)
407 * $(TD no) $(TD no) )
408 * $(TR $(TD $(PLUSMN)1.0) $(TD $(PLUSMN)$(INFIN)) $(TD -$(NAN))
409 * $(TD no) $(TD yes) )
410 * $(TR $(TD $(LT) 0.0) $(TD finite, nonintegral) $(TD $(NAN))
411 * $(TD no) $(TD yes))
412 * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(LT) 0.0) $(TD $(PLUSMNINF))
413 * $(TD yes) $(TD no) )
414 * $(TR $(TD $(PLUSMN)0.0) $(TD $(LT) 0.0, not odd integer) $(TD +$(INFIN))
415 * $(TD yes) $(TD no))
416 * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(GT) 0.0) $(TD $(PLUSMN)0.0)
417 * $(TD no) $(TD no) )
418 * $(TR $(TD $(PLUSMN)0.0) $(TD $(GT) 0.0, not odd integer) $(TD +0.0)
419 * $(TD no) $(TD no) )
420 * )
421 */
422 Unqual!(Largest!(F, G)) pow(F, G)(F x, G y) @nogc @trusted pure nothrow
423 if (isFloatingPoint!(F) && isFloatingPoint!(G))
424 {
425 return _powImpl(x, y);
426 }
427
428 ///
429 @safe pure nothrow @nogc unittest
430 {
431 import std.math.operations : isClose;
432
433 assert(isClose(pow(2.0, 3.0), 8.0));
434 assert(isClose(pow(1.5, 10.0), 57.6650390625));
435
436 // square root of 9
437 assert(isClose(pow(9.0, 0.5), 3.0));
438 // 10th root of 1024
439 assert(isClose(pow(1024.0, 0.1), 2.0));
440
441 assert(isClose(pow(-4.0, 3.0), -64.0));
442
443 // reciprocal of 4 ^^ 2
444 assert(isClose(pow(4.0, -2.0), 0.0625));
445 // reciprocal of (-2) ^^ 3
446 assert(isClose(pow(-2.0, -3.0), -0.125));
447
448 assert(isClose(pow(-2.5, 3.0), -15.625));
449 // reciprocal of 2.5 ^^ 3
450 assert(isClose(pow(2.5, -3.0), 0.064));
451 // reciprocal of (-2.5) ^^ 3
452 assert(isClose(pow(-2.5, -3.0), -0.064));
453
454 // reciprocal of square root of 4
455 assert(isClose(pow(4.0, -0.5), 0.5));
456
457 // per definition
458 assert(isClose(pow(0.0, 0.0), 1.0));
459 }
460
461 ///
462 @safe pure nothrow @nogc unittest
463 {
464 import std.math.operations : isClose;
465
466 // the result is a complex number
467 // which cannot be represented as floating point number
468 import std.math.traits : isNaN;
469 assert(isNaN(pow(-2.5, -1.5)));
470
471 // use the ^^-operator of std.complex instead
472 import std.complex : complex;
473 auto c1 = complex(-2.5, 0.0);
474 auto c2 = complex(-1.5, 0.0);
475 auto result = c1 ^^ c2;
476 // exact result apparently depends on `real` precision => increased tolerance
477 assert(isClose(result.re, -4.64705438e-17, 2e-4));
478 assert(isClose(result.im, 2.52982e-1, 2e-4));
479 }
480
481 @safe pure nothrow @nogc unittest
482 {
483 import std.math.traits : isNaN;
484
485 assert(pow(1.5, real.infinity) == real.infinity);
486 assert(pow(0.5, real.infinity) == 0.0);
487 assert(pow(1.5, -real.infinity) == 0.0);
488 assert(pow(0.5, -real.infinity) == real.infinity);
489 assert(pow(real.infinity, 1.0) == real.infinity);
490 assert(pow(real.infinity, -1.0) == 0.0);
491 assert(pow(real.infinity, real.infinity) == real.infinity);
492 assert(pow(-real.infinity, 1.0) == -real.infinity);
493 assert(pow(-real.infinity, 2.0) == real.infinity);
494 assert(pow(-real.infinity, -1.0) == -0.0);
495 assert(pow(-real.infinity, -2.0) == 0.0);
496 assert(isNaN(pow(1.0, real.infinity)));
497 assert(pow(0.0, -1.0) == real.infinity);
498 assert(pow(real.nan, 0.0) == 1.0);
499 assert(isNaN(pow(real.nan, 3.0)));
500 assert(isNaN(pow(3.0, real.nan)));
501 }
502
503 @safe @nogc nothrow unittest
504 {
505 import std.math.operations : isClose;
506
507 assert(isClose(pow(2.0L, 10.0L), 1024, 1e-18));
508 }
509
510 @safe pure nothrow @nogc unittest
511 {
512 import std.math.operations : isClose;
513 import std.math.traits : isIdentical, isNaN;
514 import std.math.constants : PI;
515
516 // Test all the special values. These unittests can be run on Windows
517 // by temporarily changing the version (linux) to version (all).
518 immutable float zero = 0;
519 immutable real one = 1;
520 immutable double two = 2;
521 immutable float three = 3;
522 immutable float fnan = float.nan;
523 immutable double dnan = double.nan;
524 immutable real rnan = real.nan;
525 immutable dinf = double.infinity;
526 immutable rninf = -real.infinity;
527
528 assert(pow(fnan, zero) == 1);
529 assert(pow(dnan, zero) == 1);
530 assert(pow(rnan, zero) == 1);
531
532 assert(pow(two, dinf) == double.infinity);
533 assert(isIdentical(pow(0.2f, dinf), +0.0));
534 assert(pow(0.99999999L, rninf) == real.infinity);
535 assert(isIdentical(pow(1.000000001, rninf), +0.0));
536 assert(pow(dinf, 0.001) == dinf);
537 assert(isIdentical(pow(dinf, -0.001), +0.0));
538 assert(pow(rninf, 3.0L) == rninf);
539 assert(pow(rninf, 2.0L) == real.infinity);
540 assert(isIdentical(pow(rninf, -3.0), -0.0));
541 assert(isIdentical(pow(rninf, -2.0), +0.0));
542
543 // @@@BUG@@@ somewhere
544 version (OSX) {} else assert(isNaN(pow(one, dinf)));
545 version (OSX) {} else assert(isNaN(pow(-one, dinf)));
546 assert(isNaN(pow(-0.2, PI)));
547 // boundary cases. Note that epsilon == 2^^-n for some n,
548 // so 1/epsilon == 2^^n is always even.
549 assert(pow(-1.0L, 1/real.epsilon - 1.0L) == -1.0L);
550 assert(pow(-1.0L, 1/real.epsilon) == 1.0L);
551 assert(isNaN(pow(-1.0L, 1/real.epsilon-0.5L)));
552 assert(isNaN(pow(-1.0L, -1/real.epsilon+0.5L)));
553
554 assert(pow(0.0, -3.0) == double.infinity);
555 assert(pow(-0.0, -3.0) == -double.infinity);
556 assert(pow(0.0, -PI) == double.infinity);
557 assert(pow(-0.0, -PI) == double.infinity);
558 assert(isIdentical(pow(0.0, 5.0), 0.0));
559 assert(isIdentical(pow(-0.0, 5.0), -0.0));
560 assert(isIdentical(pow(0.0, 6.0), 0.0));
561 assert(isIdentical(pow(-0.0, 6.0), 0.0));
562
563 // https://issues.dlang.org/show_bug.cgi?id=14786 fixed
564 immutable real maxOdd = pow(2.0L, real.mant_dig) - 1.0L;
565 assert(pow(-1.0L, maxOdd) == -1.0L);
566 assert(pow(-1.0L, -maxOdd) == -1.0L);
567 assert(pow(-1.0L, maxOdd + 1.0L) == 1.0L);
568 assert(pow(-1.0L, -maxOdd + 1.0L) == 1.0L);
569 assert(pow(-1.0L, maxOdd - 1.0L) == 1.0L);
570 assert(pow(-1.0L, -maxOdd - 1.0L) == 1.0L);
571
572 // Now, actual numbers.
573 assert(isClose(pow(two, three), 8.0));
574 assert(isClose(pow(two, -2.5), 0.1767766953));
575
576 // Test integer to float power.
577 immutable uint twoI = 2;
578 assert(isClose(pow(twoI, three), 8.0));
579 }
580
581 // https://issues.dlang.org/show_bug.cgi?id=20508
582 @safe pure nothrow @nogc unittest
583 {
584 import std.math.traits : isNaN;
585
586 assert(isNaN(pow(-double.infinity, 0.5)));
587
588 assert(isNaN(pow(-real.infinity, real.infinity)));
589 assert(isNaN(pow(-real.infinity, -real.infinity)));
590 assert(isNaN(pow(-real.infinity, 1.234)));
591 assert(isNaN(pow(-real.infinity, -0.751)));
592 assert(pow(-real.infinity, 0.0) == 1.0);
593 }
594
595 private real _powImpl(real x, real y) @safe @nogc pure nothrow
596 {
597 alias F = real;
598 import core.math : fabs, sqrt;
599 import std.math.traits : isInfinity, isNaN, signbit;
600
601 // Special cases.
602 if (isNaN(y))
603 return y;
604 if (isNaN(x) && y != 0.0)
605 return x;
606
607 // Even if x is NaN.
608 if (y == 0.0)
609 return 1.0;
610 if (y == 1.0)
611 return x;
612
613 if (isInfinity(y))
614 {
615 if (isInfinity(x))
616 {
617 if (!signbit(y) && !signbit(x))
618 return F.infinity;
619 else
620 return F.nan;
621 }
622 else if (fabs(x) > 1)
623 {
624 if (signbit(y))
625 return +0.0;
626 else
627 return F.infinity;
628 }
629 else if (fabs(x) == 1)
630 {
631 return F.nan;
632 }
633 else // < 1
634 {
635 if (signbit(y))
636 return F.infinity;
637 else
638 return +0.0;
639 }
640 }
641 if (isInfinity(x))
642 {
643 if (signbit(x))
644 {
645 long i = cast(long) y;
646 if (y > 0.0)
647 {
648 if (i == y && i & 1)
649 return -F.infinity;
650 else if (i == y)
651 return F.infinity;
652 else
653 return -F.nan;
654 }
655 else if (y < 0.0)
656 {
657 if (i == y && i & 1)
658 return -0.0;
659 else if (i == y)
660 return +0.0;
661 else
662 return F.nan;
663 }
664 }
665 else
666 {
667 if (y > 0.0)
668 return F.infinity;
669 else if (y < 0.0)
670 return +0.0;
671 }
672 }
673
674 if (x == 0.0)
675 {
676 if (signbit(x))
677 {
678 long i = cast(long) y;
679 if (y > 0.0)
680 {
681 if (i == y && i & 1)
682 return -0.0;
683 else
684 return +0.0;
685 }
686 else if (y < 0.0)
687 {
688 if (i == y && i & 1)
689 return -F.infinity;
690 else
691 return F.infinity;
692 }
693 }
694 else
695 {
696 if (y > 0.0)
697 return +0.0;
698 else if (y < 0.0)
699 return F.infinity;
700 }
701 }
702 if (x == 1.0)
703 return 1.0;
704
705 if (y >= F.max)
706 {
707 if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0))
708 return 0.0;
709 if (x > 1.0 || x < -1.0)
710 return F.infinity;
711 }
712 if (y <= -F.max)
713 {
714 if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0))
715 return F.infinity;
716 if (x > 1.0 || x < -1.0)
717 return 0.0;
718 }
719
720 if (x >= F.max)
721 {
722 if (y > 0.0)
723 return F.infinity;
724 else
725 return 0.0;
726 }
727 if (x <= -F.max)
728 {
729 long i = cast(long) y;
730 if (y > 0.0)
731 {
732 if (i == y && i & 1)
733 return -F.infinity;
734 else
735 return F.infinity;
736 }
737 else if (y < 0.0)
738 {
739 if (i == y && i & 1)
740 return -0.0;
741 else
742 return +0.0;
743 }
744 }
745
746 // Integer power of x.
747 long iy = cast(long) y;
748 if (iy == y && fabs(y) < 32_768.0)
749 return pow(x, iy);
750
751 real sign = 1.0;
752 if (x < 0)
753 {
754 // Result is real only if y is an integer
755 // Check for a non-zero fractional part
756 enum maxOdd = pow(2.0L, real.mant_dig) - 1.0L;
757 static if (maxOdd > ulong.max)
758 {
759 // Generic method, for any FP type
760 import std.math.rounding : floor;
761 if (floor(y) != y)
762 return sqrt(x); // Complex result -- create a NaN
763
764 const hy = 0.5 * y;
765 if (floor(hy) != hy)
766 sign = -1.0;
767 }
768 else
769 {
770 // Much faster, if ulong has enough precision
771 const absY = fabs(y);
772 if (absY <= maxOdd)
773 {
774 const uy = cast(ulong) absY;
775 if (uy != absY)
776 return sqrt(x); // Complex result -- create a NaN
777
778 if (uy & 1)
779 sign = -1.0;
780 }
781 }
782 x = -x;
783 }
784 version (INLINE_YL2X)
785 {
786 // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) )
787 // TODO: This is not accurate in practice. A fast and accurate
788 // (though complicated) method is described in:
789 // "An efficient rounding boundary test for pow(x, y)
790 // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007).
791 return sign * exp2( core.math.yl2x(x, y) );
792 }
793 else
794 {
795 // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) )
796 // TODO: This is not accurate in practice. A fast and accurate
797 // (though complicated) method is described in:
798 // "An efficient rounding boundary test for pow(x, y)
799 // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007).
800 auto w = exp2(y * log2(x));
801 return sign * w;
802 }
803 }
804
805 /** Computes the value of a positive integer `x`, raised to the power `n`, modulo `m`.
806 *
807 * Params:
808 * x = base
809 * n = exponent
810 * m = modulus
811 *
812 * Returns:
813 * `x` to the power `n`, modulo `m`.
814 * The return type is the largest of `x`'s and `m`'s type.
815 *
816 * The function requires that all values have unsigned types.
817 */
818 Unqual!(Largest!(F, H)) powmod(F, G, H)(F x, G n, H m)
819 if (isUnsigned!F && isUnsigned!G && isUnsigned!H)
820 {
821 import std.meta : AliasSeq;
822
823 alias T = Unqual!(Largest!(F, H));
824 static if (T.sizeof <= 4)
825 {
826 alias DoubleT = AliasSeq!(void, ushort, uint, void, ulong)[T.sizeof];
827 }
828
829 static T mulmod(T a, T b, T c)
830 {
831 static if (T.sizeof == 8)
832 {
833 static T addmod(T a, T b, T c)
834 {
835 b = c - b;
836 if (a >= b)
837 return a - b;
838 else
839 return c - b + a;
840 }
841
842 T result = 0, tmp;
843
844 b %= c;
845 while (a > 0)
846 {
847 if (a & 1)
848 result = addmod(result, b, c);
849
850 a >>= 1;
851 b = addmod(b, b, c);
852 }
853
854 return result;
855 }
856 else
857 {
858 DoubleT result = cast(DoubleT) (cast(DoubleT) a * cast(DoubleT) b);
859 return result % c;
860 }
861 }
862
863 T base = x, result = 1, modulus = m;
864 Unqual!G exponent = n;
865
866 while (exponent > 0)
867 {
868 if (exponent & 1)
869 result = mulmod(result, base, modulus);
870
871 base = mulmod(base, base, modulus);
872 exponent >>= 1;
873 }
874
875 return result;
876 }
877
878 ///
879 @safe pure nothrow @nogc unittest
880 {
881 assert(powmod(1U, 10U, 3U) == 1);
882 assert(powmod(3U, 2U, 6U) == 3);
883 assert(powmod(5U, 5U, 15U) == 5);
884 assert(powmod(2U, 3U, 5U) == 3);
885 assert(powmod(2U, 4U, 5U) == 1);
886 assert(powmod(2U, 5U, 5U) == 2);
887 }
888
889 @safe pure nothrow @nogc unittest
890 {
891 ulong a = 18446744073709551615u, b = 20u, c = 18446744073709551610u;
892 assert(powmod(a, b, c) == 95367431640625u);
893 a = 100; b = 7919; c = 18446744073709551557u;
894 assert(powmod(a, b, c) == 18223853583554725198u);
895 a = 117; b = 7919; c = 18446744073709551557u;
896 assert(powmod(a, b, c) == 11493139548346411394u);
897 a = 134; b = 7919; c = 18446744073709551557u;
898 assert(powmod(a, b, c) == 10979163786734356774u);
899 a = 151; b = 7919; c = 18446744073709551557u;
900 assert(powmod(a, b, c) == 7023018419737782840u);
901 a = 168; b = 7919; c = 18446744073709551557u;
902 assert(powmod(a, b, c) == 58082701842386811u);
903 a = 185; b = 7919; c = 18446744073709551557u;
904 assert(powmod(a, b, c) == 17423478386299876798u);
905 a = 202; b = 7919; c = 18446744073709551557u;
906 assert(powmod(a, b, c) == 5522733478579799075u);
907 a = 219; b = 7919; c = 18446744073709551557u;
908 assert(powmod(a, b, c) == 15230218982491623487u);
909 a = 236; b = 7919; c = 18446744073709551557u;
910 assert(powmod(a, b, c) == 5198328724976436000u);
911
912 a = 0; b = 7919; c = 18446744073709551557u;
913 assert(powmod(a, b, c) == 0);
914 a = 123; b = 0; c = 18446744073709551557u;
915 assert(powmod(a, b, c) == 1);
916
917 immutable ulong a1 = 253, b1 = 7919, c1 = 18446744073709551557u;
918 assert(powmod(a1, b1, c1) == 3883707345459248860u);
919
920 uint x = 100 ,y = 7919, z = 1844674407u;
921 assert(powmod(x, y, z) == 1613100340u);
922 x = 134; y = 7919; z = 1844674407u;
923 assert(powmod(x, y, z) == 734956622u);
924 x = 151; y = 7919; z = 1844674407u;
925 assert(powmod(x, y, z) == 1738696945u);
926 x = 168; y = 7919; z = 1844674407u;
927 assert(powmod(x, y, z) == 1247580927u);
928 x = 185; y = 7919; z = 1844674407u;
929 assert(powmod(x, y, z) == 1293855176u);
930 x = 202; y = 7919; z = 1844674407u;
931 assert(powmod(x, y, z) == 1566963682u);
932 x = 219; y = 7919; z = 1844674407u;
933 assert(powmod(x, y, z) == 181227807u);
934 x = 236; y = 7919; z = 1844674407u;
935 assert(powmod(x, y, z) == 217988321u);
936 x = 253; y = 7919; z = 1844674407u;
937 assert(powmod(x, y, z) == 1588843243u);
938
939 x = 0; y = 7919; z = 184467u;
940 assert(powmod(x, y, z) == 0);
941 x = 123; y = 0; z = 1844674u;
942 assert(powmod(x, y, z) == 1);
943
944 immutable ubyte x1 = 117;
945 immutable uint y1 = 7919;
946 immutable uint z1 = 1844674407u;
947 auto res = powmod(x1, y1, z1);
948 assert(is(typeof(res) == uint));
949 assert(res == 9479781u);
950
951 immutable ushort x2 = 123;
952 immutable uint y2 = 203;
953 immutable ubyte z2 = 113;
954 auto res2 = powmod(x2, y2, z2);
955 assert(is(typeof(res2) == ushort));
956 assert(res2 == 42u);
957 }
958
959 /**
960 * Calculates e$(SUPERSCRIPT x).
961 *
962 * $(TABLE_SV
963 * $(TR $(TH x) $(TH e$(SUPERSCRIPT x)) )
964 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
965 * $(TR $(TD -$(INFIN)) $(TD +0.0) )
966 * $(TR $(TD $(NAN)) $(TD $(NAN)) )
967 * )
968 */
969 pragma(inline, true)
970 real exp(real x) @trusted pure nothrow @nogc // TODO: @safe
971 {
972 version (InlineAsm_X87)
973 {
974 import std.math.constants : LOG2E;
975
976 // e^^x = 2^^(LOG2E*x)
977 // (This is valid because the overflow & underflow limits for exp
978 // and exp2 are so similar).
979 if (!__ctfe)
980 return exp2Asm(LOG2E*x);
981 }
982 return expImpl(x);
983 }
984
985 /// ditto
986 pragma(inline, true)
987 double exp(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) exp(cast(real) x) : expImpl(x); }
988
989 /// ditto
990 pragma(inline, true)
991 float exp(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) exp(cast(real) x) : expImpl(x); }
992
993 ///
994 @safe unittest
995 {
996 import std.math.operations : feqrel;
997 import std.math.constants : E;
998
999 assert(exp(0.0) == 1.0);
1000 assert(exp(3.0).feqrel(E * E * E) > 16);
1001 }
1002
1003 private T expImpl(T)(T x) @safe pure nothrow @nogc
1004 {
1005 import std.math.traits : floatTraits, RealFormat, isNaN;
1006 import std.math.rounding : floor;
1007 import std.math.algebraic : poly;
1008 import std.math.constants : LOG2E;
1009
1010 alias F = floatTraits!T;
1011 static if (F.realFormat == RealFormat.ieeeSingle)
1012 {
1013 static immutable T[6] P = [
1014 5.0000001201E-1,
1015 1.6666665459E-1,
1016 4.1665795894E-2,
1017 8.3334519073E-3,
1018 1.3981999507E-3,
1019 1.9875691500E-4,
1020 ];
1021
1022 enum T C1 = 0.693359375;
1023 enum T C2 = -2.12194440e-4;
1024
1025 // Overflow and Underflow limits.
1026 enum T OF = 88.72283905206835;
1027 enum T UF = -103.278929903431851103; // ln(2^-149)
1028 }
1029 else static if (F.realFormat == RealFormat.ieeeDouble)
1030 {
1031 // Coefficients for exp(x)
1032 static immutable T[3] P = [
1033 9.99999999999999999910E-1L,
1034 3.02994407707441961300E-2L,
1035 1.26177193074810590878E-4L,
1036 ];
1037 static immutable T[4] Q = [
1038 2.00000000000000000009E0L,
1039 2.27265548208155028766E-1L,
1040 2.52448340349684104192E-3L,
1041 3.00198505138664455042E-6L,
1042 ];
1043
1044 // C1 + C2 = LN2.
1045 enum T C1 = 6.93145751953125E-1;
1046 enum T C2 = 1.42860682030941723212E-6;
1047
1048 // Overflow and Underflow limits.
1049 enum T OF = 7.09782712893383996732E2; // ln((1-2^-53) * 2^1024)
1050 enum T UF = -7.451332191019412076235E2; // ln(2^-1075)
1051 }
1052 else static if (F.realFormat == RealFormat.ieeeExtended ||
1053 F.realFormat == RealFormat.ieeeExtended53)
1054 {
1055 // Coefficients for exp(x)
1056 static immutable T[3] P = [
1057 9.9999999999999999991025E-1L,
1058 3.0299440770744196129956E-2L,
1059 1.2617719307481059087798E-4L,
1060 ];
1061 static immutable T[4] Q = [
1062 2.0000000000000000000897E0L,
1063 2.2726554820815502876593E-1L,
1064 2.5244834034968410419224E-3L,
1065 3.0019850513866445504159E-6L,
1066 ];
1067
1068 // C1 + C2 = LN2.
1069 enum T C1 = 6.9314575195312500000000E-1L;
1070 enum T C2 = 1.4286068203094172321215E-6L;
1071
1072 // Overflow and Underflow limits.
1073 enum T OF = 1.1356523406294143949492E4L; // ln((1-2^-64) * 2^16384)
1074 enum T UF = -1.13994985314888605586758E4L; // ln(2^-16446)
1075 }
1076 else static if (F.realFormat == RealFormat.ieeeQuadruple)
1077 {
1078 // Coefficients for exp(x) - 1
1079 static immutable T[5] P = [
1080 9.999999999999999999999999999999999998502E-1L,
1081 3.508710990737834361215404761139478627390E-2L,
1082 2.708775201978218837374512615596512792224E-4L,
1083 6.141506007208645008909088812338454698548E-7L,
1084 3.279723985560247033712687707263393506266E-10L
1085 ];
1086 static immutable T[6] Q = [
1087 2.000000000000000000000000000000000000150E0,
1088 2.368408864814233538909747618894558968880E-1L,
1089 3.611828913847589925056132680618007270344E-3L,
1090 1.504792651814944826817779302637284053660E-5L,
1091 1.771372078166251484503904874657985291164E-8L,
1092 2.980756652081995192255342779918052538681E-12L
1093 ];
1094
1095 // C1 + C2 = LN2.
1096 enum T C1 = 6.93145751953125E-1L;
1097 enum T C2 = 1.428606820309417232121458176568075500134E-6L;
1098
1099 // Overflow and Underflow limits.
1100 enum T OF = 1.135583025911358400418251384584930671458833e4L;
1101 enum T UF = -1.143276959615573793352782661133116431383730e4L;
1102 }
1103 else
1104 static assert(0, "Not implemented for this architecture");
1105
1106 // Special cases.
1107 if (isNaN(x))
1108 return x;
1109 if (x > OF)
1110 return real.infinity;
1111 if (x < UF)
1112 return 0.0;
1113
1114 // Express: e^^x = e^^g * 2^^n
1115 // = e^^g * e^^(n * LOG2E)
1116 // = e^^(g + n * LOG2E)
1117 T xx = floor((cast(T) LOG2E) * x + cast(T) 0.5);
1118 const int n = cast(int) xx;
1119 x -= xx * C1;
1120 x -= xx * C2;
1121
1122 static if (F.realFormat == RealFormat.ieeeSingle)
1123 {
1124 xx = x * x;
1125 x = poly(x, P) * xx + x + 1.0f;
1126 }
1127 else
1128 {
1129 // Rational approximation for exponential of the fractional part:
1130 // e^^x = 1 + 2x P(x^^2) / (Q(x^^2) - P(x^^2))
1131 xx = x * x;
1132 const T px = x * poly(xx, P);
1133 x = px / (poly(xx, Q) - px);
1134 x = (cast(T) 1.0) + (cast(T) 2.0) * x;
1135 }
1136
1137 // Scale by power of 2.
1138 x = core.math.ldexp(x, n);
1139
1140 return x;
1141 }
1142
1143 @safe @nogc nothrow unittest
1144 {
1145 import std.math.traits : floatTraits, RealFormat;
1146 import std.math.operations : NaN, feqrel, isClose;
1147 import std.math.constants : E;
1148 import std.math.traits : isIdentical;
1149 import std.math.algebraic : abs;
1150
1151 version (IeeeFlagsSupport) import std.math.hardware : IeeeFlags, resetIeeeFlags, ieeeFlags;
1152 version (FloatingPointControlSupport)
1153 {
1154 import std.math.hardware : FloatingPointControl;
1155
1156 FloatingPointControl ctrl;
1157 if (FloatingPointControl.hasExceptionTraps)
1158 ctrl.disableExceptions(FloatingPointControl.allExceptions);
1159 ctrl.rounding = FloatingPointControl.roundToNearest;
1160 }
1161
1162 static void testExp(T)()
1163 {
1164 enum realFormat = floatTraits!T.realFormat;
1165 static if (realFormat == RealFormat.ieeeQuadruple)
1166 {
1167 static immutable T[2][] exptestpoints =
1168 [ // x exp(x)
1169 [ 1.0L, E ],
1170 [ 0.5L, 0x1.a61298e1e069bc972dfefab6df34p+0L ],
1171 [ 3.0L, E*E*E ],
1172 [ 0x1.6p+13L, 0x1.6e509d45728655cdb4840542acb5p+16250L ], // near overflow
1173 [ 0x1.7p+13L, T.infinity ], // close overflow
1174 [ 0x1p+80L, T.infinity ], // far overflow
1175 [ T.infinity, T.infinity ],
1176 [-0x1.18p+13L, 0x1.5e4bf54b4807034ea97fef0059a6p-12927L ], // near underflow
1177 [-0x1.625p+13L, 0x1.a6bd68a39d11fec3a250cd97f524p-16358L ], // ditto
1178 [-0x1.62dafp+13L, 0x0.cb629e9813b80ed4d639e875be6cp-16382L ], // near underflow - subnormal
1179 [-0x1.6549p+13L, 0x0.0000000000000000000000000001p-16382L ], // ditto
1180 [-0x1.655p+13L, 0 ], // close underflow
1181 [-0x1p+30L, 0 ], // far underflow
1182 ];
1183 }
1184 else static if (realFormat == RealFormat.ieeeExtended ||
1185 realFormat == RealFormat.ieeeExtended53)
1186 {
1187 static immutable T[2][] exptestpoints =
1188 [ // x exp(x)
1189 [ 1.0L, E ],
1190 [ 0.5L, 0x1.a61298e1e069bc97p+0L ],
1191 [ 3.0L, E*E*E ],
1192 [ 0x1.1p+13L, 0x1.29aeffefc8ec645p+12557L ], // near overflow
1193 [ 0x1.7p+13L, T.infinity ], // close overflow
1194 [ 0x1p+80L, T.infinity ], // far overflow
1195 [ T.infinity, T.infinity ],
1196 [-0x1.18p+13L, 0x1.5e4bf54b4806db9p-12927L ], // near underflow
1197 [-0x1.625p+13L, 0x1.a6bd68a39d11f35cp-16358L ], // ditto
1198 [-0x1.62dafp+13L, 0x1.96c53d30277021dp-16383L ], // near underflow - subnormal
1199 [-0x1.643p+13L, 0x1p-16444L ], // ditto
1200 [-0x1.645p+13L, 0 ], // close underflow
1201 [-0x1p+30L, 0 ], // far underflow
1202 ];
1203 }
1204 else static if (realFormat == RealFormat.ieeeDouble)
1205 {
1206 static immutable T[2][] exptestpoints =
1207 [ // x, exp(x)
1208 [ 1.0L, E ],
1209 [ 0.5L, 0x1.a61298e1e069cp+0L ],
1210 [ 3.0L, E*E*E ],
1211 [ 0x1.6p+9L, 0x1.93bf4ec282efbp+1015L ], // near overflow
1212 [ 0x1.7p+9L, T.infinity ], // close overflow
1213 [ 0x1p+80L, T.infinity ], // far overflow
1214 [ T.infinity, T.infinity ],
1215 [-0x1.6p+9L, 0x1.44a3824e5285fp-1016L ], // near underflow
1216 [-0x1.64p+9L, 0x0.06f84920bb2d4p-1022L ], // near underflow - subnormal
1217 [-0x1.743p+9L, 0x0.0000000000001p-1022L ], // ditto
1218 [-0x1.8p+9L, 0 ], // close underflow
1219 [-0x1p+30L, 0 ], // far underflow
1220 ];
1221 }
1222 else static if (realFormat == RealFormat.ieeeSingle)
1223 {
1224 static immutable T[2][] exptestpoints =
1225 [ // x, exp(x)
1226 [ 1.0L, E ],
1227 [ 0.5L, 0x1.a61299p+0L ],
1228 [ 3.0L, E*E*E ],
1229 [ 0x1.62p+6L, 0x1.99b988p+127L ], // near overflow
1230 [ 0x1.7p+6L, T.infinity ], // close overflow
1231 [ 0x1p+80L, T.infinity ], // far overflow
1232 [ T.infinity, T.infinity ],
1233 [-0x1.5cp+6L, 0x1.666d0ep-126L ], // near underflow
1234 [-0x1.7p+6L, 0x0.026a42p-126L ], // near underflow - subnormal
1235 [-0x1.9cp+6L, 0x0.000002p-126L ], // ditto
1236 [-0x1.ap+6L, 0 ], // close underflow
1237 [-0x1p+30L, 0 ], // far underflow
1238 ];
1239 }
1240 else
1241 static assert(0, "No exp() tests for real type!");
1242
1243 const minEqualMantissaBits = T.mant_dig - 13;
1244 T x;
1245 version (IeeeFlagsSupport) IeeeFlags f;
1246 foreach (ref pair; exptestpoints)
1247 {
1248 version (IeeeFlagsSupport) resetIeeeFlags();
1249 x = exp(pair[0]);
1250 //printf("exp(%La) = %La, should be %La\n", cast(real) pair[0], cast(real) x, cast(real) pair[1]);
1251 assert(feqrel(x, pair[1]) >= minEqualMantissaBits);
1252 }
1253
1254 // Ideally, exp(0) would not set the inexact flag.
1255 // Unfortunately, fldl2e sets it!
1256 // So it's not realistic to avoid setting it.
1257 assert(exp(cast(T) 0.0) == 1.0);
1258
1259 // NaN propagation. Doesn't set flags, bcos was already NaN.
1260 version (IeeeFlagsSupport)
1261 {
1262 resetIeeeFlags();
1263 x = exp(T.nan);
1264 f = ieeeFlags;
1265 assert(isIdentical(abs(x), T.nan));
1266 assert(f.flags == 0);
1267
1268 resetIeeeFlags();
1269 x = exp(-T.nan);
1270 f = ieeeFlags;
1271 assert(isIdentical(abs(x), T.nan));
1272 assert(f.flags == 0);
1273 }
1274 else
1275 {
1276 x = exp(T.nan);
1277 assert(isIdentical(abs(x), T.nan));
1278
1279 x = exp(-T.nan);
1280 assert(isIdentical(abs(x), T.nan));
1281 }
1282
1283 x = exp(NaN(0x123));
1284 assert(isIdentical(x, NaN(0x123)));
1285 }
1286
1287 import std.meta : AliasSeq;
1288 foreach (T; AliasSeq!(real, double, float))
1289 testExp!T();
1290
1291 // High resolution test (verified against GNU MPFR/Mathematica).
1292 assert(exp(0.5L) == 0x1.A612_98E1_E069_BC97_2DFE_FAB6_DF34p+0L);
1293
1294 assert(isClose(exp(3.0L), E * E * E, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1295 }
1296
1297 /**
1298 * Calculates the value of the natural logarithm base (e)
1299 * raised to the power of x, minus 1.
1300 *
1301 * For very small x, expm1(x) is more accurate
1302 * than exp(x)-1.
1303 *
1304 * $(TABLE_SV
1305 * $(TR $(TH x) $(TH e$(SUPERSCRIPT x)-1) )
1306 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) )
1307 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
1308 * $(TR $(TD -$(INFIN)) $(TD -1.0) )
1309 * $(TR $(TD $(NAN)) $(TD $(NAN)) )
1310 * )
1311 */
1312 pragma(inline, true)
1313 real expm1(real x) @trusted pure nothrow @nogc // TODO: @safe
1314 {
1315 version (InlineAsm_X87)
1316 {
1317 if (!__ctfe)
1318 return expm1Asm(x);
1319 }
1320 return expm1Impl(x);
1321 }
1322
1323 /// ditto
1324 pragma(inline, true)
1325 double expm1(double x) @safe pure nothrow @nogc
1326 {
1327 return __ctfe ? cast(double) expm1(cast(real) x) : expm1Impl(x);
1328 }
1329
1330 /// ditto
1331 pragma(inline, true)
1332 float expm1(float x) @safe pure nothrow @nogc
1333 {
1334 // no single-precision version in Cephes => use double precision
1335 return __ctfe ? cast(float) expm1(cast(real) x) : cast(float) expm1Impl(cast(double) x);
1336 }
1337
1338 ///
1339 @safe unittest
1340 {
1341 import std.math.traits : isIdentical;
1342 import std.math.operations : feqrel;
1343
1344 assert(isIdentical(expm1(0.0), 0.0));
1345 assert(expm1(1.0).feqrel(1.71828) > 16);
1346 assert(expm1(2.0).feqrel(6.3890) > 16);
1347 }
1348
1349 version (InlineAsm_X87)
1350 private real expm1Asm(real x) @trusted pure nothrow @nogc
1351 {
1352 version (X86)
1353 {
1354 enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
1355 asm pure nothrow @nogc
1356 {
1357 /* expm1() for x87 80-bit reals, IEEE754-2008 conformant.
1358 * Author: Don Clugston.
1359 *
1360 * expm1(x) = 2^^(rndint(y))* 2^^(y-rndint(y)) - 1 where y = LN2*x.
1361 * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^^(rndint(y))
1362 * and 2ym1 = (2^^(y-rndint(y))-1).
1363 * If 2rndy < 0.5*real.epsilon, result is -1.
1364 * Implementation is otherwise the same as for exp2()
1365 */
1366 naked;
1367 fld real ptr [ESP+4] ; // x
1368 mov AX, [ESP+4+8]; // AX = exponent and sign
1369 sub ESP, 12+8; // Create scratch space on the stack
1370 // [ESP,ESP+2] = scratchint
1371 // [ESP+4..+6, +8..+10, +10] = scratchreal
1372 // set scratchreal mantissa = 1.0
1373 mov dword ptr [ESP+8], 0;
1374 mov dword ptr [ESP+8+4], 0x80000000;
1375 and AX, 0x7FFF; // drop sign bit
1376 cmp AX, 0x401D; // avoid InvalidException in fist
1377 jae L_extreme;
1378 fldl2e;
1379 fmulp ST(1), ST; // y = x*log2(e)
1380 fist dword ptr [ESP]; // scratchint = rndint(y)
1381 fisub dword ptr [ESP]; // y - rndint(y)
1382 // and now set scratchreal exponent
1383 mov EAX, [ESP];
1384 add EAX, 0x3fff;
1385 jle short L_largenegative;
1386 cmp EAX,0x8000;
1387 jge short L_largepositive;
1388 mov [ESP+8+8],AX;
1389 f2xm1; // 2ym1 = 2^^(y-rndint(y)) -1
1390 fld real ptr [ESP+8] ; // 2rndy = 2^^rndint(y)
1391 fmul ST(1), ST; // ST=2rndy, ST(1)=2rndy*2ym1
1392 fld1;
1393 fsubp ST(1), ST; // ST = 2rndy-1, ST(1) = 2rndy * 2ym1 - 1
1394 faddp ST(1), ST; // ST = 2rndy * 2ym1 + 2rndy - 1
1395 add ESP,12+8;
1396 ret PARAMSIZE;
1397
1398 L_extreme: // Extreme exponent. X is very large positive, very
1399 // large negative, infinity, or NaN.
1400 fxam;
1401 fstsw AX;
1402 test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1403 jz L_was_nan; // if x is NaN, returns x
1404 test AX, 0x0200;
1405 jnz L_largenegative;
1406 L_largepositive:
1407 // Set scratchreal = real.max.
1408 // squaring it will create infinity, and set overflow flag.
1409 mov word ptr [ESP+8+8], 0x7FFE;
1410 fstp ST(0);
1411 fld real ptr [ESP+8]; // load scratchreal
1412 fmul ST(0), ST; // square it, to create havoc!
1413 L_was_nan:
1414 add ESP,12+8;
1415 ret PARAMSIZE;
1416 L_largenegative:
1417 fstp ST(0);
1418 fld1;
1419 fchs; // return -1. Underflow flag is not set.
1420 add ESP,12+8;
1421 ret PARAMSIZE;
1422 }
1423 }
1424 else version (X86_64)
1425 {
1426 asm pure nothrow @nogc
1427 {
1428 naked;
1429 }
1430 version (Win64)
1431 {
1432 asm pure nothrow @nogc
1433 {
1434 fld real ptr [RCX]; // x
1435 mov AX,[RCX+8]; // AX = exponent and sign
1436 }
1437 }
1438 else
1439 {
1440 asm pure nothrow @nogc
1441 {
1442 fld real ptr [RSP+8]; // x
1443 mov AX,[RSP+8+8]; // AX = exponent and sign
1444 }
1445 }
1446 asm pure nothrow @nogc
1447 {
1448 /* expm1() for x87 80-bit reals, IEEE754-2008 conformant.
1449 * Author: Don Clugston.
1450 *
1451 * expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x.
1452 * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y))
1453 * and 2ym1 = (2^(y-rndint(y))-1).
1454 * If 2rndy < 0.5*real.epsilon, result is -1.
1455 * Implementation is otherwise the same as for exp2()
1456 */
1457 sub RSP, 24; // Create scratch space on the stack
1458 // [RSP,RSP+2] = scratchint
1459 // [RSP+4..+6, +8..+10, +10] = scratchreal
1460 // set scratchreal mantissa = 1.0
1461 mov dword ptr [RSP+8], 0;
1462 mov dword ptr [RSP+8+4], 0x80000000;
1463 and AX, 0x7FFF; // drop sign bit
1464 cmp AX, 0x401D; // avoid InvalidException in fist
1465 jae L_extreme;
1466 fldl2e;
1467 fmul ; // y = x*log2(e)
1468 fist dword ptr [RSP]; // scratchint = rndint(y)
1469 fisub dword ptr [RSP]; // y - rndint(y)
1470 // and now set scratchreal exponent
1471 mov EAX, [RSP];
1472 add EAX, 0x3fff;
1473 jle short L_largenegative;
1474 cmp EAX,0x8000;
1475 jge short L_largepositive;
1476 mov [RSP+8+8],AX;
1477 f2xm1; // 2^(y-rndint(y)) -1
1478 fld real ptr [RSP+8] ; // 2^rndint(y)
1479 fmul ST(1), ST;
1480 fld1;
1481 fsubp ST(1), ST;
1482 fadd;
1483 add RSP,24;
1484 ret;
1485
1486 L_extreme: // Extreme exponent. X is very large positive, very
1487 // large negative, infinity, or NaN.
1488 fxam;
1489 fstsw AX;
1490 test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1491 jz L_was_nan; // if x is NaN, returns x
1492 test AX, 0x0200;
1493 jnz L_largenegative;
1494 L_largepositive:
1495 // Set scratchreal = real.max.
1496 // squaring it will create infinity, and set overflow flag.
1497 mov word ptr [RSP+8+8], 0x7FFE;
1498 fstp ST(0);
1499 fld real ptr [RSP+8]; // load scratchreal
1500 fmul ST(0), ST; // square it, to create havoc!
1501 L_was_nan:
1502 add RSP,24;
1503 ret;
1504
1505 L_largenegative:
1506 fstp ST(0);
1507 fld1;
1508 fchs; // return -1. Underflow flag is not set.
1509 add RSP,24;
1510 ret;
1511 }
1512 }
1513 else
1514 static assert(0);
1515 }
1516
1517 private T expm1Impl(T)(T x) @safe pure nothrow @nogc
1518 {
1519 import std.math.traits : floatTraits, RealFormat;
1520 import std.math.rounding : floor;
1521 import std.math.algebraic : poly;
1522 import std.math.constants : LN2;
1523
1524 // Coefficients for exp(x) - 1 and overflow/underflow limits.
1525 enum realFormat = floatTraits!T.realFormat;
1526 static if (realFormat == RealFormat.ieeeQuadruple)
1527 {
1528 static immutable T[8] P = [
1529 2.943520915569954073888921213330863757240E8L,
1530 -5.722847283900608941516165725053359168840E7L,
1531 8.944630806357575461578107295909719817253E6L,
1532 -7.212432713558031519943281748462837065308E5L,
1533 4.578962475841642634225390068461943438441E4L,
1534 -1.716772506388927649032068540558788106762E3L,
1535 4.401308817383362136048032038528753151144E1L,
1536 -4.888737542888633647784737721812546636240E-1L
1537 ];
1538
1539 static immutable T[9] Q = [
1540 1.766112549341972444333352727998584753865E9L,
1541 -7.848989743695296475743081255027098295771E8L,
1542 1.615869009634292424463780387327037251069E8L,
1543 -2.019684072836541751428967854947019415698E7L,
1544 1.682912729190313538934190635536631941751E6L,
1545 -9.615511549171441430850103489315371768998E4L,
1546 3.697714952261803935521187272204485251835E3L,
1547 -8.802340681794263968892934703309274564037E1L,
1548 1.0
1549 ];
1550
1551 enum T OF = 1.1356523406294143949491931077970764891253E4L;
1552 enum T UF = -1.143276959615573793352782661133116431383730e4L;
1553 }
1554 else static if (realFormat == RealFormat.ieeeExtended)
1555 {
1556 static immutable T[5] P = [
1557 -1.586135578666346600772998894928250240826E4L,
1558 2.642771505685952966904660652518429479531E3L,
1559 -3.423199068835684263987132888286791620673E2L,
1560 1.800826371455042224581246202420972737840E1L,
1561 -5.238523121205561042771939008061958820811E-1L,
1562 ];
1563 static immutable T[6] Q = [
1564 -9.516813471998079611319047060563358064497E4L,
1565 3.964866271411091674556850458227710004570E4L,
1566 -7.207678383830091850230366618190187434796E3L,
1567 7.206038318724600171970199625081491823079E2L,
1568 -4.002027679107076077238836622982900945173E1L,
1569 1.0
1570 ];
1571
1572 enum T OF = 1.1356523406294143949492E4L;
1573 enum T UF = -4.5054566736396445112120088E1L;
1574 }
1575 else static if (realFormat == RealFormat.ieeeDouble)
1576 {
1577 static immutable T[3] P = [
1578 9.9999999999999999991025E-1,
1579 3.0299440770744196129956E-2,
1580 1.2617719307481059087798E-4,
1581 ];
1582 static immutable T[4] Q = [
1583 2.0000000000000000000897E0,
1584 2.2726554820815502876593E-1,
1585 2.5244834034968410419224E-3,
1586 3.0019850513866445504159E-6,
1587 ];
1588 }
1589 else
1590 static assert(0, "no coefficients for expm1()");
1591
1592 static if (realFormat == RealFormat.ieeeDouble) // special case for double precision
1593 {
1594 if (x < -0.5 || x > 0.5)
1595 return exp(x) - 1.0;
1596 if (x == 0.0)
1597 return x;
1598
1599 const T xx = x * x;
1600 x = x * poly(xx, P);
1601 x = x / (poly(xx, Q) - x);
1602 return x + x;
1603 }
1604 else
1605 {
1606 // C1 + C2 = LN2.
1607 enum T C1 = 6.9314575195312500000000E-1L;
1608 enum T C2 = 1.428606820309417232121458176568075500134E-6L;
1609
1610 // Special cases.
1611 if (x > OF)
1612 return real.infinity;
1613 if (x == cast(T) 0.0)
1614 return x;
1615 if (x < UF)
1616 return -1.0;
1617
1618 // Express x = LN2 (n + remainder), remainder not exceeding 1/2.
1619 int n = cast(int) floor((cast(T) 0.5) + x / cast(T) LN2);
1620 x -= n * C1;
1621 x -= n * C2;
1622
1623 // Rational approximation:
1624 // exp(x) - 1 = x + 0.5 x^^2 + x^^3 P(x) / Q(x)
1625 T px = x * poly(x, P);
1626 T qx = poly(x, Q);
1627 const T xx = x * x;
1628 qx = x + ((cast(T) 0.5) * xx + xx * px / qx);
1629
1630 // We have qx = exp(remainder LN2) - 1, so:
1631 // exp(x) - 1 = 2^^n (qx + 1) - 1 = 2^^n qx + 2^^n - 1.
1632 px = core.math.ldexp(cast(T) 1.0, n);
1633 x = px * qx + (px - cast(T) 1.0);
1634
1635 return x;
1636 }
1637 }
1638
1639 @safe @nogc nothrow unittest
1640 {
1641 import std.math.traits : isNaN;
1642 import std.math.operations : isClose, CommonDefaultFor;
1643
1644 static void testExpm1(T)()
1645 {
1646 // NaN
1647 assert(isNaN(expm1(cast(T) T.nan)));
1648
1649 static immutable T[] xs = [ -2, -0.75, -0.3, 0.0, 0.1, 0.2, 0.5, 1.0 ];
1650 foreach (x; xs)
1651 {
1652 const T e = expm1(x);
1653 const T r = exp(x) - 1;
1654
1655 //printf("expm1(%Lg) = %Lg, should approximately be %Lg\n", cast(real) x, cast(real) e, cast(real) r);
1656 assert(isClose(r, e, CommonDefaultFor!(T,T), CommonDefaultFor!(T,T)));
1657 }
1658 }
1659
1660 import std.meta : AliasSeq;
1661 foreach (T; AliasSeq!(real, double))
1662 testExpm1!T();
1663 }
1664
1665 /**
1666 * Calculates 2$(SUPERSCRIPT x).
1667 *
1668 * $(TABLE_SV
1669 * $(TR $(TH x) $(TH exp2(x)) )
1670 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
1671 * $(TR $(TD -$(INFIN)) $(TD +0.0) )
1672 * $(TR $(TD $(NAN)) $(TD $(NAN)) )
1673 * )
1674 */
1675 pragma(inline, true)
1676 real exp2(real x) @nogc @trusted pure nothrow // TODO: @safe
1677 {
1678 version (InlineAsm_X87)
1679 {
1680 if (!__ctfe)
1681 return exp2Asm(x);
1682 }
1683 return exp2Impl(x);
1684 }
1685
1686 /// ditto
1687 pragma(inline, true)
1688 double exp2(double x) @nogc @safe pure nothrow { return __ctfe ? cast(double) exp2(cast(real) x) : exp2Impl(x); }
1689
1690 /// ditto
1691 pragma(inline, true)
1692 float exp2(float x) @nogc @safe pure nothrow { return __ctfe ? cast(float) exp2(cast(real) x) : exp2Impl(x); }
1693
1694 ///
1695 @safe unittest
1696 {
1697 import std.math.traits : isIdentical;
1698 import std.math.operations : feqrel;
1699
1700 assert(isIdentical(exp2(0.0), 1.0));
1701 assert(exp2(2.0).feqrel(4.0) > 16);
1702 assert(exp2(8.0).feqrel(256.0) > 16);
1703 }
1704
1705 @safe unittest
1706 {
1707 version (CRuntime_Microsoft) {} else // aexp2/exp2f/exp2l not implemented
1708 {
1709 assert( core.stdc.math.exp2f(0.0f) == 1 );
1710 assert( core.stdc.math.exp2 (0.0) == 1 );
1711 assert( core.stdc.math.exp2l(0.0L) == 1 );
1712 }
1713 }
1714
1715 version (InlineAsm_X87)
1716 private real exp2Asm(real x) @nogc @trusted pure nothrow
1717 {
1718 version (X86)
1719 {
1720 enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
1721
1722 asm pure nothrow @nogc
1723 {
1724 /* exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1725 * Author: Don Clugston.
1726 *
1727 * exp2(x) = 2^^(rndint(x))* 2^^(y-rndint(x))
1728 * The trick for high performance is to avoid the fscale(28cycles on core2),
1729 * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1730 *
1731 * We can do frndint by using fist. BUT we can't use it for huge numbers,
1732 * because it will set the Invalid Operation flag if overflow or NaN occurs.
1733 * Fortunately, whenever this happens the result would be zero or infinity.
1734 *
1735 * We can perform fscale by directly poking into the exponent. BUT this doesn't
1736 * work for the (very rare) cases where the result is subnormal. So we fall back
1737 * to the slow method in that case.
1738 */
1739 naked;
1740 fld real ptr [ESP+4] ; // x
1741 mov AX, [ESP+4+8]; // AX = exponent and sign
1742 sub ESP, 12+8; // Create scratch space on the stack
1743 // [ESP,ESP+2] = scratchint
1744 // [ESP+4..+6, +8..+10, +10] = scratchreal
1745 // set scratchreal mantissa = 1.0
1746 mov dword ptr [ESP+8], 0;
1747 mov dword ptr [ESP+8+4], 0x80000000;
1748 and AX, 0x7FFF; // drop sign bit
1749 cmp AX, 0x401D; // avoid InvalidException in fist
1750 jae L_extreme;
1751 fist dword ptr [ESP]; // scratchint = rndint(x)
1752 fisub dword ptr [ESP]; // x - rndint(x)
1753 // and now set scratchreal exponent
1754 mov EAX, [ESP];
1755 add EAX, 0x3fff;
1756 jle short L_subnormal;
1757 cmp EAX,0x8000;
1758 jge short L_overflow;
1759 mov [ESP+8+8],AX;
1760 L_normal:
1761 f2xm1;
1762 fld1;
1763 faddp ST(1), ST; // 2^^(x-rndint(x))
1764 fld real ptr [ESP+8] ; // 2^^rndint(x)
1765 add ESP,12+8;
1766 fmulp ST(1), ST;
1767 ret PARAMSIZE;
1768
1769 L_subnormal:
1770 // Result will be subnormal.
1771 // In this rare case, the simple poking method doesn't work.
1772 // The speed doesn't matter, so use the slow fscale method.
1773 fild dword ptr [ESP]; // scratchint
1774 fld1;
1775 fscale;
1776 fstp real ptr [ESP+8]; // scratchreal = 2^^scratchint
1777 fstp ST(0); // drop scratchint
1778 jmp L_normal;
1779
1780 L_extreme: // Extreme exponent. X is very large positive, very
1781 // large negative, infinity, or NaN.
1782 fxam;
1783 fstsw AX;
1784 test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1785 jz L_was_nan; // if x is NaN, returns x
1786 // set scratchreal = real.min_normal
1787 // squaring it will return 0, setting underflow flag
1788 mov word ptr [ESP+8+8], 1;
1789 test AX, 0x0200;
1790 jnz L_waslargenegative;
1791 L_overflow:
1792 // Set scratchreal = real.max.
1793 // squaring it will create infinity, and set overflow flag.
1794 mov word ptr [ESP+8+8], 0x7FFE;
1795 L_waslargenegative:
1796 fstp ST(0);
1797 fld real ptr [ESP+8]; // load scratchreal
1798 fmul ST(0), ST; // square it, to create havoc!
1799 L_was_nan:
1800 add ESP,12+8;
1801 ret PARAMSIZE;
1802 }
1803 }
1804 else version (X86_64)
1805 {
1806 asm pure nothrow @nogc
1807 {
1808 naked;
1809 }
1810 version (Win64)
1811 {
1812 asm pure nothrow @nogc
1813 {
1814 fld real ptr [RCX]; // x
1815 mov AX,[RCX+8]; // AX = exponent and sign
1816 }
1817 }
1818 else
1819 {
1820 asm pure nothrow @nogc
1821 {
1822 fld real ptr [RSP+8]; // x
1823 mov AX,[RSP+8+8]; // AX = exponent and sign
1824 }
1825 }
1826 asm pure nothrow @nogc
1827 {
1828 /* exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1829 * Author: Don Clugston.
1830 *
1831 * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x))
1832 * The trick for high performance is to avoid the fscale(28cycles on core2),
1833 * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1834 *
1835 * We can do frndint by using fist. BUT we can't use it for huge numbers,
1836 * because it will set the Invalid Operation flag is overflow or NaN occurs.
1837 * Fortunately, whenever this happens the result would be zero or infinity.
1838 *
1839 * We can perform fscale by directly poking into the exponent. BUT this doesn't
1840 * work for the (very rare) cases where the result is subnormal. So we fall back
1841 * to the slow method in that case.
1842 */
1843 sub RSP, 24; // Create scratch space on the stack
1844 // [RSP,RSP+2] = scratchint
1845 // [RSP+4..+6, +8..+10, +10] = scratchreal
1846 // set scratchreal mantissa = 1.0
1847 mov dword ptr [RSP+8], 0;
1848 mov dword ptr [RSP+8+4], 0x80000000;
1849 and AX, 0x7FFF; // drop sign bit
1850 cmp AX, 0x401D; // avoid InvalidException in fist
1851 jae L_extreme;
1852 fist dword ptr [RSP]; // scratchint = rndint(x)
1853 fisub dword ptr [RSP]; // x - rndint(x)
1854 // and now set scratchreal exponent
1855 mov EAX, [RSP];
1856 add EAX, 0x3fff;
1857 jle short L_subnormal;
1858 cmp EAX,0x8000;
1859 jge short L_overflow;
1860 mov [RSP+8+8],AX;
1861 L_normal:
1862 f2xm1;
1863 fld1;
1864 fadd; // 2^(x-rndint(x))
1865 fld real ptr [RSP+8] ; // 2^rndint(x)
1866 add RSP,24;
1867 fmulp ST(1), ST;
1868 ret;
1869
1870 L_subnormal:
1871 // Result will be subnormal.
1872 // In this rare case, the simple poking method doesn't work.
1873 // The speed doesn't matter, so use the slow fscale method.
1874 fild dword ptr [RSP]; // scratchint
1875 fld1;
1876 fscale;
1877 fstp real ptr [RSP+8]; // scratchreal = 2^scratchint
1878 fstp ST(0); // drop scratchint
1879 jmp L_normal;
1880
1881 L_extreme: // Extreme exponent. X is very large positive, very
1882 // large negative, infinity, or NaN.
1883 fxam;
1884 fstsw AX;
1885 test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1886 jz L_was_nan; // if x is NaN, returns x
1887 // set scratchreal = real.min
1888 // squaring it will return 0, setting underflow flag
1889 mov word ptr [RSP+8+8], 1;
1890 test AX, 0x0200;
1891 jnz L_waslargenegative;
1892 L_overflow:
1893 // Set scratchreal = real.max.
1894 // squaring it will create infinity, and set overflow flag.
1895 mov word ptr [RSP+8+8], 0x7FFE;
1896 L_waslargenegative:
1897 fstp ST(0);
1898 fld real ptr [RSP+8]; // load scratchreal
1899 fmul ST(0), ST; // square it, to create havoc!
1900 L_was_nan:
1901 add RSP,24;
1902 ret;
1903 }
1904 }
1905 else
1906 static assert(0);
1907 }
1908
1909 private T exp2Impl(T)(T x) @nogc @safe pure nothrow
1910 {
1911 import std.math.traits : floatTraits, RealFormat, isNaN;
1912 import std.math.rounding : floor;
1913 import std.math.algebraic : poly;
1914
1915 // Coefficients for exp2(x)
1916 enum realFormat = floatTraits!T.realFormat;
1917 static if (realFormat == RealFormat.ieeeQuadruple)
1918 {
1919 static immutable T[5] P = [
1920 9.079594442980146270952372234833529694788E12L,
1921 1.530625323728429161131811299626419117557E11L,
1922 5.677513871931844661829755443994214173883E8L,
1923 6.185032670011643762127954396427045467506E5L,
1924 1.587171580015525194694938306936721666031E2L
1925 ];
1926
1927 static immutable T[6] Q = [
1928 2.619817175234089411411070339065679229869E13L,
1929 1.490560994263653042761789432690793026977E12L,
1930 1.092141473886177435056423606755843616331E10L,
1931 2.186249607051644894762167991800811827835E7L,
1932 1.236602014442099053716561665053645270207E4L,
1933 1.0
1934 ];
1935 }
1936 else static if (realFormat == RealFormat.ieeeExtended)
1937 {
1938 static immutable T[3] P = [
1939 2.0803843631901852422887E6L,
1940 3.0286971917562792508623E4L,
1941 6.0614853552242266094567E1L,
1942 ];
1943 static immutable T[4] Q = [
1944 6.0027204078348487957118E6L,
1945 3.2772515434906797273099E5L,
1946 1.7492876999891839021063E3L,
1947 1.0000000000000000000000E0L,
1948 ];
1949 }
1950 else static if (realFormat == RealFormat.ieeeDouble)
1951 {
1952 static immutable T[3] P = [
1953 1.51390680115615096133E3L,
1954 2.02020656693165307700E1L,
1955 2.30933477057345225087E-2L,
1956 ];
1957 static immutable T[3] Q = [
1958 4.36821166879210612817E3L,
1959 2.33184211722314911771E2L,
1960 1.00000000000000000000E0L,
1961 ];
1962 }
1963 else static if (realFormat == RealFormat.ieeeSingle)
1964 {
1965 static immutable T[6] P = [
1966 6.931472028550421E-001L,
1967 2.402264791363012E-001L,
1968 5.550332471162809E-002L,
1969 9.618437357674640E-003L,
1970 1.339887440266574E-003L,
1971 1.535336188319500E-004L,
1972 ];
1973 }
1974 else
1975 static assert(0, "no coefficients for exp2()");
1976
1977 // Overflow and Underflow limits.
1978 enum T OF = T.max_exp;
1979 enum T UF = T.min_exp - 1;
1980
1981 // Special cases.
1982 if (isNaN(x))
1983 return x;
1984 if (x > OF)
1985 return real.infinity;
1986 if (x < UF)
1987 return 0.0;
1988
1989 static if (realFormat == RealFormat.ieeeSingle) // special case for single precision
1990 {
1991 // The following is necessary because range reduction blows up.
1992 if (x == 0.0f)
1993 return 1.0f;
1994
1995 // Separate into integer and fractional parts.
1996 const T i = floor(x);
1997 int n = cast(int) i;
1998 x -= i;
1999 if (x > 0.5f)
2000 {
2001 n += 1;
2002 x -= 1.0f;
2003 }
2004
2005 // Rational approximation:
2006 // exp2(x) = 1.0 + x P(x)
2007 x = 1.0f + x * poly(x, P);
2008 }
2009 else
2010 {
2011 // Separate into integer and fractional parts.
2012 const T i = floor(x + cast(T) 0.5);
2013 int n = cast(int) i;
2014 x -= i;
2015
2016 // Rational approximation:
2017 // exp2(x) = 1.0 + 2x P(x^^2) / (Q(x^^2) - P(x^^2))
2018 const T xx = x * x;
2019 const T px = x * poly(xx, P);
2020 x = px / (poly(xx, Q) - px);
2021 x = (cast(T) 1.0) + (cast(T) 2.0) * x;
2022 }
2023
2024 // Scale by power of 2.
2025 x = core.math.ldexp(x, n);
2026
2027 return x;
2028 }
2029
2030 @safe @nogc nothrow unittest
2031 {
2032 import std.math.operations : feqrel, NaN, isClose;
2033 import std.math.traits : isIdentical;
2034 import std.math.constants : SQRT2;
2035
2036 assert(feqrel(exp2(0.5L), SQRT2) >= real.mant_dig -1);
2037 assert(exp2(8.0L) == 256.0);
2038 assert(exp2(-9.0L)== 1.0L/512.0);
2039
2040 static void testExp2(T)()
2041 {
2042 // NaN
2043 const T specialNaN = NaN(0x0123L);
2044 assert(isIdentical(exp2(specialNaN), specialNaN));
2045
2046 // over-/underflow
2047 enum T OF = T.max_exp;
2048 enum T UF = T.min_exp - T.mant_dig;
2049 assert(isIdentical(exp2(OF + 1), cast(T) T.infinity));
2050 assert(isIdentical(exp2(UF - 1), cast(T) 0.0));
2051
2052 static immutable T[2][] vals =
2053 [
2054 // x, exp2(x)
2055 [ 0.0, 1.0 ],
2056 [ -0.0, 1.0 ],
2057 [ 0.5, SQRT2 ],
2058 [ 8.0, 256.0 ],
2059 [ -9.0, 1.0 / 512 ],
2060 ];
2061
2062 foreach (ref val; vals)
2063 {
2064 const T x = val[0];
2065 const T r = val[1];
2066 const T e = exp2(x);
2067
2068 //printf("exp2(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) e, cast(real) r);
2069 assert(isClose(r, e));
2070 }
2071 }
2072
2073 import std.meta : AliasSeq;
2074 foreach (T; AliasSeq!(real, double, float))
2075 testExp2!T();
2076 }
2077
2078 /*********************************************************************
2079 * Separate floating point value into significand and exponent.
2080 *
2081 * Returns:
2082 * Calculate and return $(I x) and $(I exp) such that
2083 * value =$(I x)*2$(SUPERSCRIPT exp) and
2084 * .5 $(LT)= |$(I x)| $(LT) 1.0
2085 *
2086 * $(I x) has same sign as value.
2087 *
2088 * $(TABLE_SV
2089 * $(TR $(TH value) $(TH returns) $(TH exp))
2090 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD 0))
2091 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD int.max))
2092 * $(TR $(TD -$(INFIN)) $(TD -$(INFIN)) $(TD int.min))
2093 * $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min))
2094 * )
2095 */
2096 T frexp(T)(const T value, out int exp) @trusted pure nothrow @nogc
2097 if (isFloatingPoint!T)
2098 {
2099 import std.math.traits : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB, isSubnormal;
2100
2101 if (__ctfe)
2102 {
2103 // Handle special cases.
2104 if (value == 0) { exp = 0; return value; }
2105 if (value == T.infinity) { exp = int.max; return value; }
2106 if (value == -T.infinity || value != value) { exp = int.min; return value; }
2107 // Handle ordinary cases.
2108 // In CTFE there is no performance advantage for having separate
2109 // paths for different floating point types.
2110 T absValue = value < 0 ? -value : value;
2111 int expCount;
2112 static if (T.mant_dig > double.mant_dig)
2113 {
2114 for (; absValue >= 0x1.0p+1024L; absValue *= 0x1.0p-1024L)
2115 expCount += 1024;
2116 for (; absValue < 0x1.0p-1021L; absValue *= 0x1.0p+1021L)
2117 expCount -= 1021;
2118 }
2119 const double dval = cast(double) absValue;
2120 int dexp = cast(int) (((*cast(const long*) &dval) >>> 52) & 0x7FF) + double.min_exp - 2;
2121 dexp++;
2122 expCount += dexp;
2123 absValue *= 2.0 ^^ -dexp;
2124 // If the original value was subnormal or if it was a real
2125 // then absValue can still be outside the [0.5, 1.0) range.
2126 if (absValue < 0.5)
2127 {
2128 assert(T.mant_dig > double.mant_dig || isSubnormal(value));
2129 do
2130 {
2131 absValue += absValue;
2132 expCount--;
2133 } while (absValue < 0.5);
2134 }
2135 else
2136 {
2137 assert(absValue < 1 || T.mant_dig > double.mant_dig);
2138 for (; absValue >= 1; absValue *= T(0.5))
2139 expCount++;
2140 }
2141 exp = expCount;
2142 return value < 0 ? -absValue : absValue;
2143 }
2144
2145 Unqual!T vf = value;
2146 ushort* vu = cast(ushort*)&vf;
2147 static if (is(immutable T == immutable float))
2148 int* vi = cast(int*)&vf;
2149 else
2150 long* vl = cast(long*)&vf;
2151 int ex;
2152 alias F = floatTraits!T;
2153
2154 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2155 static if (F.realFormat == RealFormat.ieeeExtended ||
2156 F.realFormat == RealFormat.ieeeExtended53)
2157 {
2158 if (ex)
2159 { // If exponent is non-zero
2160 if (ex == F.EXPMASK) // infinity or NaN
2161 {
2162 if (*vl & 0x7FFF_FFFF_FFFF_FFFF) // NaN
2163 {
2164 *vl |= 0xC000_0000_0000_0000; // convert NaNS to NaNQ
2165 exp = int.min;
2166 }
2167 else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity
2168 exp = int.min;
2169 else // positive infinity
2170 exp = int.max;
2171
2172 }
2173 else
2174 {
2175 exp = ex - F.EXPBIAS;
2176 vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE;
2177 }
2178 }
2179 else if (!*vl)
2180 {
2181 // vf is +-0.0
2182 exp = 0;
2183 }
2184 else
2185 {
2186 // subnormal
2187
2188 vf *= F.RECIP_EPSILON;
2189 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2190 exp = ex - F.EXPBIAS - T.mant_dig + 1;
2191 vu[F.EXPPOS_SHORT] = ((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FFE;
2192 }
2193 return vf;
2194 }
2195 else static if (F.realFormat == RealFormat.ieeeQuadruple)
2196 {
2197 if (ex) // If exponent is non-zero
2198 {
2199 if (ex == F.EXPMASK)
2200 {
2201 // infinity or NaN
2202 if (vl[MANTISSA_LSB] |
2203 (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN
2204 {
2205 // convert NaNS to NaNQ
2206 vl[MANTISSA_MSB] |= 0x0000_8000_0000_0000;
2207 exp = int.min;
2208 }
2209 else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity
2210 exp = int.min;
2211 else // positive infinity
2212 exp = int.max;
2213 }
2214 else
2215 {
2216 exp = ex - F.EXPBIAS;
2217 vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]);
2218 }
2219 }
2220 else if ((vl[MANTISSA_LSB] |
2221 (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0)
2222 {
2223 // vf is +-0.0
2224 exp = 0;
2225 }
2226 else
2227 {
2228 // subnormal
2229 vf *= F.RECIP_EPSILON;
2230 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2231 exp = ex - F.EXPBIAS - T.mant_dig + 1;
2232 vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]);
2233 }
2234 return vf;
2235 }
2236 else static if (F.realFormat == RealFormat.ieeeDouble)
2237 {
2238 if (ex) // If exponent is non-zero
2239 {
2240 if (ex == F.EXPMASK) // infinity or NaN
2241 {
2242 if (*vl == 0x7FF0_0000_0000_0000) // positive infinity
2243 {
2244 exp = int.max;
2245 }
2246 else if (*vl == 0xFFF0_0000_0000_0000) // negative infinity
2247 exp = int.min;
2248 else
2249 { // NaN
2250 *vl |= 0x0008_0000_0000_0000; // convert NaNS to NaNQ
2251 exp = int.min;
2252 }
2253 }
2254 else
2255 {
2256 exp = (ex - F.EXPBIAS) >> 4;
2257 vu[F.EXPPOS_SHORT] = cast(ushort)((0x800F & vu[F.EXPPOS_SHORT]) | 0x3FE0);
2258 }
2259 }
2260 else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF))
2261 {
2262 // vf is +-0.0
2263 exp = 0;
2264 }
2265 else
2266 {
2267 // subnormal
2268 vf *= F.RECIP_EPSILON;
2269 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2270 exp = ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1;
2271 vu[F.EXPPOS_SHORT] =
2272 cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FE0);
2273 }
2274 return vf;
2275 }
2276 else static if (F.realFormat == RealFormat.ieeeSingle)
2277 {
2278 if (ex) // If exponent is non-zero
2279 {
2280 if (ex == F.EXPMASK) // infinity or NaN
2281 {
2282 if (*vi == 0x7F80_0000) // positive infinity
2283 {
2284 exp = int.max;
2285 }
2286 else if (*vi == 0xFF80_0000) // negative infinity
2287 exp = int.min;
2288 else
2289 { // NaN
2290 *vi |= 0x0040_0000; // convert NaNS to NaNQ
2291 exp = int.min;
2292 }
2293 }
2294 else
2295 {
2296 exp = (ex - F.EXPBIAS) >> 7;
2297 vu[F.EXPPOS_SHORT] = cast(ushort)((0x807F & vu[F.EXPPOS_SHORT]) | 0x3F00);
2298 }
2299 }
2300 else if (!(*vi & 0x7FFF_FFFF))
2301 {
2302 // vf is +-0.0
2303 exp = 0;
2304 }
2305 else
2306 {
2307 // subnormal
2308 vf *= F.RECIP_EPSILON;
2309 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2310 exp = ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1;
2311 vu[F.EXPPOS_SHORT] =
2312 cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3F00);
2313 }
2314 return vf;
2315 }
2316 else // static if (F.realFormat == RealFormat.ibmExtended)
2317 {
2318 assert(0, "frexp not implemented");
2319 }
2320 }
2321
2322 ///
2323 @safe unittest
2324 {
2325 import std.math.operations : isClose;
2326
2327 int exp;
2328 real mantissa = frexp(123.456L, exp);
2329
2330 assert(isClose(mantissa * pow(2.0L, cast(real) exp), 123.456L));
2331
2332 assert(frexp(-real.nan, exp) && exp == int.min);
2333 assert(frexp(real.nan, exp) && exp == int.min);
2334 assert(frexp(-real.infinity, exp) == -real.infinity && exp == int.min);
2335 assert(frexp(real.infinity, exp) == real.infinity && exp == int.max);
2336 assert(frexp(-0.0, exp) == -0.0 && exp == 0);
2337 assert(frexp(0.0, exp) == 0.0 && exp == 0);
2338 }
2339
2340 @safe @nogc nothrow unittest
2341 {
2342 import std.math.operations : isClose;
2343
2344 int exp;
2345 real mantissa = frexp(123.456L, exp);
2346
2347 // check if values are equal to 19 decimal digits of precision
2348 assert(isClose(mantissa * pow(2.0L, cast(real) exp), 123.456L, 1e-18));
2349 }
2350
2351 @safe unittest
2352 {
2353 import std.math.traits : floatTraits, RealFormat, isIdentical;
2354 import std.meta : AliasSeq;
2355 import std.typecons : tuple, Tuple;
2356
2357 static foreach (T; AliasSeq!(real, double, float))
2358 {{
2359 Tuple!(T, T, int)[] vals = [ // x,frexp,exp
2360 tuple(T(0.0), T( 0.0 ), 0),
2361 tuple(T(-0.0), T( -0.0), 0),
2362 tuple(T(1.0), T( .5 ), 1),
2363 tuple(T(-1.0), T( -.5 ), 1),
2364 tuple(T(2.0), T( .5 ), 2),
2365 tuple(T(float.min_normal/2.0f), T(.5), -126),
2366 tuple(T.infinity, T.infinity, int.max),
2367 tuple(-T.infinity, -T.infinity, int.min),
2368 tuple(T.nan, T.nan, int.min),
2369 tuple(-T.nan, -T.nan, int.min),
2370
2371 // https://issues.dlang.org/show_bug.cgi?id=16026:
2372 tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2)
2373 ];
2374
2375 foreach (elem; vals)
2376 {
2377 T x = elem[0];
2378 T e = elem[1];
2379 int exp = elem[2];
2380 int eptr;
2381 T v = frexp(x, eptr);
2382 assert(isIdentical(e, v));
2383 assert(exp == eptr);
2384 }
2385
2386 static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended)
2387 {
2388 static T[3][] extendedvals = [ // x,frexp,exp
2389 [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal
2390 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063],
2391 [T.min_normal, .5, -16381],
2392 [T.min_normal/2.0L, .5, -16382] // subnormal
2393 ];
2394 foreach (elem; extendedvals)
2395 {
2396 T x = elem[0];
2397 T e = elem[1];
2398 int exp = cast(int) elem[2];
2399 int eptr;
2400 T v = frexp(x, eptr);
2401 assert(isIdentical(e, v));
2402 assert(exp == eptr);
2403 }
2404 }
2405 }}
2406
2407 // CTFE
2408 alias CtfeFrexpResult= Tuple!(real, int);
2409 static CtfeFrexpResult ctfeFrexp(T)(const T value)
2410 {
2411 int exp;
2412 auto significand = frexp(value, exp);
2413 return CtfeFrexpResult(significand, exp);
2414 }
2415 static foreach (T; AliasSeq!(real, double, float))
2416 {{
2417 enum Tuple!(T, T, int)[] vals = [ // x,frexp,exp
2418 tuple(T(0.0), T( 0.0 ), 0),
2419 tuple(T(-0.0), T( -0.0), 0),
2420 tuple(T(1.0), T( .5 ), 1),
2421 tuple(T(-1.0), T( -.5 ), 1),
2422 tuple(T(2.0), T( .5 ), 2),
2423 tuple(T(float.min_normal/2.0f), T(.5), -126),
2424 tuple(T.infinity, T.infinity, int.max),
2425 tuple(-T.infinity, -T.infinity, int.min),
2426 tuple(T.nan, T.nan, int.min),
2427 tuple(-T.nan, -T.nan, int.min),
2428
2429 // https://issues.dlang.org/show_bug.cgi?id=16026:
2430 tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2)
2431 ];
2432
2433 static foreach (elem; vals)
2434 {
2435 static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], elem[2]));
2436 }
2437
2438 static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended)
2439 {
2440 enum T[3][] extendedvals = [ // x,frexp,exp
2441 [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal
2442 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063],
2443 [T.min_normal, .5, -16381],
2444 [T.min_normal/2.0L, .5, -16382] // subnormal
2445 ];
2446 static foreach (elem; extendedvals)
2447 {
2448 static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], cast(int) elem[2]));
2449 }
2450 }
2451 }}
2452 }
2453
2454 @safe unittest
2455 {
2456 import std.meta : AliasSeq;
2457 void foo() {
2458 static foreach (T; AliasSeq!(real, double, float))
2459 {{
2460 int exp;
2461 const T a = 1;
2462 immutable T b = 2;
2463 auto c = frexp(a, exp);
2464 auto d = frexp(b, exp);
2465 }}
2466 }
2467 }
2468
2469 /******************************************
2470 * Extracts the exponent of x as a signed integral value.
2471 *
2472 * If x is not a special value, the result is the same as
2473 * $(D cast(int) logb(x)).
2474 *
2475 * $(TABLE_SV
2476 * $(TR $(TH x) $(TH ilogb(x)) $(TH Range error?))
2477 * $(TR $(TD 0) $(TD FP_ILOGB0) $(TD yes))
2478 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD int.max) $(TD no))
2479 * $(TR $(TD $(NAN)) $(TD FP_ILOGBNAN) $(TD no))
2480 * )
2481 */
2482 int ilogb(T)(const T x) @trusted pure nothrow @nogc
2483 if (isFloatingPoint!T)
2484 {
2485 import std.math.traits : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB;
2486
2487 import core.bitop : bsr;
2488 alias F = floatTraits!T;
2489
2490 union floatBits
2491 {
2492 T rv;
2493 ushort[T.sizeof/2] vu;
2494 uint[T.sizeof/4] vui;
2495 static if (T.sizeof >= 8)
2496 ulong[T.sizeof/8] vul;
2497 }
2498 floatBits y = void;
2499 y.rv = x;
2500
2501 int ex = y.vu[F.EXPPOS_SHORT] & F.EXPMASK;
2502 static if (F.realFormat == RealFormat.ieeeExtended ||
2503 F.realFormat == RealFormat.ieeeExtended53)
2504 {
2505 if (ex)
2506 {
2507 // If exponent is non-zero
2508 if (ex == F.EXPMASK) // infinity or NaN
2509 {
2510 if (y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) // NaN
2511 return FP_ILOGBNAN;
2512 else // +-infinity
2513 return int.max;
2514 }
2515 else
2516 {
2517 return ex - F.EXPBIAS - 1;
2518 }
2519 }
2520 else if (!y.vul[0])
2521 {
2522 // vf is +-0.0
2523 return FP_ILOGB0;
2524 }
2525 else
2526 {
2527 // subnormal
2528 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(y.vul[0]);
2529 }
2530 }
2531 else static if (F.realFormat == RealFormat.ieeeQuadruple)
2532 {
2533 if (ex) // If exponent is non-zero
2534 {
2535 if (ex == F.EXPMASK)
2536 {
2537 // infinity or NaN
2538 if (y.vul[MANTISSA_LSB] | ( y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN
2539 return FP_ILOGBNAN;
2540 else // +- infinity
2541 return int.max;
2542 }
2543 else
2544 {
2545 return ex - F.EXPBIAS - 1;
2546 }
2547 }
2548 else if ((y.vul[MANTISSA_LSB] | (y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0)
2549 {
2550 // vf is +-0.0
2551 return FP_ILOGB0;
2552 }
2553 else
2554 {
2555 // subnormal
2556 const ulong msb = y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF;
2557 const ulong lsb = y.vul[MANTISSA_LSB];
2558 if (msb)
2559 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(msb) + 64;
2560 else
2561 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(lsb);
2562 }
2563 }
2564 else static if (F.realFormat == RealFormat.ieeeDouble)
2565 {
2566 if (ex) // If exponent is non-zero
2567 {
2568 if (ex == F.EXPMASK) // infinity or NaN
2569 {
2570 if ((y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FF0_0000_0000_0000) // +- infinity
2571 return int.max;
2572 else // NaN
2573 return FP_ILOGBNAN;
2574 }
2575 else
2576 {
2577 return ((ex - F.EXPBIAS) >> 4) - 1;
2578 }
2579 }
2580 else if (!(y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF))
2581 {
2582 // vf is +-0.0
2583 return FP_ILOGB0;
2584 }
2585 else
2586 {
2587 // subnormal
2588 enum MANTISSAMASK_64 = ((cast(ulong) F.MANTISSAMASK_INT) << 32) | 0xFFFF_FFFF;
2589 return ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1 + bsr(y.vul[0] & MANTISSAMASK_64);
2590 }
2591 }
2592 else static if (F.realFormat == RealFormat.ieeeSingle)
2593 {
2594 if (ex) // If exponent is non-zero
2595 {
2596 if (ex == F.EXPMASK) // infinity or NaN
2597 {
2598 if ((y.vui[0] & 0x7FFF_FFFF) == 0x7F80_0000) // +- infinity
2599 return int.max;
2600 else // NaN
2601 return FP_ILOGBNAN;
2602 }
2603 else
2604 {
2605 return ((ex - F.EXPBIAS) >> 7) - 1;
2606 }
2607 }
2608 else if (!(y.vui[0] & 0x7FFF_FFFF))
2609 {
2610 // vf is +-0.0
2611 return FP_ILOGB0;
2612 }
2613 else
2614 {
2615 // subnormal
2616 const uint mantissa = y.vui[0] & F.MANTISSAMASK_INT;
2617 return ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1 + bsr(mantissa);
2618 }
2619 }
2620 else // static if (F.realFormat == RealFormat.ibmExtended)
2621 {
2622 assert(0, "ilogb not implemented");
2623 }
2624 }
2625 /// ditto
2626 int ilogb(T)(const T x) @safe pure nothrow @nogc
2627 if (isIntegral!T && isUnsigned!T)
2628 {
2629 import core.bitop : bsr;
2630 if (x == 0)
2631 return FP_ILOGB0;
2632 else
2633 {
2634 static assert(T.sizeof <= ulong.sizeof, "integer size too large for the current ilogb implementation");
2635 return bsr(x);
2636 }
2637 }
2638 /// ditto
2639 int ilogb(T)(const T x) @safe pure nothrow @nogc
2640 if (isIntegral!T && isSigned!T)
2641 {
2642 import std.traits : Unsigned;
2643 // Note: abs(x) can not be used because the return type is not Unsigned and
2644 // the return value would be wrong for x == int.min
2645 Unsigned!T absx = x >= 0 ? x : -x;
2646 return ilogb(absx);
2647 }
2648
2649 ///
2650 @safe pure unittest
2651 {
2652 assert(ilogb(1) == 0);
2653 assert(ilogb(3) == 1);
2654 assert(ilogb(3.0) == 1);
2655 assert(ilogb(100_000_000) == 26);
2656
2657 assert(ilogb(0) == FP_ILOGB0);
2658 assert(ilogb(0.0) == FP_ILOGB0);
2659 assert(ilogb(double.nan) == FP_ILOGBNAN);
2660 assert(ilogb(double.infinity) == int.max);
2661 }
2662
2663 /**
2664 Special return values of $(LREF ilogb).
2665 */
2666 alias FP_ILOGB0 = core.stdc.math.FP_ILOGB0;
2667 /// ditto
2668 alias FP_ILOGBNAN = core.stdc.math.FP_ILOGBNAN;
2669
2670 ///
2671 @safe pure unittest
2672 {
2673 assert(ilogb(0) == FP_ILOGB0);
2674 assert(ilogb(0.0) == FP_ILOGB0);
2675 assert(ilogb(double.nan) == FP_ILOGBNAN);
2676 }
2677
2678 @safe nothrow @nogc unittest
2679 {
2680 import std.math.traits : floatTraits, RealFormat;
2681 import std.math.operations : nextUp;
2682 import std.meta : AliasSeq;
2683 import std.typecons : Tuple;
2684 static foreach (F; AliasSeq!(float, double, real))
2685 {{
2686 alias T = Tuple!(F, int);
2687 T[13] vals = // x, ilogb(x)
2688 [
2689 T( F.nan , FP_ILOGBNAN ),
2690 T( -F.nan , FP_ILOGBNAN ),
2691 T( F.infinity, int.max ),
2692 T( -F.infinity, int.max ),
2693 T( 0.0 , FP_ILOGB0 ),
2694 T( -0.0 , FP_ILOGB0 ),
2695 T( 2.0 , 1 ),
2696 T( 2.0001 , 1 ),
2697 T( 1.9999 , 0 ),
2698 T( 0.5 , -1 ),
2699 T( 123.123 , 6 ),
2700 T( -123.123 , 6 ),
2701 T( 0.123 , -4 ),
2702 ];
2703
2704 foreach (elem; vals)
2705 {
2706 assert(ilogb(elem[0]) == elem[1]);
2707 }
2708 }}
2709
2710 // min_normal and subnormals
2711 assert(ilogb(-float.min_normal) == -126);
2712 assert(ilogb(nextUp(-float.min_normal)) == -127);
2713 assert(ilogb(nextUp(-float(0.0))) == -149);
2714 assert(ilogb(-double.min_normal) == -1022);
2715 assert(ilogb(nextUp(-double.min_normal)) == -1023);
2716 assert(ilogb(nextUp(-double(0.0))) == -1074);
2717 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
2718 {
2719 assert(ilogb(-real.min_normal) == -16382);
2720 assert(ilogb(nextUp(-real.min_normal)) == -16383);
2721 assert(ilogb(nextUp(-real(0.0))) == -16445);
2722 }
2723 else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
2724 {
2725 assert(ilogb(-real.min_normal) == -1022);
2726 assert(ilogb(nextUp(-real.min_normal)) == -1023);
2727 assert(ilogb(nextUp(-real(0.0))) == -1074);
2728 }
2729
2730 // test integer types
2731 assert(ilogb(0) == FP_ILOGB0);
2732 assert(ilogb(int.max) == 30);
2733 assert(ilogb(int.min) == 31);
2734 assert(ilogb(uint.max) == 31);
2735 assert(ilogb(long.max) == 62);
2736 assert(ilogb(long.min) == 63);
2737 assert(ilogb(ulong.max) == 63);
2738 }
2739
2740 /*******************************************
2741 * Compute n * 2$(SUPERSCRIPT exp)
2742 * References: frexp
2743 */
2744
2745 pragma(inline, true)
2746 real ldexp(real n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
2747 ///ditto
2748 pragma(inline, true)
2749 double ldexp(double n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
2750 ///ditto
2751 pragma(inline, true)
2752 float ldexp(float n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
2753
2754 ///
2755 @nogc @safe pure nothrow unittest
2756 {
2757 import std.meta : AliasSeq;
2758 static foreach (T; AliasSeq!(float, double, real))
2759 {{
2760 T r;
2761
2762 r = ldexp(3.0L, 3);
2763 assert(r == 24);
2764
2765 r = ldexp(cast(T) 3.0, cast(int) 3);
2766 assert(r == 24);
2767
2768 T n = 3.0;
2769 int exp = 3;
2770 r = ldexp(n, exp);
2771 assert(r == 24);
2772 }}
2773 }
2774
2775 @safe pure nothrow @nogc unittest
2776 {
2777 import std.math.traits : floatTraits, RealFormat;
2778
2779 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended ||
2780 floatTraits!(real).realFormat == RealFormat.ieeeExtended53 ||
2781 floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
2782 {
2783 assert(ldexp(1.0L, -16384) == 0x1p-16384L);
2784 assert(ldexp(1.0L, -16382) == 0x1p-16382L);
2785 int x;
2786 real n = frexp(0x1p-16384L, x);
2787 assert(n == 0.5L);
2788 assert(x==-16383);
2789 assert(ldexp(n, x)==0x1p-16384L);
2790 }
2791 else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
2792 {
2793 assert(ldexp(1.0L, -1024) == 0x1p-1024L);
2794 assert(ldexp(1.0L, -1022) == 0x1p-1022L);
2795 int x;
2796 real n = frexp(0x1p-1024L, x);
2797 assert(n == 0.5L);
2798 assert(x==-1023);
2799 assert(ldexp(n, x)==0x1p-1024L);
2800 }
2801 else static assert(false, "Floating point type real not supported");
2802 }
2803
2804 /* workaround https://issues.dlang.org/show_bug.cgi?id=14718
2805 float parsing depends on platform strtold
2806 @safe pure nothrow @nogc unittest
2807 {
2808 assert(ldexp(1.0, -1024) == 0x1p-1024);
2809 assert(ldexp(1.0, -1022) == 0x1p-1022);
2810 int x;
2811 double n = frexp(0x1p-1024, x);
2812 assert(n == 0.5);
2813 assert(x==-1023);
2814 assert(ldexp(n, x)==0x1p-1024);
2815 }
2816
2817 @safe pure nothrow @nogc unittest
2818 {
2819 assert(ldexp(1.0f, -128) == 0x1p-128f);
2820 assert(ldexp(1.0f, -126) == 0x1p-126f);
2821 int x;
2822 float n = frexp(0x1p-128f, x);
2823 assert(n == 0.5f);
2824 assert(x==-127);
2825 assert(ldexp(n, x)==0x1p-128f);
2826 }
2827 */
2828
2829 @safe @nogc nothrow unittest
2830 {
2831 import std.math.operations : isClose;
2832
2833 static real[3][] vals = // value,exp,ldexp
2834 [
2835 [ 0, 0, 0],
2836 [ 1, 0, 1],
2837 [ -1, 0, -1],
2838 [ 1, 1, 2],
2839 [ 123, 10, 125952],
2840 [ real.max, int.max, real.infinity],
2841 [ real.max, -int.max, 0],
2842 [ real.min_normal, -int.max, 0],
2843 ];
2844 int i;
2845
2846 for (i = 0; i < vals.length; i++)
2847 {
2848 real x = vals[i][0];
2849 int exp = cast(int) vals[i][1];
2850 real z = vals[i][2];
2851 real l = ldexp(x, exp);
2852
2853 assert(isClose(z, l, 1e-6));
2854 }
2855
2856 real function(real, int) pldexp = &ldexp;
2857 assert(pldexp != null);
2858 }
2859
2860 private
2861 {
2862 // Coefficients shared across log(), log2(), log10(), log1p().
2863 template LogCoeffs(T)
2864 {
2865 import std.math.traits : floatTraits, RealFormat;
2866
2867 static if (floatTraits!T.realFormat == RealFormat.ieeeQuadruple)
2868 {
2869 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
2870 // Theoretical peak relative error = 5.3e-37
2871 static immutable real[13] logP = [
2872 1.313572404063446165910279910527789794488E4L,
2873 7.771154681358524243729929227226708890930E4L,
2874 2.014652742082537582487669938141683759923E5L,
2875 3.007007295140399532324943111654767187848E5L,
2876 2.854829159639697837788887080758954924001E5L,
2877 1.797628303815655343403735250238293741397E5L,
2878 7.594356839258970405033155585486712125861E4L,
2879 2.128857716871515081352991964243375186031E4L,
2880 3.824952356185897735160588078446136783779E3L,
2881 4.114517881637811823002128927449878962058E2L,
2882 2.321125933898420063925789532045674660756E1L,
2883 4.998469661968096229986658302195402690910E-1L,
2884 1.538612243596254322971797716843006400388E-6L
2885 ];
2886 static immutable real[13] logQ = [
2887 3.940717212190338497730839731583397586124E4L,
2888 2.626900195321832660448791748036714883242E5L,
2889 7.777690340007566932935753241556479363645E5L,
2890 1.347518538384329112529391120390701166528E6L,
2891 1.514882452993549494932585972882995548426E6L,
2892 1.158019977462989115839826904108208787040E6L,
2893 6.132189329546557743179177159925690841200E5L,
2894 2.248234257620569139969141618556349415120E5L,
2895 5.605842085972455027590989944010492125825E4L,
2896 9.147150349299596453976674231612674085381E3L,
2897 9.104928120962988414618126155557301584078E2L,
2898 4.839208193348159620282142911143429644326E1L,
2899 1.0
2900 ];
2901
2902 // log2 uses the same coefficients as log.
2903 alias log2P = logP;
2904 alias log2Q = logQ;
2905
2906 // log10 uses the same coefficients as log.
2907 alias log10P = logP;
2908 alias log10Q = logQ;
2909
2910 // Coefficients for log(x) = z + z^^3 R(z^^2)/S(z^^2)
2911 // where z = 2(x-1)/(x+1)
2912 // Theoretical peak relative error = 1.1e-35
2913 static immutable real[6] logR = [
2914 1.418134209872192732479751274970992665513E5L,
2915 -8.977257995689735303686582344659576526998E4L,
2916 2.048819892795278657810231591630928516206E4L,
2917 -2.024301798136027039250415126250455056397E3L,
2918 8.057002716646055371965756206836056074715E1L,
2919 -8.828896441624934385266096344596648080902E-1L
2920 ];
2921 static immutable real[7] logS = [
2922 1.701761051846631278975701529965589676574E6L,
2923 -1.332535117259762928288745111081235577029E6L,
2924 4.001557694070773974936904547424676279307E5L,
2925 -5.748542087379434595104154610899551484314E4L,
2926 3.998526750980007367835804959888064681098E3L,
2927 -1.186359407982897997337150403816839480438E2L,
2928 1.0
2929 ];
2930 }
2931 else static if (floatTraits!T.realFormat == RealFormat.ieeeExtended ||
2932 floatTraits!T.realFormat == RealFormat.ieeeExtended53)
2933 {
2934 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
2935 // Theoretical peak relative error = 2.32e-20
2936 static immutable real[7] logP = [
2937 2.0039553499201281259648E1L,
2938 5.7112963590585538103336E1L,
2939 6.0949667980987787057556E1L,
2940 2.9911919328553073277375E1L,
2941 6.5787325942061044846969E0L,
2942 4.9854102823193375972212E-1L,
2943 4.5270000862445199635215E-5L,
2944 ];
2945 static immutable real[7] logQ = [
2946 6.0118660497603843919306E1L,
2947 2.1642788614495947685003E2L,
2948 3.0909872225312059774938E2L,
2949 2.2176239823732856465394E2L,
2950 8.3047565967967209469434E1L,
2951 1.5062909083469192043167E1L,
2952 1.0000000000000000000000E0L,
2953 ];
2954
2955 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
2956 // Theoretical peak relative error = 6.2e-22
2957 static immutable real[7] log2P = [
2958 1.0747524399916215149070E2L,
2959 3.4258224542413922935104E2L,
2960 4.2401812743503691187826E2L,
2961 2.5620629828144409632571E2L,
2962 7.7671073698359539859595E1L,
2963 1.0767376367209449010438E1L,
2964 4.9962495940332550844739E-1L,
2965 ];
2966 static immutable real[8] log2Q = [
2967 3.2242573199748645407652E2L,
2968 1.2695660352705325274404E3L,
2969 2.0307734695595183428202E3L,
2970 1.6911722418503949084863E3L,
2971 7.7952888181207260646090E2L,
2972 1.9444210022760132894510E2L,
2973 2.3479774160285863271658E1L,
2974 1.0000000000000000000000E0,
2975 ];
2976
2977 // log10 uses the same coefficients as log2.
2978 alias log10P = log2P;
2979 alias log10Q = log2Q;
2980
2981 // Coefficients for log(x) = z + z^^3 R(z^^2)/S(z^^2)
2982 // where z = 2(x-1)/(x+1)
2983 // Theoretical peak relative error = 6.16e-22
2984 static immutable real[4] logR = [
2985 -3.5717684488096787370998E1L,
2986 1.0777257190312272158094E1L,
2987 -7.1990767473014147232598E-1L,
2988 1.9757429581415468984296E-3L,
2989 ];
2990 static immutable real[4] logS = [
2991 -4.2861221385716144629696E2L,
2992 1.9361891836232102174846E2L,
2993 -2.6201045551331104417768E1L,
2994 1.0000000000000000000000E0L,
2995 ];
2996 }
2997 else static if (floatTraits!T.realFormat == RealFormat.ieeeDouble)
2998 {
2999 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
3000 static immutable double[6] logP = [
3001 7.70838733755885391666E0,
3002 1.79368678507819816313E1,
3003 1.44989225341610930846E1,
3004 4.70579119878881725854E0,
3005 4.97494994976747001425E-1,
3006 1.01875663804580931796E-4,
3007 ];
3008 static immutable double[6] logQ = [
3009 2.31251620126765340583E1,
3010 7.11544750618563894466E1,
3011 8.29875266912776603211E1,
3012 4.52279145837532221105E1,
3013 1.12873587189167450590E1,
3014 1.00000000000000000000E0,
3015 ];
3016
3017 // log2 uses the same coefficients as log.
3018 alias log2P = logP;
3019 alias log2Q = logQ;
3020
3021 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
3022 static immutable double[7] logp1P = [
3023 2.0039553499201281259648E1,
3024 5.7112963590585538103336E1,
3025 6.0949667980987787057556E1,
3026 2.9911919328553073277375E1,
3027 6.5787325942061044846969E0,
3028 4.9854102823193375972212E-1,
3029 4.5270000862445199635215E-5,
3030 ];
3031 static immutable double[7] logp1Q = [
3032 6.0118660497603843919306E1,
3033 2.1642788614495947685003E2,
3034 3.0909872225312059774938E2,
3035 2.2176239823732856465394E2,
3036 8.3047565967967209469434E1,
3037 1.5062909083469192043167E1,
3038 1.0000000000000000000000E0,
3039 ];
3040
3041 static immutable double[7] log10P = [
3042 1.98892446572874072159E1,
3043 5.67349287391754285487E1,
3044 6.06127134467767258030E1,
3045 2.97877425097986925891E1,
3046 6.56312093769992875930E0,
3047 4.98531067254050724270E-1,
3048 4.58482948458143443514E-5,
3049 ];
3050 static immutable double[7] log10Q = [
3051 5.96677339718622216300E1,
3052 2.14955586696422947765E2,
3053 3.07254189979530058263E2,
3054 2.20664384982121929218E2,
3055 8.27410449222435217021E1,
3056 1.50314182634250003249E1,
3057 1.00000000000000000000E0,
3058 ];
3059
3060 // Coefficients for log(x) = z + z^^3 R(z)/S(z)
3061 // where z = 2(x-1)/(x+1)
3062 static immutable double[3] logR = [
3063 -6.41409952958715622951E1,
3064 1.63866645699558079767E1,
3065 -7.89580278884799154124E-1,
3066 ];
3067 static immutable double[4] logS = [
3068 -7.69691943550460008604E2,
3069 3.12093766372244180303E2,
3070 -3.56722798256324312549E1,
3071 1.00000000000000000000E0,
3072 ];
3073 }
3074 else static if (floatTraits!T.realFormat == RealFormat.ieeeSingle)
3075 {
3076 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)
3077 static immutable float[9] logP = [
3078 3.3333331174E-1,
3079 -2.4999993993E-1,
3080 2.0000714765E-1,
3081 -1.6668057665E-1,
3082 1.4249322787E-1,
3083 -1.2420140846E-1,
3084 1.1676998740E-1,
3085 -1.1514610310E-1,
3086 7.0376836292E-2,
3087 ];
3088
3089 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
3090 static immutable float[7] logp1P = [
3091 2.0039553499E1,
3092 5.7112963590E1,
3093 6.0949667980E1,
3094 2.9911919328E1,
3095 6.5787325942E0,
3096 4.9854102823E-1,
3097 4.5270000862E-5,
3098 ];
3099 static immutable float[7] logp1Q = [
3100 6.01186604976E1,
3101 2.16427886144E2,
3102 3.09098722253E2,
3103 2.21762398237E2,
3104 8.30475659679E1,
3105 1.50629090834E1,
3106 1.00000000000E0,
3107 ];
3108
3109 // log2 and log10 uses the same coefficients as log.
3110 alias log2P = logP;
3111 alias log10P = logP;
3112 }
3113 else
3114 static assert(0, "no coefficients for log()");
3115 }
3116 }
3117
3118 /**************************************
3119 * Calculate the natural logarithm of x.
3120 *
3121 * $(TABLE_SV
3122 * $(TR $(TH x) $(TH log(x)) $(TH divide by 0?) $(TH invalid?))
3123 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
3124 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes))
3125 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no))
3126 * )
3127 */
3128 pragma(inline, true)
3129 real log(real x) @safe pure nothrow @nogc
3130 {
3131 version (INLINE_YL2X)
3132 {
3133 import std.math.constants : LN2;
3134 return core.math.yl2x(x, LN2);
3135 }
3136 else
3137 return logImpl(x);
3138 }
3139
3140 /// ditto
3141 pragma(inline, true)
3142 double log(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log(cast(real) x) : logImpl(x); }
3143
3144 /// ditto
3145 pragma(inline, true)
3146 float log(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log(cast(real) x) : logImpl(x); }
3147
3148 // @@@DEPRECATED_[2.112.0]@@@
3149 deprecated("`std.math.exponential.log` called with argument types `(int)` matches both "
3150 ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.")
3151 real log(int x) @safe pure nothrow @nogc { return log(cast(real) x); }
3152 // @@@DEPRECATED_[2.112.0]@@@
3153 deprecated("`std.math.exponential.log` called with argument types `(uint)` matches both "
3154 ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.")
3155 real log(uint x) @safe pure nothrow @nogc { return log(cast(real) x); }
3156 // @@@DEPRECATED_[2.112.0]@@@
3157 deprecated("`std.math.exponential.log` called with argument types `(long)` matches both "
3158 ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.")
3159 real log(long x) @safe pure nothrow @nogc { return log(cast(real) x); }
3160 // @@@DEPRECATED_[2.112.0]@@@
3161 deprecated("`std.math.exponential.log` called with argument types `(ulong)` matches both "
3162 ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.")
3163 real log(ulong x) @safe pure nothrow @nogc { return log(cast(real) x); }
3164
3165 ///
3166 @safe pure nothrow @nogc unittest
3167 {
3168 import std.math.operations : feqrel;
3169 import std.math.constants : E;
3170
3171 assert(feqrel(log(E), 1) >= real.mant_dig - 1);
3172 }
3173
3174 private T logImpl(T, bool LOG1P = false)(T x) @safe pure nothrow @nogc
3175 {
3176 import std.math.constants : SQRT1_2;
3177 import std.math.algebraic : poly;
3178 import std.math.traits : isInfinity, isNaN, signbit, floatTraits, RealFormat;
3179
3180 alias coeffs = LogCoeffs!T;
3181 alias F = floatTraits!T;
3182
3183 static if (LOG1P)
3184 {
3185 const T xm1 = x;
3186 x = x + 1.0;
3187 }
3188
3189 static if (F.realFormat == RealFormat.ieeeExtended ||
3190 F.realFormat == RealFormat.ieeeExtended53 ||
3191 F.realFormat == RealFormat.ieeeQuadruple)
3192 {
3193 // C1 + C2 = LN2.
3194 enum T C1 = 6.93145751953125E-1L;
3195 enum T C2 = 1.428606820309417232121458176568075500134E-6L;
3196 }
3197 else static if (F.realFormat == RealFormat.ieeeDouble)
3198 {
3199 enum T C1 = 0.693359375;
3200 enum T C2 = -2.121944400546905827679e-4;
3201 }
3202 else static if (F.realFormat == RealFormat.ieeeSingle)
3203 {
3204 enum T C1 = 0.693359375;
3205 enum T C2 = -2.12194440e-4;
3206 }
3207 else
3208 static assert(0, "Not implemented for this architecture");
3209
3210 // Special cases.
3211 if (isNaN(x))
3212 return x;
3213 if (isInfinity(x) && !signbit(x))
3214 return x;
3215 if (x == 0.0)
3216 return -T.infinity;
3217 if (x < 0.0)
3218 return T.nan;
3219
3220 // Separate mantissa from exponent.
3221 // Note, frexp is used so that denormal numbers will be handled properly.
3222 T y, z;
3223 int exp;
3224
3225 x = frexp(x, exp);
3226
3227 static if (F.realFormat == RealFormat.ieeeDouble ||
3228 F.realFormat == RealFormat.ieeeExtended ||
3229 F.realFormat == RealFormat.ieeeExtended53 ||
3230 F.realFormat == RealFormat.ieeeQuadruple)
3231 {
3232 // Logarithm using log(x) = z + z^^3 R(z) / S(z),
3233 // where z = 2(x - 1)/(x + 1)
3234 if ((exp > 2) || (exp < -2))
3235 {
3236 if (x < SQRT1_2)
3237 { // 2(2x - 1)/(2x + 1)
3238 exp -= 1;
3239 z = x - 0.5;
3240 y = 0.5 * z + 0.5;
3241 }
3242 else
3243 { // 2(x - 1)/(x + 1)
3244 z = x - 0.5;
3245 z -= 0.5;
3246 y = 0.5 * x + 0.5;
3247 }
3248 x = z / y;
3249 z = x * x;
3250 z = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS));
3251 z += exp * C2;
3252 z += x;
3253 z += exp * C1;
3254
3255 return z;
3256 }
3257 }
3258
3259 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
3260 if (x < SQRT1_2)
3261 {
3262 exp -= 1;
3263 static if (LOG1P)
3264 {
3265 if (exp != 0)
3266 x = 2.0 * x - 1.0;
3267 else
3268 x = xm1;
3269 }
3270 else
3271 x = 2.0 * x - 1.0;
3272
3273 }
3274 else
3275 {
3276 static if (LOG1P)
3277 {
3278 if (exp != 0)
3279 x = x - 1.0;
3280 else
3281 x = xm1;
3282 }
3283 else
3284 x = x - 1.0;
3285 }
3286 z = x * x;
3287 static if (F.realFormat == RealFormat.ieeeSingle)
3288 y = x * (z * poly(x, coeffs.logP));
3289 else
3290 y = x * (z * poly(x, coeffs.logP) / poly(x, coeffs.logQ));
3291 y += exp * C2;
3292 z = y - 0.5 * z;
3293
3294 // Note, the sum of above terms does not exceed x/4,
3295 // so it contributes at most about 1/4 lsb to the error.
3296 z += x;
3297 z += exp * C1;
3298
3299 return z;
3300 }
3301
3302 @safe @nogc nothrow unittest
3303 {
3304 import std.math.traits : floatTraits, RealFormat;
3305 import std.meta : AliasSeq;
3306
3307 static void testLog(T)(T[2][] vals)
3308 {
3309 import std.math.operations : isClose;
3310 import std.math.traits : isNaN;
3311 foreach (ref pair; vals)
3312 {
3313 if (isNaN(pair[1]))
3314 assert(isNaN(log(pair[0])));
3315 else
3316 assert(isClose(log(pair[0]), pair[1]));
3317 }
3318 }
3319 static foreach (F; AliasSeq!(float, double, real))
3320 {{
3321 scope F[2][] vals = [
3322 [F(1), F(0x0p+0)], [F(2), F(0x1.62e42fefa39ef358p-1)],
3323 [F(4), F(0x1.62e42fefa39ef358p+0)], [F(8), F(0x1.0a2b23f3bab73682p+1)],
3324 [F(16), F(0x1.62e42fefa39ef358p+1)], [F(32), F(0x1.bb9d3beb8c86b02ep+1)],
3325 [F(64), F(0x1.0a2b23f3bab73682p+2)], [F(128), F(0x1.3687a9f1af2b14ecp+2)],
3326 [F(256), F(0x1.62e42fefa39ef358p+2)], [F(512), F(0x1.8f40b5ed9812d1c2p+2)],
3327 [F(1024), F(0x1.bb9d3beb8c86b02ep+2)], [F(2048), F(0x1.e7f9c1e980fa8e98p+2)],
3328 [F(3), F(0x1.193ea7aad030a976p+0)], [F(5), F(0x1.9c041f7ed8d336bp+0)],
3329 [F(7), F(0x1.f2272ae325a57546p+0)], [F(15), F(0x1.5aa16394d481f014p+1)],
3330 [F(17), F(0x1.6aa6bc1fa7f79cfp+1)], [F(31), F(0x1.b78ce48912b59f12p+1)],
3331 [F(33), F(0x1.bf8d8f4d5b8d1038p+1)], [F(63), F(0x1.09291e8e3181b20ep+2)],
3332 [F(65), F(0x1.0b292939429755ap+2)], [F(-0), -F.infinity], [F(0), -F.infinity],
3333 [F(0.1), F(-0x1.26bb1bbb5551582ep+1)], [F(0.25), F(-0x1.62e42fefa39ef358p+0)],
3334 [F(0.75), F(-0x1.269621134db92784p-2)], [F(0.875), F(-0x1.1178e8227e47bde4p-3)],
3335 [F(10), F(0x1.26bb1bbb5551582ep+1)], [F(100), F(0x1.26bb1bbb5551582ep+2)],
3336 [F(10000), F(0x1.26bb1bbb5551582ep+3)],
3337 ];
3338 testLog(vals);
3339 }}
3340 {
3341 float[2][16] vals = [
3342 [float.nan, float.nan],[-float.nan, float.nan],
3343 [float.infinity, float.infinity], [-float.infinity, float.nan],
3344 [float.min_normal, -0x1.5d58ap+6f], [-float.min_normal, float.nan],
3345 [float.max, 0x1.62e43p+6f], [-float.max, float.nan],
3346 [float.min_normal / 2, -0x1.601e68p+6f], [-float.min_normal / 2, float.nan],
3347 [float.max / 2, 0x1.601e68p+6f], [-float.max / 2, float.nan],
3348 [float.min_normal / 3, -0x1.61bd9ap+6f], [-float.min_normal / 3, float.nan],
3349 [float.max / 3, 0x1.5e7f36p+6f], [-float.max / 3, float.nan],
3350 ];
3351 testLog(vals);
3352 }
3353 {
3354 double[2][16] vals = [
3355 [double.nan, double.nan],[-double.nan, double.nan],
3356 [double.infinity, double.infinity], [-double.infinity, double.nan],
3357 [double.min_normal, -0x1.6232bdd7abcd2p+9], [-double.min_normal, double.nan],
3358 [double.max, 0x1.62e42fefa39efp+9], [-double.max, double.nan],
3359 [double.min_normal / 2, -0x1.628b76e3a7b61p+9], [-double.min_normal / 2, double.nan],
3360 [double.max / 2, 0x1.628b76e3a7b61p+9], [-double.max / 2, double.nan],
3361 [double.min_normal / 3, -0x1.62bf5d2b81354p+9], [-double.min_normal / 3, double.nan],
3362 [double.max / 3, 0x1.6257909bce36ep+9], [-double.max / 3, double.nan],
3363 ];
3364 testLog(vals);
3365 }
3366 alias F = floatTraits!real;
3367 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
3368 {{
3369 real[2][16] vals = [
3370 [real.nan, real.nan],[-real.nan, real.nan],
3371 [real.infinity, real.infinity], [-real.infinity, real.nan],
3372 [real.min_normal, -0x1.62d918ce2421d66p+13L], [-real.min_normal, real.nan],
3373 [real.max, 0x1.62e42fefa39ef358p+13L], [-real.max, real.nan],
3374 [real.min_normal / 2, -0x1.62dea45ee3e064dcp+13L], [-real.min_normal / 2, real.nan],
3375 [real.max / 2, 0x1.62dea45ee3e064dcp+13L], [-real.max / 2, real.nan],
3376 [real.min_normal / 3, -0x1.62e1e2c3617857e6p+13L], [-real.min_normal / 3, real.nan],
3377 [real.max / 3, 0x1.62db65fa664871d2p+13L], [-real.max / 3, real.nan],
3378 ];
3379 testLog(vals);
3380 }}
3381 }
3382
3383 /**************************************
3384 * Calculate the base-10 logarithm of x.
3385 *
3386 * $(TABLE_SV
3387 * $(TR $(TH x) $(TH log10(x)) $(TH divide by 0?) $(TH invalid?))
3388 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
3389 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes))
3390 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no))
3391 * )
3392 */
3393 pragma(inline, true)
3394 real log10(real x) @safe pure nothrow @nogc
3395 {
3396 version (INLINE_YL2X)
3397 {
3398 import std.math.constants : LOG2;
3399 return core.math.yl2x(x, LOG2);
3400 }
3401 else
3402 return log10Impl(x);
3403 }
3404
3405 /// ditto
3406 pragma(inline, true)
3407 double log10(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log10(cast(real) x) : log10Impl(x); }
3408
3409 /// ditto
3410 pragma(inline, true)
3411 float log10(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log10(cast(real) x) : log10Impl(x); }
3412
3413 // @@@DEPRECATED_[2.112.0]@@@
3414 deprecated("`std.math.exponential.log10` called with argument types `(int)` matches both "
3415 ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.")
3416 real log10(int x) @safe pure nothrow @nogc { return log10(cast(real) x); }
3417 // @@@DEPRECATED_[2.112.0]@@@
3418 deprecated("`std.math.exponential.log10` called with argument types `(uint)` matches both "
3419 ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.")
3420 real log10(uint x) @safe pure nothrow @nogc { return log10(cast(real) x); }
3421 // @@@DEPRECATED_[2.112.0]@@@
3422 deprecated("`std.math.exponential.log10` called with argument types `(long)` matches both "
3423 ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.")
3424 real log10(long x) @safe pure nothrow @nogc { return log10(cast(real) x); }
3425 // @@@DEPRECATED_[2.112.0]@@@
3426 deprecated("`std.math.exponential.log10` called with argument types `(ulong)` matches both "
3427 ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.")
3428 real log10(ulong x) @safe pure nothrow @nogc { return log10(cast(real) x); }
3429
3430 ///
3431 @safe pure nothrow @nogc unittest
3432 {
3433 import std.math.algebraic : fabs;
3434
3435 assert(fabs(log10(1000.0L) - 3) < .000001);
3436 }
3437
3438 @safe pure nothrow @nogc unittest
3439 {
3440 import std.math.algebraic : fabs;
3441
3442 assert(fabs(log10(1000.0) - 3) < .000001);
3443 assert(fabs(log10(1000.0f) - 3) < .000001);
3444 }
3445
3446 private T log10Impl(T)(T x) @safe pure nothrow @nogc
3447 {
3448 import std.math.constants : SQRT1_2;
3449 import std.math.algebraic : poly;
3450 import std.math.traits : isNaN, isInfinity, signbit, floatTraits, RealFormat;
3451
3452 alias coeffs = LogCoeffs!T;
3453 alias F = floatTraits!T;
3454
3455 static if (F.realFormat == RealFormat.ieeeExtended ||
3456 F.realFormat == RealFormat.ieeeExtended53 ||
3457 F.realFormat == RealFormat.ieeeQuadruple)
3458 {
3459 // log10(2) split into two parts.
3460 enum T L102A = 0.3125L;
3461 enum T L102B = -1.14700043360188047862611052755069732318101185E-2L;
3462
3463 // log10(e) split into two parts.
3464 enum T L10EA = 0.5L;
3465 enum T L10EB = -6.570551809674817234887108108339491770560299E-2L;
3466 }
3467 else static if (F.realFormat == RealFormat.ieeeDouble ||
3468 F.realFormat == RealFormat.ieeeSingle)
3469 {
3470 enum T L102A = 3.0078125E-1;
3471 enum T L102B = 2.48745663981195213739E-4;
3472
3473 enum T L10EA = 4.3359375E-1;
3474 enum T L10EB = 7.00731903251827651129E-4;
3475 }
3476 else
3477 static assert(0, "Not implemented for this architecture");
3478
3479 // Special cases are the same as for log.
3480 if (isNaN(x))
3481 return x;
3482 if (isInfinity(x) && !signbit(x))
3483 return x;
3484 if (x == 0.0)
3485 return -T.infinity;
3486 if (x < 0.0)
3487 return T.nan;
3488
3489 // Separate mantissa from exponent.
3490 // Note, frexp is used so that denormal numbers will be handled properly.
3491 T y, z;
3492 int exp;
3493
3494 x = frexp(x, exp);
3495
3496 static if (F.realFormat == RealFormat.ieeeExtended ||
3497 F.realFormat == RealFormat.ieeeExtended53 ||
3498 F.realFormat == RealFormat.ieeeQuadruple)
3499 {
3500 // Logarithm using log(x) = z + z^^3 R(z) / S(z),
3501 // where z = 2(x - 1)/(x + 1)
3502 if ((exp > 2) || (exp < -2))
3503 {
3504 if (x < SQRT1_2)
3505 { // 2(2x - 1)/(2x + 1)
3506 exp -= 1;
3507 z = x - 0.5;
3508 y = 0.5 * z + 0.5;
3509 }
3510 else
3511 { // 2(x - 1)/(x + 1)
3512 z = x - 0.5;
3513 z -= 0.5;
3514 y = 0.5 * x + 0.5;
3515 }
3516 x = z / y;
3517 z = x * x;
3518 y = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS));
3519 goto Ldone;
3520 }
3521 }
3522
3523 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
3524 if (x < SQRT1_2)
3525 {
3526 exp -= 1;
3527 x = 2.0 * x - 1.0;
3528 }
3529 else
3530 x = x - 1.0;
3531
3532 z = x * x;
3533 static if (F.realFormat == RealFormat.ieeeSingle)
3534 y = x * (z * poly(x, coeffs.log10P));
3535 else
3536 y = x * (z * poly(x, coeffs.log10P) / poly(x, coeffs.log10Q));
3537 y = y - 0.5 * z;
3538
3539 // Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
3540 // This sequence of operations is critical and it may be horribly
3541 // defeated by some compiler optimizers.
3542 Ldone:
3543 z = y * L10EB;
3544 z += x * L10EB;
3545 z += exp * L102B;
3546 z += y * L10EA;
3547 z += x * L10EA;
3548 z += exp * L102A;
3549
3550 return z;
3551 }
3552
3553 @safe @nogc nothrow unittest
3554 {
3555 import std.math.traits : floatTraits, RealFormat;
3556 import std.meta : AliasSeq;
3557
3558 static void testLog10(T)(T[2][] vals)
3559 {
3560 import std.math.operations : isClose;
3561 import std.math.traits : isNaN;
3562 foreach (ref pair; vals)
3563 {
3564 if (isNaN(pair[1]))
3565 assert(isNaN(log10(pair[0])));
3566 else
3567 assert(isClose(log10(pair[0]), pair[1]));
3568 }
3569 }
3570 static foreach (F; AliasSeq!(float, double, real))
3571 {{
3572 scope F[2][] vals = [
3573 [F(1), F(0x0p+0)], [F(2), F(0x1.34413509f79fef32p-2)],
3574 [F(4), F(0x1.34413509f79fef32p-1)], [F(8), F(0x1.ce61cf8ef36fe6cap-1)],
3575 [F(16), F(0x1.34413509f79fef32p+0)], [F(32), F(0x1.8151824c7587eafep+0)],
3576 [F(64), F(0x1.ce61cf8ef36fe6cap+0)], [F(128), F(0x1.0db90e68b8abf14cp+1)],
3577 [F(256), F(0x1.34413509f79fef32p+1)], [F(512), F(0x1.5ac95bab3693ed18p+1)],
3578 [F(1024), F(0x1.8151824c7587eafep+1)], [F(2048), F(0x1.a7d9a8edb47be8e4p+1)],
3579 [F(3), F(0x1.e8927964fd5fd08cp-2)], [F(5), F(0x1.65df657b04300868p-1)],
3580 [F(7), F(0x1.b0b0b0b78cc3f296p-1)], [F(15), F(0x1.2d145116c16ff856p+0)],
3581 [F(17), F(0x1.3afeb354b7d9731ap+0)], [F(31), F(0x1.7dc9e145867e62eap+0)],
3582 [F(33), F(0x1.84bd545e4baeddp+0)], [F(63), F(0x1.cca1950e4511e192p+0)],
3583 [F(65), F(0x1.d01b16f9433cf7b8p+0)], [F(-0), -F.infinity], [F(0), -F.infinity],
3584 [F(0.1), F(-0x1p+0)], [F(0.25), F(-0x1.34413509f79fef32p-1)],
3585 [F(0.75), F(-0x1.ffbfc2bbc7803758p-4)], [F(0.875), F(-0x1.db11ed766abf432ep-5)],
3586 [F(10), F(0x1p+0)], [F(100), F(0x1p+1)], [F(10000), F(0x1p+2)],
3587 ];
3588 testLog10(vals);
3589 }}
3590 {
3591 float[2][16] vals = [
3592 [float.nan, float.nan],[-float.nan, float.nan],
3593 [float.infinity, float.infinity], [-float.infinity, float.nan],
3594 [float.min_normal, -0x1.2f703p+5f], [-float.min_normal, float.nan],
3595 [float.max, 0x1.344136p+5f], [-float.max, float.nan],
3596 [float.min_normal / 2, -0x1.31d8b2p+5f], [-float.min_normal / 2, float.nan],
3597 [float.max / 2, 0x1.31d8b2p+5f], [-float.max / 2, float.nan],
3598 [float.min_normal / 3, -0x1.334156p+5f], [-float.min_normal / 3, float.nan],
3599 [float.max / 3, 0x1.30701p+5f], [-float.max / 3, float.nan],
3600 ];
3601 testLog10(vals);
3602 }
3603 {
3604 double[2][16] vals = [
3605 [double.nan, double.nan],[-double.nan, double.nan],
3606 [double.infinity, double.infinity], [-double.infinity, double.nan],
3607 [double.min_normal, -0x1.33a7146f72a42p+8], [-double.min_normal, double.nan],
3608 [double.max, 0x1.34413509f79ffp+8], [-double.max, double.nan],
3609 [double.min_normal / 2, -0x1.33f424bcb522p+8], [-double.min_normal / 2, double.nan],
3610 [double.max / 2, 0x1.33f424bcb522p+8], [-double.max / 2, double.nan],
3611 [double.min_normal / 3, -0x1.3421390dcbe37p+8], [-double.min_normal / 3, double.nan],
3612 [double.max / 3, 0x1.33c7106b9e609p+8], [-double.max / 3, double.nan],
3613 ];
3614 testLog10(vals);
3615 }
3616 alias F = floatTraits!real;
3617 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
3618 {{
3619 real[2][16] vals = [
3620 [real.nan, real.nan],[-real.nan, real.nan],
3621 [real.infinity, real.infinity], [-real.infinity, real.nan],
3622 [real.min_normal, -0x1.343793004f503232p+12L], [-real.min_normal, real.nan],
3623 [real.max, 0x1.34413509f79fef32p+12L], [-real.max, real.nan],
3624 [real.min_normal / 2, -0x1.343c6405237810b2p+12L], [-real.min_normal / 2, real.nan],
3625 [real.max / 2, 0x1.343c6405237810b2p+12L], [-real.max / 2, real.nan],
3626 [real.min_normal / 3, -0x1.343f354a34e427bp+12L], [-real.min_normal / 3, real.nan],
3627 [real.max / 3, 0x1.343992c0120bf9b2p+12L], [-real.max / 3, real.nan],
3628 ];
3629 testLog10(vals);
3630 }}
3631 }
3632
3633 /**
3634 * Calculates the natural logarithm of 1 + x.
3635 *
3636 * For very small x, log1p(x) will be more accurate than
3637 * log(1 + x).
3638 *
3639 * $(TABLE_SV
3640 * $(TR $(TH x) $(TH log1p(x)) $(TH divide by 0?) $(TH invalid?))
3641 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) $(TD no))
3642 * $(TR $(TD -1.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
3643 * $(TR $(TD $(LT)-1.0) $(TD -$(NAN)) $(TD no) $(TD yes))
3644 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no))
3645 * )
3646 */
3647 pragma(inline, true)
3648 real log1p(real x) @safe pure nothrow @nogc
3649 {
3650 version (INLINE_YL2X)
3651 {
3652 // On x87, yl2xp1 is valid if and only if -0.5 <= lg(x) <= 0.5,
3653 // ie if -0.29 <= x <= 0.414
3654 import std.math.constants : LN2;
3655 return (core.math.fabs(x) <= 0.25) ? core.math.yl2xp1(x, LN2) : core.math.yl2x(x+1, LN2);
3656 }
3657 else
3658 return log1pImpl(x);
3659 }
3660
3661 /// ditto
3662 pragma(inline, true)
3663 double log1p(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log1p(cast(real) x) : log1pImpl(x); }
3664
3665 /// ditto
3666 pragma(inline, true)
3667 float log1p(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log1p(cast(real) x) : log1pImpl(x); }
3668
3669 // @@@DEPRECATED_[2.112.0]@@@
3670 deprecated("`std.math.exponential.log1p` called with argument types `(int)` matches both "
3671 ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.")
3672 real log1p(int x) @safe pure nothrow @nogc { return log1p(cast(real) x); }
3673 // @@@DEPRECATED_[2.112.0]@@@
3674 deprecated("`std.math.exponential.log1p` called with argument types `(uint)` matches both "
3675 ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.")
3676 real log1p(uint x) @safe pure nothrow @nogc { return log1p(cast(real) x); }
3677 // @@@DEPRECATED_[2.112.0]@@@
3678 deprecated("`std.math.exponential.log1p` called with argument types `(long)` matches both "
3679 ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.")
3680 real log1p(long x) @safe pure nothrow @nogc { return log1p(cast(real) x); }
3681 // @@@DEPRECATED_[2.112.0]@@@
3682 deprecated("`std.math.exponential.log1p` called with argument types `(ulong)` matches both "
3683 ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.")
3684 real log1p(ulong x) @safe pure nothrow @nogc { return log1p(cast(real) x); }
3685
3686 ///
3687 @safe pure unittest
3688 {
3689 import std.math.traits : isIdentical, isNaN;
3690 import std.math.operations : feqrel;
3691
3692 assert(isIdentical(log1p(0.0), 0.0));
3693 assert(log1p(1.0).feqrel(0.69314) > 16);
3694
3695 assert(log1p(-1.0) == -real.infinity);
3696 assert(isNaN(log1p(-2.0)));
3697 assert(log1p(real.nan) is real.nan);
3698 assert(log1p(-real.nan) is -real.nan);
3699 assert(log1p(real.infinity) == real.infinity);
3700 }
3701
3702 private T log1pImpl(T)(T x) @safe pure nothrow @nogc
3703 {
3704 import std.math.traits : isNaN, isInfinity, signbit;
3705 import std.math.algebraic : poly;
3706 import std.math.constants : SQRT1_2, SQRT2;
3707 import std.math.traits : floatTraits, RealFormat;
3708
3709 // Special cases.
3710 if (isNaN(x) || x == 0.0)
3711 return x;
3712 if (isInfinity(x) && !signbit(x))
3713 return x;
3714 if (x == -1.0)
3715 return -T.infinity;
3716 if (x < -1.0)
3717 return T.nan;
3718
3719 alias F = floatTraits!T;
3720 static if (F.realFormat == RealFormat.ieeeSingle ||
3721 F.realFormat == RealFormat.ieeeDouble)
3722 {
3723 // When the input is within the range 1/sqrt(2) <= x+1 <= sqrt(2), compute
3724 // log1p inline. Forwarding to log() would otherwise result in inaccuracies.
3725 const T xp1 = x + 1.0;
3726 if (xp1 >= SQRT1_2 && xp1 <= SQRT2)
3727 {
3728 alias coeffs = LogCoeffs!T;
3729
3730 T px = poly(x, coeffs.logp1P);
3731 T qx = poly(x, coeffs.logp1Q);
3732 const T xx = x * x;
3733 qx = x + ((cast(T) -0.5) * xx + x * (xx * px / qx));
3734 return qx;
3735 }
3736 }
3737
3738 return logImpl!(T, true)(x);
3739 }
3740
3741 @safe @nogc nothrow unittest
3742 {
3743 import std.math.traits : floatTraits, RealFormat;
3744 import std.meta : AliasSeq;
3745
3746 static void testLog1p(T)(T[2][] vals)
3747 {
3748 import std.math.operations : isClose;
3749 import std.math.traits : isNaN;
3750 foreach (ref pair; vals)
3751 {
3752 if (isNaN(pair[1]))
3753 assert(isNaN(log1p(pair[0])));
3754 else
3755 assert(isClose(log1p(pair[0]), pair[1]));
3756 }
3757 }
3758 static foreach (F; AliasSeq!(float, double, real))
3759 {{
3760 scope F[2][] vals = [
3761 [F(1), F(0x1.62e42fefa39ef358p-1)], [F(2), F(0x1.193ea7aad030a976p+0)],
3762 [F(4), F(0x1.9c041f7ed8d336bp+0)], [F(8), F(0x1.193ea7aad030a976p+1)],
3763 [F(16), F(0x1.6aa6bc1fa7f79cfp+1)], [F(32), F(0x1.bf8d8f4d5b8d1038p+1)],
3764 [F(64), F(0x1.0b292939429755ap+2)], [F(128), F(0x1.37072a9b5b6cb31p+2)],
3765 [F(256), F(0x1.63241004e9010ad8p+2)], [F(512), F(0x1.8f60adf041bde2a8p+2)],
3766 [F(1024), F(0x1.bbad39ebe1cc08b6p+2)], [F(2048), F(0x1.e801c1698ba4395cp+2)],
3767 [F(3), F(0x1.62e42fefa39ef358p+0)], [F(5), F(0x1.cab0bfa2a2002322p+0)],
3768 [F(7), F(0x1.0a2b23f3bab73682p+1)], [F(15), F(0x1.62e42fefa39ef358p+1)],
3769 [F(17), F(0x1.71f7b3a6b918664cp+1)], [F(31), F(0x1.bb9d3beb8c86b02ep+1)],
3770 [F(33), F(0x1.c35fc81b90df59c6p+1)], [F(63), F(0x1.0a2b23f3bab73682p+2)],
3771 [F(65), F(0x1.0c234da4a23a6686p+2)], [F(-0), F(-0x0p+0)], [F(0), F(0x0p+0)],
3772 [F(0.1), F(0x1.8663f793c46c69cp-4)], [F(0.25), F(0x1.c8ff7c79a9a21ac4p-3)],
3773 [F(0.75), F(0x1.1e85f5e7040d03ep-1)], [F(0.875), F(0x1.41d8fe84672ae646p-1)],
3774 [F(10), F(0x1.32ee3b77f374bb7cp+1)], [F(100), F(0x1.275e2271bba30be4p+2)],
3775 [F(10000), F(0x1.26bbed6fbd84182ep+3)],
3776 ];
3777 testLog1p(vals);
3778 }}
3779 {
3780 float[2][16] vals = [
3781 [float.nan, float.nan],[-float.nan, float.nan],
3782 [float.infinity, float.infinity], [-float.infinity, float.nan],
3783 [float.min_normal, 0x1p-126f], [-float.min_normal, -0x1p-126f],
3784 [float.max, 0x1.62e43p+6f], [-float.max, float.nan],
3785 [float.min_normal / 2, 0x0.8p-126f], [-float.min_normal / 2, -0x0.8p-126f],
3786 [float.max / 2, 0x1.601e68p+6f], [-float.max / 2, float.nan],
3787 [float.min_normal / 3, 0x0.555556p-126f], [-float.min_normal / 3, -0x0.555556p-126f],
3788 [float.max / 3, 0x1.5e7f36p+6f], [-float.max / 3, float.nan],
3789 ];
3790 testLog1p(vals);
3791 }
3792 {
3793 double[2][16] vals = [
3794 [double.nan, double.nan],[-double.nan, double.nan],
3795 [double.infinity, double.infinity], [-double.infinity, double.nan],
3796 [double.min_normal, 0x1p-1022], [-double.min_normal, -0x1p-1022],
3797 [double.max, 0x1.62e42fefa39efp+9], [-double.max, double.nan],
3798 [double.min_normal / 2, 0x0.8p-1022], [-double.min_normal / 2, -0x0.8p-1022],
3799 [double.max / 2, 0x1.628b76e3a7b61p+9], [-double.max / 2, double.nan],
3800 [double.min_normal / 3, 0x0.5555555555555p-1022], [-double.min_normal / 3, -0x0.5555555555555p-1022],
3801 [double.max / 3, 0x1.6257909bce36ep+9], [-double.max / 3, double.nan],
3802 ];
3803 testLog1p(vals);
3804 }
3805 alias F = floatTraits!real;
3806 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
3807 {{
3808 real[2][16] vals = [
3809 [real.nan, real.nan],[-real.nan, real.nan],
3810 [real.infinity, real.infinity], [-real.infinity, real.nan],
3811 [real.min_normal, 0x1p-16382L], [-real.min_normal, -0x1p-16382L],
3812 [real.max, 0x1.62e42fefa39ef358p+13L], [-real.max, real.nan],
3813 [real.min_normal / 2, 0x0.8p-16382L], [-real.min_normal / 2, -0x0.8p-16382L],
3814 [real.max / 2, 0x1.62dea45ee3e064dcp+13L], [-real.max / 2, real.nan],
3815 [real.min_normal / 3, 0x0.5555555555555556p-16382L], [-real.min_normal / 3, -0x0.5555555555555556p-16382L],
3816 [real.max / 3, 0x1.62db65fa664871d2p+13L], [-real.max / 3, real.nan],
3817 ];
3818 testLog1p(vals);
3819 }}
3820 }
3821
3822 /***************************************
3823 * Calculates the base-2 logarithm of x:
3824 * $(SUB log, 2)x
3825 *
3826 * $(TABLE_SV
3827 * $(TR $(TH x) $(TH log2(x)) $(TH divide by 0?) $(TH invalid?))
3828 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no) )
3829 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes) )
3830 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no) )
3831 * )
3832 */
3833 pragma(inline, true)
3834 real log2(real x) @safe pure nothrow @nogc
3835 {
3836 version (INLINE_YL2X)
3837 return core.math.yl2x(x, 1.0L);
3838 else
3839 return log2Impl(x);
3840 }
3841
3842 /// ditto
3843 pragma(inline, true)
3844 double log2(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log2(cast(real) x) : log2Impl(x); }
3845
3846 /// ditto
3847 pragma(inline, true)
3848 float log2(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log2(cast(real) x) : log2Impl(x); }
3849
3850 // @@@DEPRECATED_[2.112.0]@@@
3851 deprecated("`std.math.exponential.log2` called with argument types `(int)` matches both "
3852 ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.")
3853 real log2(int x) @safe pure nothrow @nogc { return log2(cast(real) x); }
3854 // @@@DEPRECATED_[2.112.0]@@@
3855 deprecated("`std.math.exponential.log2` called with argument types `(uint)` matches both "
3856 ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.")
3857 real log2(uint x) @safe pure nothrow @nogc { return log2(cast(real) x); }
3858 // @@@DEPRECATED_[2.112.0]@@@
3859 deprecated("`std.math.exponential.log2` called with argument types `(long)` matches both "
3860 ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.")
3861 real log2(long x) @safe pure nothrow @nogc { return log2(cast(real) x); }
3862 // @@@DEPRECATED_[2.112.0]@@@
3863 deprecated("`std.math.exponential.log2` called with argument types `(ulong)` matches both "
3864 ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.")
3865 real log2(ulong x) @safe pure nothrow @nogc { return log2(cast(real) x); }
3866
3867 ///
3868 @safe unittest
3869 {
3870 import std.math.operations : isClose;
3871
3872 assert(isClose(log2(1024.0L), 10));
3873 }
3874
3875 @safe @nogc nothrow unittest
3876 {
3877 import std.math.operations : isClose;
3878
3879 // check if values are equal to 19 decimal digits of precision
3880 assert(isClose(log2(1024.0L), 10, 1e-18));
3881 }
3882
3883 private T log2Impl(T)(T x) @safe pure nothrow @nogc
3884 {
3885 import std.math.traits : isNaN, isInfinity, signbit;
3886 import std.math.constants : SQRT1_2, LOG2E;
3887 import std.math.algebraic : poly;
3888 import std.math.traits : floatTraits, RealFormat;
3889
3890 alias coeffs = LogCoeffs!T;
3891 alias F = floatTraits!T;
3892
3893 // Special cases are the same as for log.
3894 if (isNaN(x))
3895 return x;
3896 if (isInfinity(x) && !signbit(x))
3897 return x;
3898 if (x == 0.0)
3899 return -T.infinity;
3900 if (x < 0.0)
3901 return T.nan;
3902
3903 // Separate mantissa from exponent.
3904 // Note, frexp is used so that denormal numbers will be handled properly.
3905 T y, z;
3906 int exp;
3907
3908 x = frexp(x, exp);
3909
3910 static if (F.realFormat == RealFormat.ieeeDouble ||
3911 F.realFormat == RealFormat.ieeeExtended ||
3912 F.realFormat == RealFormat.ieeeExtended53 ||
3913 F.realFormat == RealFormat.ieeeQuadruple)
3914 {
3915 // Logarithm using log(x) = z + z^^3 R(z) / S(z),
3916 // where z = 2(x - 1)/(x + 1)
3917 if ((exp > 2) || (exp < -2))
3918 {
3919 if (x < SQRT1_2)
3920 { // 2(2x - 1)/(2x + 1)
3921 exp -= 1;
3922 z = x - 0.5;
3923 y = 0.5 * z + 0.5;
3924 }
3925 else
3926 { // 2(x - 1)/(x + 1)
3927 z = x - 0.5;
3928 z -= 0.5;
3929 y = 0.5 * x + 0.5;
3930 }
3931 x = z / y;
3932 z = x * x;
3933 y = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS));
3934 goto Ldone;
3935 }
3936 }
3937
3938 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
3939 if (x < SQRT1_2)
3940 {
3941 exp -= 1;
3942 x = 2.0 * x - 1.0;
3943 }
3944 else
3945 x = x - 1.0;
3946
3947 z = x * x;
3948 static if (F.realFormat == RealFormat.ieeeSingle)
3949 y = x * (z * poly(x, coeffs.log2P));
3950 else
3951 y = x * (z * poly(x, coeffs.log2P) / poly(x, coeffs.log2Q));
3952 y = y - 0.5 * z;
3953
3954 // Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
3955 // This sequence of operations is critical and it may be horribly
3956 // defeated by some compiler optimizers.
3957 Ldone:
3958 z = y * (LOG2E - 1.0);
3959 z += x * (LOG2E - 1.0);
3960 z += y;
3961 z += x;
3962 z += exp;
3963
3964 return z;
3965 }
3966
3967 @safe @nogc nothrow unittest
3968 {
3969 import std.math.traits : floatTraits, RealFormat;
3970 import std.meta : AliasSeq;
3971
3972 static void testLog2(T)(T[2][] vals)
3973 {
3974 import std.math.operations : isClose;
3975 import std.math.traits : isNaN;
3976 foreach (ref pair; vals)
3977 {
3978 if (isNaN(pair[1]))
3979 assert(isNaN(log2(pair[0])));
3980 else
3981 assert(isClose(log2(pair[0]), pair[1]));
3982 }
3983 }
3984 static foreach (F; AliasSeq!(float, double, real))
3985 {{
3986 scope F[2][] vals = [
3987 [F(1), F(0x0p+0)], [F(2), F(0x1p+0)],
3988 [F(4), F(0x1p+1)], [F(8), F(0x1.8p+1)],
3989 [F(16), F(0x1p+2)], [F(32), F(0x1.4p+2)],
3990 [F(64), F(0x1.8p+2)], [F(128), F(0x1.cp+2)],
3991 [F(256), F(0x1p+3)], [F(512), F(0x1.2p+3)],
3992 [F(1024), F(0x1.4p+3)], [F(2048), F(0x1.6p+3)],
3993 [F(3), F(0x1.95c01a39fbd687ap+0)], [F(5), F(0x1.2934f0979a3715fcp+1)],
3994 [F(7), F(0x1.675767f54042cd9ap+1)], [F(15), F(0x1.f414fdb4982259ccp+1)],
3995 [F(17), F(0x1.0598fdbeb244c5ap+2)], [F(31), F(0x1.3d118d66c4d4e554p+2)],
3996 [F(33), F(0x1.42d75a6eb1dfb0e6p+2)], [F(63), F(0x1.7e8bc1179e0caa9cp+2)],
3997 [F(65), F(0x1.816e79685c2d2298p+2)], [F(-0), -F.infinity], [F(0), -F.infinity],
3998 [F(0.1), F(-0x1.a934f0979a3715fcp+1)], [F(0.25), F(-0x1p+1)],
3999 [F(0.75), F(-0x1.a8ff971810a5e182p-2)], [F(0.875), F(-0x1.8a8980abfbd32668p-3)],
4000 [F(10), F(0x1.a934f0979a3715fcp+1)], [F(100), F(0x1.a934f0979a3715fcp+2)],
4001 [F(10000), F(0x1.a934f0979a3715fcp+3)],
4002 ];
4003 testLog2(vals);
4004 }}
4005 {
4006 float[2][16] vals = [
4007 [float.nan, float.nan],[-float.nan, float.nan],
4008 [float.infinity, float.infinity], [-float.infinity, float.nan],
4009 [float.min_normal, -0x1.f8p+6f], [-float.min_normal, float.nan],
4010 [float.max, 0x1p+7f], [-float.max, float.nan],
4011 [float.min_normal / 2, -0x1.fcp+6f], [-float.min_normal / 2, float.nan],
4012 [float.max / 2, 0x1.fcp+6f], [-float.max / 2, float.nan],
4013 [float.min_normal / 3, -0x1.fe57p+6f], [-float.min_normal / 3, float.nan],
4014 [float.max / 3, 0x1.f9a9p+6f], [-float.max / 3, float.nan],
4015 ];
4016 testLog2(vals);
4017 }
4018 {
4019 double[2][16] vals = [
4020 [double.nan, double.nan],[-double.nan, double.nan],
4021 [double.infinity, double.infinity], [-double.infinity, double.nan],
4022 [double.min_normal, -0x1.ffp+9], [-double.min_normal, double.nan],
4023 [double.max, 0x1p+10], [-double.max, double.nan],
4024 [double.min_normal / 2, -0x1.ff8p+9], [-double.min_normal / 2, double.nan],
4025 [double.max / 2, 0x1.ff8p+9], [-double.max / 2, double.nan],
4026 [double.min_normal / 3, -0x1.ffcae00d1cfdfp+9], [-double.min_normal / 3, double.nan],
4027 [double.max / 3, 0x1.ff351ff2e3021p+9], [-double.max / 3, double.nan],
4028 ];
4029 testLog2(vals);
4030 }
4031 alias F = floatTraits!real;
4032 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
4033 {{
4034 real[2][16] vals = [
4035 [real.nan, real.nan],[-real.nan, real.nan],
4036 [real.infinity, real.infinity], [-real.infinity, real.nan],
4037 [real.min_normal, -0x1.fffp+13L], [-real.min_normal, real.nan],
4038 [real.max, 0x1p+14L], [-real.max, real.nan],
4039 [real.min_normal / 2, -0x1.fff8p+13L], [-real.min_normal / 2, real.nan],
4040 [real.max / 2, 0x1.fff8p+13L], [-real.max / 2, real.nan],
4041 [real.min_normal / 3, -0x1.fffcae00d1cfdeb4p+13L], [-real.min_normal / 3, real.nan],
4042 [real.max / 3, 0x1.fff351ff2e30214cp+13L], [-real.max / 3, real.nan],
4043 ];
4044 testLog2(vals);
4045 }}
4046 }
4047
4048 /*****************************************
4049 * Extracts the exponent of x as a signed integral value.
4050 *
4051 * If x is subnormal, it is treated as if it were normalized.
4052 * For a positive, finite x:
4053 *
4054 * 1 $(LT)= $(I x) * FLT_RADIX$(SUPERSCRIPT -logb(x)) $(LT) FLT_RADIX
4055 *
4056 * $(TABLE_SV
4057 * $(TR $(TH x) $(TH logb(x)) $(TH divide by 0?) )
4058 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no))
4059 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) )
4060 * )
4061 */
4062 pragma(inline, true)
4063 real logb(real x) @trusted pure nothrow @nogc
4064 {
4065 version (InlineAsm_X87_MSVC)
4066 return logbAsm(x);
4067 else
4068 return logbImpl(x);
4069 }
4070
4071 /// ditto
4072 pragma(inline, true)
4073 double logb(double x) @trusted pure nothrow @nogc { return logbImpl(x); }
4074
4075 /// ditto
4076 pragma(inline, true)
4077 float logb(float x) @trusted pure nothrow @nogc { return logbImpl(x); }
4078
4079 ///
4080 @safe @nogc nothrow unittest
4081 {
4082 assert(logb(1.0) == 0);
4083 assert(logb(100.0) == 6);
4084
4085 assert(logb(0.0) == -real.infinity);
4086 assert(logb(real.infinity) == real.infinity);
4087 assert(logb(-real.infinity) == real.infinity);
4088 }
4089
4090 @safe @nogc nothrow unittest
4091 {
4092 import std.meta : AliasSeq;
4093 import std.typecons : Tuple;
4094 import std.math.traits : isNaN;
4095 static foreach (F; AliasSeq!(float, double, real))
4096 {{
4097 alias T = Tuple!(F, F);
4098 T[17] vals = // x, logb(x)
4099 [
4100 T(1.0 , 0 ),
4101 T(100.0 , 6 ),
4102 T(0.0 , -F.infinity),
4103 T(-0.0 , -F.infinity),
4104 T(1024 , 10 ),
4105 T(-2000 , 10 ),
4106 T(0x0.1p-127 , -131 ),
4107 T(0x0.01p-127 , -135 ),
4108 T(0x0.011p-127 , -135 ),
4109 T(F.nan , F.nan ),
4110 T(-F.nan , F.nan ),
4111 T(F.infinity , F.infinity ),
4112 T(-F.infinity , F.infinity ),
4113 T(F.min_normal , F.min_exp-1),
4114 T(-F.min_normal, F.min_exp-1),
4115 T(F.max , F.max_exp-1),
4116 T(-F.max , F.max_exp-1),
4117 ];
4118
4119 foreach (elem; vals)
4120 {
4121 if (isNaN(elem[1]))
4122 assert(isNaN(logb(elem[1])));
4123 else
4124 assert(logb(elem[0]) == elem[1]);
4125 }
4126 }}
4127 }
4128
4129 version (InlineAsm_X87_MSVC)
4130 private T logbAsm(T)(T x) @trusted pure nothrow @nogc
4131 {
4132 version (X86_64)
4133 {
4134 asm pure nothrow @nogc
4135 {
4136 naked ;
4137 fld real ptr [RCX] ;
4138 fxtract ;
4139 fstp ST(0) ;
4140 ret ;
4141 }
4142 }
4143 else
4144 {
4145 asm pure nothrow @nogc
4146 {
4147 fld x ;
4148 fxtract ;
4149 fstp ST(0) ;
4150 }
4151 }
4152 }
4153
4154 private T logbImpl(T)(T x) @trusted pure nothrow @nogc
4155 {
4156 import std.math.traits : isFinite;
4157
4158 // Handle special cases.
4159 if (!isFinite(x))
4160 return x * x;
4161 if (x == 0)
4162 return -1 / (x * x);
4163
4164 return ilogb(x);
4165 }
4166
4167 @safe @nogc nothrow unittest
4168 {
4169 import std.math.traits : floatTraits, RealFormat;
4170 import std.meta : AliasSeq;
4171
4172 static void testLogb(T)(T[2][] vals)
4173 {
4174 import std.math.operations : isClose;
4175 import std.math.traits : isNaN;
4176 foreach (ref pair; vals)
4177 {
4178 if (isNaN(pair[1]))
4179 assert(isNaN(logb(pair[0])));
4180 else
4181 assert(isClose(logb(pair[0]), pair[1]));
4182 }
4183 }
4184 static foreach (F; AliasSeq!(float, double, real))
4185 {{
4186 scope F[2][] vals = [
4187 [F(1), F(0x0p+0)], [F(2), F(0x1p+0)],
4188 [F(4), F(0x1p+1)], [F(8), F(0x1.8p+1)],
4189 [F(16), F(0x1p+2)], [F(32), F(0x1.4p+2)],
4190 [F(64), F(0x1.8p+2)], [F(128), F(0x1.cp+2)],
4191 [F(256), F(0x1p+3)], [F(512), F(0x1.2p+3)],
4192 [F(1024), F(0x1.4p+3)], [F(2048), F(0x1.6p+3)],
4193 [F(3), F(0x1p+0)], [F(5), F(0x1p+1)],
4194 [F(7), F(0x1p+1)], [F(15), F(0x1.8p+1)],
4195 [F(17), F(0x1p+2)], [F(31), F(0x1p+2)],
4196 [F(33), F(0x1.4p+2)], [F(63), F(0x1.4p+2)],
4197 [F(65), F(0x1.8p+2)], [F(-0), -F.infinity], [F(0), -F.infinity],
4198 [F(0.1), F(-0x1p+2)], [F(0.25), F(-0x1p+1)],
4199 [F(0.75), F(-0x1p+0)], [F(0.875), F(-0x1p+0)],
4200 [F(10), F(0x1.8p+1)], [F(100), F(0x1.8p+2)],
4201 [F(10000), F(0x1.ap+3)],
4202 ];
4203 testLogb(vals);
4204 }}
4205 {
4206 float[2][16] vals = [
4207 [float.nan, float.nan],[-float.nan, float.nan],
4208 [float.infinity, float.infinity], [-float.infinity, float.infinity],
4209 [float.min_normal, -0x1.f8p+6f], [-float.min_normal, -0x1.f8p+6f],
4210 [float.max, 0x1.fcp+6f], [-float.max, 0x1.fcp+6f],
4211 [float.min_normal / 2, -0x1.fcp+6f], [-float.min_normal / 2, -0x1.fcp+6f],
4212 [float.max / 2, 0x1.f8p+6f], [-float.max / 2, 0x1.f8p+6f],
4213 [float.min_normal / 3, -0x1p+7f], [-float.min_normal / 3, -0x1p+7f],
4214 [float.max / 3, 0x1.f8p+6f], [-float.max / 3, 0x1.f8p+6f],
4215 ];
4216 testLogb(vals);
4217 }
4218 {
4219 double[2][16] vals = [
4220 [double.nan, double.nan],[-double.nan, double.nan],
4221 [double.infinity, double.infinity], [-double.infinity, double.infinity],
4222 [double.min_normal, -0x1.ffp+9], [-double.min_normal, -0x1.ffp+9],
4223 [double.max, 0x1.ff8p+9], [-double.max, 0x1.ff8p+9],
4224 [double.min_normal / 2, -0x1.ff8p+9], [-double.min_normal / 2, -0x1.ff8p+9],
4225 [double.max / 2, 0x1.ffp+9], [-double.max / 2, 0x1.ffp+9],
4226 [double.min_normal / 3, -0x1p+10], [-double.min_normal / 3, -0x1p+10],
4227 [double.max / 3, 0x1.ffp+9], [-double.max / 3, 0x1.ffp+9],
4228 ];
4229 testLogb(vals);
4230 }
4231 alias F = floatTraits!real;
4232 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
4233 {{
4234 real[2][16] vals = [
4235 [real.nan, real.nan],[-real.nan, real.nan],
4236 [real.infinity, real.infinity], [-real.infinity, real.infinity],
4237 [real.min_normal, -0x1.fffp+13L], [-real.min_normal, -0x1.fffp+13L],
4238 [real.max, 0x1.fff8p+13L], [-real.max, 0x1.fff8p+13L],
4239 [real.min_normal / 2, -0x1.fff8p+13L], [-real.min_normal / 2, -0x1.fff8p+13L],
4240 [real.max / 2, 0x1.fffp+13L], [-real.max / 2, 0x1.fffp+13L],
4241 [real.min_normal / 3, -0x1p+14L], [-real.min_normal / 3, -0x1p+14L],
4242 [real.max / 3, 0x1.fffp+13L], [-real.max / 3, 0x1.fffp+13L],
4243 ];
4244 testLogb(vals);
4245 }}
4246 }
4247
4248 /*************************************
4249 * Efficiently calculates x * 2$(SUPERSCRIPT n).
4250 *
4251 * scalbn handles underflow and overflow in
4252 * the same fashion as the basic arithmetic operators.
4253 *
4254 * $(TABLE_SV
4255 * $(TR $(TH x) $(TH scalb(x)))
4256 * $(TR $(TD $(PLUSMNINF)) $(TD $(PLUSMNINF)) )
4257 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) )
4258 * )
4259 */
4260 pragma(inline, true)
4261 real scalbn(real x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
4262
4263 /// ditto
4264 pragma(inline, true)
4265 double scalbn(double x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
4266
4267 /// ditto
4268 pragma(inline, true)
4269 float scalbn(float x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
4270
4271 ///
4272 @safe pure nothrow @nogc unittest
4273 {
4274 assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L);
4275 assert(scalbn(-real.infinity, 5) == -real.infinity);
4276 assert(scalbn(2.0,10) == 2048.0);
4277 assert(scalbn(2048.0f,-10) == 2.0f);
4278 }
4279
4280 pragma(inline, true)
4281 private F _scalbn(F)(F x, int n)
4282 {
4283 import std.math.traits : isInfinity;
4284
4285 if (__ctfe)
4286 {
4287 // Handle special cases.
4288 if (x == F(0.0) || isInfinity(x))
4289 return x;
4290 }
4291 return core.math.ldexp(x, n);
4292 }
4293
4294 @safe pure nothrow @nogc unittest
4295 {
4296 // CTFE-able test
4297 static assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L);
4298 static assert(scalbn(-real.infinity, 5) == -real.infinity);
4299 // Test with large exponent delta n where the result is in bounds but 2.0L ^^ n is not.
4300 enum initialExponent = real.min_exp + 2, resultExponent = real.max_exp - 2;
4301 enum int n = resultExponent - initialExponent;
4302 enum real x = 0x1.2345678abcdefp0L * (2.0L ^^ initialExponent);
4303 enum staticResult = scalbn(x, n);
4304 static assert(staticResult == 0x1.2345678abcdefp0L * (2.0L ^^ resultExponent));
4305 assert(scalbn(x, n) == staticResult);
4306 }
4307
Mirror of all GCC SVN branches and tags