DateTime-Precise
view release on metacpan or search on metacpan
lib/DateTime/Math/bigfloat.pl view on Meta::CPAN
require 'DateTime/Math/bigint.pl';
# Arbitrary length float math package
#
# by Mark Biggar
#
# number format
# canonical strings have the form /[+-]\d+E[+-]\d+/
# Input values can have inbedded whitespace
# Error returns
# 'NaN' An input parameter was "Not a Number" or
# divide by zero or sqrt of negative number
# Division is computed to
# max($div_scale,length(dividend)+length(divisor))
# digits by default.
# Also used for default sqrt scale
$div_scale = 40;
# Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
lib/DateTime/Math/bigfloat.pl view on Meta::CPAN
# fround(NSTR, SCALE) return NSTR round to SCALE digits
# ffround(NSTR, SCALE) return NSTR round at SCALEth place
# fnorm(NSTR) return (NSTR) normalize
# fsqrt(NSTR[, SCALE]) return NSTR sqrt to SCALE places
# Convert a number to canonical string form.
# Takes something that looks like a number and converts it to
# the form /^[+-]\d+E[+-]\d+$/.
sub DateTime::Math::fnorm { #(string) return fnum_str
local($_) = @_;
defined($_) or return 'NaN';
s/\s+//g; # strip white space
if (/^([+-]?)(\d*)(\.(\d*))?([Ee]([+-]?\d+))?$/ &&
($2 ne '' || defined($4))) {
my $x = defined($4) ? $4 : '';
my $y = defined($6) ? $6 : 0;
&norm(($1 ? "$1$2$x" : "+$2$x"), (($x ne '') ? $y-length($x) : $y));
} else {
'NaN';
}
}
# normalize number -- for internal use
sub norm { #(mantissa, exponent) return fnum_str
my ($mantissa, $exp) = @_;
defined($exp) or $exp = 0;
if ($mantissa eq 'NaN') {
'NaN';
} else {
# strip leading zeros
$mantissa =~ s/^([+-])0+/$1/;
if (length($mantissa) == 1) {
'+0E+0';
} else {
# strip trailing zeros
$exp += length($1) if ($mantissa =~ s/(0+)$//);
sprintf("%sE%+ld", $mantissa, $exp);
}
lib/DateTime/Math/bigfloat.pl view on Meta::CPAN
# absolute value
sub DateTime::Math::fabs { #(fnum_str) return fnum_str
local($_) = &DateTime::Math::fnorm($_[0]);
s/^-/+/; # mash sign
$_;
}
# multiplication
sub DateTime::Math::fmul { #(fnum_str, fnum_str) return fnum_str
my ($x,$y) = (&DateTime::Math::fnorm($_[0]),&DateTime::Math::fnorm($_[1]));
if ($x eq 'NaN' || $y eq 'NaN') {
'NaN';
} else {
my ($xm,$xe) = split('E',$x);
my ($ym,$ye) = split('E',$y);
&norm(&DateTime::Math::bmul($xm,$ym),$xe+$ye);
}
}
# addition
sub DateTime::Math::fadd { #(fnum_str, fnum_str) return fnum_str
my ($x,$y) = (&DateTime::Math::fnorm($_[0]),&DateTime::Math::fnorm($_[1]));
if ($x eq 'NaN' || $y eq 'NaN') {
'NaN';
} else {
my ($xm,$xe) = split('E',$x);
my ($ym,$ye) = split('E',$y);
($xm,$xe,$ym,$ye) = ($ym,$ye,$xm,$xe) if ($xe < $ye);
&norm(&DateTime::Math::badd($ym,$xm.('0' x ($xe-$ye))),$ye);
}
}
# subtraction
sub DateTime::Math::fsub { #(fnum_str, fnum_str) return fnum_str
&DateTime::Math::fadd($_[0],&DateTime::Math::fneg($_[1]));
}
# division
# args are dividend, divisor, scale (optional)
# result has at most max(scale, length(dividend), length(divisor)) digits
sub DateTime::Math::fdiv #(fnum_str, fnum_str[,scale]) return fnum_str
{
my ($x,$y,$scale) = (&DateTime::Math::fnorm($_[0]),&DateTime::Math::fnorm($_[1]),$_[2]);
if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0E+0') {
'NaN';
} else {
my ($xm,$xe) = split('E',$x);
my ($ym,$ye) = split('E',$y);
$scale = $div_scale if (!$scale);
$scale = length($xm)-1 if (length($xm)-1 > $scale);
$scale = length($ym)-1 if (length($ym)-1 > $scale);
$scale = $scale + length($ym) - length($xm);
&norm(&round(&DateTime::Math::bdiv($xm.('0' x $scale),$ym),$ym),
$xe-$ye-$scale);
}
}
# round int $q based on fraction $r/$base using $rnd_mode
sub round { #(int_str, int_str, int_str) return int_str
my ($q,$r,$base) = @_;
if ($q eq 'NaN' || $r eq 'NaN') {
'NaN';
} elsif ($rnd_mode eq 'trunc') {
$q; # just truncate
} else {
my $cmp = &DateTime::Math::bcmp(&DateTime::Math::bmul($r,'+2'),$base);
if ( $cmp < 0 ||
($cmp == 0 &&
( $rnd_mode eq 'zero' ||
($rnd_mode eq '-inf' && (substr($q,0,1) eq '+')) ||
($rnd_mode eq '+inf' && (substr($q,0,1) eq '-')) ||
($rnd_mode eq 'even' && $q =~ /[24680]$/) ||
lib/DateTime/Math/bigfloat.pl view on Meta::CPAN
} else {
&DateTime::Math::badd($q, ((substr($q,0,1) eq '-') ? '-1' : '+1'));
# round up
}
}
}
# round the mantissa of $x to $scale digits
sub DateTime::Math::fround { #(fnum_str, scale) return fnum_str
my ($x,$scale) = (&DateTime::Math::fnorm($_[0]),$_[1]);
if ($x eq 'NaN' || $scale <= 0) {
$x;
} else {
my ($xm,$xe) = split('E',$x);
if (length($xm)-1 <= $scale) {
$x;
} else {
&norm(&round(substr($xm,0,$scale+1),
"+0".substr($xm,$scale+1,1),"+10"),
$xe+length($xm)-$scale-1);
}
}
}
# round $x at the 10 to the $scale digit place
sub DateTime::Math::ffround { #(fnum_str, scale) return fnum_str
my ($x,$scale) = (&DateTime::Math::fnorm($_[0]),$_[1]);
if ($x eq 'NaN') {
'NaN';
} else {
my ($xm,$xe) = split('E',$x);
if ($xe >= $scale) {
$x;
} else {
$xe = length($xm)+$xe-$scale;
if ($xe < 1) {
'+0E+0';
} elsif ($xe == 1) {
&norm(&round('+0',"+0".substr($xm,1,1),"+10"), $scale);
lib/DateTime/Math/bigfloat.pl view on Meta::CPAN
}
}
}
}
# compare 2 values returns one of undef, <0, =0, >0
# returns undef if either or both input value are not numbers
sub DateTime::Math::fcmp #(fnum_str, fnum_str) return cond_code
{
my ($x, $y) = (&DateTime::Math::fnorm($_[0]),&DateTime::Math::fnorm($_[1]));
if ($x eq "NaN" || $y eq "NaN") {
return;
} else {
# Compare signs between the two numbers.
my $ret = (ord($y) <=> ord($x));
$ret and return $ret;
# Compare the numbers by making both of them integer and using the
# integer compare routines. Make the numbers into integers by
# taking the number with the larger exponent and adding either
# abs($xe - $ye) to the end of it so that the two numbers have the
# same exponent.
my ($xm,$xe,$ym,$ye) = split('E', $x."E$y");
my $diff = abs($xe - $ye);
(($xe > $ye) ? $xm : $ym) .= '0' x $diff;
&DateTime::Math::bcmp($xm,$ym) <=> 0;
}
}
# square root by Newtons method.
sub DateTime::Math::fsqrt { #(fnum_str[, scale]) return fnum_str
my ($x, $scale) = (&DateTime::Math::fnorm($_[0]), $_[1]);
if ($x eq 'NaN' || $x =~ /^-/) {
'NaN';
} elsif ($x eq '+0E+0') {
'+0E+0';
} else {
my ($xm, $xe) = split('E',$x);
$scale = $div_scale if (!$scale);
$scale = length($xm)-1 if ($scale < length($xm)-1);
my ($gs, $guess) = (1, sprintf("1E%+d", (length($xm)+$xe-1)/2));
while ($gs < 2*$scale) {
$guess = &DateTime::Math::fmul(&DateTime::Math::fadd($guess,&DateTime::Math::fdiv($x,$guess,$gs*2)),".5");
$gs *= 2;
lib/DateTime/Math/bigint.pl view on Meta::CPAN
# /^\s*[+-]?[\d\s]+$/.
# Examples:
# '+0' canonical zero value
# ' -123 123 123' canonical value '-123123123'
# '1 23 456 7890' canonical value '+1234567890'
# Output values always always in canonical form
#
# Actual math is done in an internal format consisting of an array
# whose first element is the sign (/^[+-]$/) and whose remaining
# elements are base 100000 digits with the least significant digit first.
# The string 'NaN' is used to represent the result when input arguments
# are not numbers, as well as the result of dividing by zero
#
# routines provided are:
#
# bneg(BINT) return BINT negation
# babs(BINT) return BINT absolute value
# bcmp(BINT,BINT) return CODE compare numbers (undef,<0,=0,>0)
# badd(BINT,BINT) return BINT addition
# bsub(BINT,BINT) return BINT subtraction
# bmul(BINT,BINT) return BINT multiplication
lib/DateTime/Math/bigint.pl view on Meta::CPAN
# bmod(BINT,BINT) return BINT modulus
# bgcd(BINT,BINT) return BINT greatest common divisor
# bnorm(BINT) return BINT normalization
#
$zero = 0;
# normalize string form of number. Strip leading zeros. Strip any
# white space and add a sign, if missing.
# Strings that are not numbers result the value 'NaN'.
sub DateTime::Math::bnorm { #(num_str) return num_str
local($_) = @_;
defined($_) or return 'NaN';
s/\s+//g; # strip white space
if (s/^([+-]?)0*(\d+)$/$1$2/) { # test if number
substr($_,0,0) = '+' unless $1; # Add missing sign
s/^-0/+0/;
$_;
} else {
'NaN';
}
}
# Convert a number from string format to internal base 100000 format.
# Assumes normalized value as input.
sub internal { #(num_str) return int_num_array
my $d = shift;
my ($is,$il) = (substr($d,0,1),length($d)-2);
substr($d,0,1) = '';
($is, reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $d)));
lib/DateTime/Math/bigint.pl view on Meta::CPAN
sub abs { # post-normalized abs for internal use
my $num = shift;
$num =~ s/^-/+/;
$num;
}
# Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
sub DateTime::Math::bcmp { #(num_str, num_str) return cond_code
my ($x,$y) = (&DateTime::Math::bnorm($_[0]),&DateTime::Math::bnorm($_[1]));
if ($x eq 'NaN') {
return;
} elsif ($y eq 'NaN') {
return;
} else {
&cmp($x,$y);
}
}
sub cmp { # post-normalized compare for internal use
my ($cx, $cy) = @_;
return 0 if ($cx eq $cy);
lib/DateTime/Math/bigint.pl view on Meta::CPAN
return -1 if ($sy eq '+');
$ld = length($cy) - length($cx);
return $ld if ($ld);
return $cy cmp $cx;
}
}
sub DateTime::Math::badd { #(num_str, num_str) return num_str
my ($x, $y) = (&DateTime::Math::bnorm($_[0]),&DateTime::Math::bnorm($_[1]));
if ($x eq 'NaN') {
'NaN';
} elsif ($y eq 'NaN') {
'NaN';
} else {
my @x = &internal($x); # convert to internal form
my @y = &internal($y);
my ($sx, $sy) = (shift @x, shift @y); # get signs
if ($sx eq $sy) {
&external($sx, &add(\@x, \@y)); # if same sign add
} else {
($x, $y) = (&abs($x),&abs($y)); # make abs
if (&cmp($y,$x) > 0) {
&external($sy, &sub(\@y, \@x));
lib/DateTime/Math/bigint.pl view on Meta::CPAN
}
}
sub DateTime::Math::bsub { #(num_str, num_str) return num_str
&DateTime::Math::badd($_[0],&DateTime::Math::bneg($_[1]));
}
# GCD -- Euclids algorithm Knuth Vol 2 pg 296
sub DateTime::Math::bgcd { #(num_str, num_str) return num_str
my ($x,$y) = (&DateTime::Math::bnorm($_[0]),&DateTime::Math::bnorm($_[1]));
if ($x eq 'NaN' || $y eq 'NaN') {
'NaN';
} else {
($x, $y) = ($y,&DateTime::Math::bmod($x,$y)) while $y ne '+0';
$x;
}
}
# routine to add two base 1e5 numbers
# stolen from Knuth Vol 2 Algorithm A pg 231
# there are separate routines to add and sub as per Kunth pg 233
sub add { #(int_num_array, int_num_array) return int_num_array
lib/DateTime/Math/bigint.pl view on Meta::CPAN
for my $sx (@$sx_array) {
last unless @$sy_array || $bar;
$sx += 1e5 if $bar = (($sx -= (@$sy_array ? shift(@$sy_array) : 0) + $bar) < 0);
}
@$sx_array;
}
# multiply two numbers -- stolen from Knuth Vol 2 pg 233
sub DateTime::Math::bmul { #(num_str, num_str) return num_str
my ($x, $y) = (&DateTime::Math::bnorm($_[0]), &DateTime::Math::bnorm($_[1]));
if ($x eq 'NaN') {
'NaN';
} elsif ($y eq 'NaN') {
'NaN';
} else {
my @x = &internal($x);
my @y = &internal($y);
&external(&mul(\@x, \@y));
}
}
# multiply two numbers in internal representation
# destroys the arguments, supposes that two arguments are different
sub mul { #(*int_num_array, *int_num_array) return int_num_array
lib/DateTime/Math/bigint.pl view on Meta::CPAN
}
# modulus
sub DateTime::Math::bmod { #(num_str, num_str) return num_str
(&DateTime::Math::bdiv(@_))[1];
}
sub DateTime::Math::bdiv { #(dividend: num_str, divisor: num_str) return num_str
my ($x, $y) = (&DateTime::Math::bnorm($_[0]), &DateTime::Math::bnorm($_[1]));
return wantarray ? ('NaN','NaN') : 'NaN'
if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0');
return wantarray ? ('+0',$x) : '+0' if (&cmp(&abs($x),&abs($y)) < 0);
my @x = &internal($x);
my @y = &internal($y);
my $srem = $y[0];
my $sr = (shift @x ne shift @y) ? '-' : '+';
my ($car, $bar, $prd, $dd) = (0, 0, 0, 0);
if (($dd = int(1e5/($y[$#y]+1))) != 1) {
for $x (@x) {
$x = $x * $dd + $car;
$x -= ($car = int($x * 1e-5)) * 1e5;
lib/DateTime/Math/bigint.pl view on Meta::CPAN
}
(&external($sr, @q), &external($srem, @d, $zero));
} else {
&external($sr, @q);
}
}
# compute power of two numbers -- stolen from Knuth Vol 2 pg 233
sub DateTime::Math::bpow { #(num_str, num_str) return num_str
my ($x, $y) = (&DateTime::Math::bnorm($_[0]), &DateTime::Math::bnorm($_[1]));
if ($x eq 'NaN') {
'NaN';
} elsif ($y eq 'NaN') {
'NaN';
} elsif ($x eq '+1') {
'+1';
} elsif ($x eq '-1') {
&bmod($x,2) ? '-1': '+1';
} elsif ($y =~ /^-/) {
'NaN';
} elsif ($x eq '+0' && $y eq '+0') {
'NaN';
} else {
my @x = &internal($x);
my @pow2 = @x;
my @pow = &internal("+1");
my ($y1, $res, @tmp1, @tmp2) = (1); # need tmp to send to mul
while ($y ne '+0') {
($y,$res)=&DateTime::Math::bdiv($y,2);
if ($res ne '+0') {my @tmp=@pow2; @pow =&mul(\@pow, \@tmp);}
if ($y ne '+0') {my @tmp=@pow2; @pow2=&mul(\@pow2, \@tmp);}
}
( run in 0.738 second using v1.01-cache-2.11-cpan-4d50c553e7e )