File Coverage

blib/lib/Math/NV.pm

Criterion	Covered	Total	%
statement	30	108	27.7
branch	4	50	8.0
condition	0	3	0.0
subroutine	8	12	66.6
pod	0	7	0.0
total	42	180	23.3

## This file generated by InlineX::C2XS (version 0.22) using Inline::C (version 0.53)
package Math::NV;
use warnings;
use strict;
use 5.010;

require Exporter;
*import = \&Exporter::import;
require DynaLoader;

$Math::NV::VERSION = '1.02';

DynaLoader::bootstrap Math::NV $Math::NV::VERSION;

@Math::NV::EXPORT = ();
@Math::NV::EXPORT_OK = qw(
    nv nv_type mant_dig ld2binary ld_str2binary is_eq
    bin2val Cprintf Csprintf nv_mpfr is_eq_mpfr
    );

%Math::NV::EXPORT_TAGS = (all => [qw(
    nv nv_type mant_dig ld2binary ld_str2binary is_eq
    bin2val Cprintf Csprintf nv_mpfr is_eq_mpfr
    )]);

eval {require Math::MPFR;};

if($@) {$Math::NV::no_mpfr = $@}
else   {$Math::NV_no_mpfr = 0 }

$Math::NV::no_warn = 0; # set to true to disable warning about non-string argument

# %_itsa is utilised in the formulation of the diagnostic message
# when it's detected that the provided arg is not a string.

my %_itsa = (
  1 => 'UV',
  2 => 'IV',
  3 => 'NV',
  4 => 'string',
  0 => 'unknown',
);

sub dl_load_flags {0} # Prevent DynaLoader from complaining and croaking

sub ld2binary {
  my @ret = _ld2binary($_[0]);
  my $prec = pop(@ret);
  my $exp = pop(@ret);
  my $mantissa = join '', @ret;
  return ($mantissa, $exp, $prec);
}

sub ld_str2binary {
  my @ret = _ld_str2binary($_[0]);
  my $prec = pop(@ret);
  my $exp = pop(@ret);
  my $mantissa = join '', @ret;
  return ($mantissa, $exp, $prec);
}

sub bin2val {
  my($mantissa, $exp, $prec) = (shift, shift, shift);
  my $sign = $mantissa =~ /^\-/ ? '-' : '';
  # Remove everything upto and including the radix point
  # as it now contains no useful information.
  $mantissa =~ s/.+\.//;
  # For our purposes the values $prec and $exp need
  # to be reduced by 1.
  $exp--;

  # Perl bugs make the following (commented out) code unreliable,
  # so we now hand the calculations over to C.
  # (And there's no need to decrement $prec.)
  #$prec--;
  #for(0..$prec) {
  #  if(substr($mantissa, $_, 1)) {$val += 2**$exp}
  #  $exp--;
  #}
  my @mantissa = split //, $mantissa;
  my $val = _bin2val($prec, $exp, \@mantissa);
  $sign eq '-' ? return -$val : return $val;
}

sub is_eq {
  unless($Math::NV::no_warn) {
    my $itsa = $_[0];
    $itsa = _itsa($itsa); # make sure that $_[0] IOK flag is not set
    warn "Argument given to is_eq() is $_itsa{$itsa}, not a string - probably not what you want"
    if $itsa != 4;
  }
  my $nv = $_[0];
  my $check = nv($_[0]);
  return 1 if $nv == $check;
  if(mant_dig() == 64) {
    # pack/unpack like to deliver irrelevant (ie ignored) leading bytes
    # if NV is 80-bit long double
    my $first = scalar(reverse(unpack("h*", pack("F<", $nv))));
    $first = substr($first, length($first) - 20, 20);
    my $second = scalar(reverse(unpack("h*", pack("F<", $check))));
    $second = substr($second, length($second) - 20, 20);
    warn "In is_eq:\nperl: $first vs C: $second\n\n";
  }
  else {
    warn "In is_eq:\nperl: ",
      scalar(reverse(unpack("h*", pack("F<", $nv)))), " vs C: ",
      scalar(reverse(unpack("h*", pack("F<", $check)))), "\n\n";
  }

  return 0;
}

sub is_eq_mpfr {

  if($Math::NV::no_mpfr) {die "In is_eq_mpfr(): $Math::NV::no_mpfr"}

  unless($Math::NV::no_warn) {
    my $itsa = $_[0];
    $itsa = _itsa($itsa); # make sure that $_[0] IOK flag is not set
    warn "Argument given to is_eq() is $_itsa{$itsa}, not a string - probably not what you want"
    if $itsa != 4;
  }

  my $ret = 0;
  my $nv = $_[0];
  my $bits = mant_dig();
  $bits = 2098 if $bits == 106; # double double

  my $fr = Math::MPFR::Rmpfr_init2($bits);
  Math::MPFR::Rmpfr_set_str($fr, $nv, 10, 0);

  if($nv == Math::MPFR::Rmpfr_get_NV($fr, 0)) {$ret = 1}
  else {
    if($bits == 64) {
      # pack/unpack like to deliver irrelevant (ie ignored) leading bytes
      # if NV is 80-bit long double
      my $first = scalar(reverse(unpack("h*", pack("F<", $nv))));
      $first = substr($first, length($first) - 20, 20);
      my $second = scalar(reverse(unpack("h*", pack("F<", Math::MPFR::Rmpfr_get_NV($fr, 0)))));
      $second = substr($second, length($second) - 20, 20);
      warn "In is_eq_mpfr:\nperl: $first vs mpfr: $second\n\n";
    }
    else {
      warn "In is_eq_mpfr:\nperl: ",
        scalar(reverse(unpack("h*", pack("F<", $nv)))), " vs mpfr: ",
        scalar(reverse(unpack("h*", pack("F<", Math::MPFR::Rmpfr_get_NV($fr, 0))))), "\n\n";
    }
  }

  return $ret;

}

sub nv_mpfr {

  if($Math::NV::no_mpfr) {die "In nv_mpfr(): $Math::NV::no_mpfr"}

  unless($Math::NV::no_warn) {
    my $itsa = $_[0];
    $itsa = _itsa($itsa); # make sure that $_[0] IOK flag is not set
    warn "Argument given to is_eq() is $_itsa{$itsa}, not a string - probably not what you want"
    if $itsa != 4;
  }

  my $bits;

  if(!defined($_[1])) {$bits = mant_dig()}
  else {$bits = $_[1]}

  return _double_double($_[0]) if $bits == 106; # doubledouble

  my $val = Math::MPFR::Rmpfr_init2($bits);
  Math::MPFR::Rmpfr_set_str($val, $_[0], 10, 0);

  if($bits == mant_dig() ) { # 53, 64, 106 or 113 bits
    my $nv = Math::MPFR::Rmpfr_get_NV($val, 0);
    my $ret = scalar(reverse(unpack("h*", pack("F<", $nv))));

    # pack/unpack like to deliver irrelevant (ie ignored) leading bytes
    # if NV is 80-bit long double
    if($bits == 64) {
      return substr($ret, length($ret) - 20, 20);
    }
    return $ret;
  }

  if($bits == 53) {
    my $nv = Math::MPFR::Rmpfr_get_d($val, 0);
    return scalar(reverse(unpack("h*", pack("d<", $nv))));
  }

  if($bits == 64) {
    die "No _ld_bytes with this version ($Math::MPFR::VERSION) - need at least 3.27"
      unless $Math::MPFR::VERSION > '3.26';
    my @bytes = Math::MPFR::_ld_bytes($_[0], 64);
    return join('', @bytes);
  }

  if($bits == 113) {
    my $t;
    eval{$t = Math::MPFR::_have_IEEE_754_long_double();}; # needs Math-MPFR-3.33, perl-5.22.
    if(!$@ && $t) {
      my @bytes = Math::MPFR::_ld_bytes($_[0], 113);
      return join('', @bytes);
    }
    else { # assume __float128 (though that might not be the case)
      die "No _f128_bytes with this version ($Math::MPFR::VERSION) - need at least 3.26"
        unless $Math::MPFR::VERSION > 3.25;
      my @bytes = Math::MPFR::_f128_bytes($_[0], 113);
      return join('', @bytes);
    }
  }

  die "Unrecognized value for bits ($bits)";
}


sub _double_double {
  my $val = Math::MPFR::Rmpfr_init2(2098);
  Math::MPFR::Rmpfr_set_str($val, shift, 10, 0);
  my @val = _dd_obj($val);
  return [scalar(reverse(unpack("h*", pack("d<", $val[0])))),
          scalar(reverse(unpack("h*", pack("d<", $val[1]))))];
}

sub _dd_obj {
  my $obj = shift;
  my $msd = Math::MPFR::Rmpfr_get_d($obj, 0);
  $obj -= $msd;
  return ($msd, Math::MPFR::Rmpfr_get_d($obj, 0));
}


1;

__END__