Cross Reference: bigfloat.pl
xref
: /
osnet-11
/
usr
/
src
/
cmd
/
perl
/
5.8.4
/
distrib
/
lib
/
bigfloat.pl
Home
History
Annotate
Line#
Navigate
Download
Search
only in
./
1
N/A
package
bigfloat
;
1
N/A
require
"
bigint.pl
"
;
1
N/A
#
1
N/A
# This library is no longer being maintained, and is included for backward
1
N/A
# compatibility with Perl 4 programs which may require it.
1
N/A
#
1
N/A
# In particular, this should not be used as an example of modern Perl
1
N/A
# programming techniques.
1
N/A
#
1
N/A
# Suggested alternative: Math::BigFloat
1
N/A
#
1
N/A
# Arbitrary length float math package
1
N/A
#
1
N/A
# by Mark Biggar
1
N/A
#
1
N/A
# number format
1
N/A
# canonical strings have the form /[+-]\d+E[+-]\d+/
1
N/A
# Input values can have embedded whitespace
1
N/A
# Error returns
1
N/A
# 'NaN' An input parameter was "Not a Number" or
1
N/A
# divide by zero or sqrt of negative number
1
N/A
# Division is computed to
1
N/A
# max($div_scale,length(dividend)+length(divisor))
1
N/A
# digits by default.
1
N/A
# Also used for default sqrt scale
1
N/A
1
N/A
$
div_scale
=
40
;
1
N/A
1
N/A
# Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
1
N/A
1
N/A
$
rnd_mode
=
'even'
;
1
N/A
1
N/A
# bigfloat routines
1
N/A
#
1
N/A
# fadd(NSTR, NSTR) return NSTR addition
1
N/A
# fsub(NSTR, NSTR) return NSTR subtraction
1
N/A
# fmul(NSTR, NSTR) return NSTR multiplication
1
N/A
# fdiv(NSTR, NSTR[,SCALE]) returns NSTR division to SCALE places
1
N/A
# fneg(NSTR) return NSTR negation
1
N/A
# fabs(NSTR) return NSTR absolute value
1
N/A
# fcmp(NSTR,NSTR) return CODE compare undef,<0,=0,>0
1
N/A
# fround(NSTR, SCALE) return NSTR round to SCALE digits
1
N/A
# ffround(NSTR, SCALE) return NSTR round at SCALEth place
1
N/A
# fnorm(NSTR) return (NSTR) normalize
1
N/A
# fsqrt(NSTR[, SCALE]) return NSTR sqrt to SCALE places
1
N/A
1
N/A
# Convert a number to canonical string form.
1
N/A
# Takes something that looks like a number and converts it to
1
N/A
# the form /^[+-]\d+E[+-]\d+$/.
1
N/A
sub
main
'fnorm { #(string) return fnum_str
1
N/A
local
($_) = @_;
1
N/A
s/\s+//g;
# strip white space
1
N/A
if
(/^([+-]?)(\d*)(\.(\d*))?([
Ee
]([+-]?\d+))?$/
1
N/A
&& ($
2
ne
''
||
defined
($
4
))) {
1
N/A
my
$x =
defined
($
4
) ? $
4
:
''
;
1
N/A
&
norm
(($
1
?
"$1$2$x"
:
"+$2$x"
), (($x
ne
''
) ? $
6
-
length
($x) : $
6
));
1
N/A
}
else
{
1
N/A
'NaN'
;
1
N/A
}
1
N/A
}
1
N/A
1
N/A
# normalize number -- for internal use
1
N/A
sub
norm
{
#(mantissa, exponent) return fnum_str
1
N/A
local
($_, $
exp
) = @_;
1
N/A
if
($_
eq
'NaN'
) {
1
N/A
'NaN'
;
1
N/A
}
else
{
1
N/A
s/^([+-])
0
+/$
1
/;
# strip leading zeros
1
N/A
if
(
length
($_) ==
1
) {
1
N/A
'+0E+0'
;
1
N/A
}
else
{
1
N/A
$
exp
+=
length
($
1
)
if
(s/(
0
+)$//);
# strip trailing zeros
1
N/A
sprintf
(
"%sE%+ld"
, $_, $
exp
);
1
N/A
}
1
N/A
}
1
N/A
}
1
N/A
1
N/A
# negation
1
N/A
sub
main
'fneg { #(fnum_str) return fnum_str
1
N/A
local
($_) = &
'fnorm($_[$[]);
1
N/A
vec
($_,
0
,
8
) ^=
ord
(
'+'
) ^
ord
(
'-'
)
unless
$_
eq
'+0E+0'
;
# flip sign
1
N/A
if
(
ord
(
"\t"
) ==
9
) {
# ascii
1
N/A
s/^H/N/;
1
N/A
}
1
N/A
else
{
# ebcdic character set
1
N/A
s/\
373
/N/;
1
N/A
}
1
N/A
$_;
1
N/A
}
1
N/A
1
N/A
# absolute value
1
N/A
sub
main
'fabs { #(fnum_str) return fnum_str
1
N/A
local
($_) = &
'fnorm($_[$[]);
1
N/A
s/^-/+/;
# mash sign
1
N/A
$_;
1
N/A
}
1
N/A
1
N/A
# multiplication
1
N/A
sub
main
'fmul { #(fnum_str, fnum_str) return fnum_str
1
N/A
local
($x,$y) = (&
'fnorm($_[$[]),&'
fnorm
($_[$[+
1
]));
1
N/A
if
($x
eq
'NaN'
|| $y
eq
'NaN'
) {
1
N/A
'NaN'
;
1
N/A
}
else
{
1
N/A
local
($
xm
,$
xe
) =
split
(
'E'
,$x);
1
N/A
local
($
ym
,$
ye
) =
split
(
'E'
,$y);
1
N/A
&
norm
(&
'bmul($xm,$ym),$xe+$ye);
1
N/A
}
1
N/A
}
1
N/A
1
N/A
# addition
1
N/A
sub
main
'fadd { #(fnum_str, fnum_str) return fnum_str
1
N/A
local
($x,$y) = (&
'fnorm($_[$[]),&'
fnorm
($_[$[+
1
]));
1
N/A
if
($x
eq
'NaN'
|| $y
eq
'NaN'
) {
1
N/A
'NaN'
;
1
N/A
}
else
{
1
N/A
local
($
xm
,$
xe
) =
split
(
'E'
,$x);
1
N/A
local
($
ym
,$
ye
) =
split
(
'E'
,$y);
1
N/A
($
xm
,$
xe
,$
ym
,$
ye
) = ($
ym
,$
ye
,$
xm
,$
xe
)
if
($
xe
< $
ye
);
1
N/A
&
norm
(&
'badd($ym,$xm.('
0
' x ($xe-$ye))),$ye);
1
N/A
}
1
N/A
}
1
N/A
1
N/A
# subtraction
1
N/A
sub
main
'fsub { #(fnum_str, fnum_str) return fnum_str
1
N/A
&
'fadd($_[$[],&'
fneg
($_[$[+
1
]));
1
N/A
}
1
N/A
1
N/A
# division
1
N/A
# args are dividend, divisor, scale (optional)
1
N/A
# result has at most max(scale, length(dividend), length(divisor)) digits
1
N/A
sub
main
'fdiv #(fnum_str, fnum_str[,scale]) return fnum_str
1
N/A
{
1
N/A
local
($x,$y,$
scale
) = (&
'fnorm($_[$[]),&'
fnorm
($_[$[+
1
]),$_[$[+
2
]);
1
N/A
if
($x
eq
'NaN'
|| $y
eq
'NaN'
|| $y
eq
'+0E+0'
) {
1
N/A
'NaN'
;
1
N/A
}
else
{
1
N/A
local
($
xm
,$
xe
) =
split
(
'E'
,$x);
1
N/A
local
($
ym
,$
ye
) =
split
(
'E'
,$y);
1
N/A
$
scale
= $
div_scale
if
(!$
scale
);
1
N/A
$
scale
=
length
($
xm
)-
1
if
(
length
($
xm
)-
1
> $
scale
);
1
N/A
$
scale
=
length
($
ym
)-
1
if
(
length
($
ym
)-
1
> $
scale
);
1
N/A
$
scale
= $
scale
+
length
($
ym
) -
length
($
xm
);
1
N/A
&
norm
(&
round
(&
'bdiv($xm.('
0
' x $scale),$ym),&'
babs
($
ym
)),
1
N/A
$
xe
-$
ye
-$
scale
);
1
N/A
}
1
N/A
}
1
N/A
1
N/A
# round int $q based on fraction $r/$base using $rnd_mode
1
N/A
sub
round
{
#(int_str, int_str, int_str) return int_str
1
N/A
local
($q,$r,$
base
) = @_;
1
N/A
if
($q
eq
'NaN'
|| $r
eq
'NaN'
) {
1
N/A
'NaN'
;
1
N/A
}
elsif
($
rnd_mode
eq
'trunc'
) {
1
N/A
$q;
# just truncate
1
N/A
}
else
{
1
N/A
local
($
cmp
) = &
'bcmp(&'
bmul
($r,
'+2'
),$
base
);
1
N/A
if
( $
cmp
<
0
||
1
N/A
($
cmp
==
0
&&
1
N/A
( $
rnd_mode
eq
'zero'
||
1
N/A
($
rnd_mode
eq
'-inf'
&& (
substr
($q,$[,
1
)
eq
'+'
)) ||
1
N/A
($
rnd_mode
eq
'+inf'
&& (
substr
($q,$[,
1
)
eq
'-'
)) ||
1
N/A
($
rnd_mode
eq
'even'
&& $q =~ /[
24680
]$/) ||
1
N/A
($
rnd_mode
eq
'odd'
&& $q =~ /[
13579
]$/) )) ) {
1
N/A
$q;
# round down
1
N/A
}
else
{
1
N/A
&
'badd($q, ((substr($q,$[,1) eq '
-
') ? '
-
1
' : '
+
1
'));
1
N/A
# round up
1
N/A
}
1
N/A
}
1
N/A
}
1
N/A
1
N/A
# round the mantissa of $x to $scale digits
1
N/A
sub
main
'fround { #(fnum_str, scale) return fnum_str
1
N/A
local
($x,$
scale
) = (&
'fnorm($_[$[]),$_[$[+1]);
1
N/A
if
($x
eq
'NaN'
|| $
scale
<=
0
) {
1
N/A
$x;
1
N/A
}
else
{
1
N/A
local
($
xm
,$
xe
) =
split
(
'E'
,$x);
1
N/A
if
(
length
($
xm
)-
1
<= $
scale
) {
1
N/A
$x;
1
N/A
}
else
{
1
N/A
&
norm
(&
round
(
substr
($
xm
,$[,$
scale
+
1
),
1
N/A
"+0"
.
substr
($
xm
,$[+$
scale
+
1
,
1
),
"+10"
),
1
N/A
$
xe
+
length
($
xm
)-$
scale
-
1
);
1
N/A
}
1
N/A
}
1
N/A
}
1
N/A
1
N/A
# round $x at the 10 to the $scale digit place
1
N/A
sub
main
'ffround { #(fnum_str, scale) return fnum_str
1
N/A
local
($x,$
scale
) = (&
'fnorm($_[$[]),$_[$[+1]);
1
N/A
if
($x
eq
'NaN'
) {
1
N/A
'NaN'
;
1
N/A
}
else
{
1
N/A
local
($
xm
,$
xe
) =
split
(
'E'
,$x);
1
N/A
if
($
xe
>= $
scale
) {
1
N/A
$x;
1
N/A
}
else
{
1
N/A
$
xe
=
length
($
xm
)+$
xe
-$
scale
;
1
N/A
if
($
xe
<
1
) {
1
N/A
'+0E+0'
;
1
N/A
}
elsif
($
xe
==
1
) {
1
N/A
# The first substr preserves the sign, which means that
1
N/A
# we'll pass a non-normalized "-0" to &round when rounding
1
N/A
# -0.006 (for example), purely so that &round won't lose
1
N/A
# the sign.
1
N/A
&
norm
(&
round
(
substr
($
xm
,$[,
1
).
'0'
,
1
N/A
"+0"
.
substr
($
xm
,$[+
1
,
1
),
"+10"
), $
scale
);
1
N/A
}
else
{
1
N/A
&
norm
(&
round
(
substr
($
xm
,$[,$
xe
),
1
N/A
"+0"
.
substr
($
xm
,$[+$
xe
,
1
),
"+10"
), $
scale
);
1
N/A
}
1
N/A
}
1
N/A
}
1
N/A
}
1
N/A
1
N/A
# compare 2 values returns one of undef, <0, =0, >0
1
N/A
# returns undef if either or both input value are not numbers
1
N/A
sub
main
'fcmp #(fnum_str, fnum_str) return cond_code
1
N/A
{
1
N/A
local
($x, $y) = (&
'fnorm($_[$[]),&'
fnorm
($_[$[+
1
]));
1
N/A
if
($x
eq
"NaN"
|| $y
eq
"NaN"
) {
1
N/A
undef
;
1
N/A
}
else
{
1
N/A
ord
($y) <=>
ord
($x)
1
N/A
||
1
N/A
(
local
($
xm
,$
xe
,$
ym
,$
ye
) =
split
(
'E'
, $x.
"E$y"
),
1
N/A
(($
xe
<=> $
ye
) * (
substr
($x,$[,
1
).
'1'
)
1
N/A
|| &
bigint
'cmp($xm,$ym))
1
N/A
);
1
N/A
}
1
N/A
}
1
N/A
1
N/A
# square root by Newtons method.
1
N/A
sub
main
'fsqrt { #(fnum_str[, scale]) return fnum_str
1
N/A
local
($x, $
scale
) = (&
'fnorm($_[$[]), $_[$[+1]);
1
N/A
if
($x
eq
'NaN'
|| $x =~ /^-/) {
1
N/A
'NaN'
;
1
N/A
}
elsif
($x
eq
'+0E+0'
) {
1
N/A
'+0E+0'
;
1
N/A
}
else
{
1
N/A
local
($
xm
, $
xe
) =
split
(
'E'
,$x);
1
N/A
$
scale
= $
div_scale
if
(!$
scale
);
1
N/A
$
scale
=
length
($
xm
)-
1
if
($
scale
<
length
($
xm
)-
1
);
1
N/A
local
($
gs
, $
guess
) = (
1
,
sprintf
(
"1E%+d"
, (
length
($
xm
)+$
xe
-
1
)/
2
));
1
N/A
while
($
gs
<
2
*$
scale
) {
1
N/A
$
guess
= &
'fmul(&'
fadd
($
guess
,&
'fdiv($x,$guess,$gs*2)),".5");
1
N/A
$
gs
*=
2
;
1
N/A
}
1
N/A
&
'fround($guess, $scale);
1
N/A
}
1
N/A
}
1
N/A
1
N/A
1
;