1 warn "Legacy library @{[(caller(0))[6]]} will be removed from the Perl core distribution in the next major release. Please install it from the CPAN distribution Perl4::CoreLibs. It is being used at @{[(caller)[1]]}, line @{[(caller)[2]]}.\n";
6 # This library is no longer being maintained, and is included for backward
7 # compatibility with Perl 4 programs which may require it.
8 # This legacy library is deprecated and will be removed in a future
11 # In particular, this should not be used as an example of modern Perl
12 # programming techniques.
14 # Suggested alternative: Math::BigFloat
16 # Arbitrary length float math package
21 # canonical strings have the form /[+-]\d+E[+-]\d+/
22 # Input values can have embedded whitespace
24 # 'NaN' An input parameter was "Not a Number" or
25 # divide by zero or sqrt of negative number
26 # Division is computed to
27 # max($div_scale,length(dividend)+length(divisor))
29 # Also used for default sqrt scale
33 # Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
39 # fadd(NSTR, NSTR) return NSTR addition
40 # fsub(NSTR, NSTR) return NSTR subtraction
41 # fmul(NSTR, NSTR) return NSTR multiplication
42 # fdiv(NSTR, NSTR[,SCALE]) returns NSTR division to SCALE places
43 # fneg(NSTR) return NSTR negation
44 # fabs(NSTR) return NSTR absolute value
45 # fcmp(NSTR,NSTR) return CODE compare undef,<0,=0,>0
46 # fround(NSTR, SCALE) return NSTR round to SCALE digits
47 # ffround(NSTR, SCALE) return NSTR round at SCALEth place
48 # fnorm(NSTR) return (NSTR) normalize
49 # fsqrt(NSTR[, SCALE]) return NSTR sqrt to SCALE places
51 # Convert a number to canonical string form.
52 # Takes something that looks like a number and converts it to
53 # the form /^[+-]\d+E[+-]\d+$/.
54 sub main'fnorm { #(string) return fnum_str
56 s/\s+//g; # strip white space
57 if (/^([+-]?)(\d*)(\.(\d*))?([Ee]([+-]?\d+))?$/
58 && ($2 ne '' || defined($4))) {
59 my $x = defined($4) ? $4 : '';
60 &norm(($1 ? "$1$2$x" : "+$2$x"), (($x ne '') ? $6-length($x) : $6));
66 # normalize number -- for internal use
67 sub norm { #(mantissa, exponent) return fnum_str
72 s/^([+-])0+/$1/; # strip leading zeros
73 if (length($_) == 1) {
76 $exp += length($1) if (s/(0+)$//); # strip trailing zeros
77 sprintf("%sE%+ld", $_, $exp);
83 sub main'fneg { #(fnum_str) return fnum_str
84 local($_) = &'fnorm($_[0]);
85 vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0E+0'; # flip sign
86 if ( ord("\t") == 9 ) { # ascii
89 else { # ebcdic character set
96 sub main'fabs { #(fnum_str) return fnum_str
97 local($_) = &'fnorm($_[0]);
103 sub main'fmul { #(fnum_str, fnum_str) return fnum_str
104 local($x,$y) = (&'fnorm($_[0]),&'fnorm($_[1]));
105 if ($x eq 'NaN' || $y eq 'NaN') {
108 local($xm,$xe) = split('E',$x);
109 local($ym,$ye) = split('E',$y);
110 &norm(&'bmul($xm,$ym),$xe+$ye);
115 sub main'fadd { #(fnum_str, fnum_str) return fnum_str
116 local($x,$y) = (&'fnorm($_[0]),&'fnorm($_[1]));
117 if ($x eq 'NaN' || $y eq 'NaN') {
120 local($xm,$xe) = split('E',$x);
121 local($ym,$ye) = split('E',$y);
122 ($xm,$xe,$ym,$ye) = ($ym,$ye,$xm,$xe) if ($xe < $ye);
123 &norm(&'badd($ym,$xm.('0' x ($xe-$ye))),$ye);
128 sub main'fsub { #(fnum_str, fnum_str) return fnum_str
129 &'fadd($_[0],&'fneg($_[1]));
133 # args are dividend, divisor, scale (optional)
134 # result has at most max(scale, length(dividend), length(divisor)) digits
135 sub main'fdiv #(fnum_str, fnum_str[,scale]) return fnum_str
137 local($x,$y,$scale) = (&'fnorm($_[0]),&'fnorm($_[1]),$_[2]);
138 if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0E+0') {
141 local($xm,$xe) = split('E',$x);
142 local($ym,$ye) = split('E',$y);
143 $scale = $div_scale if (!$scale);
144 $scale = length($xm)-1 if (length($xm)-1 > $scale);
145 $scale = length($ym)-1 if (length($ym)-1 > $scale);
146 $scale = $scale + length($ym) - length($xm);
147 &norm(&round(&'bdiv($xm.('0' x $scale),$ym),&'babs($ym)),
152 # round int $q based on fraction $r/$base using $rnd_mode
153 sub round { #(int_str, int_str, int_str) return int_str
154 local($q,$r,$base) = @_;
155 if ($q eq 'NaN' || $r eq 'NaN') {
157 } elsif ($rnd_mode eq 'trunc') {
160 local($cmp) = &'bcmp(&'bmul($r,'+2'),$base);
163 ( $rnd_mode eq 'zero' ||
164 ($rnd_mode eq '-inf' && (substr($q,0,1) eq '+')) ||
165 ($rnd_mode eq '+inf' && (substr($q,0,1) eq '-')) ||
166 ($rnd_mode eq 'even' && $q =~ /[24680]$/) ||
167 ($rnd_mode eq 'odd' && $q =~ /[13579]$/) )) ) {
170 &'badd($q, ((substr($q,0,1) eq '-') ? '-1' : '+1'));
176 # round the mantissa of $x to $scale digits
177 sub main'fround { #(fnum_str, scale) return fnum_str
178 local($x,$scale) = (&'fnorm($_[0]),$_[1]);
179 if ($x eq 'NaN' || $scale <= 0) {
182 local($xm,$xe) = split('E',$x);
183 if (length($xm)-1 <= $scale) {
186 &norm(&round(substr($xm,0,$scale+1),
187 "+0".substr($xm,$scale+1,1),"+10"),
188 $xe+length($xm)-$scale-1);
193 # round $x at the 10 to the $scale digit place
194 sub main'ffround { #(fnum_str, scale) return fnum_str
195 local($x,$scale) = (&'fnorm($_[0]),$_[1]);
199 local($xm,$xe) = split('E',$x);
203 $xe = length($xm)+$xe-$scale;
207 # The first substr preserves the sign, which means that
208 # we'll pass a non-normalized "-0" to &round when rounding
209 # -0.006 (for example), purely so that &round won't lose
211 &norm(&round(substr($xm,0,1).'0',
212 "+0".substr($xm,1,1),"+10"), $scale);
214 &norm(&round(substr($xm,0,$xe),
215 "+0".substr($xm,$xe,1),"+10"), $scale);
221 # compare 2 values returns one of undef, <0, =0, >0
222 # returns undef if either or both input value are not numbers
223 sub main'fcmp #(fnum_str, fnum_str) return cond_code
225 local($x, $y) = (&'fnorm($_[0]),&'fnorm($_[1]));
226 if ($x eq "NaN" || $y eq "NaN") {
231 ( local($xm,$xe,$ym,$ye) = split('E', $x."E$y"),
232 (($xe <=> $ye) * (substr($x,0,1).'1')
233 || &bigint'cmp($xm,$ym))
238 # square root by Newtons method.
239 sub main'fsqrt { #(fnum_str[, scale]) return fnum_str
240 local($x, $scale) = (&'fnorm($_[0]), $_[1]);
241 if ($x eq 'NaN' || $x =~ /^-/) {
243 } elsif ($x eq '+0E+0') {
246 local($xm, $xe) = split('E',$x);
247 $scale = $div_scale if (!$scale);
248 $scale = length($xm)-1 if ($scale < length($xm)-1);
249 local($gs, $guess) = (1, sprintf("1E%+d", (length($xm)+$xe-1)/2));
250 while ($gs < 2*$scale) {
251 $guess = &'fmul(&'fadd($guess,&'fdiv($x,$guess,$gs*2)),".5");
254 &'fround($guess, $scale);