decimal.floats source code

1 module decimal.floats;
2 
3 
4 /* ****************************************************************************************************************** */
5 /* FLOAT UTILITY FUNCTIONS                                                                                            */
6 /* ****************************************************************************************************************** */
7 
8 private import decimal.integrals: uint128, clz, divrem;
9 
10 package:
11 
12 union FU
13 {
14     float f;
15     uint u;
16 }
17 
18 @safe pure nothrow @nogc
19 float fpack(const bool sign, int exp, uint mantissa)
20 {
21     if (mantissa == 0)
22         return sign ? -0.0f : +0.0f;
23     auto shift = clz(mantissa) - 8;
24     if (shift < 0)
25     {
26         if (exp > int.max + shift)
27             return sign ? -float.infinity : +float.infinity;
28         else       
29             mantissa >>= -shift;
30     }
31     else
32     {
33         if (exp < int.min + shift)
34             return sign ? -0.0f : +0.0f;
35         mantissa <<= shift;
36     }
37     exp -= shift;
38 
39     if (exp > int.max - 150)
40         return sign ? -float.infinity : +float.infinity;
41     exp += 150;
42 
43     if (exp >= 0xFF)
44         return sign ? -float.infinity : +float.infinity;
45 
46     if (exp <= 0)
47     {
48         --exp;
49         if (exp < -23)
50             return sign ? -0.0f : +0.0f;
51         mantissa >>= -exp;
52         exp = 0;
53     }
54 
55     FU fu;
56     fu.u = (mantissa & 0x7FFFFF) | (exp << 23);
57     if (sign)
58         fu.u |= 0x80000000U;
59     return fu.f;
60 }
61 
62 @safe pure nothrow @nogc
63 bool funpack(const float f, out int exp, out uint mantissa, out bool inf, out bool nan)
64 {
65     FU fu;
66     fu.f = f;
67 
68     exp = (fu.u >> 23) & 0xFF;
69     mantissa = fu.u & 0x7FFFFFU;
70     if (exp == 0)
71     {
72         inf = false; nan = false;
73         if (mantissa)
74             exp -= 149;
75     }
76     else if (exp == 0xFF)
77     {
78         inf = mantissa == 0;
79         nan = !inf;
80     }
81     else
82     {
83         inf = false; nan = false;
84         mantissa |= 0x00800000;
85         exp -= 150;
86     }
87 
88     return (fu.u & 0x80000000U) != 0;
89 }
90 
91 union DU
92 {
93     double d;
94     ulong u;
95 }
96 
97 @safe pure nothrow @nogc
98 double dpack(const bool sign, int exp, ulong mantissa)
99 {
100     if (mantissa == 0)
101         return sign ? -0.0 : +0.0;
102 
103     auto shift = clz(mantissa) - 11;
104     if (shift < 0)
105     {
106         if (exp > int.max + shift)
107             return sign ? -double.infinity : +double.infinity;
108         else       
109             mantissa >>= -shift;
110     }
111     else
112     {
113         if (exp < int.min + shift)
114             return sign ? -0.0 : +0.0;
115         mantissa <<= shift;
116     }
117     exp -= shift;
118 
119     if (exp > int.max - 1075)
120         return sign ? -double.infinity : +double.infinity;
121     exp += 1075;
122 
123     if (exp >= 0x7FF)
124         return sign ? -double.infinity : +double.infinity;
125 
126     if (exp <= 0)
127     {
128         --exp;
129         if (exp < -52)
130             return sign ? -0.0 : +0.0;
131         mantissa >>= -exp;
132         exp = 0;
133     }
134 
135     DU du;
136     du.u = (mantissa & 0x000FFFFFFFFFFFFFUL) | (cast(ulong)exp << 52);
137     if (sign)
138         du.u |= 0x8000000000000000UL;
139     return du.d;
140 }
141 
142 @safe pure nothrow @nogc
143 bool dunpack(const double d, out int exp, out ulong mantissa, out bool inf, out bool nan)
144 {
145     DU du;
146     du.d = d;
147 
148     exp = (du.u >> 52) & 0x7FF;
149     mantissa = du.u & 0xFFFFFFFFFFFFF;
150 
151     if (exp == 0)
152     {
153         inf = false; nan = false;
154         if (mantissa)
155             exp -= 1074;
156     }
157     else if (exp == 0x7FF)
158     {
159         inf = mantissa == 0;
160         nan = !inf;
161     }
162     else
163     {
164         inf = false; nan = false;
165         mantissa |= 0x10000000000000;
166         exp -= 1075;
167     }
168 
169     return (du.u & 0x8000000000000000UL) != 0;
170 }
171 
172 union RU
173 {
174     real r;
175     struct
176     {   //align(1):
177         version(LittleEndian)
178         {
179             ulong m;
180             ushort e;    
181         }
182         else
183         {
184             ushort e;
185             ulong m;
186         }
187     }
188 }
189 
190 @safe pure nothrow @nogc
191 real rpack(const bool sign, int exp, ulong mantissa)
192 {
193     if (mantissa == 0)
194         return sign ? -0.0L : +0.0L;
195 
196     auto shift = clz(mantissa);
197     if (exp < int.min + shift)
198         return sign ? -0.0L : +0.0L;
199     mantissa <<= shift;
200     exp -= shift;
201 
202     if (exp > int.max - 16447) //16383 + 64
203         return sign ? -real.infinity : +real.infinity;
204     exp += 16447;
205 
206     if (exp >= 0x7FFF)
207         return sign ? -real.infinity : +real.infinity;
208 
209     if (exp <= 0)
210     {
211         --exp;
212         if (exp < -64)
213             return sign ? -0.0L : +0.0L;
214         mantissa >>= -exp;
215         exp = 0;
216     }
217 
218     RU ru;
219     ru.m = mantissa;
220     ru.e = cast(ushort)exp;
221     if (sign)
222         ru.e |= 0x8000U;
223     return ru.r;
224 }
225 
226 @safe pure nothrow @nogc
227 bool runpack(const real r, out int exp, out ulong mantissa, out bool inf, out bool nan)
228 {
229     RU ru;
230     ru.r = r;
231 
232     exp = ru.e & 0x7FFF;
233     mantissa = ru.m;
234 
235     if (exp == 0)
236     {
237         inf = false; nan = false;
238         if (mantissa)
239             exp -= 16445;
240     }
241     else if (exp == 0x7FFF)
242     {
243         inf = (mantissa & 0x7FFF_FFFF_FFFF_FFFF) == 0;
244         nan = !inf;
245     }
246     else
247     {
248         inf = false; nan = false;
249         exp -= 16446;
250     }
251 
252     return (ru.e & 0x8000) != 0;
253 }
254 
255 
256 
257 
258 void floatExtract(float f, out uint coefficient, out int exponent)
259 {
260     // x * 2^n = y * 10^m -> 2^n = n * log2(10)
261 }
262 
263 @nogc @safe pure nothrow
264 bool exp2to10(ref uint coefficient, ref int exponent)
265 {
266     enum maxMultiplicable = uint.max / 5U;
267 
268     bool inexact;
269 
270     auto e5 = -exponent;
271 
272     
273 
274     while (e5 > 0)
275     {
276         if (!(coefficient & 1U))
277         {
278             ++exponent;
279             --e5;
280             coefficient >>>= 1;
281         }
282         else
283         {
284             --e5;
285             if (coefficient <= maxMultiplicable)
286                 coefficient *= 5U;
287             else
288             {
289                 inexact = true;
290                 ++exponent;
291                 coefficient >>>= 1;
292             }
293         }
294     }
295 
296     while (e5 < 0)
297     {
298         if ((coefficient & 0x80000000U) != 0x80000000U)
299         {
300             --exponent;
301             ++e5;
302             coefficient <<= 1;
303         }
304         else
305         {
306             ++e5;
307             if (coefficient % 5U != 0)
308                 inexact = true;
309             coefficient /= 5U;
310         }
311     }
312 
313     return inexact;
314 }
315 
316 public
317 @nogc @safe pure nothrow
318 bool exp2to10(ref ulong coefficient, ref int exponent)
319 {
320     enum maxMultiplicable = ulong.max / 5UL;
321 
322     bool inexact;
323 
324     auto e5 = -exponent;
325 
326     while (e5 > 0)
327     {
328         if (!(coefficient & 1UL))
329         {
330             ++exponent;
331             --e5;
332             coefficient >>>= 1;
333         }
334         else
335         {
336             --e5;
337             if (coefficient <= maxMultiplicable)
338                 coefficient *= 5UL;
339             else
340             {
341                 inexact = true;
342                 ++exponent;
343                 coefficient >>>= 1;
344             }
345         }
346     }
347 
348     while (e5 < 0)
349     {
350         if ((coefficient & 0x8000000000000000UL) != 0x8000000000000000UL)
351         {
352             --exponent;
353             ++e5;
354             coefficient <<= 1;
355         }
356         else
357         {
358             ++e5;
359             if (coefficient % 5UL != 0)
360                 inexact = true;
361             coefficient /= 5UL;
362         }
363     }
364 
365     return inexact;
366 }
367 
368 @nogc @safe pure nothrow
369 bool exp2to10(ref uint128 coefficient, ref int exponent)
370 {
371     enum maxMultiplicable = uint128.max / 5UL;
372 
373     bool inexact;
374 
375     auto e5 = -exponent;
376 
377     while (e5 > 0)
378     {
379         if (!(coefficient & 1UL))
380         {
381             ++exponent;
382             --e5;
383             coefficient >>>= 1;
384         }
385         else
386         {
387             --e5;
388             if (coefficient <= maxMultiplicable)
389                 coefficient *= 5U;
390             else
391             {
392                 inexact = true;
393                 ++exponent;
394                 coefficient >>>= 1;
395             }
396         }
397     }
398 
399     while (e5 < 0)
400     {
401         if ((coefficient.hi & 0x8000000000000000UL) != 0x8000000000000000UL)
402         {
403             --exponent;
404             ++e5;
405             coefficient <<= 1;
406         }
407         else
408         {
409             ++e5;
410             if (divrem(coefficient, 5U))
411                 inexact = true;
412         }
413     }
414 
415     return inexact;
416 }
417 
418