1 // Written in the D programming language.
2
3 /**
4 This is a submodule of $(MREF std, math).
5
6 It contains several functions for work with floating point numbers.
7
8 Copyright: Copyright The D Language Foundation 2000 - 2011.
9 License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
10 Authors: $(HTTP digitalmars.com, Walter Bright), Don Clugston,
11 Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger
12 Source: $(PHOBOSSRC std/math/operations.d)
13
14 Macros:
15 TABLE_SV = <table border="1" cellpadding="4" cellspacing="0">
16 <caption>Special Values</caption>
17 $0</table>
18 SVH = $(TR $(TH $1) $(TH $2))
19 SV = $(TR $(TD $1) $(TD $2))
20 NAN = $(RED NAN)
21 PLUSMN = ±
22 INFIN = ∞
23 LT = <
24 GT = >
25 */
26
27 module std.math.operations;
28
29 import std.traits : CommonType, isFloatingPoint, isIntegral, Unqual;
30
31 // Functions for NaN payloads
32 /*
33 * A 'payload' can be stored in the significand of a $(NAN). One bit is required
34 * to distinguish between a quiet and a signalling $(NAN). This leaves 22 bits
35 * of payload for a float; 51 bits for a double; 62 bits for an 80-bit real;
36 * and 111 bits for a 128-bit quad.
37 */
38 /**
39 * Create a quiet $(NAN), storing an integer inside the payload.
40 *
41 * For floats, the largest possible payload is 0x3F_FFFF.
42 * For doubles, it is 0x3_FFFF_FFFF_FFFF.
43 * For 80-bit or 128-bit reals, it is 0x3FFF_FFFF_FFFF_FFFF.
44 */
NaN(ulong payload)45 real NaN(ulong payload) @trusted pure nothrow @nogc
46 {
47 import std.math : floatTraits, RealFormat;
48
49 alias F = floatTraits!(real);
50 static if (F.realFormat == RealFormat.ieeeExtended ||
51 F.realFormat == RealFormat.ieeeExtended53)
52 {
53 // real80 (in x86 real format, the implied bit is actually
54 // not implied but a real bit which is stored in the real)
55 ulong v = 3; // implied bit = 1, quiet bit = 1
56 }
57 else
58 {
59 ulong v = 1; // no implied bit. quiet bit = 1
60 }
61 if (__ctfe)
62 {
63 v = 1; // We use a double in CTFE.
64 assert(payload >>> 51 == 0,
65 "Cannot set more than 51 bits of NaN payload in CTFE.");
66 }
67
68
69 ulong a = payload;
70
71 // 22 Float bits
72 ulong w = a & 0x3F_FFFF;
73 a -= w;
74
75 v <<=22;
76 v |= w;
77 a >>=22;
78
79 // 29 Double bits
80 v <<=29;
81 w = a & 0xFFF_FFFF;
82 v |= w;
83 a -= w;
84 a >>=29;
85
86 if (__ctfe)
87 {
88 v |= 0x7FF0_0000_0000_0000;
89 return *cast(double*) &v;
90 }
91 else static if (F.realFormat == RealFormat.ieeeDouble)
92 {
93 v |= 0x7FF0_0000_0000_0000;
94 real x;
95 * cast(ulong *)(&x) = v;
96 return x;
97 }
98 else
99 {
100 v <<=11;
101 a &= 0x7FF;
102 v |= a;
103 real x = real.nan;
104
105 // Extended real bits
106 static if (F.realFormat == RealFormat.ieeeQuadruple)
107 {
108 v <<= 1; // there's no implicit bit
109
110 version (LittleEndian)
111 {
112 *cast(ulong*)(6+cast(ubyte*)(&x)) = v;
113 }
114 else
115 {
116 *cast(ulong*)(2+cast(ubyte*)(&x)) = v;
117 }
118 }
119 else
120 {
121 *cast(ulong *)(&x) = v;
122 }
123 return x;
124 }
125 }
126
127 ///
128 @safe @nogc pure nothrow unittest
129 {
130 import std.math.traits : isNaN;
131
132 real a = NaN(1_000_000);
133 assert(isNaN(a));
134 assert(getNaNPayload(a) == 1_000_000);
135 }
136
137 @system pure nothrow @nogc unittest // not @safe because taking address of local.
138 {
139 import std.math : floatTraits, RealFormat;
140
141 static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
142 {
143 auto x = NaN(1);
144 auto xl = *cast(ulong*)&x;
145 assert(xl & 0x8_0000_0000_0000UL); //non-signaling bit, bit 52
146 assert((xl & 0x7FF0_0000_0000_0000UL) == 0x7FF0_0000_0000_0000UL); //all exp bits set
147 }
148 }
149
150 /**
151 * Extract an integral payload from a $(NAN).
152 *
153 * Returns:
154 * the integer payload as a ulong.
155 *
156 * For floats, the largest possible payload is 0x3F_FFFF.
157 * For doubles, it is 0x3_FFFF_FFFF_FFFF.
158 * For 80-bit or 128-bit reals, it is 0x3FFF_FFFF_FFFF_FFFF.
159 */
getNaNPayload(real x)160 ulong getNaNPayload(real x) @trusted pure nothrow @nogc
161 {
162 import std.math : floatTraits, RealFormat;
163
164 // assert(isNaN(x));
165 alias F = floatTraits!(real);
166 ulong m = void;
167 if (__ctfe)
168 {
169 double y = x;
170 m = *cast(ulong*) &y;
171 // Make it look like an 80-bit significand.
172 // Skip exponent, and quiet bit
173 m &= 0x0007_FFFF_FFFF_FFFF;
174 m <<= 11;
175 }
176 else static if (F.realFormat == RealFormat.ieeeDouble)
177 {
178 m = *cast(ulong*)(&x);
179 // Make it look like an 80-bit significand.
180 // Skip exponent, and quiet bit
181 m &= 0x0007_FFFF_FFFF_FFFF;
182 m <<= 11;
183 }
184 else static if (F.realFormat == RealFormat.ieeeQuadruple)
185 {
186 version (LittleEndian)
187 {
188 m = *cast(ulong*)(6+cast(ubyte*)(&x));
189 }
190 else
191 {
192 m = *cast(ulong*)(2+cast(ubyte*)(&x));
193 }
194
195 m >>= 1; // there's no implicit bit
196 }
197 else
198 {
199 m = *cast(ulong*)(&x);
200 }
201
202 // ignore implicit bit and quiet bit
203
204 const ulong f = m & 0x3FFF_FF00_0000_0000L;
205
206 ulong w = f >>> 40;
207 w |= (m & 0x00FF_FFFF_F800L) << (22 - 11);
208 w |= (m & 0x7FF) << 51;
209 return w;
210 }
211
212 ///
213 @safe @nogc pure nothrow unittest
214 {
215 import std.math.traits : isNaN;
216
217 real a = NaN(1_000_000);
218 assert(isNaN(a));
219 assert(getNaNPayload(a) == 1_000_000);
220 }
221
222 @safe @nogc pure nothrow unittest
223 {
224 import std.math.traits : isIdentical, isNaN;
225
226 enum real a = NaN(1_000_000);
227 static assert(isNaN(a));
228 static assert(getNaNPayload(a) == 1_000_000);
229 real b = NaN(1_000_000);
230 assert(isIdentical(b, a));
231 // The CTFE version of getNaNPayload relies on it being impossible
232 // for a CTFE-constructed NaN to have more than 51 bits of payload.
233 enum nanNaN = NaN(getNaNPayload(real.nan));
234 assert(isIdentical(real.nan, nanNaN));
235 static if (real.init != real.init)
236 {
237 enum initNaN = NaN(getNaNPayload(real.init));
238 assert(isIdentical(real.init, initNaN));
239 }
240 }
241
debug(UnitTest)242 debug(UnitTest)
243 {
244 @safe pure nothrow @nogc unittest
245 {
246 real nan4 = NaN(0x789_ABCD_EF12_3456);
247 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended
248 || floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
249 {
250 assert(getNaNPayload(nan4) == 0x789_ABCD_EF12_3456);
251 }
252 else
253 {
254 assert(getNaNPayload(nan4) == 0x1_ABCD_EF12_3456);
255 }
256 double nan5 = nan4;
257 assert(getNaNPayload(nan5) == 0x1_ABCD_EF12_3456);
258 float nan6 = nan4;
259 assert(getNaNPayload(nan6) == 0x12_3456);
260 nan4 = NaN(0xFABCD);
261 assert(getNaNPayload(nan4) == 0xFABCD);
262 nan6 = nan4;
263 assert(getNaNPayload(nan6) == 0xFABCD);
264 nan5 = NaN(0x100_0000_0000_3456);
265 assert(getNaNPayload(nan5) == 0x0000_0000_3456);
266 }
267 }
268
269 /**
270 * Calculate the next largest floating point value after x.
271 *
272 * Return the least number greater than x that is representable as a real;
273 * thus, it gives the next point on the IEEE number line.
274 *
275 * $(TABLE_SV
276 * $(SVH x, nextUp(x) )
277 * $(SV -$(INFIN), -real.max )
278 * $(SV $(PLUSMN)0.0, real.min_normal*real.epsilon )
279 * $(SV real.max, $(INFIN) )
280 * $(SV $(INFIN), $(INFIN) )
281 * $(SV $(NAN), $(NAN) )
282 * )
283 */
nextUp(real x)284 real nextUp(real x) @trusted pure nothrow @nogc
285 {
286 import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB;
287
288 alias F = floatTraits!(real);
289 static if (F.realFormat != RealFormat.ieeeDouble)
290 {
291 if (__ctfe)
292 {
293 if (x == -real.infinity)
294 return -real.max;
295 if (!(x < real.infinity)) // Infinity or NaN.
296 return x;
297 real delta;
298 // Start with a decent estimate of delta.
299 if (x <= 0x1.ffffffffffffep+1023 && x >= -double.max)
300 {
301 const double d = cast(double) x;
302 delta = (cast(real) nextUp(d) - cast(real) d) * 0x1p-11L;
303 while (x + (delta * 0x1p-100L) > x)
304 delta *= 0x1p-100L;
305 }
306 else
307 {
308 delta = 0x1p960L;
309 while (!(x + delta > x) && delta < real.max * 0x1p-100L)
310 delta *= 0x1p100L;
311 }
312 if (x + delta > x)
313 {
314 while (x + (delta / 2) > x)
315 delta /= 2;
316 }
317 else
318 {
319 do { delta += delta; } while (!(x + delta > x));
320 }
321 if (x < 0 && x + delta == 0)
322 return -0.0L;
323 return x + delta;
324 }
325 }
326 static if (F.realFormat == RealFormat.ieeeDouble)
327 {
328 return nextUp(cast(double) x);
329 }
330 else static if (F.realFormat == RealFormat.ieeeQuadruple)
331 {
332 ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
333 if (e == F.EXPMASK)
334 {
335 // NaN or Infinity
336 if (x == -real.infinity) return -real.max;
337 return x; // +Inf and NaN are unchanged.
338 }
339
340 auto ps = cast(ulong *)&x;
341 if (ps[MANTISSA_MSB] & 0x8000_0000_0000_0000)
342 {
343 // Negative number
344 if (ps[MANTISSA_LSB] == 0 && ps[MANTISSA_MSB] == 0x8000_0000_0000_0000)
345 {
346 // it was negative zero, change to smallest subnormal
347 ps[MANTISSA_LSB] = 1;
348 ps[MANTISSA_MSB] = 0;
349 return x;
350 }
351 if (ps[MANTISSA_LSB] == 0) --ps[MANTISSA_MSB];
352 --ps[MANTISSA_LSB];
353 }
354 else
355 {
356 // Positive number
357 ++ps[MANTISSA_LSB];
358 if (ps[MANTISSA_LSB] == 0) ++ps[MANTISSA_MSB];
359 }
360 return x;
361 }
362 else static if (F.realFormat == RealFormat.ieeeExtended ||
363 F.realFormat == RealFormat.ieeeExtended53)
364 {
365 // For 80-bit reals, the "implied bit" is a nuisance...
366 ushort *pe = cast(ushort *)&x;
367 ulong *ps = cast(ulong *)&x;
368 // EPSILON is 1 for 64-bit, and 2048 for 53-bit precision reals.
369 enum ulong EPSILON = 2UL ^^ (64 - real.mant_dig);
370
371 if ((pe[F.EXPPOS_SHORT] & F.EXPMASK) == F.EXPMASK)
372 {
373 // First, deal with NANs and infinity
374 if (x == -real.infinity) return -real.max;
375 return x; // +Inf and NaN are unchanged.
376 }
377 if (pe[F.EXPPOS_SHORT] & 0x8000)
378 {
379 // Negative number -- need to decrease the significand
380 *ps -= EPSILON;
381 // Need to mask with 0x7FFF... so subnormals are treated correctly.
382 if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_FFFF_FFFF_FFFF)
383 {
384 if (pe[F.EXPPOS_SHORT] == 0x8000) // it was negative zero
385 {
386 *ps = 1;
387 pe[F.EXPPOS_SHORT] = 0; // smallest subnormal.
388 return x;
389 }
390
391 --pe[F.EXPPOS_SHORT];
392
393 if (pe[F.EXPPOS_SHORT] == 0x8000)
394 return x; // it's become a subnormal, implied bit stays low.
395
396 *ps = 0xFFFF_FFFF_FFFF_FFFF; // set the implied bit
397 return x;
398 }
399 return x;
400 }
401 else
402 {
403 // Positive number -- need to increase the significand.
404 // Works automatically for positive zero.
405 *ps += EPSILON;
406 if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0)
407 {
408 // change in exponent
409 ++pe[F.EXPPOS_SHORT];
410 *ps = 0x8000_0000_0000_0000; // set the high bit
411 }
412 }
413 return x;
414 }
415 else // static if (F.realFormat == RealFormat.ibmExtended)
416 {
417 assert(0, "nextUp not implemented");
418 }
419 }
420
421 /** ditto */
nextUp(double x)422 double nextUp(double x) @trusted pure nothrow @nogc
423 {
424 ulong s = *cast(ulong *)&x;
425
426 if ((s & 0x7FF0_0000_0000_0000) == 0x7FF0_0000_0000_0000)
427 {
428 // First, deal with NANs and infinity
429 if (x == -x.infinity) return -x.max;
430 return x; // +INF and NAN are unchanged.
431 }
432 if (s & 0x8000_0000_0000_0000) // Negative number
433 {
434 if (s == 0x8000_0000_0000_0000) // it was negative zero
435 {
436 s = 0x0000_0000_0000_0001; // change to smallest subnormal
437 return *cast(double*) &s;
438 }
439 --s;
440 }
441 else
442 { // Positive number
443 ++s;
444 }
445 return *cast(double*) &s;
446 }
447
448 /** ditto */
nextUp(float x)449 float nextUp(float x) @trusted pure nothrow @nogc
450 {
451 uint s = *cast(uint *)&x;
452
453 if ((s & 0x7F80_0000) == 0x7F80_0000)
454 {
455 // First, deal with NANs and infinity
456 if (x == -x.infinity) return -x.max;
457
458 return x; // +INF and NAN are unchanged.
459 }
460 if (s & 0x8000_0000) // Negative number
461 {
462 if (s == 0x8000_0000) // it was negative zero
463 {
464 s = 0x0000_0001; // change to smallest subnormal
465 return *cast(float*) &s;
466 }
467
468 --s;
469 }
470 else
471 {
472 // Positive number
473 ++s;
474 }
475 return *cast(float*) &s;
476 }
477
478 ///
479 @safe @nogc pure nothrow unittest
480 {
481 assert(nextUp(1.0 - 1.0e-6).feqrel(0.999999) > 16);
482 assert(nextUp(1.0 - real.epsilon).feqrel(1.0) > 16);
483 }
484
485 /**
486 * Calculate the next smallest floating point value before x.
487 *
488 * Return the greatest number less than x that is representable as a real;
489 * thus, it gives the previous point on the IEEE number line.
490 *
491 * $(TABLE_SV
492 * $(SVH x, nextDown(x) )
493 * $(SV $(INFIN), real.max )
494 * $(SV $(PLUSMN)0.0, -real.min_normal*real.epsilon )
495 * $(SV -real.max, -$(INFIN) )
496 * $(SV -$(INFIN), -$(INFIN) )
497 * $(SV $(NAN), $(NAN) )
498 * )
499 */
nextDown(real x)500 real nextDown(real x) @safe pure nothrow @nogc
501 {
502 return -nextUp(-x);
503 }
504
505 /** ditto */
nextDown(double x)506 double nextDown(double x) @safe pure nothrow @nogc
507 {
508 return -nextUp(-x);
509 }
510
511 /** ditto */
nextDown(float x)512 float nextDown(float x) @safe pure nothrow @nogc
513 {
514 return -nextUp(-x);
515 }
516
517 ///
518 @safe pure nothrow @nogc unittest
519 {
520 assert( nextDown(1.0 + real.epsilon) == 1.0);
521 }
522
523 @safe pure nothrow @nogc unittest
524 {
525 import std.math : floatTraits, RealFormat;
526 import std.math.traits : isIdentical;
527
528 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended ||
529 floatTraits!(real).realFormat == RealFormat.ieeeDouble ||
530 floatTraits!(real).realFormat == RealFormat.ieeeExtended53 ||
531 floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
532 {
533 // Tests for reals
534 assert(isIdentical(nextUp(NaN(0xABC)), NaN(0xABC)));
535 //static assert(isIdentical(nextUp(NaN(0xABC)), NaN(0xABC)));
536 // negative numbers
537 assert( nextUp(-real.infinity) == -real.max );
538 assert( nextUp(-1.0L-real.epsilon) == -1.0 );
539 assert( nextUp(-2.0L) == -2.0 + real.epsilon);
540 static assert( nextUp(-real.infinity) == -real.max );
541 static assert( nextUp(-1.0L-real.epsilon) == -1.0 );
542 static assert( nextUp(-2.0L) == -2.0 + real.epsilon);
543 // subnormals and zero
544 assert( nextUp(-real.min_normal) == -real.min_normal*(1-real.epsilon) );
545 assert( nextUp(-real.min_normal*(1-real.epsilon)) == -real.min_normal*(1-2*real.epsilon) );
546 assert( isIdentical(-0.0L, nextUp(-real.min_normal*real.epsilon)) );
547 assert( nextUp(-0.0L) == real.min_normal*real.epsilon );
548 assert( nextUp(0.0L) == real.min_normal*real.epsilon );
549 assert( nextUp(real.min_normal*(1-real.epsilon)) == real.min_normal );
550 assert( nextUp(real.min_normal) == real.min_normal*(1+real.epsilon) );
551 static assert( nextUp(-real.min_normal) == -real.min_normal*(1-real.epsilon) );
552 static assert( nextUp(-real.min_normal*(1-real.epsilon)) == -real.min_normal*(1-2*real.epsilon) );
553 static assert( -0.0L is nextUp(-real.min_normal*real.epsilon) );
554 static assert( nextUp(-0.0L) == real.min_normal*real.epsilon );
555 static assert( nextUp(0.0L) == real.min_normal*real.epsilon );
556 static assert( nextUp(real.min_normal*(1-real.epsilon)) == real.min_normal );
557 static assert( nextUp(real.min_normal) == real.min_normal*(1+real.epsilon) );
558 // positive numbers
559 assert( nextUp(1.0L) == 1.0 + real.epsilon );
560 assert( nextUp(2.0L-real.epsilon) == 2.0 );
561 assert( nextUp(real.max) == real.infinity );
562 assert( nextUp(real.infinity)==real.infinity );
563 static assert( nextUp(1.0L) == 1.0 + real.epsilon );
564 static assert( nextUp(2.0L-real.epsilon) == 2.0 );
565 static assert( nextUp(real.max) == real.infinity );
566 static assert( nextUp(real.infinity)==real.infinity );
567 // ctfe near double.max boundary
568 static assert(nextUp(nextDown(cast(real) double.max)) == cast(real) double.max);
569 }
570
571 double n = NaN(0xABC);
572 assert(isIdentical(nextUp(n), n));
573 // negative numbers
574 assert( nextUp(-double.infinity) == -double.max );
575 assert( nextUp(-1-double.epsilon) == -1.0 );
576 assert( nextUp(-2.0) == -2.0 + double.epsilon);
577 // subnormals and zero
578
579 assert( nextUp(-double.min_normal) == -double.min_normal*(1-double.epsilon) );
580 assert( nextUp(-double.min_normal*(1-double.epsilon)) == -double.min_normal*(1-2*double.epsilon) );
581 assert( isIdentical(-0.0, nextUp(-double.min_normal*double.epsilon)) );
582 assert( nextUp(0.0) == double.min_normal*double.epsilon );
583 assert( nextUp(-0.0) == double.min_normal*double.epsilon );
584 assert( nextUp(double.min_normal*(1-double.epsilon)) == double.min_normal );
585 assert( nextUp(double.min_normal) == double.min_normal*(1+double.epsilon) );
586 // positive numbers
587 assert( nextUp(1.0) == 1.0 + double.epsilon );
588 assert( nextUp(2.0-double.epsilon) == 2.0 );
589 assert( nextUp(double.max) == double.infinity );
590
591 float fn = NaN(0xABC);
592 assert(isIdentical(nextUp(fn), fn));
593 float f = -float.min_normal*(1-float.epsilon);
594 float f1 = -float.min_normal;
595 assert( nextUp(f1) == f);
596 f = 1.0f+float.epsilon;
597 f1 = 1.0f;
598 assert( nextUp(f1) == f );
599 f1 = -0.0f;
600 assert( nextUp(f1) == float.min_normal*float.epsilon);
601 assert( nextUp(float.infinity)==float.infinity );
602
603 assert(nextDown(1.0L+real.epsilon)==1.0);
604 assert(nextDown(1.0+double.epsilon)==1.0);
605 f = 1.0f+float.epsilon;
606 assert(nextDown(f)==1.0);
607 assert(nextafter(1.0+real.epsilon, -real.infinity)==1.0);
608
609 // CTFE
610
611 enum double ctfe_n = NaN(0xABC);
612 //static assert(isIdentical(nextUp(ctfe_n), ctfe_n)); // FIXME: https://issues.dlang.org/show_bug.cgi?id=20197
613 static assert(nextUp(double.nan) is double.nan);
614 // negative numbers
615 static assert( nextUp(-double.infinity) == -double.max );
616 static assert( nextUp(-1-double.epsilon) == -1.0 );
617 static assert( nextUp(-2.0) == -2.0 + double.epsilon);
618 // subnormals and zero
619
620 static assert( nextUp(-double.min_normal) == -double.min_normal*(1-double.epsilon) );
621 static assert( nextUp(-double.min_normal*(1-double.epsilon)) == -double.min_normal*(1-2*double.epsilon) );
622 static assert( -0.0 is nextUp(-double.min_normal*double.epsilon) );
623 static assert( nextUp(0.0) == double.min_normal*double.epsilon );
624 static assert( nextUp(-0.0) == double.min_normal*double.epsilon );
625 static assert( nextUp(double.min_normal*(1-double.epsilon)) == double.min_normal );
626 static assert( nextUp(double.min_normal) == double.min_normal*(1+double.epsilon) );
627 // positive numbers
628 static assert( nextUp(1.0) == 1.0 + double.epsilon );
629 static assert( nextUp(2.0-double.epsilon) == 2.0 );
630 static assert( nextUp(double.max) == double.infinity );
631
632 enum float ctfe_fn = NaN(0xABC);
633 //static assert(isIdentical(nextUp(ctfe_fn), ctfe_fn)); // FIXME: https://issues.dlang.org/show_bug.cgi?id=20197
634 static assert(nextUp(float.nan) is float.nan);
635 static assert(nextUp(-float.min_normal) == -float.min_normal*(1-float.epsilon));
636 static assert(nextUp(1.0f) == 1.0f+float.epsilon);
637 static assert(nextUp(-0.0f) == float.min_normal*float.epsilon);
638 static assert(nextUp(float.infinity)==float.infinity);
639 static assert(nextDown(1.0L+real.epsilon)==1.0);
640 static assert(nextDown(1.0+double.epsilon)==1.0);
641 static assert(nextDown(1.0f+float.epsilon)==1.0);
642 static assert(nextafter(1.0+real.epsilon, -real.infinity)==1.0);
643 }
644
645
646
647 /******************************************
648 * Calculates the next representable value after x in the direction of y.
649 *
650 * If y > x, the result will be the next largest floating-point value;
651 * if y < x, the result will be the next smallest value.
652 * If x == y, the result is y.
653 * If x or y is a NaN, the result is a NaN.
654 *
655 * Remarks:
656 * This function is not generally very useful; it's almost always better to use
657 * the faster functions nextUp() or nextDown() instead.
658 *
659 * The FE_INEXACT and FE_OVERFLOW exceptions will be raised if x is finite and
660 * the function result is infinite. The FE_INEXACT and FE_UNDERFLOW
661 * exceptions will be raised if the function value is subnormal, and x is
662 * not equal to y.
663 */
nextafter(T)664 T nextafter(T)(const T x, const T y) @safe pure nothrow @nogc
665 {
666 import std.math.traits : isNaN;
667
668 if (x == y || isNaN(y))
669 {
670 return y;
671 }
672
673 if (isNaN(x))
674 {
675 return x;
676 }
677
678 return ((y>x) ? nextUp(x) : nextDown(x));
679 }
680
681 ///
682 @safe pure nothrow @nogc unittest
683 {
684 import std.math.traits : isNaN;
685
686 float a = 1;
687 assert(is(typeof(nextafter(a, a)) == float));
688 assert(nextafter(a, a.infinity) > a);
689 assert(isNaN(nextafter(a, a.nan)));
690 assert(isNaN(nextafter(a.nan, a)));
691
692 double b = 2;
693 assert(is(typeof(nextafter(b, b)) == double));
694 assert(nextafter(b, b.infinity) > b);
695 assert(isNaN(nextafter(b, b.nan)));
696 assert(isNaN(nextafter(b.nan, b)));
697
698 real c = 3;
699 assert(is(typeof(nextafter(c, c)) == real));
700 assert(nextafter(c, c.infinity) > c);
701 assert(isNaN(nextafter(c, c.nan)));
702 assert(isNaN(nextafter(c.nan, c)));
703 }
704
705 @safe pure nothrow @nogc unittest
706 {
707 import std.math.traits : isNaN, signbit;
708
709 // CTFE
710 enum float a = 1;
711 static assert(is(typeof(nextafter(a, a)) == float));
712 static assert(nextafter(a, a.infinity) > a);
713 static assert(isNaN(nextafter(a, a.nan)));
714 static assert(isNaN(nextafter(a.nan, a)));
715
716 enum double b = 2;
717 static assert(is(typeof(nextafter(b, b)) == double));
718 static assert(nextafter(b, b.infinity) > b);
719 static assert(isNaN(nextafter(b, b.nan)));
720 static assert(isNaN(nextafter(b.nan, b)));
721
722 enum real c = 3;
723 static assert(is(typeof(nextafter(c, c)) == real));
724 static assert(nextafter(c, c.infinity) > c);
725 static assert(isNaN(nextafter(c, c.nan)));
726 static assert(isNaN(nextafter(c.nan, c)));
727
728 enum real negZero = nextafter(+0.0L, -0.0L);
729 static assert(negZero == -0.0L);
730 static assert(signbit(negZero));
731
732 static assert(nextafter(c, c) == c);
733 }
734
735 //real nexttoward(real x, real y) { return core.stdc.math.nexttowardl(x, y); }
736
737 /**
738 * Returns the positive difference between x and y.
739 *
740 * Equivalent to `fmax(x-y, 0)`.
741 *
742 * Returns:
743 * $(TABLE_SV
744 * $(TR $(TH x, y) $(TH fdim(x, y)))
745 * $(TR $(TD x $(GT) y) $(TD x - y))
746 * $(TR $(TD x $(LT)= y) $(TD +0.0))
747 * )
748 */
fdim(real x,real y)749 real fdim(real x, real y) @safe pure nothrow @nogc
750 {
751 return (x < y) ? +0.0 : x - y;
752 }
753
754 ///
755 @safe pure nothrow @nogc unittest
756 {
757 import std.math.traits : isNaN;
758
759 assert(fdim(2.0, 0.0) == 2.0);
760 assert(fdim(-2.0, 0.0) == 0.0);
761 assert(fdim(real.infinity, 2.0) == real.infinity);
762 assert(isNaN(fdim(real.nan, 2.0)));
763 assert(isNaN(fdim(2.0, real.nan)));
764 assert(isNaN(fdim(real.nan, real.nan)));
765 }
766
767 /**
768 * Returns the larger of `x` and `y`.
769 *
770 * If one of the arguments is a `NaN`, the other is returned.
771 *
772 * See_Also: $(REF max, std,algorithm,comparison) is faster because it does not perform the `isNaN` test.
773 */
774 F fmax(F)(const F x, const F y) @safe pure nothrow @nogc
775 if (__traits(isFloating, F))
776 {
777 import std.math.traits : isNaN;
778
779 // Do the more predictable test first. Generates 0 branches with ldc and 1 branch with gdc.
780 // See https://godbolt.org/z/erxrW9
781 if (isNaN(x)) return y;
782 return y > x ? y : x;
783 }
784
785 ///
786 @safe pure nothrow @nogc unittest
787 {
788 import std.meta : AliasSeq;
789 static foreach (F; AliasSeq!(float, double, real))
790 {
791 assert(fmax(F(0.0), F(2.0)) == 2.0);
792 assert(fmax(F(-2.0), 0.0) == F(0.0));
793 assert(fmax(F.infinity, F(2.0)) == F.infinity);
794 assert(fmax(F.nan, F(2.0)) == F(2.0));
795 assert(fmax(F(2.0), F.nan) == F(2.0));
796 }
797 }
798
799 /**
800 * Returns the smaller of `x` and `y`.
801 *
802 * If one of the arguments is a `NaN`, the other is returned.
803 *
804 * See_Also: $(REF min, std,algorithm,comparison) is faster because it does not perform the `isNaN` test.
805 */
806 F fmin(F)(const F x, const F y) @safe pure nothrow @nogc
807 if (__traits(isFloating, F))
808 {
809 import std.math.traits : isNaN;
810
811 // Do the more predictable test first. Generates 0 branches with ldc and 1 branch with gdc.
812 // See https://godbolt.org/z/erxrW9
813 if (isNaN(x)) return y;
814 return y < x ? y : x;
815 }
816
817 ///
818 @safe pure nothrow @nogc unittest
819 {
820 import std.meta : AliasSeq;
821 static foreach (F; AliasSeq!(float, double, real))
822 {
823 assert(fmin(F(0.0), F(2.0)) == 0.0);
824 assert(fmin(F(-2.0), F(0.0)) == -2.0);
825 assert(fmin(F.infinity, F(2.0)) == 2.0);
826 assert(fmin(F.nan, F(2.0)) == 2.0);
827 assert(fmin(F(2.0), F.nan) == 2.0);
828 }
829 }
830
831 /**************************************
832 * Returns (x * y) + z, rounding only once according to the
833 * current rounding mode.
834 *
835 * BUGS: Not currently implemented - rounds twice.
836 */
pragma(inline,true)837 pragma(inline, true)
838 real fma(real x, real y, real z) @safe pure nothrow @nogc { return (x * y) + z; }
839
840 ///
841 @safe pure nothrow @nogc unittest
842 {
843 assert(fma(0.0, 2.0, 2.0) == 2.0);
844 assert(fma(2.0, 2.0, 2.0) == 6.0);
845 assert(fma(real.infinity, 2.0, 2.0) == real.infinity);
846 assert(fma(real.nan, 2.0, 2.0) is real.nan);
847 assert(fma(2.0, 2.0, real.nan) is real.nan);
848 }
849
850 /**************************************
851 * To what precision is x equal to y?
852 *
853 * Returns: the number of mantissa bits which are equal in x and y.
854 * eg, 0x1.F8p+60 and 0x1.F1p+60 are equal to 5 bits of precision.
855 *
856 * $(TABLE_SV
857 * $(TR $(TH x) $(TH y) $(TH feqrel(x, y)))
858 * $(TR $(TD x) $(TD x) $(TD real.mant_dig))
859 * $(TR $(TD x) $(TD $(GT)= 2*x) $(TD 0))
860 * $(TR $(TD x) $(TD $(LT)= x/2) $(TD 0))
861 * $(TR $(TD $(NAN)) $(TD any) $(TD 0))
862 * $(TR $(TD any) $(TD $(NAN)) $(TD 0))
863 * )
864 */
865 int feqrel(X)(const X x, const X y) @trusted pure nothrow @nogc
866 if (isFloatingPoint!(X))
867 {
868 import std.math : floatTraits, RealFormat;
869 import core.math : fabs;
870
871 /* Public Domain. Author: Don Clugston, 18 Aug 2005.
872 */
873 alias F = floatTraits!(X);
874 static if (F.realFormat == RealFormat.ieeeSingle
875 || F.realFormat == RealFormat.ieeeDouble
876 || F.realFormat == RealFormat.ieeeExtended
877 || F.realFormat == RealFormat.ieeeExtended53
878 || F.realFormat == RealFormat.ieeeQuadruple)
879 {
880 if (x == y)
881 return X.mant_dig; // ensure diff != 0, cope with INF.
882
883 Unqual!X diff = fabs(x - y);
884
885 ushort *pa = cast(ushort *)(&x);
886 ushort *pb = cast(ushort *)(&y);
887 ushort *pd = cast(ushort *)(&diff);
888
889
890 // The difference in abs(exponent) between x or y and abs(x-y)
891 // is equal to the number of significand bits of x which are
892 // equal to y. If negative, x and y have different exponents.
893 // If positive, x and y are equal to 'bitsdiff' bits.
894 // AND with 0x7FFF to form the absolute value.
895 // To avoid out-by-1 errors, we subtract 1 so it rounds down
896 // if the exponents were different. This means 'bitsdiff' is
897 // always 1 lower than we want, except that if bitsdiff == 0,
898 // they could have 0 or 1 bits in common.
899
900 int bitsdiff = ((( (pa[F.EXPPOS_SHORT] & F.EXPMASK)
901 + (pb[F.EXPPOS_SHORT] & F.EXPMASK)
902 - (1 << F.EXPSHIFT)) >> 1)
903 - (pd[F.EXPPOS_SHORT] & F.EXPMASK)) >> F.EXPSHIFT;
904 if ( (pd[F.EXPPOS_SHORT] & F.EXPMASK) == 0)
905 { // Difference is subnormal
906 // For subnormals, we need to add the number of zeros that
907 // lie at the start of diff's significand.
908 // We do this by multiplying by 2^^real.mant_dig
909 diff *= F.RECIP_EPSILON;
910 return bitsdiff + X.mant_dig - ((pd[F.EXPPOS_SHORT] & F.EXPMASK) >> F.EXPSHIFT);
911 }
912
913 if (bitsdiff > 0)
914 return bitsdiff + 1; // add the 1 we subtracted before
915
916 // Avoid out-by-1 errors when factor is almost 2.
917 if (bitsdiff == 0
918 && ((pa[F.EXPPOS_SHORT] ^ pb[F.EXPPOS_SHORT]) & F.EXPMASK) == 0)
919 {
920 return 1;
921 } else return 0;
922 }
923 else
924 {
925 static assert(false, "Not implemented for this architecture");
926 }
927 }
928
929 ///
930 @safe pure unittest
931 {
932 assert(feqrel(2.0, 2.0) == 53);
933 assert(feqrel(2.0f, 2.0f) == 24);
934 assert(feqrel(2.0, double.nan) == 0);
935
936 // Test that numbers are within n digits of each
937 // other by testing if feqrel > n * log2(10)
938
939 // five digits
940 assert(feqrel(2.0, 2.00001) > 16);
941 // ten digits
942 assert(feqrel(2.0, 2.00000000001) > 33);
943 }
944
945 @safe pure nothrow @nogc unittest
946 {
testFeqrel(F)947 void testFeqrel(F)()
948 {
949 // Exact equality
950 assert(feqrel(F.max, F.max) == F.mant_dig);
951 assert(feqrel!(F)(0.0, 0.0) == F.mant_dig);
952 assert(feqrel(F.infinity, F.infinity) == F.mant_dig);
953
954 // a few bits away from exact equality
955 F w=1;
956 for (int i = 1; i < F.mant_dig - 1; ++i)
957 {
958 assert(feqrel!(F)(1.0 + w * F.epsilon, 1.0) == F.mant_dig-i);
959 assert(feqrel!(F)(1.0 - w * F.epsilon, 1.0) == F.mant_dig-i);
960 assert(feqrel!(F)(1.0, 1 + (w-1) * F.epsilon) == F.mant_dig - i + 1);
961 w*=2;
962 }
963
964 assert(feqrel!(F)(1.5+F.epsilon, 1.5) == F.mant_dig-1);
965 assert(feqrel!(F)(1.5-F.epsilon, 1.5) == F.mant_dig-1);
966 assert(feqrel!(F)(1.5-F.epsilon, 1.5+F.epsilon) == F.mant_dig-2);
967
968
969 // Numbers that are close
970 assert(feqrel!(F)(0x1.Bp+84, 0x1.B8p+84) == 5);
971 assert(feqrel!(F)(0x1.8p+10, 0x1.Cp+10) == 2);
972 assert(feqrel!(F)(1.5 * (1 - F.epsilon), 1.0L) == 2);
973 assert(feqrel!(F)(1.5, 1.0) == 1);
974 assert(feqrel!(F)(2 * (1 - F.epsilon), 1.0L) == 1);
975
976 // Factors of 2
977 assert(feqrel(F.max, F.infinity) == 0);
978 assert(feqrel!(F)(2 * (1 - F.epsilon), 1.0L) == 1);
979 assert(feqrel!(F)(1.0, 2.0) == 0);
980 assert(feqrel!(F)(4.0, 1.0) == 0);
981
982 // Extreme inequality
983 assert(feqrel(F.nan, F.nan) == 0);
984 assert(feqrel!(F)(0.0L, -F.nan) == 0);
985 assert(feqrel(F.nan, F.infinity) == 0);
986 assert(feqrel(F.infinity, -F.infinity) == 0);
987 assert(feqrel(F.max, -F.max) == 0);
988
989 assert(feqrel(F.min_normal / 8, F.min_normal / 17) == 3);
990
991 const F Const = 2;
992 immutable F Immutable = 2;
993 auto Compiles = feqrel(Const, Immutable);
994 }
995
996 assert(feqrel(7.1824L, 7.1824L) == real.mant_dig);
997
998 testFeqrel!(real)();
999 testFeqrel!(double)();
1000 testFeqrel!(float)();
1001 }
1002
1003 /**
1004 Computes whether a values is approximately equal to a reference value,
1005 admitting a maximum relative difference, and a maximum absolute difference.
1006
1007 Warning:
1008 This template is considered out-dated. It will be removed from
1009 Phobos in 2.106.0. Please use $(LREF isClose) instead. To achieve
1010 a similar behaviour to `approxEqual(a, b)` use
1011 `isClose(a, b, 1e-2, 1e-5)`. In case of comparing to 0.0,
1012 `isClose(a, b, 0.0, eps)` should be used, where `eps`
1013 represents the accepted deviation from 0.0."
1014
1015 Params:
1016 value = Value to compare.
1017 reference = Reference value.
1018 maxRelDiff = Maximum allowable difference relative to `reference`.
1019 Setting to 0.0 disables this check. Defaults to `1e-2`.
1020 maxAbsDiff = Maximum absolute difference. This is mainly usefull
1021 for comparing values to zero. Setting to 0.0 disables this check.
1022 Defaults to `1e-5`.
1023
1024 Returns:
1025 `true` if `value` is approximately equal to `reference` under
1026 either criterium. It is sufficient, when `value ` satisfies
1027 one of the two criteria.
1028
1029 If one item is a range, and the other is a single value, then
1030 the result is the logical and-ing of calling `approxEqual` on
1031 each element of the ranged item against the single item. If
1032 both items are ranges, then `approxEqual` returns `true` if
1033 and only if the ranges have the same number of elements and if
1034 `approxEqual` evaluates to `true` for each pair of elements.
1035
1036 See_Also:
1037 Use $(LREF feqrel) to get the number of equal bits in the mantissa.
1038 */
1039 deprecated("approxEqual will be removed in 2.106.0. Please use isClose instead.")
approxEqual(T,U,V)1040 bool approxEqual(T, U, V)(T value, U reference, V maxRelDiff = 1e-2, V maxAbsDiff = 1e-5)
1041 {
1042 import core.math : fabs;
1043 import std.range.primitives : empty, front, isInputRange, popFront;
1044 static if (isInputRange!T)
1045 {
1046 static if (isInputRange!U)
1047 {
1048 // Two ranges
1049 for (;; value.popFront(), reference.popFront())
1050 {
1051 if (value.empty) return reference.empty;
1052 if (reference.empty) return value.empty;
1053 if (!approxEqual(value.front, reference.front, maxRelDiff, maxAbsDiff))
1054 return false;
1055 }
1056 }
1057 else static if (isIntegral!U)
1058 {
1059 // convert reference to real
1060 return approxEqual(value, real(reference), maxRelDiff, maxAbsDiff);
1061 }
1062 else
1063 {
1064 // value is range, reference is number
1065 for (; !value.empty; value.popFront())
1066 {
1067 if (!approxEqual(value.front, reference, maxRelDiff, maxAbsDiff))
1068 return false;
1069 }
1070 return true;
1071 }
1072 }
1073 else
1074 {
1075 static if (isInputRange!U)
1076 {
1077 // value is number, reference is range
1078 for (; !reference.empty; reference.popFront())
1079 {
1080 if (!approxEqual(value, reference.front, maxRelDiff, maxAbsDiff))
1081 return false;
1082 }
1083 return true;
1084 }
1085 else static if (isIntegral!T || isIntegral!U)
1086 {
1087 // convert both value and reference to real
1088 return approxEqual(real(value), real(reference), maxRelDiff, maxAbsDiff);
1089 }
1090 else
1091 {
1092 // two numbers
1093 //static assert(is(T : real) && is(U : real));
1094 if (reference == 0)
1095 {
1096 return fabs(value) <= maxAbsDiff;
1097 }
1098 static if (is(typeof(value.infinity)) && is(typeof(reference.infinity)))
1099 {
1100 if (value == value.infinity && reference == reference.infinity ||
1101 value == -value.infinity && reference == -reference.infinity) return true;
1102 }
1103 return fabs((value - reference) / reference) <= maxRelDiff
1104 || maxAbsDiff != 0 && fabs(value - reference) <= maxAbsDiff;
1105 }
1106 }
1107 }
1108
1109 deprecated @safe pure nothrow unittest
1110 {
1111 assert(approxEqual(1.0, 1.0099));
1112 assert(!approxEqual(1.0, 1.011));
1113 assert(approxEqual(0.00001, 0.0));
1114 assert(!approxEqual(0.00002, 0.0));
1115
1116 assert(approxEqual(3.0, [3, 3.01, 2.99])); // several reference values is strange
1117 assert(approxEqual([3, 3.01, 2.99], 3.0)); // better
1118
1119 float[] arr1 = [ 1.0, 2.0, 3.0 ];
1120 double[] arr2 = [ 1.001, 1.999, 3 ];
1121 assert(approxEqual(arr1, arr2));
1122 }
1123
1124 deprecated @safe pure nothrow unittest
1125 {
1126 // relative comparison depends on reference, make sure proper
1127 // side is used when comparing range to single value. Based on
1128 // https://issues.dlang.org/show_bug.cgi?id=15763
1129 auto a = [2e-3 - 1e-5];
1130 auto b = 2e-3 + 1e-5;
1131 assert(a[0].approxEqual(b));
1132 assert(!b.approxEqual(a[0]));
1133 assert(a.approxEqual(b));
1134 assert(!b.approxEqual(a));
1135 }
1136
1137 deprecated @safe pure nothrow @nogc unittest
1138 {
1139 assert(!approxEqual(0.0,1e-15,1e-9,0.0));
1140 assert(approxEqual(0.0,1e-15,1e-9,1e-9));
1141 assert(!approxEqual(1.0,3.0,0.0,1.0));
1142
1143 assert(approxEqual(1.00000000099,1.0,1e-9,0.0));
1144 assert(!approxEqual(1.0000000011,1.0,1e-9,0.0));
1145 }
1146
1147 deprecated @safe pure nothrow @nogc unittest
1148 {
1149 // maybe unintuitive behavior
1150 assert(approxEqual(1000.0,1010.0));
1151 assert(approxEqual(9_090_000_000.0,9_000_000_000.0));
1152 assert(approxEqual(0.0,1e30,1.0));
1153 assert(approxEqual(0.00001,1e-30));
1154 assert(!approxEqual(-1e-30,1e-30,1e-2,0.0));
1155 }
1156
1157 deprecated @safe pure nothrow @nogc unittest
1158 {
1159 int a = 10;
1160 assert(approxEqual(10, a));
1161
1162 assert(!approxEqual(3, 0));
1163 assert(approxEqual(3, 3));
1164 assert(approxEqual(3.0, 3));
1165 assert(approxEqual(3, 3.0));
1166
1167 assert(approxEqual(0.0,0.0));
1168 assert(approxEqual(-0.0,0.0));
1169 assert(approxEqual(0.0f,0.0));
1170 }
1171
1172 deprecated @safe pure nothrow @nogc unittest
1173 {
1174 real num = real.infinity;
1175 assert(num == real.infinity);
1176 assert(approxEqual(num, real.infinity));
1177 num = -real.infinity;
1178 assert(num == -real.infinity);
1179 assert(approxEqual(num, -real.infinity));
1180
1181 assert(!approxEqual(1,real.nan));
1182 assert(!approxEqual(real.nan,real.max));
1183 assert(!approxEqual(real.nan,real.nan));
1184 }
1185
1186 deprecated @safe pure nothrow unittest
1187 {
1188 assert(!approxEqual([1.0,2.0,3.0],[1.0,2.0]));
1189 assert(!approxEqual([1.0,2.0],[1.0,2.0,3.0]));
1190
1191 assert(approxEqual!(real[],real[])([],[]));
1192 assert(approxEqual(cast(real[])[],cast(real[])[]));
1193 }
1194
1195
1196 /**
1197 Computes whether two values are approximately equal, admitting a maximum
1198 relative difference, and a maximum absolute difference.
1199
1200 Params:
1201 lhs = First item to compare.
1202 rhs = Second item to compare.
1203 maxRelDiff = Maximum allowable relative difference.
1204 Setting to 0.0 disables this check. Default depends on the type of
1205 `lhs` and `rhs`: It is approximately half the number of decimal digits of
1206 precision of the smaller type.
1207 maxAbsDiff = Maximum absolute difference. This is mainly usefull
1208 for comparing values to zero. Setting to 0.0 disables this check.
1209 Defaults to `0.0`.
1210
1211 Returns:
1212 `true` if the two items are approximately equal under either criterium.
1213 It is sufficient, when `value ` satisfies one of the two criteria.
1214
1215 If one item is a range, and the other is a single value, then
1216 the result is the logical and-ing of calling `isClose` on
1217 each element of the ranged item against the single item. If
1218 both items are ranges, then `isClose` returns `true` if
1219 and only if the ranges have the same number of elements and if
1220 `isClose` evaluates to `true` for each pair of elements.
1221
1222 See_Also:
1223 Use $(LREF feqrel) to get the number of equal bits in the mantissa.
1224 */
1225 bool isClose(T, U, V = CommonType!(FloatingPointBaseType!T,FloatingPointBaseType!U))
1226 (T lhs, U rhs, V maxRelDiff = CommonDefaultFor!(T,U), V maxAbsDiff = 0.0)
1227 {
1228 import std.range.primitives : empty, front, isInputRange, popFront;
1229 import std.complex : Complex;
1230 static if (isInputRange!T)
1231 {
1232 static if (isInputRange!U)
1233 {
1234 // Two ranges
1235 for (;; lhs.popFront(), rhs.popFront())
1236 {
1237 if (lhs.empty) return rhs.empty;
1238 if (rhs.empty) return lhs.empty;
1239 if (!isClose(lhs.front, rhs.front, maxRelDiff, maxAbsDiff))
1240 return false;
1241 }
1242 }
1243 else
1244 {
1245 // lhs is range, rhs is number
1246 for (; !lhs.empty; lhs.popFront())
1247 {
1248 if (!isClose(lhs.front, rhs, maxRelDiff, maxAbsDiff))
1249 return false;
1250 }
1251 return true;
1252 }
1253 }
1254 else static if (isInputRange!U)
1255 {
1256 // lhs is number, rhs is range
1257 for (; !rhs.empty; rhs.popFront())
1258 {
1259 if (!isClose(lhs, rhs.front, maxRelDiff, maxAbsDiff))
1260 return false;
1261 }
1262 return true;
1263 }
1264 else static if (is(T TE == Complex!TE))
1265 {
1266 static if (is(U UE == Complex!UE))
1267 {
1268 // Two complex numbers
1269 return isClose(lhs.re, rhs.re, maxRelDiff, maxAbsDiff)
1270 && isClose(lhs.im, rhs.im, maxRelDiff, maxAbsDiff);
1271 }
1272 else
1273 {
1274 // lhs is complex, rhs is number
1275 return isClose(lhs.re, rhs, maxRelDiff, maxAbsDiff)
1276 && isClose(lhs.im, 0.0, maxRelDiff, maxAbsDiff);
1277 }
1278 }
1279 else static if (is(U UE == Complex!UE))
1280 {
1281 // lhs is number, rhs is complex
1282 return isClose(lhs, rhs.re, maxRelDiff, maxAbsDiff)
1283 && isClose(0.0, rhs.im, maxRelDiff, maxAbsDiff);
1284 }
1285 else
1286 {
1287 // two numbers
1288 if (lhs == rhs) return true;
1289
1290 static if (is(typeof(lhs.infinity)) && is(typeof(rhs.infinity)))
1291 {
1292 if (lhs == lhs.infinity || rhs == rhs.infinity ||
1293 lhs == -lhs.infinity || rhs == -rhs.infinity) return false;
1294 }
1295
1296 import std.math.algebraic : abs;
1297
1298 auto diff = abs(lhs - rhs);
1299
1300 return diff <= maxRelDiff*abs(lhs)
1301 || diff <= maxRelDiff*abs(rhs)
1302 || diff <= maxAbsDiff;
1303 }
1304 }
1305
1306 ///
1307 @safe pure nothrow @nogc unittest
1308 {
1309 assert(isClose(1.0,0.999_999_999));
1310 assert(isClose(0.001, 0.000_999_999_999));
1311 assert(isClose(1_000_000_000.0,999_999_999.0));
1312
1313 assert(isClose(17.123_456_789, 17.123_456_78));
1314 assert(!isClose(17.123_456_789, 17.123_45));
1315
1316 // use explicit 3rd parameter for less (or more) accuracy
1317 assert(isClose(17.123_456_789, 17.123_45, 1e-6));
1318 assert(!isClose(17.123_456_789, 17.123_45, 1e-7));
1319
1320 // use 4th parameter when comparing close to zero
1321 assert(!isClose(1e-100, 0.0));
1322 assert(isClose(1e-100, 0.0, 0.0, 1e-90));
1323 assert(!isClose(1e-10, -1e-10));
1324 assert(isClose(1e-10, -1e-10, 0.0, 1e-9));
1325 assert(!isClose(1e-300, 1e-298));
1326 assert(isClose(1e-300, 1e-298, 0.0, 1e-200));
1327
1328 // different default limits for different floating point types
1329 assert(isClose(1.0f, 0.999_99f));
1330 assert(!isClose(1.0, 0.999_99));
1331 static if (real.sizeof > double.sizeof)
1332 assert(!isClose(1.0L, 0.999_999_999L));
1333 }
1334
1335 ///
1336 @safe pure nothrow unittest
1337 {
1338 assert(isClose([1.0, 2.0, 3.0], [0.999_999_999, 2.000_000_001, 3.0]));
1339 assert(!isClose([1.0, 2.0], [0.999_999_999, 2.000_000_001, 3.0]));
1340 assert(!isClose([1.0, 2.0, 3.0], [0.999_999_999, 2.000_000_001]));
1341
1342 assert(isClose([2.0, 1.999_999_999, 2.000_000_001], 2.0));
1343 assert(isClose(2.0, [2.0, 1.999_999_999, 2.000_000_001]));
1344 }
1345
1346 @safe pure nothrow unittest
1347 {
1348 assert(!isClose([1.0, 2.0, 3.0], [0.999_999_999, 3.0, 3.0]));
1349 assert(!isClose([2.0, 1.999_999, 2.000_000_001], 2.0));
1350 assert(!isClose(2.0, [2.0, 1.999_999_999, 2.000_000_999]));
1351 }
1352
1353 @safe pure nothrow @nogc unittest
1354 {
1355 immutable a = 1.00001f;
1356 const b = 1.000019;
1357 assert(isClose(a,b));
1358
1359 assert(isClose(1.00001f,1.000019f));
1360 assert(isClose(1.00001f,1.000019));
1361 assert(isClose(1.00001,1.000019f));
1362 assert(!isClose(1.00001,1.000019));
1363
1364 real a1 = 1e-300L;
1365 real a2 = a1.nextUp;
1366 assert(isClose(a1,a2));
1367 }
1368
1369 @safe pure nothrow unittest
1370 {
1371 float[] arr1 = [ 1.0, 2.0, 3.0 ];
1372 double[] arr2 = [ 1.00001, 1.99999, 3 ];
1373 assert(isClose(arr1, arr2));
1374 }
1375
1376 @safe pure nothrow @nogc unittest
1377 {
1378 assert(!isClose(1000.0,1010.0));
1379 assert(!isClose(9_090_000_000.0,9_000_000_000.0));
1380 assert(isClose(0.0,1e30,1.0));
1381 assert(!isClose(0.00001,1e-30));
1382 assert(!isClose(-1e-30,1e-30,1e-2,0.0));
1383 }
1384
1385 @safe pure nothrow @nogc unittest
1386 {
1387 assert(!isClose(3, 0));
1388 assert(isClose(3, 3));
1389 assert(isClose(3.0, 3));
1390 assert(isClose(3, 3.0));
1391
1392 assert(isClose(0.0,0.0));
1393 assert(isClose(-0.0,0.0));
1394 assert(isClose(0.0f,0.0));
1395 }
1396
1397 @safe pure nothrow @nogc unittest
1398 {
1399 real num = real.infinity;
1400 assert(num == real.infinity);
1401 assert(isClose(num, real.infinity));
1402 num = -real.infinity;
1403 assert(num == -real.infinity);
1404 assert(isClose(num, -real.infinity));
1405
1406 assert(!isClose(1,real.nan));
1407 assert(!isClose(real.nan,real.max));
1408 assert(!isClose(real.nan,real.nan));
1409 }
1410
1411 @safe pure nothrow @nogc unittest
1412 {
1413 assert(isClose!(real[],real[],real)([],[]));
1414 assert(isClose(cast(real[])[],cast(real[])[]));
1415 }
1416
1417 @safe pure nothrow @nogc unittest
1418 {
1419 import std.conv : to;
1420
1421 float f = 31.79f;
1422 double d = 31.79;
1423 double f2d = f.to!double;
1424
1425 assert(isClose(f,f2d));
1426 assert(!isClose(d,f2d));
1427 }
1428
1429 @safe pure nothrow @nogc unittest
1430 {
1431 import std.conv : to;
1432
1433 double d = 31.79;
1434 float f = d.to!float;
1435 double f2d = f.to!double;
1436
1437 assert(isClose(f,f2d));
1438 assert(!isClose(d,f2d));
1439 assert(isClose(d,f2d,1e-4));
1440 }
1441
CommonDefaultFor(T,U)1442 package(std.math) template CommonDefaultFor(T,U)
1443 {
1444 import std.algorithm.comparison : min;
1445
1446 alias baseT = FloatingPointBaseType!T;
1447 alias baseU = FloatingPointBaseType!U;
1448
1449 enum CommonType!(baseT, baseU) CommonDefaultFor = 10.0L ^^ -((min(baseT.dig, baseU.dig) + 1) / 2 + 1);
1450 }
1451
FloatingPointBaseType(T)1452 private template FloatingPointBaseType(T)
1453 {
1454 import std.range.primitives : ElementType;
1455 static if (isFloatingPoint!T)
1456 {
1457 alias FloatingPointBaseType = Unqual!T;
1458 }
1459 else static if (isFloatingPoint!(ElementType!(Unqual!T)))
1460 {
1461 alias FloatingPointBaseType = Unqual!(ElementType!(Unqual!T));
1462 }
1463 else
1464 {
1465 alias FloatingPointBaseType = real;
1466 }
1467 }
1468
1469 /***********************************
1470 * Defines a total order on all floating-point numbers.
1471 *
1472 * The order is defined as follows:
1473 * $(UL
1474 * $(LI All numbers in [-$(INFIN), +$(INFIN)] are ordered
1475 * the same way as by built-in comparison, with the exception of
1476 * -0.0, which is less than +0.0;)
1477 * $(LI If the sign bit is set (that is, it's 'negative'), $(NAN) is less
1478 * than any number; if the sign bit is not set (it is 'positive'),
1479 * $(NAN) is greater than any number;)
1480 * $(LI $(NAN)s of the same sign are ordered by the payload ('negative'
1481 * ones - in reverse order).)
1482 * )
1483 *
1484 * Returns:
1485 * negative value if `x` precedes `y` in the order specified above;
1486 * 0 if `x` and `y` are identical, and positive value otherwise.
1487 *
1488 * See_Also:
1489 * $(MYREF isIdentical)
1490 * Standards: Conforms to IEEE 754-2008
1491 */
1492 int cmp(T)(const(T) x, const(T) y) @nogc @trusted pure nothrow
1493 if (isFloatingPoint!T)
1494 {
1495 import std.math : floatTraits, RealFormat;
1496
1497 alias F = floatTraits!T;
1498
1499 static if (F.realFormat == RealFormat.ieeeSingle
1500 || F.realFormat == RealFormat.ieeeDouble)
1501 {
1502 static if (T.sizeof == 4)
1503 alias UInt = uint;
1504 else
1505 alias UInt = ulong;
1506
1507 union Repainter
1508 {
1509 T number;
1510 UInt bits;
1511 }
1512
1513 enum msb = ~(UInt.max >>> 1);
1514
1515 import std.typecons : Tuple;
1516 Tuple!(Repainter, Repainter) vars = void;
1517 vars[0].number = x;
1518 vars[1].number = y;
1519
1520 foreach (ref var; vars)
1521 if (var.bits & msb)
1522 var.bits = ~var.bits;
1523 else
1524 var.bits |= msb;
1525
1526 if (vars[0].bits < vars[1].bits)
1527 return -1;
1528 else if (vars[0].bits > vars[1].bits)
1529 return 1;
1530 else
1531 return 0;
1532 }
1533 else static if (F.realFormat == RealFormat.ieeeExtended53
1534 || F.realFormat == RealFormat.ieeeExtended
1535 || F.realFormat == RealFormat.ieeeQuadruple)
1536 {
1537 static if (F.realFormat == RealFormat.ieeeQuadruple)
1538 alias RemT = ulong;
1539 else
1540 alias RemT = ushort;
1541
1542 struct Bits
1543 {
1544 ulong bulk;
1545 RemT rem;
1546 }
1547
1548 union Repainter
1549 {
1550 T number;
1551 Bits bits;
1552 ubyte[T.sizeof] bytes;
1553 }
1554
1555 import std.typecons : Tuple;
1556 Tuple!(Repainter, Repainter) vars = void;
1557 vars[0].number = x;
1558 vars[1].number = y;
1559
1560 foreach (ref var; vars)
1561 if (var.bytes[F.SIGNPOS_BYTE] & 0x80)
1562 {
1563 var.bits.bulk = ~var.bits.bulk;
1564 var.bits.rem = cast(typeof(var.bits.rem))(-1 - var.bits.rem); // ~var.bits.rem
1565 }
1566 else
1567 {
1568 var.bytes[F.SIGNPOS_BYTE] |= 0x80;
1569 }
1570
version(LittleEndian)1571 version (LittleEndian)
1572 {
1573 if (vars[0].bits.rem < vars[1].bits.rem)
1574 return -1;
1575 else if (vars[0].bits.rem > vars[1].bits.rem)
1576 return 1;
1577 else if (vars[0].bits.bulk < vars[1].bits.bulk)
1578 return -1;
1579 else if (vars[0].bits.bulk > vars[1].bits.bulk)
1580 return 1;
1581 else
1582 return 0;
1583 }
1584 else
1585 {
1586 if (vars[0].bits.bulk < vars[1].bits.bulk)
1587 return -1;
1588 else if (vars[0].bits.bulk > vars[1].bits.bulk)
1589 return 1;
1590 else if (vars[0].bits.rem < vars[1].bits.rem)
1591 return -1;
1592 else if (vars[0].bits.rem > vars[1].bits.rem)
1593 return 1;
1594 else
1595 return 0;
1596 }
1597 }
1598 else
1599 {
1600 // IBM Extended doubledouble does not follow the general
1601 // sign-exponent-significand layout, so has to be handled generically
1602
1603 import std.math.traits : signbit, isNaN;
1604
1605 const int xSign = signbit(x),
1606 ySign = signbit(y);
1607
1608 if (xSign == 1 && ySign == 1)
1609 return cmp(-y, -x);
1610 else if (xSign == 1)
1611 return -1;
1612 else if (ySign == 1)
1613 return 1;
1614 else if (x < y)
1615 return -1;
1616 else if (x == y)
1617 return 0;
1618 else if (x > y)
1619 return 1;
1620 else if (isNaN(x) && !isNaN(y))
1621 return 1;
1622 else if (isNaN(y) && !isNaN(x))
1623 return -1;
1624 else if (getNaNPayload(x) < getNaNPayload(y))
1625 return -1;
1626 else if (getNaNPayload(x) > getNaNPayload(y))
1627 return 1;
1628 else
1629 return 0;
1630 }
1631 }
1632
1633 /// Most numbers are ordered naturally.
1634 @safe unittest
1635 {
1636 assert(cmp(-double.infinity, -double.max) < 0);
1637 assert(cmp(-double.max, -100.0) < 0);
1638 assert(cmp(-100.0, -0.5) < 0);
1639 assert(cmp(-0.5, 0.0) < 0);
1640 assert(cmp(0.0, 0.5) < 0);
1641 assert(cmp(0.5, 100.0) < 0);
1642 assert(cmp(100.0, double.max) < 0);
1643 assert(cmp(double.max, double.infinity) < 0);
1644
1645 assert(cmp(1.0, 1.0) == 0);
1646 }
1647
1648 /// Positive and negative zeroes are distinct.
1649 @safe unittest
1650 {
1651 assert(cmp(-0.0, +0.0) < 0);
1652 assert(cmp(+0.0, -0.0) > 0);
1653 }
1654
1655 /// Depending on the sign, $(NAN)s go to either end of the spectrum.
1656 @safe unittest
1657 {
1658 assert(cmp(-double.nan, -double.infinity) < 0);
1659 assert(cmp(double.infinity, double.nan) < 0);
1660 assert(cmp(-double.nan, double.nan) < 0);
1661 }
1662
1663 /// $(NAN)s of the same sign are ordered by the payload.
1664 @safe unittest
1665 {
1666 assert(cmp(NaN(10), NaN(20)) < 0);
1667 assert(cmp(-NaN(20), -NaN(10)) < 0);
1668 }
1669
1670 @safe unittest
1671 {
1672 import std.meta : AliasSeq;
1673 static foreach (T; AliasSeq!(float, double, real))
1674 {{
1675 T[] values = [-cast(T) NaN(20), -cast(T) NaN(10), -T.nan, -T.infinity,
1676 -T.max, -T.max / 2, T(-16.0), T(-1.0).nextDown,
1677 T(-1.0), T(-1.0).nextUp,
1678 T(-0.5), -T.min_normal, (-T.min_normal).nextUp,
1679 -2 * T.min_normal * T.epsilon,
1680 -T.min_normal * T.epsilon,
1681 T(-0.0), T(0.0),
1682 T.min_normal * T.epsilon,
1683 2 * T.min_normal * T.epsilon,
1684 T.min_normal.nextDown, T.min_normal, T(0.5),
1685 T(1.0).nextDown, T(1.0),
1686 T(1.0).nextUp, T(16.0), T.max / 2, T.max,
1687 T.infinity, T.nan, cast(T) NaN(10), cast(T) NaN(20)];
1688
foreach(i,x;values)1689 foreach (i, x; values)
1690 {
1691 foreach (y; values[i + 1 .. $])
1692 {
1693 assert(cmp(x, y) < 0);
1694 assert(cmp(y, x) > 0);
1695 }
1696 assert(cmp(x, x) == 0);
1697 }
1698 }}
1699 }
1700
1701 package(std): // not yet public
1702
1703 struct FloatingPointBitpattern(T)
1704 if (isFloatingPoint!T)
1705 {
1706 static if (T.mant_dig <= 64)
1707 {
1708 ulong mantissa;
1709 }
1710 else
1711 {
1712 ulong mantissa_lsb;
1713 ulong mantissa_msb;
1714 }
1715
1716 int exponent;
1717 bool negative;
1718 }
1719
1720 FloatingPointBitpattern!T extractBitpattern(T)(T val) @trusted
1721 if (isFloatingPoint!T)
1722 {
1723 import std.math : floatTraits, RealFormat;
1724
1725 FloatingPointBitpattern!T ret;
1726
1727 alias F = floatTraits!T;
1728 static if (F.realFormat == RealFormat.ieeeExtended)
1729 {
1730 if (__ctfe)
1731 {
1732 import core.math : fabs, ldexp;
1733 import std.math.rounding : floor;
1734 import std.math.traits : isInfinity, isNaN, signbit;
1735 import std.math.exponential : log2;
1736
1737 if (isNaN(val) || isInfinity(val))
1738 ret.exponent = 32767;
1739 else if (fabs(val) < real.min_normal)
1740 ret.exponent = 0;
1741 else if (fabs(val) >= nextUp(real.max / 2))
1742 ret.exponent = 32766;
1743 else
1744 ret.exponent = cast(int) (val.fabs.log2.floor() + 16383);
1745
1746 if (ret.exponent == 32767)
1747 {
1748 // NaN or infinity
1749 ret.mantissa = isNaN(val) ? ((1L << 63) - 1) : 0;
1750 }
1751 else
1752 {
1753 auto delta = 16382 + 64 // bias + bits of ulong
1754 - (ret.exponent == 0 ? 1 : ret.exponent); // -1 in case of subnormals
1755 val = ldexp(val, delta); // val *= 2^^delta
1756
1757 ulong tmp = cast(ulong) fabs(val);
1758 if (ret.exponent != 32767 && ret.exponent > 0 && tmp <= ulong.max / 2)
1759 {
1760 // correction, due to log2(val) being rounded up:
1761 ret.exponent--;
1762 val *= 2;
1763 tmp = cast(ulong) fabs(val);
1764 }
1765
1766 ret.mantissa = tmp & long.max;
1767 }
1768
1769 ret.negative = (signbit(val) == 1);
1770 }
1771 else
1772 {
1773 ushort* vs = cast(ushort*) &val;
1774 ret.mantissa = (cast(ulong*) vs)[0] & long.max;
1775 ret.exponent = vs[4] & short.max;
1776 ret.negative = (vs[4] >> 15) & 1;
1777 }
1778 }
1779 else
1780 {
1781 static if (F.realFormat == RealFormat.ieeeSingle)
1782 {
1783 ulong ival = *cast(uint*) &val;
1784 }
1785 else static if (F.realFormat == RealFormat.ieeeDouble)
1786 {
1787 ulong ival = *cast(ulong*) &val;
1788 }
1789 else
1790 {
1791 static assert(false, "Floating point type `" ~ F.realFormat ~ "` not supported.");
1792 }
1793
1794 import std.math.exponential : log2;
1795 enum log2_max_exp = cast(int) log2(T.max_exp);
1796
1797 ret.mantissa = ival & ((1L << (T.mant_dig - 1)) - 1);
1798 ret.exponent = (ival >> (T.mant_dig - 1)) & ((1L << (log2_max_exp + 1)) - 1);
1799 ret.negative = (ival >> (T.mant_dig + log2_max_exp)) & 1;
1800 }
1801
1802 // add leading 1 for normalized values and correct exponent for denormalied values
1803 if (ret.exponent != 0 && ret.exponent != 2 * T.max_exp - 1)
1804 ret.mantissa |= 1L << (T.mant_dig - 1);
1805 else if (ret.exponent == 0)
1806 ret.exponent = 1;
1807
1808 ret.exponent -= T.max_exp - 1;
1809
1810 return ret;
1811 }
1812
1813 @safe pure unittest
1814 {
1815 float f = 1.0f;
1816 auto bp = extractBitpattern(f);
1817 assert(bp.mantissa == 0x80_0000);
1818 assert(bp.exponent == 0);
1819 assert(bp.negative == false);
1820
1821 f = float.max;
1822 bp = extractBitpattern(f);
1823 assert(bp.mantissa == 0xff_ffff);
1824 assert(bp.exponent == 127);
1825 assert(bp.negative == false);
1826
1827 f = -1.5432e-17f;
1828 bp = extractBitpattern(f);
1829 assert(bp.mantissa == 0x8e_55c8);
1830 assert(bp.exponent == -56);
1831 assert(bp.negative == true);
1832
1833 // using double literal due to https://issues.dlang.org/show_bug.cgi?id=20361
1834 f = 2.3822073893521890206e-44;
1835 bp = extractBitpattern(f);
1836 assert(bp.mantissa == 0x00_0011);
1837 assert(bp.exponent == -126);
1838 assert(bp.negative == false);
1839
1840 f = -float.infinity;
1841 bp = extractBitpattern(f);
1842 assert(bp.mantissa == 0);
1843 assert(bp.exponent == 128);
1844 assert(bp.negative == true);
1845
1846 f = float.nan;
1847 bp = extractBitpattern(f);
1848 assert(bp.mantissa != 0); // we don't guarantee payloads
1849 assert(bp.exponent == 128);
1850 assert(bp.negative == false);
1851 }
1852
1853 @safe pure unittest
1854 {
1855 double d = 1.0;
1856 auto bp = extractBitpattern(d);
1857 assert(bp.mantissa == 0x10_0000_0000_0000L);
1858 assert(bp.exponent == 0);
1859 assert(bp.negative == false);
1860
1861 d = double.max;
1862 bp = extractBitpattern(d);
1863 assert(bp.mantissa == 0x1f_ffff_ffff_ffffL);
1864 assert(bp.exponent == 1023);
1865 assert(bp.negative == false);
1866
1867 d = -1.5432e-222;
1868 bp = extractBitpattern(d);
1869 assert(bp.mantissa == 0x11_d9b6_a401_3b04L);
1870 assert(bp.exponent == -737);
1871 assert(bp.negative == true);
1872
1873 d = 0.0.nextUp;
1874 bp = extractBitpattern(d);
1875 assert(bp.mantissa == 0x00_0000_0000_0001L);
1876 assert(bp.exponent == -1022);
1877 assert(bp.negative == false);
1878
1879 d = -double.infinity;
1880 bp = extractBitpattern(d);
1881 assert(bp.mantissa == 0);
1882 assert(bp.exponent == 1024);
1883 assert(bp.negative == true);
1884
1885 d = double.nan;
1886 bp = extractBitpattern(d);
1887 assert(bp.mantissa != 0); // we don't guarantee payloads
1888 assert(bp.exponent == 1024);
1889 assert(bp.negative == false);
1890 }
1891
1892 @safe pure unittest
1893 {
1894 import std.math : floatTraits, RealFormat;
1895
1896 alias F = floatTraits!real;
1897 static if (F.realFormat == RealFormat.ieeeExtended)
1898 {
1899 real r = 1.0L;
1900 auto bp = extractBitpattern(r);
1901 assert(bp.mantissa == 0x8000_0000_0000_0000L);
1902 assert(bp.exponent == 0);
1903 assert(bp.negative == false);
1904
1905 r = real.max;
1906 bp = extractBitpattern(r);
1907 assert(bp.mantissa == 0xffff_ffff_ffff_ffffL);
1908 assert(bp.exponent == 16383);
1909 assert(bp.negative == false);
1910
1911 r = -1.5432e-3333L;
1912 bp = extractBitpattern(r);
1913 assert(bp.mantissa == 0xc768_a2c7_a616_cc22L);
1914 assert(bp.exponent == -11072);
1915 assert(bp.negative == true);
1916
1917 r = 0.0L.nextUp;
1918 bp = extractBitpattern(r);
1919 assert(bp.mantissa == 0x0000_0000_0000_0001L);
1920 assert(bp.exponent == -16382);
1921 assert(bp.negative == false);
1922
1923 r = -float.infinity;
1924 bp = extractBitpattern(r);
1925 assert(bp.mantissa == 0);
1926 assert(bp.exponent == 16384);
1927 assert(bp.negative == true);
1928
1929 r = float.nan;
1930 bp = extractBitpattern(r);
1931 assert(bp.mantissa != 0); // we don't guarantee payloads
1932 assert(bp.exponent == 16384);
1933 assert(bp.negative == false);
1934
1935 r = nextDown(0x1p+16383L);
1936 bp = extractBitpattern(r);
1937 assert(bp.mantissa == 0xffff_ffff_ffff_ffffL);
1938 assert(bp.exponent == 16382);
1939 assert(bp.negative == false);
1940 }
1941 }
1942
1943 @safe pure unittest
1944 {
1945 import std.math : floatTraits, RealFormat;
1946 import std.math.exponential : log2;
1947
1948 alias F = floatTraits!real;
1949
1950 // log2 is broken for x87-reals on some computers in CTFE
1951 // the following test excludes these computers from the test
1952 // (issue 21757)
1953 enum test = cast(int) log2(3.05e2312L);
1954 static if (F.realFormat == RealFormat.ieeeExtended && test == 7681)
1955 {
1956 enum r1 = 1.0L;
1957 enum bp1 = extractBitpattern(r1);
1958 static assert(bp1.mantissa == 0x8000_0000_0000_0000L);
1959 static assert(bp1.exponent == 0);
1960 static assert(bp1.negative == false);
1961
1962 enum r2 = real.max;
1963 enum bp2 = extractBitpattern(r2);
1964 static assert(bp2.mantissa == 0xffff_ffff_ffff_ffffL);
1965 static assert(bp2.exponent == 16383);
1966 static assert(bp2.negative == false);
1967
1968 enum r3 = -1.5432e-3333L;
1969 enum bp3 = extractBitpattern(r3);
1970 static assert(bp3.mantissa == 0xc768_a2c7_a616_cc22L);
1971 static assert(bp3.exponent == -11072);
1972 static assert(bp3.negative == true);
1973
1974 enum r4 = 0.0L.nextUp;
1975 enum bp4 = extractBitpattern(r4);
1976 static assert(bp4.mantissa == 0x0000_0000_0000_0001L);
1977 static assert(bp4.exponent == -16382);
1978 static assert(bp4.negative == false);
1979
1980 enum r5 = -real.infinity;
1981 enum bp5 = extractBitpattern(r5);
1982 static assert(bp5.mantissa == 0);
1983 static assert(bp5.exponent == 16384);
1984 static assert(bp5.negative == true);
1985
1986 enum r6 = real.nan;
1987 enum bp6 = extractBitpattern(r6);
1988 static assert(bp6.mantissa != 0); // we don't guarantee payloads
1989 static assert(bp6.exponent == 16384);
1990 static assert(bp6.negative == false);
1991
1992 enum r7 = nextDown(0x1p+16383L);
1993 enum bp7 = extractBitpattern(r7);
1994 static assert(bp7.mantissa == 0xffff_ffff_ffff_ffffL);
1995 static assert(bp7.exponent == 16382);
1996 static assert(bp7.negative == false);
1997 }
1998 }
1999