File Coverage

blib/lib/Data/IEEE754/Tools.pm
Criterion Covered Total %
statement 211 211 100.0
branch 128 130 98.4
condition 55 60 91.6
subroutine 57 57 100.0
pod 24 46 52.1
total 475 504 94.2


line stmt bran cond sub pod time code
1             package Data::IEEE754::Tools;
2 11     11   159845 use 5.006;
  11         32  
3 11     11   38 use warnings;
  11         15  
  11         283  
4 11     11   35 use strict;
  11         21  
  11         193  
5 11     11   32 use Carp;
  11         15  
  11         841  
6 11     11   39 use Exporter 'import'; # just use the import() function, without the rest of the overhead of ISA
  11         15  
  11         347  
7            
8 11     11   4161 use version 0.77; our $VERSION = version->declare('0.013_006');
  11         14605  
  11         54  
9            
10             =pod
11            
12             =head1 NAME
13            
14             Data::IEEE754::Tools - Various tools for understanding and manipulating the underlying IEEE-754 representation of floating point values
15            
16             =head1 SYNOPSIS
17            
18             use Data::IEEE754::Tools qw/:floatingpoint :ulp/;
19            
20             # return -12.875 as decimal and hexadecimal floating point numbers
21             to_dec_floatingpoint(-12.875); # -0d1.6093750000000000p+0003
22             to_hex_floatingpoint(-12.875); # -0x1.9c00000000000p+0003
23            
24             # shows the smallest value you can add or subtract to 16.16 (ulp = "Unit in the Last Place")
25             print ulp( 16.16 ); # 3.5527136788005e-015
26            
27             # toggles the ulp: returns a float that has the ULP of 16.16 toggled
28             # (if it was a 1, it will be 0, and vice versa);
29             # running it twice should give the original value
30             print $t16 = toggle_ulp( 16.16 ); # 16.159999999999997
31             print $v16 = toggle_ulp( $t16 ); # 16.160000000000000
32            
33             =head1 DESCRIPTION
34            
35             These tools give access to the underlying IEEE 754 floating-point 64bit representation
36             used by many instances of Perl (see L). They include functions for converting
37             from the 64bit internal representation to a string that shows those bits (either as
38             hexadecimal or binary) and back, functions for converting that encoded value
39             into a more human-readable format to give insight into the meaning of the encoded
40             values, and functions to manipulate the smallest possible change for a given
41             floating-point value (which is the L or
42             "Unit in the Last Place").
43            
44             =head2 IEEE 754 Encoding
45            
46             The L standard describes
47             various floating-point encodings. The double format (`binary64') is a 64-bit base-2
48             encoding, and correpsonds to the usual Perl floating value (NV). The format includes the
49             sign (s), the power of 2 (q), and a significand (aka, mantissa; the coefficient, c):
50             C. The C<(-1)**s> term evaluates to the sign of the
51             number, where s=0 means the sign is +1 and s=1 means the sign is -1.
52            
53             For most numbers, the coefficient is an implied 1 plus an encoded fraction,
54             which is itself encoded as a 52-bit integer divided by an implied 2**52. The range of
55             valid exponents is from -1022 to +1023, which are encoded as an 11bit integer from 1
56             to 2046 (where C). With an 11bit integer,
57             there are two exponent values (C<0b000_0000_0000 = 0 - 1023 = -1023> and
58             C<0b111_1111_1111 = 2047 - 1023 = +1024>), which are used to indicate conditions outside
59             the normal range: The first special encoded-exponent, C<0b000_0000_0000>, indicates that
60             the coefficient is 0 plus the encoded fraction, at an exponent of -1022; thus, the
61             floating-point zero is encoded using an encoded-exponent of 0 and an encoded-fraction of 0
62             (C<[0 + 0/(2**52)] * [2**-1022] = 0*(2**-1022) = 0>); other numbers
63             smaller than can normally be encoded (so-called "denormals" or "subnormals"), lying
64             between 0 and 1 (non-inclusive) are encoded with the same exponent, but have a non-zero
65             encoded-fraction. The second special encoded-exponent, C<0b111_1111_1111>, indicates a
66             number that is infinite (too big to represent), or something that is not a number (NAN);
67             infinities are indicated by that special exponent and an encoded-fraction of 0; NAN
68             is indicated by that special exponent and a non-zero encoded-fraction.
69            
70             =head2 Justification for the existence of B
71            
72             L, or the equivalent L recipe L>, do a
73             good job of converting a perl floating value (NV) into the big-endian bytes
74             that encode that value, but they don't help you interpret the value.
75            
76             L has a similar suite of tools to B, but
77             uses numerical methods rather than accessing the underlying bits. It L
78             shown|http://perlmonks.org/?node_id=1167146> that its interpretation function can take
79             an order of magnitude longer than a routine that manipulates the underlying bits
80             to gather the information.
81            
82             This B module combines the two sets of functions, giving
83             access to the raw IEEE 754 encoding, or a stringification of the encoding which
84             interprets the encoding as a sign and a coefficient and a power of 2, or access to
85             the ULP and ULP-manipulating features, all using direct bit manipulation when
86             appropriate.
87            
88             =head2 Compatibility
89            
90             B works with 64bit floating-point representations.
91            
92             If you have a Perl setup which uses a larger representation (for example,
93             C; print $Config{nvsize}; # 16 =E 128bit>), values reported by
94             this module will be reduced in precision to fit the 64bit representation.
95            
96             If you have a Perl setup which uses a smaller representation (for example,
97             C; print $Config{nvsize}; # 4 =E 32bit>), the installation
98             will likely fail, because the unit tests were not set up for lower precision
99             inputs. However, forcing the installation I still allow coercion
100             from the smaller Perl NV into a true IEEE 754 double (64bit) floating-point,
101             but there is no guarantee it will work.
102            
103             =head1 INTERFACE NOT YET STABLE
104            
105             Please note: the interface to this module is not yet stable. There may be changes
106             to function naming conventions (under_scores vs. camelCase, argument order, etc).
107             Once B hits v1.000, the interface should be stable for all
108             sub-versions of v1: existing functions should keep the same calling conventions,
109             though new functions may be added; significant changes to the interface will cause
110             a transition to v2.
111            
112             =over
113            
114             =item v0.013_003: C renamed to C
115            
116             =item v0.013_003: C renamed to C
117            
118             =item v0.013_003: C renamed to C
119            
120             =back
121            
122             =head1 EXPORTABLE FUNCTIONS AND VARIABLES
123            
124             =cut
125            
126             my @EXPORT_RAW754 = qw(
127             hexstr754_from_double binstr754_from_double
128             hexstr754_to_double binstr754_to_double
129             );
130             my @EXPORT_FLOATING = qw(to_hex_floatingpoint to_dec_floatingpoint);
131             my @EXPORT_ULP = qw(ulp toggle_ulp nextUp nextDown nextAfter);
132             my @EXPORT_CONST = qw(
133             POS_ZERO
134             POS_DENORM_SMALLEST POS_DENORM_BIGGEST
135             POS_NORM_SMALLEST POS_NORM_BIGGEST
136             POS_INF
137             POS_SNAN_FIRST POS_SNAN_LAST
138             POS_IND POS_QNAN_FIRST POS_QNAN_LAST
139             NEG_ZERO
140             NEG_DENORM_SMALLEST NEG_DENORM_BIGGEST
141             NEG_NORM_SMALLEST NEG_NORM_BIGGEST
142             NEG_INF
143             NEG_SNAN_FIRST NEG_SNAN_LAST
144             NEG_IND NEG_QNAN_FIRST NEG_QNAN_LAST
145             );
146             my @EXPORT_INFO = qw(isSignMinus isNormal isFinite isZero isSubnormal
147             isInfinite isNaN isSignaling isSignalingConvertedToQuiet isCanonical
148             class radix totalOrder totalOrderMag);
149             my @EXPORT_SIGNBIT = qw(negate absolute copySign isSignMinus);
150            
151             our @EXPORT_OK = (@EXPORT_FLOATING, @EXPORT_RAW754, @EXPORT_ULP, @EXPORT_CONST, @EXPORT_INFO, @EXPORT_SIGNBIT);
152             our %EXPORT_TAGS = (
153             floating => [@EXPORT_FLOATING],
154             floatingpoint => [@EXPORT_FLOATING],
155             raw754 => [@EXPORT_RAW754],
156             ulp => [@EXPORT_ULP],
157             constants => [@EXPORT_CONST],
158             info => [@EXPORT_INFO],
159             signbit => [@EXPORT_SIGNBIT],
160             all => [@EXPORT_OK],
161             );
162            
163             =head2 :raw754
164            
165             These are the functions to do raw conversion from a floating-point value to a hexadecimal or binary
166             string of the underlying IEEE754 encoded value, and back.
167            
168             =head3 hexstr754_from_double( I )
169            
170             Converts the floating-point I into a big-endian hexadecimal representation of the underlying
171             IEEE754 encoding.
172            
173             hexstr754_from_double(12.875); # 4029C00000000000
174             # ^^^
175             # : ^^^^^^^^^^^^^
176             # : :
177             # : `- fraction
178             # :
179             # `- sign+exponent
180            
181             The first three nibbles (hexadecimal digits) encode the sign and the exponent. The sign is
182             the most significant bit of the three nibbles (so AND the first nibble with 8; if it's true,
183             the number is negative, else it's positive). The remaining 11 bits of the nibbles encode the
184             exponent: convert the 11bits to decimal, then subtract 1023. If the resulting exponent is -1023,
185             it indicates a zero or denormal value; if the exponent is +1024, it indicates an infinite (Inf) or
186             not-a-number (NaN) value, which are generally used to indicate the calculation has grown to large
187             to fit in an IEEE754 double (Inf) or has tried an performed some other undefined operation (divide
188             by zero or the logarithm of a zero or negative value) (NaN).
189            
190             The final thirteen nibbles are the encoding of the fractional value (usually C<1 + thirteennibbles /
191             16**13>, unless it's zero, denormal, infinite, or not a number).
192            
193             Of course, this is easier to decode using the L function, which interprets
194             the sign, fraction, and exponent for you. (See below for more details.)
195            
196             to_dec_floatingpoint(12.875); # +0d1.6093750000000000p+0003
197             # ^ ^^^^^^^^^^^^^^^^^^ ^^^^
198             # : : :
199             # : `- coefficient `- exponent (power of 2)
200             # :
201             # `- sign
202            
203             =head3 binstr754_from_double( I )
204            
205             Converts the floating-point I into a big-endian binary representation of the underlying
206             IEEE754 encoding.
207            
208             binstr754_from_double(12.875); # 0100000000101001110000000000000000000000000000000000000000000000
209             # ^
210             # `- sign
211             # ^^^^^^^^^^^
212             # `- exponent
213             # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
214             # `- fraction
215            
216             The first bit is the sign, the next 11 are the exponent's encoding
217            
218             =head3 hexstr754_to_double( I )
219            
220             The inverse of C: it takes a string representing the 16 nibbles
221             of the IEEE754 double value, and converts it back to a perl floating-point value.
222            
223             print hexstr754_to_double('4029C00000000000');
224             12.875
225            
226             =head3 binstr754_to_double( I )
227            
228             The inverse of C: it takes a string representing the 64 bits
229             of the IEEE754 double value, and converts it back to a perl floating-point value.
230            
231             print binstr754_to_double('0100000000101001110000000000000000000000000000000000000000000000');
232             12.875
233            
234             =cut
235             # Perl 5.10 introduced the ">" and "<" modifiers for pack which can be used to
236             # force a specific endianness.
237             if( $] lt '5.010' ) {
238             my $str = join('', unpack("H*", pack 'L' => 0x12345678));
239             if('78563412' eq $str) { # little endian, so reverse byteorder
240             *hexstr754_from_double = sub { return uc unpack('H*' => reverse pack 'd' => shift); };
241             *binstr754_from_double = sub { return uc unpack('B*' => reverse pack 'd' => shift); };
242            
243             *hexstr754_to_double = sub { return unpack('d' => reverse pack 'H*' => shift); };
244             *binstr754_to_double = sub { return unpack('d' => reverse pack 'B*' => shift); };
245            
246             *arr2x32b_from_double = sub { return unpack('N2' => reverse pack 'd' => shift); };
247             } elsif('12345678' eq $str) { # big endian, so keep default byteorder
248             *hexstr754_from_double = sub { return uc unpack('H*' => pack 'd' => shift); };
249             *binstr754_from_double = sub { return uc unpack('B*' => pack 'd' => shift); };
250            
251             *hexstr754_to_double = sub { return unpack('d' => pack 'H*' => shift); };
252             *binstr754_to_double = sub { return unpack('d' => pack 'B*' => shift); };
253            
254             *arr2x32b_from_double = sub { return unpack('N2' => pack 'd' => shift); };
255             } else {
256             # I don't handle middle-endian / mixed-endian; sorry
257             *hexstr754_from_double = sub { undef };
258             *binstr754_from_double = sub { undef };
259            
260             *hexstr754_to_double = sub { undef };
261             *binstr754_to_double = sub { undef };
262             }
263             } else {
264 4725     4725   152456 *hexstr754_from_double = sub { return uc unpack('H*' => pack 'd>' => shift); };
265 4017     4017   15244 *binstr754_from_double = sub { return uc unpack('B*' => pack 'd>' => shift); };
266            
267 1912     1912   274288 *hexstr754_to_double = sub { return unpack('d>' => pack 'H*' => shift); };
268 1967     1967   10524 *binstr754_to_double = sub { return unpack('d>' => pack 'B*' => shift); };
269            
270 318     318   707 *arr2x32b_from_double = sub { return unpack('N2' => pack 'd>' => shift); };
271             }
272            
273             =head2 :floatingpoint
274            
275             =head3 to_hex_floatingpoint( I )
276            
277             =head3 to_dec_floatingpoint( I )
278            
279             Converts value to a hexadecimal or decimal floating-point notation that indicates the sign and
280             the coefficient and the power of two, with the coefficient either in hexadecimal or decimal
281             notation.
282            
283             to_hex_floatingpoint(-3.9999999999999996) # -0x1.fffffffffffffp+0001
284             to_dec_floatingpoint(-3.9999999999999996) # -0d1.9999999999999998p+0001
285            
286             It displays the value as (sign)(0base)(implied).(fraction)p(exponent):
287            
288             =cut
289            
290             sub to_hex_floatingpoint {
291             # thanks to BrowserUK @ http://perlmonks.org/?node_id=1167146 for slighly better decision factors
292             # I tweaked it to use the two 32bit words instead of one 64bit word (which wouldn't work on some systems)
293 112     112 1 8656 my $v = shift;
294 112         145 my ($msb,$lsb) = arr2x32b_from_double($v);
295 112         112 my $sbit = ($msb & 0x80000000) >> 31;
296 112 100       166 my $sign = $sbit ? '-' : '+';
297 112         101 my $exp = (($msb & 0x7FF00000) >> 20) - 1023;
298 112         256 my $mant = sprintf '%05x%08x', $msb & 0x000FFFFF, $lsb & 0xFFFFFFFF;
299 112 100       180 if($exp == 1024) {
300 36 100       99 return $sign . "0x1.#INF000000000p+0000" if $mant eq '0000000000000';
301 30 100 100     83 return $sign . "0x1.#IND000000000p+0000" if $mant eq '8000000000000' and $sign eq '-';
302 27 100       195 return $sign . ( (($msb & 0x00080000) != 0x00080000) ? "0x1.#SNAN00000000p+0000" : "0x1.#QNAN00000000p+0000"); # v0.012 coverage note: '!=' condition only triggered on systems with SNAN; ignore Devel::Cover failures on this line on systems which quiet all SNAN to QNAN
303             }
304 76         53 my $implied = 1;
305 76 100       113 if( $exp == -1023 ) { # zero or denormal
306 18         15 $implied = 0;
307 18 100       25 $exp = $mant eq '0000000000000' ? 0 : -1022; # 0 for zero, -1022 for denormal
308             }
309 76         1151 sprintf '%s0x%1u.%13.13sp%+05d', $sign, $implied, $mant, $exp;
310             }
311            
312             sub to_dec_floatingpoint {
313             # based on to_hex_floatingpoint
314 68     68 1 65 my $v = shift;
315 68         79 my ($msb,$lsb) = arr2x32b_from_double($v);
316 68         70 my $sbit = ($msb & 0x80000000) >> 31;
317 68 100       88 my $sign = $sbit ? '-' : '+';
318 68         62 my $exp = (($msb & 0x7FF00000) >> 20) - 1023;
319 68         163 my $mant = sprintf '%05x%08x', $msb & 0x000FFFFF, $lsb & 0xFFFFFFFF;
320 68 100       115 if($exp == 1024) {
321 12 100       45 return $sign . "0d1.#INF000000000000p+0000" if $mant eq '0000000000000';
322 10 100 100     41 return $sign . "0d1.#IND000000000000p+0000" if $mant eq '8000000000000' and $sign eq '-';
323 9 100       161 return $sign . ( (($msb & 0x00080000) != 0x00080000) ? "0d1.#SNAN00000000000p+0000" : "0d1.#QNAN00000000000p+0000"); # v0.012 coverage note: '!=' condition only triggered on systems with SNAN; ignore Devel::Cover failures on this line on systems which quiet all SNAN to QNAN
324             }
325 56         42 my $implied = 1;
326 56 100       73 if( $exp == -1023 ) { # zero or denormal
327 6         4 $implied = 0;
328 6 100       11 $exp = $mant eq '0000000000000' ? 0 : -1022; # 0 for zero, -1022 for denormal
329             }
330             #$mant = (($msb & 0x000FFFFF)*4_294_967_296.0 + ($lsb & 0xFFFFFFFF)*1.0) / (2.0**52);
331             #sprintf '%s0d%1u.%.16fp%+05d', $sign, $implied, $mant, $exp;
332 56         96 my $other = abs($v) / (2.0**$exp);
333 56         1132 sprintf '%s0d%.16fp%+05d', $sign, $other, $exp;
334             }
335            
336            
337             =over
338            
339             =item sign
340            
341             The I will be + or -
342            
343             =item 0base
344            
345             The I<0base> will be C<0x> for hexadecimal, C<0d> for decimal
346            
347             =item implied.fraction
348            
349             The I indicates the hexadecimal or decimal equivalent for the coefficient
350            
351             I will be 0 for zero or denormal numbers, 1 for everything else
352            
353             I will indicate infinities (#INF), signaling not-a-numbers (#SNAN), and quiet not-a-numbers (#QNAN).
354            
355             I will range from decimal 0.0000000000000000 to 0.9999999999999998 for zero thru all the denormals,
356             and from 1.0000000000000000 to 1.9999999999999998 for normal values.
357            
358             =item p
359            
360             The I

introduces the "power" of 2. (It is analogous to the C in C<1.0e3> introducing the power of 10 in a

361             standard decimal floating-point notation, but indicates that the exponent is 2**exp instead of 10**exp.)
362            
363             =item exponent
364            
365             The I is the power of 2. Is is always a decimal number, whether the coefficient's base is hexadecimal or decimal.
366            
367             +0d1.500000000000000p+0010
368             = 1.5 * (2**10)
369             = 1.5 * 1024.0
370             = 1536.0.
371            
372             The I can range from -1022 to +1023.
373            
374             Internally, the IEEE 754 representation uses the encoding of -1023 for zero and denormals; to
375             aid in understanding the actual number, the C conversions represent
376             them as +0000 for zero, and -1022 for denormals: since denormals are C<(0+fraction)*(2**min_exp)>,
377             they are really multiples of 2**-1022, not 2**-1023.
378            
379             =back
380            
381             =head2 :constants
382            
383             These can be useful as references for the specialty values, and include the positive and negative
384             zeroes, infinities, a variety of signaling and quiet NAN values.
385            
386             POS_ZERO # +0x0.0000000000000p+0000 # signed zero (positive)
387             POS_DENORM_SMALLEST # +0x0.0000000000001p-1022 # smallest positive value that requires denormal representation in 64bit floating-point
388             POS_DENORM_BIGGEST # +0x0.fffffffffffffp-1022 # largest positive value that requires denormal representation in 64bit floating-point
389             POS_NORM_SMALLEST # +0x1.0000000000000p-1022 # smallest positive value that allows for normal representation in 64bit floating-point
390             POS_NORM_BIGGEST # +0x1.fffffffffffffp+1023 # largest positive value that allows for normal representation in 64bit floating-point
391             POS_INF # +0x1.#INF000000000p+0000 # positive infinity: indicates that the answer is out of the range of a 64bit floating-point
392             POS_SNAN_FIRST # +0x1.#SNAN00000000p+0000 # positive signaling NAN with "0x0000000000001" as the system-dependent information; note that many perl interpreters will internally convert this to a Quiet NaN (QNAN)
393             POS_SNAN_LAST # +0x1.#SNAN00000000p+0000 # positive signaling NAN with "0x7FFFFFFFFFFFF" as the system-dependent information; note that many perl interpreters will internally convert this to a Quiet NaN (QNAN)
394             POS_IND # +0x1.#QNAN00000000p+0000 # positive quiet NAN with "0x8000000000000" as the system-dependent information; some perl interpreters define the NEG_IND as an "indeterminate" value (IND), and it wouldn't surprise me if some also defined this positive version as "indeterminate" as well
395             POS_QNAN_FIRST # +0x1.#QNAN00000000p+0000 # positive quiet NAN with "0x8000000000001" as the system-dependent information
396             POS_QNAN_LAST # +0x1.#QNAN00000000p+0000 # positive quiet NAN with "0xFFFFFFFFFFFFF" as the system-dependent information
397            
398             NEG_ZERO # -0x0.0000000000000p+0000 # signed zero (negative)
399             NEG_DENORM_SMALLEST # -0x0.0000000000001p-1022 # smallest negative value that requires denormal representation in 64bit floating-point
400             NEG_DENORM_BIGGEST # -0x0.fffffffffffffp-1022 # largest negative value that requires denormal representation in 64bit floating-point
401             NEG_NORM_SMALLEST # -0x1.0000000000000p-1022 # smallest negative value that allows for normal representation in 64bit floating-point
402             NEG_NORM_BIGGEST # -0x1.fffffffffffffp+1023 # largest negative value that allows for normal representation in 64bit floating-point
403             NEG_INF # -0x1.#INF000000000p+0000 # negative infinity: indicates that the answer is out of the range of a 64bit floating-point
404             NEG_SNAN_FIRST # -0x1.#SNAN00000000p+0000 # negative signaling NAN with "0x0000000000001" as the system-dependent information; note that many perl interpreters will internally convert this to a Quiet NaN (QNAN)
405             NEG_SNAN_LAST # -0x1.#SNAN00000000p+0000 # negative signaling NAN with "0x7FFFFFFFFFFFF" as the system-dependent information; note that many perl interpreters will internally convert this to a Quiet NaN (QNAN)
406             NEG_IND # -0x1.#IND000000000p+0000 # negative quiet NAN with "0x8000000000000" as the system-dependent information; some perl interpreters define the NEG_IND as an "indeterminate" value (IND)
407             NEG_QNAN_FIRST # -0x1.#QNAN00000000p+0000 # negative quiet NAN with "0x8000000000001" as the system-dependent information
408             NEG_QNAN_LAST # -0x1.#QNAN00000000p+0000 # negative quiet NAN with "0xFFFFFFFFFFFFF" as the system-dependent information
409            
410             =cut
411            
412 10 100   10 0 68 { my $local; sub POS_ZERO () { $local = hexstr754_to_double('000'.'0000000000000') unless defined $local; return $local; } }
  10         30  
413 16 100   16 0 53 { my $local; sub POS_DENORM_SMALLEST() { $local = hexstr754_to_double('000'.'0000000000001') unless defined $local; return $local; } }
  16         35  
414 10 100   10 0 55 { my $local; sub POS_DENORM_BIGGEST () { $local = hexstr754_to_double('000'.'FFFFFFFFFFFFF') unless defined $local; return $local; } }
  10         27  
415 10 100   10 0 33 { my $local; sub POS_NORM_SMALLEST () { $local = hexstr754_to_double('001'.'0000000000000') unless defined $local; return $local; } }
  10         24  
416 10 100   10 0 30 { my $local; sub POS_NORM_BIGGEST () { $local = hexstr754_to_double('7FE'.'FFFFFFFFFFFFF') unless defined $local; return $local; } }
  10         33  
417 10 100   10 0 32 { my $local; sub POS_INF () { $local = hexstr754_to_double('7FF'.'0000000000000') unless defined $local; return $local; } }
  10         25  
418 522 100   522 0 1015 { my $local; sub POS_SNAN_FIRST () { $local = hexstr754_to_double('7FF'.'0000000000001') unless defined $local; return $local; } }
  522         709  
419 10 100   10 0 36 { my $local; sub POS_SNAN_LAST () { $local = hexstr754_to_double('7FF'.'7FFFFFFFFFFFF') unless defined $local; return $local; } }
  10         28  
420 10 100   10 0 33 { my $local; sub POS_IND () { $local = hexstr754_to_double('7FF'.'8000000000000') unless defined $local; return $local; } }
  10         29  
421 10 100   10 0 30 { my $local; sub POS_QNAN_FIRST () { $local = hexstr754_to_double('7FF'.'8000000000001') unless defined $local; return $local; } }
  10         24  
422 10 100   10 0 39 { my $local; sub POS_QNAN_LAST () { $local = hexstr754_to_double('7FF'.'FFFFFFFFFFFFF') unless defined $local; return $local; } }
  10         26  
423 6 100   6 0 37 { my $local; sub NEG_ZERO () { $local = hexstr754_to_double('800'.'0000000000000') unless defined $local; return $local; } }
  6         18  
424 6 100   6 0 28 { my $local; sub NEG_DENORM_SMALLEST() { $local = hexstr754_to_double('800'.'0000000000001') unless defined $local; return $local; } }
  6         19  
425 6 100   6 0 32 { my $local; sub NEG_DENORM_BIGGEST () { $local = hexstr754_to_double('800'.'FFFFFFFFFFFFF') unless defined $local; return $local; } }
  6         16  
426 6 100   6 0 24 { my $local; sub NEG_NORM_SMALLEST () { $local = hexstr754_to_double('801'.'0000000000000') unless defined $local; return $local; } }
  6         20  
427 6 100   6 0 28 { my $local; sub NEG_NORM_BIGGEST () { $local = hexstr754_to_double('FFE'.'FFFFFFFFFFFFF') unless defined $local; return $local; } }
  6         18  
428 6 100   6 0 28 { my $local; sub NEG_INF () { $local = hexstr754_to_double('FFF'.'0000000000000') unless defined $local; return $local; } }
  6         23  
429 6 100   6 0 25 { my $local; sub NEG_SNAN_FIRST () { $local = hexstr754_to_double('FFF'.'0000000000001') unless defined $local; return $local; } }
  6         15  
430 6 100   6 0 29 { my $local; sub NEG_SNAN_LAST () { $local = hexstr754_to_double('FFF'.'7FFFFFFFFFFFF') unless defined $local; return $local; } }
  6         18  
431 6 100   6 0 31 { my $local; sub NEG_IND () { $local = hexstr754_to_double('FFF'.'8000000000000') unless defined $local; return $local; } }
  6         23  
432 6 100   6 0 26 { my $local; sub NEG_QNAN_FIRST () { $local = hexstr754_to_double('FFF'.'8000000000001') unless defined $local; return $local; } }
  6         16  
433 6 100   6 0 56 { my $local; sub NEG_QNAN_LAST () { $local = hexstr754_to_double('FFF'.'FFFFFFFFFFFFF') unless defined $local; return $local; } }
  6         23  
434            
435             =head2 :ulp
436            
437             =head3 ulp( I )
438            
439             Returns the ULP ("Unit in the Last Place") for the given I, which is the smallest number
440             that you can add to or subtract from I and still be able to discern a difference between
441             the original and modified. Under normal (or denormal) circumstances, C $val>
442             is true.
443            
444             If the I is a zero or a denormal, C will return the smallest possible denormal.
445            
446             Since INF and NAN are not really numbers, C will just return the same I. Because
447             of the way they are handled, C $val> no longer makes sense (infinity plus
448             anything is still infinity, and adding NAN to NAN is not numerically defined, so a numerical
449             comparison is meaningless on both).
450            
451             =cut
452            
453             my $TWONEG52 = sub { 2.0**-52 };
454             my $FIFTYTWOZEROES = sub { '0'x52 };
455            
456             sub ulp { # ulp_by_div
457 30     30 1 85 my $val = shift;
458 30         39 my $rawbin = binstr754_from_double($val);
459 30         115 my ($sgn, $exp, $frac) = ($rawbin =~ /(.)(.{11})(.{52})/);
460            
461 30 100       78 return $val if $exp eq '11111111111'; # return SELF for INF or NAN
462 18 100       30 return POS_DENORM_SMALLEST if $exp eq '00000000000'; # return first positive denorm for 0 or denorm
463            
464             # this method will multiply by 2**-52 (as a constant) after
465 12         13 $sgn = '0';
466 12         17 $frac = $FIFTYTWOZEROES->();
467 12         23 $val = binstr754_to_double( $sgn . $exp . $frac );
468 12         18 $val *= $TWONEG52->();
469             }
470            
471             =head3 toggle_ulp( I )
472            
473             Returns the orginal I, but with the ULP toggled. In other words, if the ULP bit
474             was a 0, it will return a value with the ULP of 1 (equivalent to adding one ULP to a positive
475             I); if the ULP bit was a 1, it will return a value with the ULP of 0 (equivalent to
476             subtracting one ULP from a positive I). Under normal (or denormal) circumstances,
477             C is true.
478            
479             Since INF and NAN are not really numbers, C will just return the same I. Because
480             of the way they are handled, C no longer makes sense.
481            
482             =cut
483            
484             sub toggle_ulp {
485 26     26 1 56 my $val = shift;
486 26         32 my $rawbin = binstr754_from_double($val);
487 26         93 my ($sign, $exp, $fract) = ($rawbin =~ /(.)(.{11})(.{52})/);
488            
489             # INF and NAN do not have a meaningful ULP; just return SELF
490 26 100       73 if( $exp == '1'x11 ) {
491 12         21 return $val;
492             }
493            
494             # ZERO, DENORMAL, and NORMAL: toggle the last bit of fract
495 14         20 my $ulp_bit = substr $fract, -1;
496 14         23 substr $fract, -1, 1, (1-$ulp_bit);
497 14         21 $rawbin = join '', $sign, $exp, $fract;
498 14         21 return binstr754_to_double($rawbin);
499             }
500            
501             =head3 nextUp( I )
502            
503             Returns the next floating point value numerically greater than I; that is, it adds one ULP.
504             Returns infinite when I is the highest normal floating-point value.
505             Returns I when I is positive-infinite or NAN; returns the largest negative normal
506             floating-point value when I is negative-infinite.
507            
508             C is an IEEE 754r standard function (754-2008 #5.3.1).
509            
510             =cut
511            
512             sub nextUp {
513             # thanks again to BrowserUK: http://perlmonks.org/?node_id=1167146
514 176     176 1 211 my $val = shift;
515 176 100       320 return $val if $val != $val; # interestingly, NAN != NAN
516 156         237 my $h754 = hexstr754_from_double($val);
517 156 100       249 return $val if $h754 eq '7FF0000000000000'; # return self for +INF
518 154 100       206 return hexstr754_to_double('FFEFFFFFFFFFFFFF') if $h754 eq 'FFF0000000000000'; # return largest negative for -INF
519 142 100       183 return hexstr754_to_double('0000000000000001') if $h754 eq '8000000000000000'; # return +SmallestDenormal for -ZERO
520 138         168 my ($msb,$lsb) = arr2x32b_from_double($val);
521 138 100       247 $lsb += ($msb & 0x80000000) ? -1.0 : +1.0;
522 138 100       259 if($lsb == 4_294_967_296.0) {
    100          
523 22         20 $lsb = 0.0;
524 22         23 $msb += 1.0; # v0.012: LSB==4e9 only happens if you add one to LSB = 0xFFFFFFFF, so only when +msb; thus, remove extra check for msb sign here
525             } elsif ($lsb == -1.0) {
526 30         34 $lsb = 0xFFFFFFFF;
527 30         35 $msb -= 1.0; # v0.012: LSB==-1 only happens if you subtract one from LSB = 0x00000000, so only when -msb; thus, remove extra check for msb sign here
528             }
529 138         92 $msb &= 0xFFFFFFFF; # v0.011_001: bugfix: ensure 32bit MSB
530 138         99 $lsb &= 0xFFFFFFFF; # v0.011_001: bugfix: ensure 32bit MSB
531 138         466 return hexstr754_to_double( sprintf '%08X%08X', $msb, $lsb );
532             }
533            
534             =head3 nextDown( I )
535            
536             Returns the next floating point value numerically lower than I; that is, it subtracts one ULP.
537             Returns -infinity when I is the largest negative normal floating-point value.
538             Returns I when I is negative-infinite or NAN; returns the largest positive normal
539             floating-point value when I is positive-infinite.
540            
541             C is an IEEE 754r standard function (754-2008 #5.3.1).
542            
543             =cut
544            
545             sub nextDown {
546 88     88 1 187 return - nextUp( - $_[0] );
547             }
548            
549             =head3 nextAfter( I, I )
550            
551             Returns the next floating point value after I in the direction of I. If the
552             two are identical, return I; if I is numerically above I, return
553             C)>; if I is numerically below I, return C)>.
554            
555             =cut
556            
557             sub nextAfter {
558 270 100   270 1 1048 return $_[0] if $_[0] != $_[0]; # return value when value is NaN
559 160 100       294 return $_[1] if $_[1] != $_[1]; # return direction when direction is NaN
560 140 100       236 return $_[1] if $_[1] == $_[0]; # return direction when the two are equal
561 112 100       231 return nextUp($_[0]) if $_[1] > $_[0]; # return nextUp if direction > value
562 56         80 return nextDown($_[0]); # otherwise, return nextDown()
563             }
564            
565             =head2 :info
566            
567             The informational functions include various operations (defined in 754-2008 #5.7.2) that provide general
568             information about the floating-point value: most define whether a value is a special condition of
569             floating-point or not (such as normal, finite, zero, ...).
570            
571             =head3 isSignMinus( I )
572            
573             Returns 1 if I has negative sign (even applies to zeroes and NaNs); otherwise, returns 0.
574            
575             =cut
576            
577             sub isSignMinus {
578             # look at leftmost nibble, and determine whether it has the 8-bit or not (which is the sign bit)
579 1030     1030 1 237735 return (hex(substr(hexstr754_from_double(shift),0,1)) & 8) >> 3;
580             }
581            
582             =head3 isNormal( I )
583            
584             Returns 1 if I is a normal number (not zero, subnormal, infinite, or NaN); otherwise, returns 0.
585            
586             =cut
587            
588             sub isNormal {
589             # it's normal if the leftmost three nibbles (excluding sign bit) are not 000 or 7FF
590 36     36 1 5272 my $exp = hex(substr(hexstr754_from_double(shift),0,3)) & 0x7FF;
591 36   100     240 return (0 < $exp) && ($exp < 0x7FF) || 0;
592             }
593            
594             =head3 isFinite( I )
595            
596             Returns 1 if I is a finite number (zero, subnormal, or normal; not infinite or NaN); otherwise, returns 0.
597            
598             =cut
599            
600             sub isFinite {
601             # it's finite if the leftmost three nibbles (excluding sign bit) are not 7FF
602 22     22 1 5162 my $exp = hex(substr(hexstr754_from_double(shift),0,3)) & 0x7FF;
603 22   100     140 return ($exp < 0x7FF) || 0;
604             }
605            
606             =head3 isZero( I )
607            
608             Returns 1 if I is positive or negative zero; otherwise, returns 0.
609            
610             =cut
611            
612             sub isZero {
613             # it's zero if it's 0x[80]000000000000000
614 371     371 1 5520 my $str = substr(hexstr754_from_double(shift),1,15);
615 371   100     1467 return ($str eq '0'x15) || 0;
616             }
617            
618             =head3 isSubnormal( I )
619            
620             Returns 1 if I is subnormal (not zero, normal, infinite, nor NaN); otherwise, returns 0.
621            
622             =cut
623            
624             sub isSubnormal {
625             # it's subnormal if it's 0x[80]00___ and the last 13 digits are not all zero
626 36     36 1 5093 my $h = hexstr754_from_double(shift);
627 36         39 my $exp = substr($h,0,3);
628 36         33 my $frc = substr($h,3,13);
629 36   100     269 return ($exp eq '000' || $exp eq '800') && ($frc ne '0'x13) || 0;
630             }
631            
632             =head3 isInfinite( I )
633            
634             Returns 1 if I is positive or negative infinity (not zero, subnormal, normal, nor NaN); otherwise, returns 0.
635            
636             =cut
637            
638             sub isInfinite {
639             # it's infinite if it's 0x[F7]FF_0000000000000
640 35     35 1 5233 my $h = hexstr754_from_double(shift);
641 35         37 my $exp = substr($h,0,3);
642 35         28 my $frc = substr($h,3,13);
643 35   100     259 return ($exp eq '7FF' || $exp eq 'FFF') && ($frc eq '0'x13) || 0;
644             }
645            
646             =head3 isNaN( I )
647            
648             Returns 1 if I is NaN (not zero, subnormal, normal, nor infinite); otherwise, returns 0.
649            
650             =cut
651            
652             sub isNaN {
653             # it's infinite if it's 0x[F7]FF_0000000000000
654 1976     1976 1 6987 my $h = hexstr754_from_double(shift);
655 1976         1870 my $exp = substr($h,0,3);
656 1976         1594 my $frc = substr($h,3,13);
657 1976   100     9882 return ($exp eq '7FF' || $exp eq 'FFF') && ($frc ne '0'x13) || 0;
658             }
659            
660             =head3 isSignaling( I )
661            
662             Returns 1 if I is a signaling NaN (not zero, subnormal, normal, nor infinite), otherwise, returns 0.
663            
664             Note that some perl implementations convert some or all signaling NaNs to quiet NaNs, in which case,
665             C might return only 0.
666            
667             =cut
668            
669             sub isSignaling {
670             # it's signaling if isNaN and MSB of fractional portion is 1
671 556     556 1 639 my $h = hexstr754_from_double(shift);
672 556         585 my $exp = substr($h,0,3);
673 556         475 my $frc = substr($h,3,13);
674 556         654 my $qbit = (0x8 && hex(substr($h,3,1))) >> 3; # 1: quiet, 0: signaling
675 556   100     4550 return ($exp eq '7FF' || $exp eq 'FFF') && ($frc ne '0'x13) && (!$qbit) || 0;
676             }
677            
678             =head4 isSignalingConvertedToQuiet()
679            
680             Returns 1 if your implementation of perl converts a SignalingNaN to a QuietNaN, otherwise returns 0.
681            
682             This is I a standard IEEE 754 function; but this is used to determine if the C
683             function is meaningful in your implementation of perl.
684            
685             =cut
686            
687             sub isSignalingConvertedToQuiet {
688 512 50   512 1 8963 !isSignaling( POS_SNAN_FIRST ) || 0
689             }
690            
691             =head3 isCanonical( I )
692            
693             Returns 1 to indicate that I is Canonical.
694            
695             Per IEEE Std 754-2008, "Canonical" is the "preferred" encoding. Based on the
696             B's author's reading of the standard, non-canonical
697             applies to decimal floating-point encodings, not the binary floating-point
698             encodings that B handles. Since there are not multiple
699             choicesfor the representation of a binary-encoded floating-point, all
700             Is seem canonical, and thus return 1.
701            
702             =cut
703            
704 22     22 1 5224 sub isCanonical { 1 }
705            
706             =head3 class( I )
707            
708             Returns the "class" of the I:
709            
710             signalingNaN
711             quietNaN
712             negativeInfinity
713             negativeNormal
714             negativeSubnormal
715             negativeZero
716             positiveZero
717             positiveSubnormal
718             positiveNormal
719             positiveInfinity
720            
721             =cut
722            
723             sub class {
724 22 100   22 1 4231 return 'signalingNaN' if isSignaling($_[0]);
725 18 100       30 return 'quietNaN' if isNaN($_[0]);
726 12 100 100     14 return 'negativeInfinity' if isInfinite($_[0]) && isSignMinus($_[0]);
727 11 100 100     15 return 'negativeNormal' if isNormal($_[0]) && isSignMinus($_[0]);
728 9 100 100     11 return 'negativeSubnormal' if isSubnormal($_[0]) && isSignMinus($_[0]);
729 7 100 100     11 return 'negativeZero' if isZero($_[0]) && isSignMinus($_[0]);
730 6 100 66     8 return 'positiveZero' if isZero($_[0]) && !isSignMinus($_[0]);
731 5 100 66     10 return 'positiveSubnormal' if isSubnormal($_[0]) && !isSignMinus($_[0]);
732 3 100 66     5 return 'positiveNormal' if isNormal($_[0]) && !isSignMinus($_[0]);
733 1 50 33     4 return 'positiveInfinity' if isInfinite($_[0]) && !isSignMinus($_[0]);
734             }
735            
736             =head3 radix( I )
737            
738             Returns the base (radix) of the internal floating point representation.
739             This module works with the binary floating-point representations, so will always return 2.
740            
741             =cut
742            
743 22     22 1 5367 sub radix { 2 }
744            
745             =head3 totalOrder( I, I )
746            
747             Returns TRUE if I E I, FALSE if I E I.
748            
749             Special cases are ordered as below:
750            
751             -quietNaN < -signalingNaN < -infinity < ...
752             ... < -normal < -subnormal < -zero < ...
753             ... < +zero < +subnormal < +normal < ...
754             ... < +infinity < +signalingNaN < +quietNaN
755            
756             =cut
757            
758             sub totalOrder {
759 968     968 1 1982 my ($x, $y) = @_[0,1];
760 968         931 my ($bx,$by) = map { binstr754_from_double($_) } $x, $y; # convert to binary strings
  1936         1792  
761 968         3763 my @xsegs = ($bx =~ /(.)(.{11})(.{20})(.{32})/); # split into sign, exponent, MSB, LSB
762 968         2447 my @ysegs = ($by =~ /(.)(.{11})(.{20})(.{32})/); # split into sign, exponent, MSB, LSB
763 968         865 my ($xin, $yin) = map { isNaN($_) } $x, $y; # determine if NaN: used twice each, so save the values rather than running twice each during if-switch
  1936         2092  
764            
765 968 100 100     2922 if( $xin && $yin ) { # BOTH NaN
    100 100        
    100          
    100          
766             # use a trick: the rules for both-NaN treat it as if it's just another floating point,
767             # so lie about the exponent and do a normal comparison
768 200         253 ($bx, $by) = map { $_->[1] = '1' . '0'x10; join '', @$_ } \@xsegs, \@ysegs;
  400         314  
  400         840  
769 200         233 ($x, $y) = map { binstr754_to_double($_) } $bx, $by;
  400         387  
770 200   100     852 return (($x <= $y) || 0);
771             } elsif ( $xin ) { # just x NaN: TRUE if x is NEG
772 240   100     866 return ( ($xsegs[0]) || 0 );
773             } elsif ( $yin ) { # just y NaN: TRUE if y is not NEG
774 240   100     787 return ( (!$ysegs[0]) || 0 );
775             } elsif ( isZero($x) && isZero($y) ) { # both zero: TRUE if x NEG, or if x==y
776             # trick = -signbit(x) <= -signbit(y), since signbit is 1 for negative, -signbit = -1 for negative
777 8   100     42 return ( (-$xsegs[0] <= -$ysegs[0]) || 0 );
778             } else { # numeric comparison (works for inf, normal, subnormal, or only one +/-zero)
779 280   100     1041 return( ($x <= $y) || 0 );
780             }
781             }
782            
783             =head3 totalOrderMag( I, I )
784            
785             Returns TRUE if I E I, otherwise FALSE.
786             Equivalent to
787            
788             totalOrder( abs(x), abs(y) )
789            
790             Special cases are ordered as below:
791            
792             zero < subnormal < normal < infinity < signalingNaN < quietNaN
793            
794             =cut
795            
796             sub totalOrderMag {
797 484     484 1 141328 my ($x, $y) = @_[0,1];
798 484         578 my ($bx,$by) = map { binstr754_from_double($_) } $x, $y; # convert to binary strings
  968         1153  
799 484         481 ($x, $y) = map { substr $_, 0, 1, '0'; binstr754_to_double($_) } $bx, $by; # set sign bit to 0, and convert back to number
  968         981  
  968         971  
800 484         675 return totalOrder( $x, $y ); # compare normally
801             }
802            
803             # TODO = spaceship() and spaceshipMag() similar to totalOrder and totalOrderMag, but with '<' => -1, '==' => 0, '>' => +1
804            
805             =head2 :signbit
806            
807             These functions, from IEEE Std 754-2008, manipulate the sign bits
808             of the argument(s)set P.
809            
810             =head3 negate( I )
811            
812             Reverses the sign bit of I. (If the sign bit is set on I,
813             it will not be set on the output, and vice versa; this will work on
814             signed zeroes, on infinities, and on NaNs.)
815            
816             =cut
817            
818             sub negate {
819 22     22 1 7843 my $b = binstr754_from_double(shift); # convert to binary string
820 22         64 my $s = 1 - substr $b, 0, 1; # toggle sign
821 22         40 substr $b, 0, 1, $s; # replace sign
822 22         42 return binstr754_to_double($b); # convert to floating-point
823             }
824            
825             =head3 CORE::abs( I )
826             =head3 absolute( I )
827            
828             The C function behaves properly (per the IEEE 754 description)
829             for all classes of I (including signed zeroes, infinities, and NaNs),
830             per the installation test suite.
831            
832             However, C is provided as a module-based alternative.
833            
834             The test suite runs the expected input/output pairs thru both C and
835             C, so if any of those fail on your system, or
836             if you find an untested edge case where C (or C)
837             do not behave as expected, please file a bug report.
838            
839             =cut
840            
841             sub absolute {
842 44     44 1 23795 my $b = binstr754_from_double(shift); # convert to binary string
843 44         81 substr $b, 0, 1, '0'; # replace sign
844 44         83 return binstr754_to_double($b); # convert to floating-point
845             }
846            
847             =head3 copySign( I, I )
848            
849             Copies the sign from I, but uses the value from I. For example,
850            
851             $new = copySign( 1.25, -5.5); # $new is -1.25: the value of x, but the sign of y
852            
853             =cut
854            
855             sub copySign {
856 484     484 1 1486 my ($x, $y) = @_[0,1];
857 484         539 my ($bx,$by) = map { binstr754_from_double($_) } $x, $y; # convert to binary strings
  968         971  
858 484         596 substr($bx, 0, 1) = substr($by, 0, 1); # copy the sign bit from y to x
859 484         542 return binstr754_to_double($bx); # convert binary-x (with modified sign) back to double
860             }
861            
862             =head3 also exports C I C<)> (see :info)
863            
864             (:signbit also exports the C function, described in :info, above)
865            
866             =head2 :all
867            
868             Include all of the above.
869            
870             =head1 INSTALLATION
871            
872             To install this module, use your favorite CPAN client.
873            
874             For a manual install, type the following:
875            
876             perl Makefile.PL
877             make
878             make test
879             make install
880            
881             (On Windows machines, you may need to use "dmake" instead of "make".)
882            
883             =head1 SEE ALSO
884            
885             =over
886            
887             =item * L
888            
889             =item * L for
890             inspiring me to go down the IEEE754-expansion trail in perl.
891            
892             =item * L as a resource
893             for how perl interacts with the various "edge cases" (+/-infinity, L,
894             signaling and quiet L.
895            
896             =item * L: I really wanted to use this module, but it didn't get me very far down the "Tools" track,
897             and included a lot of overhead modules for its install/test that I didn't want to require for B.
898             However, I was inspired by his byteorder-dependent anonymous subs (which were in turn derived from L);
899             they were more efficient, on a per-call-to-subroutine basis, than my original inclusion of the if(byteorder) in every call to
900             the sub.
901            
902             =item * L: Similar to this module, but uses numeric manipulation.
903            
904             =back
905            
906             =head1 AUTHOR
907            
908             Peter C. Jones Cpetercj AT cpan DOT orgE>
909            
910             Please report any bugs or feature requests emailing Cbug-Data-IEEE754-Tools AT rt.cpan.orgE>
911             or thru the web interface at L.
912            
913             =head1 COPYRIGHT
914            
915             Copyright (C) 2016 Peter C. Jones
916            
917             =head1 LICENSE
918            
919             This program is free software; you can redistribute it and/or modify it
920             under the terms of either: the GNU General Public License as published
921             by the Free Software Foundation; or the Artistic License.
922            
923             See L for more information.
924            
925             =cut
926            
927             1;