Download Reference Manual
The Developer's Library for D
About Wiki Forums Source Search Contact

root/trunk/tango/math/IEEE.d

Revision 5711, 57.4 kB (checked in by lindquist, 1 year ago)

Fix some problems with (now invalid) two operand FSTP instructions in
math code.

  • Property svn:mime-type set to text/x-dsrc
  • Property svn:eol-style set to native
Line 
1 /**
2  * Low-level Mathematical Functions which take advantage of the IEEE754 ABI.
3  *
4  * Copyright: Portions Copyright (C) 2001-2005 Digital Mars.
5  * License:   BSD style: $(LICENSE), Digital Mars.
6  * Authors:   Don Clugston, Walter Bright, Sean Kelly
7  */
8 /* Portions of this code were taken from Phobos std.math, which has the following
9  * copyright notice:
10  *
11  * Author:
12  *  Walter Bright
13  * Copyright:
14  *  Copyright (c) 2001-2005 by Digital Mars,
15  *  All Rights Reserved,
16  *  www.digitalmars.com
17  * License:
18  *  This software is provided 'as-is', without any express or implied
19  *  warranty. In no event will the authors be held liable for any damages
20  *  arising from the use of this software.
21  *
22  *  Permission is granted to anyone to use this software for any purpose,
23  *  including commercial applications, and to alter it and redistribute it
24  *  freely, subject to the following restrictions:
25  *
26  *  <ul>
27  *  <li> The origin of this software must not be misrepresented; you must not
28  *       claim that you wrote the original software. If you use this software
29  *       in a product, an acknowledgment in the product documentation would be
30  *       appreciated but is not required.
31  *  </li>
32  *  <li> Altered source versions must be plainly marked as such, and must not
33  *       be misrepresented as being the original software.
34  *  </li>
35  *  <li> This notice may not be removed or altered from any source
36  *       distribution.
37  *  </li>
38  *  </ul>
39  */
40 /**
41  * Macros:
42  *
43  *  TABLE_SV = <table border=1 cellpadding=4 cellspacing=0>
44  *      <caption>Special Values</caption>
45  *      $0</table>
46  *  SVH = $(TR $(TH $1) $(TH $2))
47  *  SV  = $(TR $(TD $1) $(TD $2))
48  *  SVH3 = $(TR $(TH $1) $(TH $2) $(TH $3))
49  *  SV3  = $(TR $(TD $1) $(TD $2) $(TD $3))
50  *  NAN = $(RED NAN)
51  *  PLUSMN = &plusmn;
52  *  INFIN = &infin;
53  *  PLUSMNINF = &plusmn;&infin;
54  *  PI = &pi;
55  *  LT = &lt;
56  *  GT = &gt;
57  *  SQRT = &radix;
58  *  HALF = &frac12;
59  */
60 module tango.math.IEEE;
61
62 version(GNU){
63     // GDC is a filthy liar. It can't actually do inline asm.
64 } else version(TangoNoAsm) {
65
66 } else version(D_InlineAsm_X86) {
67     version = Naked_D_InlineAsm_X86;
68 }
69
70 version (X86){
71     version = X86_Any;
72 }
73
74 version (X86_64){
75     version = X86_Any;
76 }
77
78 version (Naked_D_InlineAsm_X86) {
79     // Don't include this extra dependency unless we need to.
80     debug(UnitTest) {
81         static import tango.stdc.math;
82     }
83 } else {
84     // Needed for cos(), sin(), tan() on GNU.
85     static import tango.stdc.math;
86 }
87
88
89 version(Windows) {
90     version(DigitalMars) {
91         version = DMDWindows;
92     }
93 }
94    
95 // Standard Tango NaN payloads.
96 // NOTE: These values may change in future Tango releases
97 // The lowest three bits indicate the cause of the NaN:
98 // 0 = error other than those listed below:
99 // 1 = domain error
100 // 2 = singularity
101 // 3 = range
102 // 4-7 = reserved.
103 enum TANGO_NAN {
104     // General errors
105     DOMAIN_ERROR = 0x0101,
106     SINGULARITY  = 0x0102,
107     RANGE_ERROR  = 0x0103,
108     // NaNs created by functions in the basic library
109     TAN_DOMAIN   = 0x1001,
110     POW_DOMAIN   = 0x1021,
111     GAMMA_DOMAIN = 0x1101,
112     GAMMA_POLE   = 0x1102,
113     SGNGAMMA     = 0x1112,
114     BETA_DOMAIN  = 0x1131,
115     // NaNs from statistical functions
116     NORMALDISTRIBUTION_INV_DOMAIN = 0x2001,
117     STUDENTSDDISTRIBUTION_DOMAIN  = 0x2011
118 }
119
120 private:
121 /* Most of the functions depend on the format of the largest IEEE floating-point type.
122  * These code will differ depending on whether 'real' is 64, 80, or 128 bits,
123  * and whether it is a big-endian or little-endian architecture.
124  * Only five 'real' ABIs are currently supported:
125  * 64 bit Big-endian  'double' (eg PowerPC)
126  * 128 bit Big-endian 'quadruple' (eg SPARC)
127  * 64 bit Little-endian 'double' (eg x86-SSE2)
128  * 80 bit Little-endian, with implied bit 'real80' (eg x87, Itanium).
129  * 128 bit Little-endian 'quadruple' (not implemented on any known processor!)
130  *
131  * There is also an unsupported ABI which does not follow IEEE; several of its functions
132  *  will generate run-time errors if used.
133  * 128 bit Big-endian 'doubledouble' (used by GDC <= 0.23 for PowerPC)
134  */
135
136 version(LittleEndian) {
137     static assert(real.mant_dig == 53 || real.mant_dig==64 || real.mant_dig == 113,
138         "Only 64-bit, 80-bit, and 128-bit reals are supported for LittleEndian CPUs");
139 } else {
140     static assert(real.mant_dig == 53 || real.mant_dig==106 || real.mant_dig == 113,
141      "Only 64-bit and 128-bit reals are supported for BigEndian CPUs. double-double reals have partial support");
142 }
143
144 // Constants used for extracting the components of the representation.
145 // They supplement the built-in floating point properties.
146 template floatTraits(T) {
147  // EXPMASK is a ushort mask to select the exponent portion (without sign)
148  // SIGNMASK is a ushort mask to select the sign bit.
149  // EXPPOS_SHORT is the index of the exponent when represented as a ushort array.
150  // SIGNPOS_BYTE is the index of the sign when represented as a ubyte array.
151  // RECIP_EPSILON is the value such that (smallest_denormal) * RECIP_EPSILON == T.min
152  const T RECIP_EPSILON = (1/T.epsilon);
153
154  static if (T.mant_dig == 24) { // float
155     enum : ushort {
156         EXPMASK = 0x7F80,
157         SIGNMASK = 0x8000,
158         EXPBIAS = 0x3F00
159     }
160     const uint EXPMASK_INT = 0x7F80_0000;
161     const uint MANTISSAMASK_INT = 0x007F_FFFF;
162     version(LittleEndian) {       
163       const EXPPOS_SHORT = 1;
164     } else {
165       const EXPPOS_SHORT = 0;
166     }
167  } else static if (T.mant_dig==53) { // double, or real==double
168      enum : ushort {
169          EXPMASK = 0x7FF0,
170          SIGNMASK = 0x8000,
171          EXPBIAS = 0x3FE0
172     }
173     const uint EXPMASK_INT = 0x7FF0_0000;
174     const uint MANTISSAMASK_INT = 0x000F_FFFF; // for the MSB only
175     version(LittleEndian) {
176       const EXPPOS_SHORT = 3;
177       const SIGNPOS_BYTE = 7;
178     } else {
179       const EXPPOS_SHORT = 0;
180       const SIGNPOS_BYTE = 0;
181     }
182  } else static if (T.mant_dig==64) { // real80
183      enum : ushort {
184          EXPMASK = 0x7FFF,
185          SIGNMASK = 0x8000,
186          EXPBIAS = 0x3FFE
187      }
188 //    const ulong QUIETNANMASK = 0xC000_0000_0000_0000; // Converts a signaling NaN to a quiet NaN.
189     version(LittleEndian) {
190       const EXPPOS_SHORT = 4;
191       const SIGNPOS_BYTE = 9;
192     } else {
193       const EXPPOS_SHORT = 0;
194       const SIGNPOS_BYTE = 0;
195     }
196  } else static if (real.mant_dig==113){ // quadruple
197      enum : ushort {
198          EXPMASK = 0x7FFF,
199          SIGNMASK = 0x8000,
200          EXPBIAS = 0x3FFE
201      }
202     version(LittleEndian) {
203       const EXPPOS_SHORT = 7;
204       const SIGNPOS_BYTE = 15;
205     } else {
206       const EXPPOS_SHORT = 0;
207       const SIGNPOS_BYTE = 0;
208     }
209  } else static if (real.mant_dig==106) { // doubledouble
210      enum : ushort {
211          EXPMASK = 0x7FF0,
212          SIGNMASK = 0x8000
213 //         EXPBIAS = 0x3FE0
214      }
215     // the exponent byte is not unique
216     version(LittleEndian) {
217       const EXPPOS_SHORT = 7; // 3 is also an exp short
218       const SIGNPOS_BYTE = 15;
219     } else {
220       const EXPPOS_SHORT = 0; // 4 is also an exp short
221       const SIGNPOS_BYTE = 0;
222     }
223  }
224 }
225
226 // These apply to all floating-point types
227 version(LittleEndian) {
228     const MANTISSA_LSB = 0;
229     const MANTISSA_MSB = 1;   
230 } else {
231     const MANTISSA_LSB = 1;
232     const MANTISSA_MSB = 0;
233 }
234
235 public:
236
237 /** IEEE exception status flags
238
239  These flags indicate that an exceptional floating-point condition has occured.
240  They indicate that a NaN or an infinity has been generated, that a result
241  is inexact, or that a signalling NaN has been encountered.
242  The return values of the properties should be treated as booleans, although
243  each is returned as an int, for speed.
244
245  Example:
246  ----
247     real a=3.5;
248     // Set all the flags to zero
249     resetIeeeFlags();
250     assert(!ieeeFlags.divByZero);
251     // Perform a division by zero.
252     a/=0.0L;
253     assert(a==real.infinity);
254     assert(ieeeFlags.divByZero);
255     // Create a NaN
256     a*=0.0L;
257     assert(ieeeFlags.invalid);
258     assert(isNaN(a));
259
260     // Check that calling func() has no effect on the
261     // status flags.
262     IeeeFlags f = ieeeFlags;
263     func();
264     assert(ieeeFlags == f);
265
266  ----
267  */
268 struct IeeeFlags
269 {
270 private:
271     // The x87 FPU status register is 16 bits.
272     // The Pentium SSE2 status register is 32 bits.
273     int m_flags;
274     version (X86_Any) {
275         // Applies to both x87 status word (16 bits) and SSE2 status word(32 bits).
276         enum : int {
277             INEXACT_MASK   = 0x20,
278             UNDERFLOW_MASK = 0x10,
279             OVERFLOW_MASK  = 0x08,
280             DIVBYZERO_MASK = 0x04,
281             INVALID_MASK   = 0x01
282         }
283         // Don't bother about denormals, they are not supported on most CPUs.
284         //  DENORMAL_MASK = 0x02;
285     } else version (PPC) {
286         // PowerPC FPSCR is a 32-bit register.
287         enum : int {
288             INEXACT_MASK   = 0x600,
289             UNDERFLOW_MASK = 0x010,
290             OVERFLOW_MASK  = 0x008,
291             DIVBYZERO_MASK = 0x020,
292             INVALID_MASK   = 0xF80
293         }
294     } else { // SPARC FSR is a 32bit register
295              //(64 bits for Sparc 7 & 8, but high 32 bits are uninteresting).
296         enum : int {
297             INEXACT_MASK   = 0x020,
298             UNDERFLOW_MASK = 0x080,
299             OVERFLOW_MASK  = 0x100,
300             DIVBYZERO_MASK = 0x040,
301             INVALID_MASK   = 0x200
302         }
303     }
304 private:
305     static IeeeFlags getIeeeFlags()
306     {
307         // This is a highly time-critical operation, and
308         // should really be an intrinsic.
309         version(D_InlineAsm_X86) {
310             version(DMDWindows) {
311              // In this case, we
312              // take advantage of the fact that for DMD-Windows
313              // a struct containing only a int is returned in EAX.
314              asm {
315                  fstsw AX;
316                  // NOTE: If compiler supports SSE2, need to OR the result with
317                  // the SSE2 status register.
318                  // Clear all irrelevant bits
319                  and EAX, 0x03D;
320              }
321            }
322            else {
323              IeeeFlags tmp1;
324              asm {
325                  fstsw AX;
326                  // NOTE: If compiler supports SSE2, need to OR the result with
327                  // the SSE2 status register.
328                  // Clear all irrelevant bits
329                  and EAX, 0x03D;
330                  mov tmp1, EAX;
331              }
332              return tmp1;
333            }
334        } else version (PPC) {
335            assert(0, "Not yet supported");
336        } else {
337            /*   SPARC:
338                int retval;
339                asm { st %fsr, retval; }
340                return retval;
341             */
342            assert(0, "Not yet supported");
343        }
344     }
345     static void resetIeeeFlags()
346     {
347        version(D_InlineAsm_X86) {
348             asm {
349                 fnclex;
350             }
351         } else {
352             /* SPARC:
353               int tmpval;
354               asm { st %fsr, tmpval; }
355               tmpval &=0xFFFF_FC00;
356               asm { ld tmpval, %fsr; }
357             */
358            assert(0, "Not yet supported");
359         }
360     }
361 public:
362     /// The result cannot be represented exactly, so rounding occured.
363     /// (example: x = sin(0.1); }
364     int inexact() { return m_flags & INEXACT_MASK; }
365     /// A zero was generated by underflow (example: x = real.min*real.epsilon/2;)
366     int underflow() { return m_flags & UNDERFLOW_MASK; }
367     /// An infinity was generated by overflow (example: x = real.max*2;)
368     int overflow() { return m_flags & OVERFLOW_MASK; }
369     /// An infinity was generated by division by zero (example: x = 3/0.0; )
370     int divByZero() { return m_flags & DIVBYZERO_MASK; }
371     /// A machine NaN was generated. (example: x = real.infinity * 0.0; )
372     int invalid() { return m_flags & INVALID_MASK; }
373 }
374
375 /// Return a snapshot of the current state of the floating-point status flags.
376 IeeeFlags ieeeFlags() { return IeeeFlags.getIeeeFlags(); }
377
378 /// Set all of the floating-point status flags to false.
379 void resetIeeeFlags() { IeeeFlags.resetIeeeFlags; }
380
381 /** IEEE rounding modes.
382  * The default mode is ROUNDTONEAREST.
383  */
384 enum RoundingMode : short {
385     ROUNDTONEAREST = 0x0000,
386     ROUNDDOWN      = 0x0400,
387     ROUNDUP        = 0x0800,
388     ROUNDTOZERO    = 0x0C00
389 };
390
391 /** Change the rounding mode used for all floating-point operations.
392  *
393  * Returns the old rounding mode.
394  *
395  * When changing the rounding mode, it is almost always necessary to restore it
396  * at the end of the function. Typical usage:
397 ---
398     auto oldrounding = setIeeeRounding(RoundingMode.ROUNDDOWN);
399     scope (exit) setIeeeRounding(oldrounding);
400 ---
401  */
402 RoundingMode setIeeeRounding(RoundingMode roundingmode) {
403    version(D_InlineAsm_X86) {
404         // TODO: For SSE/SSE2, do we also need to set the SSE rounding mode?
405         short cont;
406         asm {
407             fstcw cont;
408             mov CX, cont;
409             mov AX, cont;
410             and EAX, 0x0C00; // Form the return value
411             and CX, 0xF3FF;
412             or CX, roundingmode;
413             mov cont, CX;
414             fldcw cont;
415         }
416     } else {
417            assert(0, "Not yet supported");
418     }
419 }
420
421 /** Get the IEEE rounding mode which is in use.
422  *
423  */
424 RoundingMode getIeeeRounding() {
425    version(D_InlineAsm_X86) {
426         // TODO: For SSE/SSE2, do we also need to check the SSE rounding mode?
427         short cont;
428         asm {
429             mov EAX, 0x0C00;
430             fstcw cont;
431             and AX, cont;
432         }
433     } else {
434            assert(0, "Not yet supported");
435     }
436 }
437
438 debug(UnitTest) {
439    version(D_InlineAsm_X86) { // Won't work for anything else yet
440 unittest {
441     real a = 3.5;
442     resetIeeeFlags();
443     assert(!ieeeFlags.divByZero);
444     a /= 0.0L;
445     assert(ieeeFlags.divByZero);
446     assert(a == real.infinity);
447     a *= 0.0L;
448     assert(ieeeFlags.invalid);
449     assert(isNaN(a));
450     a = real.max;
451     a *= 2;
452     assert(ieeeFlags.overflow);
453     a = real.min * real.epsilon;
454     a /= 99;
455     assert(ieeeFlags.underflow);
456     assert(ieeeFlags.inexact);
457
458     int r = getIeeeRounding;
459     assert(r == RoundingMode.ROUNDTONEAREST);
460 }
461 }
462 }
463
464 // Note: Itanium supports more precision options than this. SSE/SSE2 does not support any.
465 enum PrecisionControl : short {
466     PRECISION80 = 0x300,
467     PRECISION64 = 0x200,
468     PRECISION32 = 0x000
469 };
470
471 /** Set the number of bits of precision used by 'real'.
472  *
473  * Returns: the old precision.
474  * This is not supported on all platforms.
475  */
476 PrecisionControl reduceRealPrecision(PrecisionControl prec) {
477    version(D_InlineAsm_X86) {
478         short cont;
479         asm {
480             fstcw cont;
481             mov CX, cont;
482             mov AX, cont;
483             and EAX, 0x0300; // Form the return value
484             and CX,  0xFCFF;
485             or  CX,  prec;
486             mov cont, CX;
487             fldcw cont;
488         }
489     } else {
490            assert(0, "Not yet supported");
491     }
492 }
493
494 /*********************************************************************
495  * Separate floating point value into significand and exponent.
496  *
497  * Returns:
498  *      Calculate and return $(I x) and $(I exp) such that
499  *      value =$(I x)*2$(SUP exp) and
500  *      .5 $(LT)= |$(I x)| $(LT) 1.0
501  *     
502  *      $(I x) has same sign as value.
503  *
504  *      $(TABLE_SV
505  *      $(TR $(TH value)           $(TH returns)         $(TH exp))
506  *      $(TR $(TD $(PLUSMN)0.0)    $(TD $(PLUSMN)0.0)    $(TD 0))
507  *      $(TR $(TD +$(INFIN))       $(TD +$(INFIN))       $(TD int.max))
508  *      $(TR $(TD -$(INFIN))       $(TD -$(INFIN))       $(TD int.min))
509  *      $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min))
510  *      )
511  */
512 real frexp(real value, out int exp)
513 {
514     ushort* vu = cast(ushort*)&value;
515     long* vl = cast(long*)&value;
516     uint ex;
517     alias floatTraits!(real) F;
518
519     ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
520   static if (real.mant_dig == 64) { // real80
521     if (ex) { // If exponent is non-zero
522         if (ex == F.EXPMASK) {   // infinity or NaN
523             if (*vl &  0x7FFF_FFFF_FFFF_FFFF) {  // NaN
524                 *vl |= 0xC000_0000_0000_0000;  // convert $(NAN)S to $(NAN)Q
525                 exp = int.min;
526             } else if (vu[F.EXPPOS_SHORT] & 0x8000) {   // negative infinity
527                 exp = int.min;
528             } else {   // positive infinity
529                 exp = int.max;
530             }
531         } else {
532             exp = ex - F.EXPBIAS;
533             vu[F.EXPPOS_SHORT] = cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE);
534         }
535     } else if (!*vl) {
536         // value is +-0.0
537         exp = 0;
538     } else {
539         // denormal
540         value *= F.RECIP_EPSILON;
541         ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
542         exp = ex - F.EXPBIAS - 63;
543         vu[F.EXPPOS_SHORT] = cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE);
544     }
545     return value;
546   } else static if (real.mant_dig == 113) { // quadruple     
547         if (ex) { // If exponent is non-zero
548             if (ex == F.EXPMASK) {   // infinity or NaN
549                 if (vl[MANTISSA_LSB] |( vl[MANTISSA_MSB]&0x0000_FFFF_FFFF_FFFF)) {  // NaN
550                     vl[MANTISSA_MSB] |= 0x0000_8000_0000_0000;  // convert $(NAN)S to $(NAN)Q
551                     exp = int.min;
552                 } else if (vu[F.EXPPOS_SHORT] & 0x8000) {   // negative infinity
553                     exp = int.min;
554                 } else {   // positive infinity
555                     exp = int.max;
556                 }
557             } else {
558                 exp = ex - F.EXPBIAS;
559                 vu[F.EXPPOS_SHORT] = cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE);
560             }
561         } else if ((vl[MANTISSA_LSB] |(vl[MANTISSA_MSB]&0x0000_FFFF_FFFF_FFFF))==0) {
562             // value is +-0.0
563             exp = 0;
564     } else {
565         // denormal
566         value *= F.RECIP_EPSILON;
567         ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
568         exp = ex - F.EXPBIAS - 113;
569         vu[F.EXPPOS_SHORT] = cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE);
570     }
571     return value;
572   } else static if (real.mant_dig==53) { // real is double
573     if (ex) { // If exponent is non-zero
574         if (ex == F.EXPMASK) {   // infinity or NaN
575             if (*vl==0x7FF0_0000_0000_0000) {  // positive infinity
576                 exp = int.max;
577             } else if (*vl==0xFFF0_0000_0000_0000) { // negative infinity
578                 exp = int.min;
579             } else { // NaN
580                 *vl |= 0x0008_0000_0000_0000;  // convert $(NAN)S to $(NAN)Q
581                 exp = int.min;
582             }
583         } else {
584             exp = (ex - F.EXPBIAS) >>> 4;
585             vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FE0;
586         }
587     } else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF)) {
588         // value is +-0.0
589         exp = 0;
590     } else {
591         // denormal
592         ushort sgn;
593         sgn = (0x8000 & vu[F.EXPPOS_SHORT])| 0x3FE0;
594         *vl &= 0x7FFF_FFFF_FFFF_FFFF;
595
596         int i = -0x3FD+11;
597         do {
598             i--;
599             *vl <<= 1;
600         } while (*vl > 0);
601         exp = i;
602         vu[F.EXPPOS_SHORT] = sgn;
603     }
604     return value;
605   }else { //static if(real.mant_dig==106) // doubledouble
606         assert(0, "Unsupported");
607   }
608 }
609
610 debug(UnitTest) {
611
612 unittest
613 {
614     static real vals[][3] = // x,frexp,exp
615     [
616         [0.0,   0.0,    0],
617         [-0.0,  -0.0,   0],
618         [1.0,   .5, 1],
619         [-1.0,  -.5,    1],
620         [2.0,   .5, 2],
621         [double.min/2.0, .5, -1022],
622         [real.infinity,real.infinity,int.max],
623         [-real.infinity,-real.infinity,int.min],
624     ];   
625    
626     int i;
627     int eptr;
628     real v = frexp(NaN(0xABC), eptr);
629     assert(isIdentical(NaN(0xABC), v));
630     assert(eptr ==int.min);
631     v = frexp(-NaN(0xABC), eptr);
632     assert(isIdentical(-NaN(0xABC), v));
633     assert(eptr ==int.min);
634
635     for (i = 0; i < vals.length; i++) {
636         real x = vals[i][0];
637         real e = vals[i][1];
638         int exp = cast(int)vals[i][2];
639         v = frexp(x, eptr);
640 //        printf("frexp(%La) = %La, should be %La, eptr = %d, should be %d\n", x, v, e, eptr, exp);
641         assert(isIdentical(e, v));
642         assert(exp == eptr);
643
644     }
645    static if (real.mant_dig == 64) {
646      static real extendedvals[][3] = [ // x,frexp,exp
647         [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L,   74],    // normal
648         [0x1.fa01712e8f0471ap-1064L,  0x1.fa01712e8f0471ap-1L,     -1063],
649         [real.min,  .5,     -16381],
650         [real.min/2.0L, .5,     -16382]    // denormal
651      ];
652
653     for (i = 0; i < extendedvals.length; i++) {
654         real x = extendedvals[i][0];
655         real e = extendedvals[i][1];
656         int exp = cast(int)extendedvals[i][2];
657         v = frexp(x, eptr);
658         assert(isIdentical(e, v));
659         assert(exp == eptr);
660
661     }
662   }
663 }
664 }
665
666 /**
667  * Compute n * 2$(SUP exp)
668  * References: frexp
669  */
670 real ldexp(real n, int exp) /* intrinsic */
671 {
672     version(Naked_D_InlineAsm_X86)
673     {
674         asm {
675             fild exp;
676             fld n;
677             fscale;
678             fstp ST(1);
679         }
680     }
681     else
682     {
683         return tango.stdc.math.ldexpl(n, exp);
684     }
685 }
686
687 /******************************************
688  * Extracts the exponent of x as a signed integral value.
689  *
690  * If x is not a special value, the result is the same as
691  * $(D cast(int)logb(x)).
692  *
693  * Remarks: This function is consistent with IEEE754R, but it
694  * differs from the C function of the same name
695  * in the return value of infinity. (in C, ilogb(real.infinity)== int.max).
696  * Note that the special return values may all be equal.
697  *
698  *      $(TABLE_SV
699  *      $(TR $(TH x)                $(TH ilogb(x))     $(TH Invalid?))
700  *      $(TR $(TD 0)                 $(TD FP_ILOGB0)   $(TD yes))
701  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD FP_ILOGBINFINITY) $(TD yes))
702  *      $(TR $(TD $(NAN))            $(TD FP_ILOGBNAN) $(TD yes))
703  *      )
704  */
705 int ilogb(real x)
706 {
707         version(Naked_D_InlineAsm_X86)
708         {
709             int y;
710             asm {
711                 fld x;
712                 fxtract;
713                 fstp ST(0); // drop significand
714                 fistp y; // and return the exponent
715             }
716             return y;
717         } else static if (real.mant_dig==64) { // 80-bit reals
718             alias floatTraits!(real) F;
719             short e = cast(short)((cast(short *)&x)[F.EXPPOS_SHORT] & F.EXPMASK);
720             if (e == F.EXPMASK) {
721                 // BUG: should also set the invalid exception
722                 ulong s = *cast(ulong *)&x;
723                 if (s == 0x8000_0000_0000_0000) {
724                     return FP_ILOGBINFINITY;
725                 }
726                 else return FP_ILOGBNAN;
727             }
728             if (e==0) {
729                 ulong s = *cast(ulong *)&x;
730                 if (s == 0x0000_0000_0000_0000) {
731                     // BUG: should also set the invalid exception
732                     return FP_ILOGB0;
733                 }
734                 // Denormals
735                 x *= F.RECIP_EPSILON;
736                 short f = (cast(short *)&x)[F.EXPPOS_SHORT];
737                 return -0x3FFF - (63-f);
738             }
739             return e - 0x3FFF;
740         } else {
741         return tango.stdc.math.ilogbl(x);
742     }
743 }
744
745 version (X86)
746 {
747     const int FP_ILOGB0        = -int.max-1;
748     const int FP_ILOGBNAN      = -int.max-1;
749     const int FP_ILOGBINFINITY = -int.max-1;
750 } else {
751     alias tango.stdc.math.FP_ILOGB0   FP_ILOGB0;
752     alias tango.stdc.math.FP_ILOGBNAN FP_ILOGBNAN;
753     const int FP_ILOGBINFINITY = int.max;
754 }
755
756 debug(UnitTest) {
757 unittest {
758     assert(ilogb(1.0) == 0);
759     assert(ilogb(65536) == 16);
760     assert(ilogb(-65536) == 16);
761     assert(ilogb(1.0 / 65536) == -16);
762     assert(ilogb(real.nan) == FP_ILOGBNAN);
763     assert(ilogb(0.0) == FP_ILOGB0);
764     assert(ilogb(-0.0) == FP_ILOGB0);
765     // denormal
766     assert(ilogb(0.125 * real.min) == real.min_exp - 4);
767     assert(ilogb(real.infinity) == FP_ILOGBINFINITY);
768 }
769 }
770
771 /*****************************************
772  * Extracts the exponent of x as a signed integral value.
773  *
774  * If x is subnormal, it is treated as if it were normalized.
775  * For a positive, finite x:
776  *
777  * 1 $(LT)= $(I x) * FLT_RADIX$(SUP -logb(x)) $(LT) FLT_RADIX
778  *
779  *      $(TABLE_SV
780  *      $(TR $(TH x)                 $(TH logb(x))   $(TH divide by 0?) )
781  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no))
782  *      $(TR $(TD $(PLUSMN)0.0)      $(TD -$(INFIN)) $(TD yes) )
783  *      )
784  */
785 real logb(real x)
786 {
787     version(Naked_D_InlineAsm_X86)
788     {
789         asm {
790             fld x;
791             fxtract;
792             fstp ST(0); // drop significand
793         }
794     } else {
795         return tango.stdc.math.logbl(x);
796     }
797 }
798
799 debug(UnitTest) {
800 unittest {
801     assert(logb(real.infinity)== real.infinity);
802     assert(isIdentical(logb(NaN(0xFCD)), NaN(0xFCD)));
803     assert(logb(1.0)== 0.0);
804     assert(logb(-65536) == 16);
805     assert(logb(0.0)== -real.infinity);
806     assert(ilogb(0.125*real.min) == real.min_exp-4);
807 }
808 }
809
810 /*************************************
811  * Efficiently calculates x * 2$(SUP n).
812  *
813  * scalbn handles underflow and overflow in
814  * the same fashion as the basic arithmetic operators.
815  *
816  *  $(TABLE_SV
817  *      $(TR $(TH x)                 $(TH scalb(x)))
818  *      $(TR $(TD $(PLUSMNINF))      $(TD $(PLUSMNINF)) )
819  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) )
820  *  )
821  */
822 real scalbn(real x, int n)
823 {
824     version(Naked_D_InlineAsm_X86)
825     {
826         asm {
827             fild n;
828             fld x;
829             fscale;
830             fstp ST(1);
831         }
832     } else {
833         // NOTE: Not implemented in DMD
834         return tango.stdc.math.scalbnl(x, n);
835     }
836 }
837
838 debug(UnitTest) {
839 unittest {
840     assert(scalbn(-real.infinity, 5) == -real.infinity);
841     assert(isIdentical(scalbn(NaN(0xABC),7), NaN(0xABC)));
842 }
843 }
844
845 /**
846  * Returns the positive difference between x and y.
847  *
848  * If either of x or y is $(NAN), it will be returned.
849  * Returns:
850  * $(TABLE_SV
851  *  $(SVH Arguments, fdim(x, y))
852  *  $(SV x $(GT) y, x - y)
853  *  $(SV x $(LT)= y, +0.0)
854  * )
855  */
856 real fdim(real x, real y)
857 {
858     return (x !<= y) ? x - y : +0.0;
859 }
860
861 debug(UnitTest) {
862 unittest {
863     assert(isIdentical(fdim(NaN(0xABC), 58.2), NaN(0xABC)));
864 }
865 }
866
867 /*******************************
868  * Returns |x|
869  *
870  *      $(TABLE_SV
871  *      $(TR $(TH x)                 $(TH fabs(x)))
872  *      $(TR $(TD $(PLUSMN)0.0)      $(TD +0.0) )
873  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) )
874  *      )
875  */
876 real fabs(real x) /* intrinsic */
877 {
878     version(D_InlineAsm_X86)
879     {
880         asm {
881             fld x;
882             fabs;
883         }
884     }
885     else
886     {
887         return tango.stdc.math.fabsl(x);
888     }
889 }
890
891 unittest {
892     assert(isIdentical(fabs(NaN(0xABC)), NaN(0xABC)));
893 }
894
895 /**
896  * Returns (x * y) + z, rounding only once according to the
897  * current rounding mode.
898  *
899  * BUGS: Not currently implemented - rounds twice.
900  */
901 real fma(float x, float y, float z)
902 {
903     return (x * y) + z;
904 }
905
906 /**
907  * Calculate cos(y) + i sin(y).
908  *
909  * On x86 CPUs, this is a very efficient operation;
910  * almost twice as fast as calculating sin(y) and cos(y)
911  * seperately, and is the preferred method when both are required.
912  */
913 creal expi(real y)
914 {
915     version(Naked_D_InlineAsm_X86)
916     {
917         asm {           
918             fld y;
919             fsincos;
920             fxch ST(1), ST(0);
921         }
922     }
923     else
924     {
925         return tango.stdc.math.cosl(y) + tango.stdc.math.sinl(y)*1i;
926     }
927 }
928
929 debug(UnitTest) {
930 unittest
931 {
932     assert(expi(1.3e5L) == tango.stdc.math.cosl(1.3e5L) + tango.stdc.math.sinl(1.3e5L) * 1i);
933     assert(expi(0.0L) == 1L + 0.0Li);
934 }
935 }
936
937 /*********************************
938  * Returns !=0 if e is a NaN.
939  */
940
941 int isNaN(real x)
942 {
943   alias floatTraits!(real) F;
944   static if (real.mant_dig==53) { // double
945         ulong*  p = cast(ulong *)&x;
946         return (*p & 0x7FF0_0000_0000_0000 == 0x7FF0_0000_0000_0000) && *p & 0x000F_FFFF_FFFF_FFFF;
947   } else static if (real.mant_dig==64) {     // real80
948         ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
949         ulong*  ps = cast(ulong *)&x;
950         return e == F.EXPMASK &&
951             *ps & 0x7FFF_FFFF_FFFF_FFFF; // not infinity
952   } else static if (real.mant_dig==113) {  // quadruple
953         ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
954         ulong*  ps = cast(ulong *)&x;
955         return e == F.EXPMASK &&
956            (ps[MANTISSA_LSB] | (ps[MANTISSA_MSB]& 0x0000_FFFF_FFFF_FFFF))!=0;
957   } else {
958       return x!=x;
959   }
960 }
961
962
963 debug(UnitTest) {
964 unittest
965 {
966     assert(isNaN(float.nan));
967     assert(isNaN(-double.nan));
968     assert(isNaN(real.nan));
969
970     assert(!isNaN(53.6));
971     assert(!isNaN(float.infinity));
972 }
973 }
974
975 /**
976  * Returns !=0 if x is normalized.
977  *
978  * (Need one for each format because subnormal
979  *  floats might be converted to normal reals)
980  */
981 int isNormal(X)(X x)
982 {
983     alias floatTraits!(X) F;
984    
985     static if(real.mant_dig==106) { // doubledouble
986     // doubledouble is normal if the least significant part is normal.
987         return isNormal((cast(double*)&x)[MANTISSA_LSB]);
988     } else {
989         ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
990         return (e != F.EXPMASK && e!=0);
991     }
992 }
993
994 debug(UnitTest) {
995 unittest
996 {
997     float f = 3;
998     double d = 500;
999     real e = 10e+48;
1000
1001     assert(isNormal(f));
1002     assert(isNormal(d));
1003     assert(isNormal(e));
1004     f=d=e=0;
1005     assert(!isNormal(f));
1006     assert(!isNormal(d));
1007     assert(!isNormal(e));
1008     assert(!isNormal(real.infinity));
1009     assert(isNormal(-real.max));
1010     assert(!isNormal(real.min/4));
1011    
1012 }
1013 }
1014
1015 /*********************************
1016  * Is the binary representation of x identical to y?
1017  *
1018  * Same as ==, except that positive and negative zero are not identical,
1019  * and two $(NAN)s are identical if they have the same 'payload'.
1020  */
1021
1022 bool isIdentical(real x, real y)
1023 {
1024     // We're doing a bitwise comparison so the endianness is irrelevant.
1025     long*   pxs = cast(long *)&x;
1026     long*   pys = cast(long *)&y;
1027   static if (real.mant_dig == 53){ //double
1028     return pxs[0] == pys[0];
1029   } else static if (real.mant_dig == 113 || real.mant_dig==106) {
1030       // quadruple or doubledouble
1031     return pxs[0] == pys[0] && pxs[1] == pys[1];
1032   } else { // real80
1033     ushort* pxe = cast(ushort *)&x;
1034     ushort* pye = cast(ushort *)&y;
1035     return pxe[4] == pye[4] && pxs[0] == pys[0];
1036   }
1037 }
1038
1039 /** ditto */
1040 bool isIdentical(ireal x, ireal y) {
1041     return isIdentical(x.im, y.im);
1042 }
1043
1044 /** ditto */
1045 bool isIdentical(creal x, creal y) {
1046     return isIdentical(x.re, y.re) && isIdentical(x.im, y.im);
1047 }
1048
1049 debug(UnitTest) {
1050 unittest {
1051     assert(isIdentical(0.0, 0.0));
1052     assert(!isIdentical(0.0, -0.0));
1053     assert(isIdentical(NaN(0xABC), NaN(0xABC)));
1054     assert(!isIdentical(NaN(0xABC), NaN(218)));
1055     assert(isIdentical(1.234e56, 1.234e56));
1056     assert(isNaN(NaN(0x12345)));
1057     assert(isIdentical(3.1 + NaN(0xDEF) * 1i, 3.1 + NaN(0xDEF)*1i));
1058     assert(!isIdentical(3.1+0.0i, 3.1-0i));
1059     assert(!isIdentical(0.0i, 2.5e58i));
1060 }
1061 }
1062
1063 /*********************************
1064  * Is number subnormal? (Also called "denormal".)
1065  * Subnormals have a 0 exponent and a 0 most significant significand bit,
1066  * but are non-zero.
1067  */
1068
1069 /* Need one for each format because subnormal floats might
1070  * be converted to normal reals.
1071  */
1072
1073 int isSubnormal(float f)
1074 {
1075     uint *p = cast(uint *)&f;
1076     return (*p & 0x7F80_0000) == 0 && *p & 0x007F_FFFF;
1077 }
1078
1079 debug(UnitTest) {
1080 unittest
1081 {
1082     float f = -float.min;
1083     assert(!isSubnormal(f));
1084     f/=4;
1085     assert(isSubnormal(f));
1086 }
1087 }
1088
1089 /// ditto
1090
1091 int isSubnormal(double d)
1092 {
1093     uint *p = cast(uint *)&d;
1094     return (p[MANTISSA_MSB] & 0x7FF0_0000) == 0 && (p[MANTISSA_LSB] || p[MANTISSA_MSB] & 0x000F_FFFF);
1095 }
1096
1097 debug(UnitTest) {
1098 unittest
1099 {
1100     double f;
1101
1102     for (f = 1; !isSubnormal(f); f /= 2)
1103     assert(f != 0);
1104 }
1105 }
1106
1107 /// ditto
1108
1109 int isSubnormal(real x)
1110 {
1111     alias floatTraits!(real) F;
1112     static if (real.mant_dig == 53) { // double
1113         return isSubnormal(cast(double)x);
1114     } else static if (real.mant_dig == 113) { // quadruple       
1115         ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
1116         long*   ps = cast(long *)&x;
1117         return (e == 0 && (((ps[MANTISSA_LSB]|(ps[MANTISSA_MSB]& 0x0000_FFFF_FFFF_FFFF))) !=0));
1118     } else static if (real.mant_dig==64) { // real80
1119         ushort* pe = cast(ushort *)&x;
1120         long*   ps = cast(long *)&x;
1121
1122         return (pe[F.EXPPOS_SHORT] & F.EXPMASK) == 0 && *ps > 0;
1123     } else { // double double
1124         return isSubnormal((cast(double*)&x)[MANTISSA_MSB]);
1125     }
1126 }
1127
1128 debug(UnitTest) {
1129 unittest
1130 {
1131     real f;
1132
1133     for (f = 1; !isSubnormal(f); f /= 2)
1134     assert(f != 0);
1135 }
1136 }
1137
1138 /*********************************
1139  * Return !=0 if x is $(PLUSMN)0.
1140  *
1141  * Does not affect any floating-point flags
1142  */
1143 int isZero(real x)
1144 {
1145     alias floatTraits!(real) F;
1146     static if (real.mant_dig == 53) { // double
1147         return ((*cast(ulong *)&x) & 0x7FFF_FFFF_FFFF_FFFF) == 0;
1148     } else static if (real.mant_dig == 113) { // quadruple   
1149         long*   ps = cast(long *)&x;
1150         return (ps[MANTISSA_LSB] | (ps[MANTISSA_MSB]& 0x7FFF_FFFF_FFFF_FFFF)) == 0;
1151     } else { // real80
1152         ushort* pe = cast(ushort *)&x;
1153         ulong*  ps = cast(ulong  *)&x;
1154         return (pe[F.EXPPOS_SHORT] & F.EXPMASK) == 0 && *ps == 0;
1155     }
1156 }
1157
1158 debug(UnitTest) {
1159 unittest
1160 {
1161     assert(isZero(0.0));
1162     assert(isZero(-0.0));
1163     assert(!isZero(2.5));
1164     assert(!isZero(real.min / 1000));
1165 }
1166 }
1167
1168 /*********************************
1169  * Return !=0 if e is $(PLUSMNINF);.
1170  */
1171
1172 int isInfinity(real x)
1173 {
1174     alias floatTraits!(real) F;
1175     static if (real.mant_dig == 53) { // double
1176         return ((*cast(ulong *)&x) & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FF8_0000_0000_0000;
1177     } else static if(real.mant_dig == 106) { //doubledouble
1178         return (((cast(ulong *)&x)[MANTISSA_MSB]) & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FF8_0000_0000_0000;   
1179     } else static if (real.mant_dig == 113) { // quadruple   
1180         long*   ps = cast(long *)&x;
1181         return (ps[MANTISSA_LSB] == 0)
1182          && (ps[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_0000_0000_0000;
1183     } else { // real80
1184         ushort e = cast(ushort)(F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]);
1185         ulong*  ps = cast(ulong *)&x;
1186
1187         return e == F.EXPMASK && *ps == 0x8000_0000_0000_0000;
1188    }
1189 }
1190
1191 debug(UnitTest) {
1192 unittest
1193 {
1194     assert(isInfinity(float.infinity));
1195     assert(!isInfinity(float.nan));
1196     assert(isInfinity(double.infinity));
1197     assert(isInfinity(-real.infinity));
1198
1199     assert(isInfinity(-1.0 / 0.0));
1200 }
1201 }
1202
1203 /**
1204  * Calculate the next largest floating point value after x.
1205  *
1206  * Return the least number greater than x that is representable as a real;
1207  * thus, it gives the next point on the IEEE number line.
1208  *
1209  *  $(TABLE_SV
1210  *    $(SVH x,            nextUp(x)   )
1211  *    $(SV  -$(INFIN),    -real.max   )
1212  *    $(SV  $(PLUSMN)0.0, real.min*real.epsilon )
1213  *    $(SV  real.max,     $(INFIN) )
1214  *    $(SV  $(INFIN),     $(INFIN) )
1215  *    $(SV  $(NAN),       $(NAN)   )
1216  * )
1217  *
1218  * Remarks:
1219  * This function is included in the IEEE 754-2008 standard.
1220  *
1221  * nextDoubleUp and nextFloatUp are the corresponding functions for
1222  * the IEEE double and IEEE float number lines.
1223  */
1224 real nextUp(real x)
1225 {
1226     alias floatTraits!(real) F;
1227     static if (real.mant_dig == 53) { // double
1228         return nextDoubleUp(x);
1229     } else static if(real.mant_dig==113) {  // quadruple
1230         ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
1231         if (e == F.EXPMASK) { // NaN or Infinity
1232              if (x == -real.infinity) return -real.max;
1233              return x; // +Inf and NaN are unchanged.
1234         }     
1235         ulong*   ps = cast(ulong *)&e;
1236         if (ps[MANTISSA_LSB] & 0x8000_0000_0000_0000)  { // Negative number
1237             if (ps[MANTISSA_LSB]==0 && ps[MANTISSA_MSB] == 0x8000_0000_0000_0000) { // it was negative zero
1238                 ps[MANTISSA_LSB] = 0x0000_0000_0000_0001; // change to smallest subnormal
1239                 ps[MANTISSA_MSB] = 0;
1240                 return x;
1241             }
1242             --*ps;
1243             if (ps[MANTISSA_LSB]==0) --ps[MANTISSA_MSB];
1244         } else { // Positive number
1245             ++ps[MANTISSA_LSB];
1246             if (ps[MANTISSA_LSB]==0) ++ps[MANTISSA_MSB];
1247         }
1248         return x;
1249          
1250     } else static if(real.mant_dig==64){ // real80
1251         // For 80-bit reals, the "implied bit" is a nuisance...
1252         ushort *pe = cast(ushort *)&x;
1253         ulong  *ps = cast(ulong  *)&x;
1254
1255         if ((pe[F.EXPPOS_SHORT] & F.EXPMASK) == F.EXPMASK) {
1256             // First, deal with NANs and infinity
1257             if (x == -real.infinity) return -real.max;
1258             return x; // +Inf and NaN are unchanged.
1259         }
1260         if (pe[F.EXPPOS_SHORT] & 0x8000)  { // Negative number -- need to decrease the significand
1261             --*ps;
1262             // Need to mask with 0x7FFF... so subnormals are treated correctly.
1263             if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_FFFF_FFFF_FFFF) {
1264                 if (pe[F.EXPPOS_SHORT] == 0x8000) { // it was negative zero
1265                     *ps = 1;
1266                     pe[F.EXPPOS_SHORT] = 0; // smallest subnormal.
1267                     return x;
1268                 }
1269                 --pe[F.EXPPOS_SHORT];
1270                 if (pe[F.EXPPOS_SHORT] == 0x8000) {
1271                     return x; // it's become a subnormal, implied bit stays low.
1272                 }
1273                 *ps = 0xFFFF_FFFF_FFFF_FFFF; // set the implied bit
1274                 return x;
1275             }
1276             return x;
1277         } else {
1278             // Positive number -- need to increase the significand.
1279             // Works automatically for positive zero.
1280             ++*ps;
1281             if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0) {
1282                 // change in exponent
1283                 ++pe[F.EXPPOS_SHORT];
1284                 *ps = 0x8000_0000_0000_0000; // set the high bit
1285             }
1286         }
1287         return x;
1288     } else { // doubledouble
1289         assert(0, "Not implemented");
1290     }
1291 }
1292
1293 /** ditto */
1294 double nextDoubleUp(double x)
1295 {
1296     ulong *ps = cast(ulong *)&x;
1297
1298     if ((*ps & 0x7FF0_0000_0000_0000) == 0x7FF0_0000_0000_0000) {
1299         // First, deal with NANs and infinity
1300         if (x == -x.infinity) return -x.max;
1301         return x; // +INF and NAN are unchanged.
1302     }
1303     if (*ps & 0x8000_0000_0000_0000)  { // Negative number
1304         if (*ps == 0x8000_0000_0000_0000) { // it was negative zero
1305             *ps = 0x0000_0000_0000_0001; // change to smallest subnormal
1306             return x;
1307         }
1308         --*ps;
1309     } else { // Positive number
1310         ++*ps;
1311     }
1312     return x;
1313 }
1314
1315 /** ditto */
1316 float nextFloatUp(float x)
1317 {
1318     uint *ps = cast(uint *)&x;
1319
1320     if ((*ps & 0x7F80_0000) == 0x7F80_0000) {
1321         // First, deal with NANs and infinity
1322         if (x == -x.infinity) return -x.max;
1323         return x; // +INF and NAN are unchanged.
1324     }
1325     if (*ps & 0x8000_0000)  { // Negative number
1326         if (*ps == 0x8000_0000) { // it was negative zero
1327             *ps = 0x0000_0001; // change to smallest subnormal
1328             return x;
1329         }
1330         --*ps;
1331     } else { // Positive number
1332         ++*ps;
1333     }
1334     return x;
1335 }
1336
1337 debug(UnitTest) {
1338 unittest {
1339  static if (real.mant_dig == 64) {
1340
1341   // Tests for 80-bit reals
1342
1343     assert(isIdentical(nextUp(NaN(0xABC)), NaN(0xABC)));
1344     // negative numbers
1345     assert( nextUp(-real.infinity) == -real.max );
1346     assert( nextUp(-1-real.epsilon) == -1.0 );
1347     assert( nextUp(-2) == -2.0 + real.epsilon);
1348     // denormals and zero
1349     assert( nextUp(-real.min) == -real.min*(1-real.epsilon) );
1350     assert( nextUp(-real.min*(1-real.epsilon) == -real.min*(1-2*real.epsilon)) );
1351     assert( isIdentical(-0.0L, nextUp(-real.min*real.epsilon)) );
1352     assert( nextUp(-0.0) == real.min*real.epsilon );
1353     assert( nextUp(0.0) == real.min*real.epsilon );
1354     assert( nextUp(real.min*(1-real.epsilon)) == real.min );
1355     assert( nextUp(real.min) == real.min*(1+real.epsilon) );
1356     // positive numbers
1357     assert( nextUp(1) == 1.0 + real.epsilon );
1358     assert( nextUp(2.0-real.epsilon) == 2.0 );
1359     assert( nextUp(real.max) == real.infinity );
1360     assert( nextUp(real.infinity)==real.infinity );
1361  }
1362
1363     assert(isIdentical(nextDoubleUp(NaN(0xABC)), NaN(0xABC)));
1364     // negative numbers
1365     assert( nextDoubleUp(-double.infinity) == -double.max );
1366     assert( nextDoubleUp(-1-double.epsilon) == -1.0 );
1367     assert( nextDoubleUp(-2) == -2.0 + double.epsilon);
1368     // denormals and zero
1369
1370     assert( nextDoubleUp(-double.min) == -double.min*(1-double.epsilon) );
1371     assert( nextDoubleUp(-double.min*(1-double.epsilon) == -double.min*(1-2*double.epsilon)) );
1372     assert( isIdentical(-0.0, nextDoubleUp(-double.min*double.epsilon)) );
1373     assert( nextDoubleUp(0.0) == double.min*double.epsilon );
1374     assert( nextDoubleUp(-0.0) == double.min*double.epsilon );
1375     assert( nextDoubleUp(double.min*(1-double.epsilon)) == double.min );
1376     assert( nextDoubleUp(double.min) == double.min*(1+double.epsilon) );
1377     // positive numbers
1378     assert( nextDoubleUp(1) == 1.0 + double.epsilon );
1379     assert( nextDoubleUp(2.0-double.epsilon) == 2.0 );
1380     assert( nextDoubleUp(double.max) == double.infinity );
1381
1382     assert(isIdentical(nextFloatUp(NaN(0xABC)), NaN(0xABC)));
1383     assert( nextFloatUp(-float.min) == -float.min*(1-float.epsilon) );
1384     assert( nextFloatUp(1.0) == 1.0+float.epsilon );
1385     assert( nextFloatUp(-0.0) == float.min*float.epsilon);
1386     assert( nextFloatUp(float.infinity)==float.infinity );
1387
1388     assert(nextDown(1.0+real.epsilon)==1.0);
1389     assert(nextDoubleDown(1.0+double.epsilon)==1.0);
1390     assert(nextFloatDown(1.0+float.epsilon)==1.0);
1391     assert(nextafter(1.0+real.epsilon, -real.infinity)==1.0);
1392 }
1393 }
1394
1395 package {
1396 /** Reduces the magnitude of x, so the bits in the lower half of its significand
1397  * are all zero. Returns the amount which needs to be added to x to restore its
1398  * initial value; this amount will also have zeros in all bits in the lower half
1399  * of its significand.
1400  */
1401 X splitSignificand(X)(ref X x)
1402 {
1403     if (fabs(x) !< X.infinity) return 0; // don't change NaN or infinity
1404     X y = x; // copy the original value
1405     static if (X.mant_dig == float.mant_dig) {
1406         uint *ps = cast(uint *)&x;
1407         (*ps) &= 0xFFFF_FC00;
1408     } else static if (X.mant_dig == 53) {
1409         ulong *ps = cast(ulong *)&x;
1410         (*ps) &= 0xFFFF_FFFF_FC00_0000L;
1411     } else static if (X.mant_dig == 64){ // 80-bit real
1412         // An x87 real80 has 63 bits, because the 'implied' bit is stored explicitly.
1413         // This is annoying, because it means the significand cannot be
1414         // precisely halved. Instead, we split it into 31+32 bits.
1415         ulong *ps = cast(ulong *)&x;
1416         (*ps) &= 0xFFFF_FFFF_0000_0000L;
1417     } else static if (X.mant_dig==113) { // quadruple
1418         ulong *ps = cast(ulong *)&x;
1419         ps[MANTISSA_LSB] &= 0xFF00_0000_0000_0000L;
1420     }
1421     //else static assert(0, "Unsupported size");
1422
1423     return y - x;
1424 }
1425
1426 unittest {
1427     double x = -0x1.234_567A_AAAA_AAp+250;
1428     double y = splitSignificand(x);
1429     assert(x == -0x1.234_5678p+250);
1430     assert(y == -0x0.000_000A_AAAA_A8p+248);
1431     assert(x + y == -0x1.234_567A_AAAA_AAp+250);
1432 }
1433 }
1434
1435 /**
1436  * Calculate the next smallest floating point value before x.
1437  *
1438  * Return the greatest number less than x that is representable as a real;
1439  * thus, it gives the previous point on the IEEE number line.
1440  *
1441  *  $(TABLE_SV
1442  *    $(SVH x,            nextDown(x)   )
1443  *    $(SV  $(INFIN),     real.max  )
1444  *    $(SV  $(PLUSMN)0.0, -real.min*real.epsilon )
1445  *    $(SV  -real.max,    -$(INFIN) )
1446  *    $(SV  -$(INFIN),    -$(INFIN) )
1447  *    $(SV  $(NAN),       $(NAN)    )
1448  * )
1449  *
1450  * Remarks:
1451  * This function is included in the IEEE 754-2008 standard.
1452  *
1453  * nextDoubleDown and nextFloatDown are the corresponding functions for
1454  * the IEEE double and IEEE float number lines.
1455  */
1456 real nextDown(real x)
1457 {
1458     return -nextUp(-x);
1459 }
1460
1461 /** ditto */
1462 double nextDoubleDown(double x)
1463 {
1464     return -nextDoubleUp(-x);
1465 }
1466
1467 /** ditto */
1468 float nextFloatDown(float x)
1469 {
1470     return -nextFloatUp(-x);
1471 }
1472
1473 debug(UnitTest) {
1474 unittest {
1475     assert( nextDown(1.0 + real.epsilon) == 1.0);
1476 }
1477 }
1478
1479 /**
1480  * Calculates the next representable value after x in the direction of y.
1481  *
1482  * If y > x, the result will be the next largest floating-point value;
1483  * if y < x, the result will be the next smallest value.
1484  * If x == y, the result is y.
1485  *
1486  * Remarks:
1487  * This function is not generally very useful; it's almost always better to use
1488  * the faster functions nextUp() or nextDown() instead.
1489  *
1490  * IEEE 754 requirements not implemented:
1491  * The FE_INEXACT and FE_OVERFLOW exceptions will be raised if x is finite and
1492  * the function result is infinite. The FE_INEXACT and FE_UNDERFLOW
1493  * exceptions will be raised if the function value is subnormal, and x is
1494  * not equal to y.
1495  */
1496 real nextafter(real x, real y)
1497 {
1498     if (x==y) return y;
1499     return (y>x) ? nextUp(x) : nextDown(x);
1500 }
1501
1502 /**************************************
1503  * To what precision is x equal to y?
1504  *
1505  * Returns: the number of significand bits which are equal in x and y.
1506  * eg, 0x1.F8p+60 and 0x1.F1p+60 are equal to 5 bits of precision.
1507  *
1508  *  $(TABLE_SV
1509  *    $(SVH3 x,      y,         feqrel(x, y)  )
1510  *    $(SV3  x,      x,         typeof(x).mant_dig )
1511  *    $(SV3  x,      $(GT)= 2*x, 0 )
1512  *    $(SV3  x,      $(LE)= x/2, 0 )
1513  *    $(SV3  $(NAN), any,       0 )
1514  *    $(SV3  any,    $(NAN),    0 )
1515  *  )
1516  *
1517  * Remarks:
1518  * This is a very fast operation, suitable for use in speed-critical code.
1519  */
1520 int feqrel(X)(X x, X y)
1521 {
1522     /* Public Domain. Author: Don Clugston, 18 Aug 2005.
1523      */
1524   static assert(is(X==real) || is(X==double) || is(X==float), "Only float, double, and real are supported by feqrel");
1525  
1526   static if (X.mant_dig == 106) { // doubledouble.
1527      int a = feqrel(cast(double*)(&x)[MANTISSA_MSB], cast(double*)(&y)[MANTISSA_MSB]);
1528      if (a != double.mant_dig) return a;
1529      return double.mant_dig + feqrel(cast(double*)(&x)[MANTISSA_LSB], cast(double*)(&y)[MANTISSA_LSB]);     
1530   } else static if (X.mant_dig==64 || X.mant_dig==113
1531                  || X.mant_dig==53 || X.mant_dig == 24) {
1532     if (x == y) return X.mant_dig; // ensure diff!=0, cope with INF.
1533
1534     X diff = fabs(x - y);
1535
1536     ushort *pa = cast(ushort *)(&x);
1537     ushort *pb = cast(ushort *)(&y);
1538     ushort *pd = cast(ushort *)(&diff);
1539
1540     alias floatTraits!(X) F;
1541
1542     // The difference in abs(exponent) between x or y and abs(x-y)
1543     // is equal to the number of significand bits of x which are
1544     // equal to y. If negative, x and y have different exponents.
1545     // If positive, x and y are equal to 'bitsdiff' bits.
1546     // AND with 0x7FFF to form the absolute value.
1547     // To avoid out-by-1 errors, we subtract 1 so it rounds down
1548     // if the exponents were different. This means 'bitsdiff' is
1549     // always 1 lower than we want, except that if bitsdiff==0,
1550     // they could have 0 or 1 bits in common.
1551
1552  static if (X.mant_dig==64 || X.mant_dig==113) { // real80 or quadruple
1553     int bitsdiff = ( ((pa[F.EXPPOS_SHORT] & F.EXPMASK)
1554                      + (pb[F.EXPPOS_SHORT]& F.EXPMASK)
1555                      - (0x8000-F.EXPMASK))>>1)
1556                 - pd[F.EXPPOS_SHORT];
1557  } else static if (X.mant_dig==53) { // double
1558     int bitsdiff = (( ((pa[F.EXPPOS_SHORT] & F.EXPMASK)
1559                      + (pb[F.EXPPOS_SHORT] & F.EXPMASK)
1560                      - (0x8000-F.EXPMASK))>>1)
1561                  - (pd[F.EXPPOS_SHORT] & F.EXPMASK))>>4;
1562  } else static if (X.mant_dig == 24) { // float
1563      int bitsdiff = (( ((pa[F.EXPPOS_SHORT] & F.EXPMASK)
1564                       + (pb[F.EXPPOS_SHORT] & F.EXPMASK)
1565                       - (0x8000-F.EXPMASK))>>1)
1566              - (pd[F.EXPPOS_SHORT] & F.EXPMASK))>>7;     
1567  }
1568     if (pd[F.EXPPOS_SHORT] == 0)
1569     {   // Difference is denormal
1570         // For denormals, we need to add the number of zeros that
1571         // lie at the start of diff's significand.
1572         // We do this by multiplying by 2^real.mant_dig
1573         diff *= F.RECIP_EPSILON;
1574         return bitsdiff + X.mant_dig - pd[F.EXPPOS_SHORT];
1575     }
1576
1577     if (bitsdiff > 0)
1578         return bitsdiff + 1; // add the 1 we subtracted before
1579        
1580     // Avoid out-by-1 errors when factor is almost 2.   
1581      static if (X.mant_dig==64 || X.mant_dig==113) { // real80 or quadruple   
1582         return (bitsdiff == 0) ? (pa[F.EXPPOS_SHORT] == pb[F.EXPPOS_SHORT]) : 0;
1583      } else static if (X.mant_dig == 53 || X.mant_dig == 24) { // double or float
1584         return (bitsdiff == 0 && !((pa[F.EXPPOS_SHORT] ^ pb[F.EXPPOS_SHORT])& F.EXPMASK)) ? 1 : 0;
1585      }
1586  } else {
1587     assert(0, "Unsupported");
1588  }
1589 }
1590
1591 debug(UnitTest) {
1592 unittest
1593 {
1594    // Exact equality
1595    assert(feqrel(real.max,real.max)==real.mant_dig);
1596    assert(feqrel(0.0L,0.0L)==real.mant_dig);
1597    assert(feqrel(7.1824L,7.1824L)==real.mant_dig);
1598    assert(feqrel(real.infinity,real.infinity)==real.mant_dig);
1599
1600    // a few bits away from exact equality
1601    real w=1;
1602    for (int i=1; i<real.mant_dig-1; ++i) {
1603       assert(feqrel(1+w*real.epsilon,1.0L)==real.mant_dig-i);
1604       assert(feqrel(1-w*real.epsilon,1.0L)==real.mant_dig-i);
1605       assert(feqrel(1.0L,1+(w-1)*real.epsilon)==real.mant_dig-i+1);
1606       w*=2;
1607    }
1608    assert(feqrel(1.5+real.epsilon,1.5L)==real.mant_dig-1);
1609    assert(feqrel(1.5-real.epsilon,1.5L)==real.mant_dig-1);
1610    assert(feqrel(1.5-real.epsilon,1.5+real.epsilon)==real.mant_dig-2);
1611    
1612    assert(feqrel(real.min/8,real.min/17)==3);;
1613    
1614    // Numbers that are close
1615    assert(feqrel(0x1.Bp+84, 0x1.B8p+84)==5);
1616    assert(feqrel(0x1.8p+10, 0x1.Cp+10)==2);
1617    assert(feqrel(1.5*(1-real.epsilon), 1.0L)==2);
1618    assert(feqrel(1.5, 1.0)==1);
1619    assert(feqrel(2*(1-real.epsilon), 1.0L)==1);
1620
1621    // Factors of 2
1622    assert(feqrel(real.max,real.infinity)==0);
1623    assert(feqrel(2*(1-real.epsilon), 1.0L)==1);
1624    assert(feqrel(1.0, 2.0)==0);
1625    assert(feqrel(4.0, 1.0)==0);
1626
1627    // Extreme inequality
1628    assert(feqrel(real.nan,real.nan)==0);
1629    assert(feqrel(0.0L,-real.nan)==0);
1630    assert(feqrel(real.nan,real.infinity)==0);
1631    assert(feqrel(real.infinity,-real.infinity)==0);
1632    assert(feqrel(-real.max,real.infinity)==0);
1633    assert(feqrel(real.max,-real.max)==0);
1634    
1635    // floats
1636    assert(feqrel(2.1f, 2.1f)==float.mant_dig);
1637    assert(feqrel(1.5f, 1.0f)==1);
1638 }
1639 }
1640
1641 /*********************************
1642  * Return 1 if sign bit of e is set, 0 if not.
1643  */
1644
1645 int signbit(real x)
1646 {
1647     return ((cast(ubyte *)&x)[floatTraits!(real).SIGNPOS_BYTE] & 0x80) != 0;
1648 }
1649
1650 debug(UnitTest) {
1651 unittest
1652 {
1653     assert(!signbit(float.nan));
1654     assert(signbit(-float.nan));
1655     assert(!signbit(168.1234));
1656     assert(signbit(-168.1234));
1657     assert(!signbit(0.0));
1658     assert(signbit(-0.0));
1659 }
1660 }
1661
1662
1663 /*********************************
1664  * Return a value composed of to with from's sign bit.
1665  */
1666
1667 real copysign(real to, real from)
1668 {
1669     ubyte* pto   = cast(ubyte *)&to;
1670     ubyte* pfrom = cast(ubyte *)&from;
1671    
1672     alias floatTraits!(real) F;
1673     pto[F.SIGNPOS_BYTE] &= 0x7F;
1674     pto[F.SIGNPOS_BYTE] |= pfrom[F.SIGNPOS_BYTE] & 0x80;
1675     return to;
1676 }
1677
1678 debug(UnitTest) {
1679 unittest
1680 {
1681     real e;
1682
1683     e = copysign(21, 23.8);
1684     assert(e == 21);
1685
1686     e = copysign(-21, 23.8);
1687     assert(e == 21);
1688
1689     e = copysign(21, -23.8);
1690     assert(e == -21);
1691
1692     e = copysign(-21, -23.8);
1693     assert(e == -21);
1694
1695     e = copysign(real.nan, -23.8);
1696     assert(isNaN(e) && signbit(e));
1697 }
1698 }
1699
1700 /** Return the value that lies halfway between x and y on the IEEE number line.
1701  *
1702  * Formally, the result is the arithmetic mean of the binary significands of x
1703  * and y, multiplied by the geometric mean of the binary exponents of x and y.
1704  * x and y must have the same sign, and must not be NaN.
1705  * Note: this function is useful for ensuring O(log n) behaviour in algorithms
1706  * involving a 'binary chop'.
1707  *
1708  * Special cases:
1709  * If x and y are within a factor of 2, (ie, feqrel(x, y) > 0), the return value
1710  * is the arithmetic mean (x + y) / 2.
1711  * If x and y are even powers of 2, the return value is the geometric mean,
1712  *   ieeeMean(x, y) = sqrt(x * y).
1713  *
1714  */
1715 T ieeeMean(T)(T x, T y)
1716 in {
1717     // both x and y must have the same sign, and must not be NaN.
1718     assert(signbit(x) == signbit(y));
1719     assert(x<>=0 && y<>=0);
1720 }
1721 body {
1722     // Runtime behaviour for contract violation:
1723     // If signs are opposite, or one is a NaN, return 0.
1724     if (!((x>=0 && y>=0) || (x<=0 && y<=0))) return 0.0;
1725
1726     // The implementation is simple: cast x and y to integers,
1727     // average them (avoiding overflow), and cast the result back to a floating-point number.
1728
1729     alias floatTraits!(real) F;
1730     T u;
1731     static if (T.mant_dig==64) { // real80
1732         // There's slight additional complexity because they are actually
1733         // 79-bit reals...
1734         ushort *ue = cast(ushort *)&u;
1735         ulong *ul = cast(ulong *)&u;
1736         ushort *xe = cast(ushort *)&x;
1737         ulong *xl = cast(ulong *)&x;
1738         ushort *ye = cast(ushort *)&y;
1739         ulong *yl = cast(ulong *)&y;
1740         // Ignore the useless implicit bit. (Bonus: this prevents overflows)
1741         ulong m = ((*xl) & 0x7FFF_FFFF_FFFF_FFFFL) + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL);
1742
1743         ushort e = cast(ushort)((xe[F.EXPPOS_SHORT] & 0x7FFF) + (ye[F.EXPPOS_SHORT] & 0x7FFF));
1744         if (m & 0x8000_0000_0000_0000L) {
1745             ++e;
1746             m &= 0x7FFF_FFFF_FFFF_FFFFL;
1747         }
1748         // Now do a multi-byte right shift
1749         uint c = e & 1; // carry
1750         e >>= 1;
1751         m >>>= 1;
1752         if (c) m |= 0x4000_0000_0000_0000L; // shift carry into significand
1753         if (e) *ul = m | 0x8000_0000_0000_0000L; // set implicit bit...
1754         else *ul = m; // ... unless exponent is 0 (denormal or zero).
1755         ue[4]=  e | (xe[F.EXPPOS_SHORT]& F.SIGNMASK); // restore sign bit
1756     } else static if(T.mant_dig == 113) { //quadruple
1757         // This would be trivial if 'ucent' were implemented...
1758         ulong *ul = cast(ulong *)&u;
1759         ulong *xl = cast(ulong *)&x;
1760         ulong *yl = cast(ulong *)&y;
1761         // Multi-byte add, then multi-byte right shift.       
1762         ulong mh = ((xl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL)
1763                   + (yl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL));
1764         // Discard the lowest bit (to avoid overflow)
1765         ulong ml = (xl[MANTISSA_LSB]>>>1) + (yl[MANTISSA_LSB]>>>1);
1766         // add the lowest bit back in, if necessary.
1767         if (xl[MANTISSA_LSB] & yl[MANTISSA_LSB] & 1) {
1768             ++ml;
1769             if (ml==0) ++mh;
1770         }
1771         mh >>>=1;
1772         ul[MANTISSA_MSB] = mh | (xl[MANTISSA_MSB] & 0x8000_0000_0000_0000);
1773         ul[MANTISSA_LSB] = ml;
1774     } else static if (T.mant_dig == double.mant_dig) {
1775         ulong *ul = cast(ulong *)&u;
1776         ulong *xl = cast(ulong *)&x;
1777         ulong *yl = cast(ulong *)&y;
1778         ulong m = (((*xl) & 0x7FFF_FFFF_FFFF_FFFFL) + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL)) >>> 1;
1779         m |= ((*xl) & 0x8000_0000_0000_0000L);
1780         *ul = m;
1781     } else static if (T.mant_dig == float.mant_dig) {
1782         uint *ul = cast(uint *)&u;
1783         uint *xl = cast(uint *)&x;
1784         uint *yl = cast(uint *)&y;
1785         uint m = (((*xl) & 0x7FFF_FFFF) + ((*yl) & 0x7FFF_FFFF)) >>> 1;
1786         m |= ((*xl) & 0x8000_0000);
1787         *ul = m;
1788     } else {
1789         assert(0, "Not implemented");
1790     }
1791     return u;
1792 }
1793
1794 debug(UnitTest) {
1795 unittest {
1796     assert(ieeeMean(-0.0,-1e-20)<0);
1797     assert(ieeeMean(0.0,1e-20)>0);
1798
1799     assert(ieeeMean(1.0L,4.0L)==2L);
1800     assert(ieeeMean(2.0*1.013,8.0*1.013)==4*1.013);
1801     assert(ieeeMean(-1.0L,-4.0L)==-2L);
1802     assert(ieeeMean(-1.0,-4.0)==-2);
1803     assert(ieeeMean(-1.0f,-4.0f)==-2f);
1804     assert(ieeeMean(-1.0,-2.0)==-1.5);
1805     assert(ieeeMean(-1*(1+8*real.epsilon),-2*(1+8*real.epsilon))==-1.5*(1+5*real.epsilon));
1806     assert(ieeeMean(0x1p60,0x1p-10)==0x1p25);
1807     static if (real.mant_dig==64) { // x87, 80-bit reals
1808       assert(ieeeMean(1.0L,real.infinity)==0x1p8192L);
1809       assert(ieeeMean(0.0L,real.infinity)==1.5);
1810     }
1811     assert(ieeeMean(0.5*real.min*(1-4*real.epsilon),0.5*real.min)==0.5*real.min*(1-2*real.epsilon));
1812 }
1813 }
1814
1815 // Functions for NaN payloads
1816 /*
1817  * A 'payload' can be stored in the significand of a $(NAN). One bit is required
1818  * to distinguish between a quiet and a signalling $(NAN). This leaves 22 bits
1819  * of payload for a float; 51 bits for a double; 62 bits for an 80-bit real;
1820  * and 111 bits for a 128-bit quad.
1821 */
1822 /**
1823  * Create a $(NAN), storing an integer inside the payload.
1824  *
1825  * For 80-bit or 128-bit reals, the largest possible payload is 0x3FFF_FFFF_FFFF_FFFF.
1826  * For doubles, it is 0x3_FFFF_FFFF_FFFF.
1827  * For floats, it is 0x3F_FFFF.
1828  */
1829 real NaN(ulong payload)
1830 {
1831     static if (real.mant_dig == 64) { //real80
1832       ulong v = 3; // implied bit = 1, quiet bit = 1
1833     } else {
1834       ulong v = 2; // no implied bit. quiet bit = 1
1835     }
1836
1837     ulong a = payload;
1838
1839     // 22 Float bits
1840     ulong w = a & 0x3F_FFFF;
1841     a -= w;
1842
1843     v <<=22;
1844     v |= w;
1845     a >>=22;
1846
1847     // 29 Double bits
1848     v <<=29;
1849     w = a & 0xFFF_FFFF;
1850     v |= w;
1851     a -= w;
1852     a >>=29;
1853
1854     static if (real.mant_dig == 53) { // double
1855         v |=0x7FF0_0000_0000_0000;
1856         real x;
1857         * cast(ulong *)(&x) = v;
1858         return x;
1859     } else {
1860         v <<=11;
1861         a &= 0x7FF;
1862         v |= a;
1863         real x = real.nan;
1864         // Extended real bits
1865         static if (real.mant_dig==113) { //quadruple
1866           v<<=1; // there's no implicit bit
1867           version(LittleEndian) {
1868             *cast(ulong*)(6+cast(ubyte*)(&x)) = v;
1869           } else {
1870             *cast(ulong*)(2+cast(ubyte*)(&x)) = v;
1871           }       
1872         } else { // real80
1873             * cast(ulong *)(&x) = v;
1874         }
1875         return x;
1876     }
1877 }
1878
1879 /**
1880  * Extract an integral payload from a $(NAN).
1881  *
1882  * Returns:
1883  * the integer payload as a ulong.
1884  *
1885  * For 80-bit or 128-bit reals, the largest possible payload is 0x3FFF_FFFF_FFFF_FFFF.
1886  * For doubles, it is 0x3_FFFF_FFFF_FFFF.
1887  * For floats, it is 0x3F_FFFF.
1888  */
1889 ulong getNaNPayload(real x)
1890 {
1891     assert(isNaN(x));
1892     static if (real.mant_dig == 53) {
1893         ulong m = *cast(ulong *)(&x);
1894         // Make it look like an 80-bit significand.
1895         // Skip exponent, and quiet bit
1896         m &= 0x0007_FFFF_FFFF_FFFF;
1897         m <<= 10;
1898     } else static if (real.mant_dig==113) { // quadruple
1899         version(LittleEndian) {
1900             ulong m = *cast(ulong*)(6+cast(ubyte*)(&x));
1901         } else {
1902             ulong m = *cast(ulong*)(2+cast(ubyte*)(&x));
1903         }
1904         m>>=1; // there's no implicit bit
1905     } else {
1906         ulong m = *cast(ulong *)(&x);
1907     }
1908     // ignore implicit bit and quiet bit
1909     ulong f = m & 0x3FFF_FF00_0000_0000L;
1910     ulong w = f >>> 40;
1911     w |= (m & 0x00FF_FFFF_F800L) << (22 - 11);
1912     w |= (m & 0x7FF) << 51;
1913     return w;
1914 }
1915
1916 debug(UnitTest) {
1917 unittest {
1918   real nan4 = NaN(0x789_ABCD_EF12_3456);
1919   static if (real.mant_dig == 64 || real.mant_dig==113) {
1920       assert (getNaNPayload(nan4) == 0x789_ABCD_EF12_3456);
1921   } else {
1922       assert (getNaNPayload(nan4) == 0x1_ABCD_EF12_3456);
1923   }
1924   double nan5 = nan4;
1925   assert (getNaNPayload(nan5) == 0x1_ABCD_EF12_3456);
1926   float nan6 = nan4;
1927   assert (getNaNPayload(nan6) == 0x12_3456);
1928   nan4 = NaN(0xFABCD);
1929   assert (getNaNPayload(nan4) == 0xFABCD);
1930   nan6 = nan4;
1931   assert (getNaNPayload(nan6) == 0xFABCD);
1932   nan5 = NaN(0x100_0000_0000_3456);
1933   assert(getNaNPayload(nan5) == 0x0000_0000_3456);
1934 }
1935 }
Note: See TracBrowser for help on using the browser.