Data-Float-DoubleDouble
view release on metacpan or search on metacpan
DoubleDouble.pm view on Meta::CPAN
}
##############################
##############################
# A function to return the exponents of the NV
sub get_exp {
my $hex = NV2H($_[0]);
my $exp1 = hex(substr($hex, 0, 3));
my $exp2 = hex(substr($hex, 16, 3));
$exp1 -= 2048 if $exp1 > 2047; # Remove sign bit
$exp2 -= 2048 if $exp2 > 2047; # Remove sign bit
$exp1++ unless $exp1; # increment if 0
$exp2++ unless $exp2; # increment if 0
return ($exp1 - 1023, $exp2 - 1023);
}
##############################
##############################
# Return the 2 doubles encapsulated in the given doubledouble.
sub get_doubles {
my $ld_hex = NV2H($_[0]);
return (H2D(substr($ld_hex, 0, 16)), H2D(substr($ld_hex, 16, 16)));
}
##############################
##############################
# Return the (hex) mantissas of the 2 doubles.
# Does not return the leading (implied) bit, but
# returns the hex representation of the next 52 bits.
sub get_mant_H {
my $hex = shift;
return (substr($hex, 3, 13), substr($hex,19, 13));
}
##############################
##############################
# Return a 3-element list for the given double-double:
# 1) sign
# 2) mantissa (in binary, implicit radix point after first digit)
# 3) exponent
# For nan/inf, the mantissa is 'nan' or 'inf' respectively unless
# 2nd arg is literally 'raw'.
sub float_B {
my $hex = NV2H($_[0]);
my $raw = @_ == 2 && $_[1] eq 'raw' ? 1 : 0;
# If the 2nd arg is 'raw' we do the calculations for the arg
# even if it is an Inf/NaN.
unless($raw) {
if($hex eq '7ff00000000000000000000000000000') { # +inf
return ('+', 'inf', 1024);
}
if($hex eq 'fff00000000000000000000000000000') { # -inf
return ('-', 'inf', 1024);
}
if($hex eq '7ff80000000000000000000000000000') { # + nan
return ('+', 'nan', 1024);
}
if($hex eq 'fff80000000000008000000000000000') { # - nan
return ('-', 'nan', 1024);
}
}
my $pre1 = hex(substr($hex, 0, 3));
my $pre2 = hex(substr($hex, 16, 3));
my $discard = 0;
my ($sign1, $sign2) = ('+', '+');
if($pre1 > 2047) {
$pre1 -= 2048;
$sign1 = '-';
}
my $single_sign = $sign1;
if($pre2 > 2047) {
$pre2 -= 2048;
$sign2 = '-';
}
my($s1, $s2) = get_mant_H($hex);
die "\$s1 is too long: ", length($s1)
if length $s1 > 13;
die "\$s2 is too long: ", length($s2)
if length $s2 > 13;
my $bin_str1 = unpack("B52", (pack "H*", $s1));
my $bin_str2 = unpack("B52", (pack "H*", $s2));
my $sign_compare;
# Check whether either string is zero and modify signs accordingly
# Check $bin_str2 *first*
$sign2 = $sign1 if ($bin_str2 !~ /1/ && !$pre2);
$sign1 = $sign2 if ($bin_str1 !~ /1/ && !$pre1);
my $bin_pre1 = $pre1 ? '1' : '0';
$pre1++ unless $pre1;
my $bin_pre2 = $pre2 ? '1' : '0';
$pre2++ unless $pre2;
my($exp1, $exp2) = ($pre1 - 1023, $pre2 - 1023);
my $inter_zero = inter_zero($exp1, $exp2);
DoubleDouble.pm view on Meta::CPAN
if($inter_zero < 0) {
$bin_pre2 = '';
$inter_zero++;
$bin_str2 = substr($bin_str2, $inter_zero * -1);
}
my($bin, $pow_adjust);
if($sign1 eq $sign2) {
$pow_adjust = 0;
$bin = $bin_pre1 . $bin_str1 . $zeroes . $bin_pre2 . $bin_str2;
}
else {
($bin, $pow_adjust) = _subtract_p($bin_pre1 . $bin_str1, $zeroes . $bin_pre2 . $bin_str2);
}
my $single_exp = $pre1 - 1023 - $pow_adjust;
$bin =~ s/0+$//; # Remove trailing zeroes
$bin .= '0' while length($bin) < 108; # Make sure the mantissa of $hex is always at least 108 bits.
$bin .= '0' while length($bin) % 4;
return($single_sign, $bin, $single_exp);
}
##############################
##############################
# Return a 3-element list for the given double-double:
# 1) sign
# 2) mantissa (in binary, implicit radix point after first digit)
# 3) exponent
# For nan/inf, the mantissa is 'nan' or 'inf' respectively.
sub NV2binary {
my @ret = _NV2binary($_[0]);
my $prec = pop(@ret);
my $exp = pop(@ret);
my $mantissa = join '', @ret;
my $sign = substr($mantissa, 0, 1, '');
$mantissa =~ s/0\.//;
$exp--; # because we've right shifted the radix point by one place.
return ($sign, $mantissa, $exp);
}
##############################
##############################
# Return the NV from the binary representation (sign, mantissa, exponent).
sub B_float {
die "Wrong number of args to B_float (", scalar @_, ")"
unless @_ == 3;
my $hex = B2float_H(@_);
return H_float($hex);
}
##############################
##############################
# Return a hex string representation as per perl Data::Float
# For NaN and Inf returns 'nan' or 'inf' (prefixed with either
# '+' or '-' as appropriate) unless an additional arg of 'raw'
# has been provided - in which case it does the calculations
# and returns the hex string it has calculated.
# This function returns a hex representation of the *actual*
# value - even if that value requires more than 106 bits.
sub float_H {
my ($sign, $mant, $exp);
if(@_ == 1) {($sign, $mant, $exp) = float_B($_[0])}
elsif(@_ == 2) {
if($_[1] eq 'raw') {
($sign, $mant, $exp) = float_B($_[0], 'raw');
}
else {
($sign, $mant, $exp) = float_B($_[0]);
}
}
else { die "Expected either 1 or 2 args to float_H() - received ", scalar @_}
if($mant eq 'nan') {
$sign eq '-' ? return '-nan'
: return '+nan';
}
if($mant eq 'inf') {
$sign eq '-' ? return '-inf'
: return '+inf';
}
my $mant_len = length $mant;
# Mantissa returned by float_B is at least 108 bits
die "Mantissa calculated by float_H() is too short ($mant_len)"
if $mant_len < 108;
# Length of mantissa returned by float_B() is always
# evenly divisible by 4
die "Mantissa calculated by float_H() is not divisible by 4 ($mant_len)"
if $mant_len % 4;
my $prefix = $sign . '0x' . substr($mant, 0, 1, '') . '.';
$mant .= '0' while length($mant) % 4;
my $middle = _bin2hex($mant);
my $suffix = "p$exp";
return $prefix . $middle . $suffix;
}
##############################
##############################
# Standardise the float_H output to match the "%La" or "%LA" format
# that C provides on my machine.
sub std_float_H {
my $str = float_H($_[0]);
$str =~ s/^\+//;
DoubleDouble.pm view on Meta::CPAN
||
($mant eq '1' . '0' x 107 && $exp == 1024)) {
$sign eq '-' ? !$raw ? return '-inf'
: return '-0x1.000000000000000000000000000p1024'
: !$raw ? return '+inf'
: return '+0x1.000000000000000000000000000p1024';
}
if($mant eq 'nan'
||
($mant eq '11' . '0' x 106 && $exp == 1024)) {
$sign eq '-' ? !$raw ? return '-nan'
: return '-0x1.800000000000000000000000000p1024'
: !$raw ? return '+nan'
: return '+0x1.800000000000000000000000000p1024';
}
my $lead = substr($mant, 0, 1, '');
$mant .= '0' while (length($mant) % 4); # _bin2hex() expects length($mant) to be divisible by 4.
$mant = _bin2hex($mant);
return $sign . '0x' . $lead . '.' . $mant . 'p' . $exp;
}
##############################
##############################
##############################
##############################
# A function to return the number of zeroes between the 2 doubles.
# Takes the 2 exponents (eg, as provided by get_exp) as args.
sub inter_zero {
die "inter_zero() takes 2 arguments"
unless @_ == 2;
my($exp1, $exp2) = (shift, shift);
return $exp1 - 53 - $exp2;
}
##############################
##############################
# Return true iff at least one argument is infinite.
sub are_inf {
for(@_) {
if($_ == 0 || $_ / $_ == 1 || $_ != $_) {
return 0;
}
}
return 1;
}
##############################
##############################
# Return true iff at least one argument is a NaN.
sub are_nan {
for(@_) {
return 0 if $_ == $_;
}
return 1;
}
##############################
##############################
##############################
##############################
##############################
##############################
##############################
# Binary subtract second arg from first arg - args must be of same length.
sub _subtract_b {
my($bin_str1, $bin_str2) = (shift, shift);
my($len1, $len2) = (length $bin_str1, length $bin_str2);
if($len1 != $len2) {
warn "\n$bin_str1\n$bin_str2\n";
die "Binary strings must be of same length - we have lengths $len1 & $len2";
}
my $ret = '';
my $borrow = 0;
for(my $i = -1; $i >= -$len1; $i--) {
my $bottom = substr($bin_str2, $i, 1);
if($borrow) {
$bottom++;
$borrow = 0;
}
my $top = substr($bin_str1, $i, 1);
if($bottom > $top) {
$top += 2;
$borrow = 1;
}
$ret = ($top - $bottom) . $ret;
}
die "_subtract_b returned wrong value: $ret"
if length $ret != $len1;
return $ret;
}
##############################
##############################
# Binary-subtract the second arg from the first arg.
# This sub written specifically for float_B() the output of which,
# is, in turn, needed for float_H().
DoubleDouble.pm view on Meta::CPAN
values, positive and zero exponents are prefixed with a +,
trailing zeroes in the mantissa are removed, and zeroes are
presented as (-)0x0p+0 or (-)0X0P+0. As for DD2HEX, the second
arg ($fmt) can be either "%La" or "%LA" (nothing else) and that
determines whether the alphabetic characters are lower case or
upper case.
Unlike float_H, this function cannot take the 'raw' argument.
Like float_H it will, however, accurately express the value
that's encapsulated in the double-double (even though that
minimum may exceed the usual 27 hex digits).
#############################################
$readable = express($nv, $opt); # $opt is an optional arg.
An alternative way of assessing the value of the double-double.
Express the double as msd + lsd, where the 2 doubles (msd and lsd)
are written in scientic notation. The doubles will be written in
decimal format unless a second arg of 'h' or 'H' is provided - in
which case they will be written in hex (respectively capitalised
hex) format.
The second arg ($opt), if provided, must be either 'h' or 'H'.
#############################################
$nv = H_float($hex);
For $hex written in the format returned by float_H(), returns
the NV that corresponds to $hex.
#############################################
@bin = float_B($nv, $opt); # Second arg isoptional
Returns the sign, the mantissa (as a base 2 string), and the
exponent of $nv. (There's an implied radix point between the
first and second digits of the mantissa).
For nan/inf, the mantissa is 'nan' or 'inf' respectively unless
2nd arg is literally 'raw' - in which case it will be a base 2
version of the nan/inf encoding.
#############################################
@bin = float_H2B($hex, $opt); # Second arg is optional
As for the above float_B() function - but takes the hex
string of the NV (as returned by float_H) as its argument,
instead of the actual NV.
For a more direct way of obtaining the array, use float_B
instead.
If a second arg is provided, it must be the string 'raw' - in
which case inf/nan mantissas will be returned in hex format
instead of as "inf"/"nan" strings.
#############################################
@bin = NV2binary($nv);
Another way of arriving at (almost) the same binary representation
of the NV -ie as an array consisting of (sign, mantissa, exponent).
The mantissa if Infs and NaNs will be returned as 'inf' or 'nan'
respectively and the sign associated with the nan will always
be '+'.
With this function, trailing zeroes are stripped from the mantissa
and exponents for 0, inf and nan might not match the other binary
representations.
This function is based on code from the mpfr library's
tests/tset_ld.c file.
#############################################
$hex = B2float_H(@bin, $opt); # $opt is an optional arg
The reverse of float_H2B. It takes the array returned by
either float_B or float_H2B as its arguments, and returns
the corresponding hex form.
If $opt is provided and is the string 'raw', the actual
hex encoding of any nan/inf will be returned - instead of
the string "inf" or "nan" respectively.
#############################################
($sign1, $sign2) = get_sign($nv);
Returns the signs of the two doubles contained in $nv.
#############################################
($exp1, $exp2) = get_exp($nv);
Returns the exponents of the two doubles contained in $nv.
#############################################
($double1, $double2) = get_doubles($nv);
Returns the two doubles contained in $nv.
#############################################
($mantissa1, $mantissa2) = get_mant_H(NV2H($nv));
Returns an array of the two 52-bit mantissa components of
the two doubles in their hex form. The values of the
implied leading (most significant) bits are not provided,
nor are the values of the two exponents.
#############################################
$intermediate_zeroes = inter_zero(get_exp($nv));
Returns the number of zeroes that need to come between the
mantissas of the 2 doubles when $nv is translated to the
representation that float_H() returns.
#############################################
$bool = are_inf(@nv); # Aliased to float_is_infinite.
Returns true if and only if all of the (NV) arguments are
infinities.
Else returns false.
#############################################
$bool = are_nan(@nv); # Aliased to float_is_nan.
Returns true if and only if all of the (NV) arguments are
NaNs. Else returns false.
#############################################
$hex = dd_bytes($nv);
Returns same as NV2H($nv).
#############################################
For Compatibility with Data::Float:
#############################################
$class = float_class($nv);
Returns one of either "NAN", "INFINITE", "ZERO", "NORMAL"
or "SUBNORMAL" - whichever is appropriate. (The NV must
belong to one (and only one) class.
#############################################
$bool = float_is_nan($nv); # Alias for are_nan()
Returns true if $nv is a NaN.
Else returns false.
#############################################
$bool = float_is_infinite($nv); # Alias for are_inf()
Returns true if $nv is infinite.
Else returns false.
#############################################
$bool = float_is_finite($nv);
Returns true if NV is neither infinite nor a NaN.
Else returns false.
#############################################
$bool = float_is_nzfinite($nv);
Returns true if NV is neither infinite, nor a NaN, nor zero.
Else returns false.
#############################################
$bool = float_is_zero($nv);
Returns true if NV is zero.
Else returns false.
#############################################
$bool = float_is_normal($nv);
Returns true if NV is finite && non-zero && the implied
leading digit in its internal representation is '1'.
Else returns false.
#############################################
$bool = float_is_subnormal($nv);
Returns true if NV is finite && non-zero && the implied
leading digit in its internal representation is '0'.
#############################################
$nv = nextafter($nv1, $nv2);
$nv1 and $nv2 must both be floating point values. Returns the
next representable floating point value adjacent to $nv1 in the
direction of $nv2, or returns $nv2 if it is numerically
equal to $nv1. Infinite values are regarded as being adjacent to
the largest representable finite values. Zero counts as one value,
even if it is signed, and it is adjacent to the positive and
negative smallest representable finite values. If a zero is returned
then it has the same sign as $nv1. Returns
NaN if either argument is a NaN.
#############################################
$nv = nextup($nv1);
$nv1 must be a floating point value. Returns the next representable
floating point value adjacent to $nv1 with a numerical value that
is strictly greater than $nv1, or returns $nv1 unchanged if there
is no such value. Infinite values are regarded as being adjacent to
the largest representable finite values. Zero counts as one value,
even if it is signed, and it is adjacent to the smallest
representable positive and negative finite values. If a zero is
returned, because $nv1 is the smallest representable negative
value, and zeroes are signed, it is a negative zero that is
returned. Returns NaN if $nv1 is a NaN.
#############################################
$nv = nextdown($nv1);
$nv1 must be a floating point value. Returns the next representable
floating point value adjacent to $nv1 with a numerical value that
is strictly less than $nv1, or returns $nv1 unchanged if there is
no such value. Infinite values are regarded as being adjacent to the
largest representable finite values. Zero counts as one value, even
if it is signed, and it is adjacent to the smallest representable
positive and negative finite values. If a zero is returned, because
$nv is the smallest representable positive value, and zeroes are
signed, it is a positive zero that is returned. Returns NaN if VALUE
is a NaN.
#############################################
#############################################
=head1 TODO
Over time, introduce the features of (and functions provided by)
Data::Float
=head1 LICENSE
This program is free software; you may redistribute it and/or
modify it under the same terms as Perl itself.
Copyright 2014 Sisyphus
=head1 AUTHOR
Sisyphus <sisyphus at(@) cpan dot (.) org>
=cut
( run in 0.809 second using v1.01-cache-2.11-cpan-39bf76dae61 )