4 # This library is no longer being maintained, and is included for backward
5 # compatibility with Perl 4 programs which may require it.
7 # In particular, this should not be used as an example of modern Perl
8 # programming techniques.
10 # Suggested alternative: Math::BigFloat
12 # Arbitrary length float math package
17 # canonical strings have the form /[+-]\d+E[+-]\d+/
18 # Input values can have embedded whitespace
20 # 'NaN' An input parameter was "Not a Number" or
21 # divide by zero or sqrt of negative number
22 # Division is computed to
23 # max($div_scale,length(dividend)+length(divisor))
25 # Also used for default sqrt scale
29 # Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
35 # fadd(NSTR, NSTR) return NSTR addition
36 # fsub(NSTR, NSTR) return NSTR subtraction
37 # fmul(NSTR, NSTR) return NSTR multiplication
38 # fdiv(NSTR, NSTR[,SCALE]) returns NSTR division to SCALE places
39 # fneg(NSTR) return NSTR negation
40 # fabs(NSTR) return NSTR absolute value
41 # fcmp(NSTR,NSTR) return CODE compare undef,<0,=0,>0
42 # fround(NSTR, SCALE) return NSTR round to SCALE digits
43 # ffround(NSTR, SCALE) return NSTR round at SCALEth place
44 # fnorm(NSTR) return (NSTR) normalize
45 # fsqrt(NSTR[, SCALE]) return NSTR sqrt to SCALE places
47 # Convert a number to canonical string form.
48 # Takes something that looks like a number and converts it to
49 # the form /^[+-]\d+E[+-]\d+$/.
50 sub main'fnorm { #(string) return fnum_str
52 s/\s+//g; # strip white space
53 if (/^([+-]?)(\d*)(\.(\d*))?([Ee]([+-]?\d+))?$/
54 && ($2 ne '' || defined($4))) {
55 my $x = defined($4) ? $4 : '';
56 &norm(($1 ? "$1$2$x" : "+$2$x"), (($x ne '') ? $6-length($x) : $6));
62 # normalize number -- for internal use
63 sub norm { #(mantissa, exponent) return fnum_str
68 s/^([+-])0+/$1/; # strip leading zeros
69 if (length($_) == 1) {
72 $exp += length($1) if (s/(0+)$//); # strip trailing zeros
73 sprintf("%sE%+ld", $_, $exp);
79 sub main'fneg { #(fnum_str) return fnum_str
80 local($_) = &'fnorm($_[$[]);
81 vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0E+0'; # flip sign
87 sub main'fabs { #(fnum_str) return fnum_str
88 local($_) = &'fnorm($_[$[]);
94 sub main'fmul { #(fnum_str, fnum_str) return fnum_str
95 local($x,$y) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]));
96 if ($x eq 'NaN' || $y eq 'NaN') {
99 local($xm,$xe) = split('E',$x);
100 local($ym,$ye) = split('E',$y);
101 &norm(&'bmul($xm,$ym),$xe+$ye);
106 sub main'fadd { #(fnum_str, fnum_str) return fnum_str
107 local($x,$y) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]));
108 if ($x eq 'NaN' || $y eq 'NaN') {
111 local($xm,$xe) = split('E',$x);
112 local($ym,$ye) = split('E',$y);
113 ($xm,$xe,$ym,$ye) = ($ym,$ye,$xm,$xe) if ($xe < $ye);
114 &norm(&'badd($ym,$xm.('0' x ($xe-$ye))),$ye);
119 sub main'fsub { #(fnum_str, fnum_str) return fnum_str
120 &'fadd($_[$[],&'fneg($_[$[+1]));
124 # args are dividend, divisor, scale (optional)
125 # result has at most max(scale, length(dividend), length(divisor)) digits
126 sub main'fdiv #(fnum_str, fnum_str[,scale]) return fnum_str
128 local($x,$y,$scale) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]),$_[$[+2]);
129 if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0E+0') {
132 local($xm,$xe) = split('E',$x);
133 local($ym,$ye) = split('E',$y);
134 $scale = $div_scale if (!$scale);
135 $scale = length($xm)-1 if (length($xm)-1 > $scale);
136 $scale = length($ym)-1 if (length($ym)-1 > $scale);
137 $scale = $scale + length($ym) - length($xm);
138 &norm(&round(&'bdiv($xm.('0' x $scale),$ym),&'babs($ym)),
143 # round int $q based on fraction $r/$base using $rnd_mode
144 sub round { #(int_str, int_str, int_str) return int_str
145 local($q,$r,$base) = @_;
146 if ($q eq 'NaN' || $r eq 'NaN') {
148 } elsif ($rnd_mode eq 'trunc') {
151 local($cmp) = &'bcmp(&'bmul($r,'+2'),$base);
154 ( $rnd_mode eq 'zero' ||
155 ($rnd_mode eq '-inf' && (substr($q,$[,1) eq '+')) ||
156 ($rnd_mode eq '+inf' && (substr($q,$[,1) eq '-')) ||
157 ($rnd_mode eq 'even' && $q =~ /[24680]$/) ||
158 ($rnd_mode eq 'odd' && $q =~ /[13579]$/) )) ) {
161 &'badd($q, ((substr($q,$[,1) eq '-') ? '-1' : '+1'));
167 # round the mantissa of $x to $scale digits
168 sub main'fround { #(fnum_str, scale) return fnum_str
169 local($x,$scale) = (&'fnorm($_[$[]),$_[$[+1]);
170 if ($x eq 'NaN' || $scale <= 0) {
173 local($xm,$xe) = split('E',$x);
174 if (length($xm)-1 <= $scale) {
177 &norm(&round(substr($xm,$[,$scale+1),
178 "+0".substr($xm,$[+$scale+1,1),"+10"),
179 $xe+length($xm)-$scale-1);
184 # round $x at the 10 to the $scale digit place
185 sub main'ffround { #(fnum_str, scale) return fnum_str
186 local($x,$scale) = (&'fnorm($_[$[]),$_[$[+1]);
190 local($xm,$xe) = split('E',$x);
194 $xe = length($xm)+$xe-$scale;
198 # The first substr preserves the sign, which means that
199 # we'll pass a non-normalized "-0" to &round when rounding
200 # -0.006 (for example), purely so that &round won't lose
202 &norm(&round(substr($xm,$[,1).'0',
203 "+0".substr($xm,$[+1,1),"+10"), $scale);
205 &norm(&round(substr($xm,$[,$xe),
206 "+0".substr($xm,$[+$xe,1),"+10"), $scale);
212 # compare 2 values returns one of undef, <0, =0, >0
213 # returns undef if either or both input value are not numbers
214 sub main'fcmp #(fnum_str, fnum_str) return cond_code
216 local($x, $y) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]));
217 if ($x eq "NaN" || $y eq "NaN") {
222 ( local($xm,$xe,$ym,$ye) = split('E', $x."E$y"),
223 (($xe <=> $ye) * (substr($x,$[,1).'1')
224 || &bigint'cmp($xm,$ym))
229 # square root by Newtons method.
230 sub main'fsqrt { #(fnum_str[, scale]) return fnum_str
231 local($x, $scale) = (&'fnorm($_[$[]), $_[$[+1]);
232 if ($x eq 'NaN' || $x =~ /^-/) {
234 } elsif ($x eq '+0E+0') {
237 local($xm, $xe) = split('E',$x);
238 $scale = $div_scale if (!$scale);
239 $scale = length($xm)-1 if ($scale < length($xm)-1);
240 local($gs, $guess) = (1, sprintf("1E%+d", (length($xm)+$xe-1)/2));
241 while ($gs < 2*$scale) {
242 $guess = &'fmul(&'fadd($guess,&'fdiv($x,$guess,$gs*2)),".5");
245 &'fround($guess, $scale);