3 # Arbitrary length float math package
8 # canonical strings have the form /[+-]\d+E[+-]\d+/
9 # Input values can have inbedded whitespace
11 # 'NaN' An input parameter was "Not a Number" or
12 # divide by zero or sqrt of negative number
13 # Division is computed to
14 # max($div_scale,length(dividend)+length(divisor))
16 # Also used for default sqrt scale
20 # Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
26 # fadd(NSTR, NSTR) return NSTR addition
27 # fsub(NSTR, NSTR) return NSTR subtraction
28 # fmul(NSTR, NSTR) return NSTR multiplication
29 # fdiv(NSTR, NSTR[,SCALE]) returns NSTR division to SCALE places
30 # fneg(NSTR) return NSTR negation
31 # fabs(NSTR) return NSTR absolute value
32 # fcmp(NSTR,NSTR) return CODE compare undef,<0,=0,>0
33 # fround(NSTR, SCALE) return NSTR round to SCALE digits
34 # ffround(NSTR, SCALE) return NSTR round at SCALEth place
35 # fnorm(NSTR) return (NSTR) normalize
36 # fsqrt(NSTR[, SCALE]) return NSTR sqrt to SCALE places
38 # Convert a number to canonical string form.
39 # Takes something that looks like a number and converts it to
40 # the form /^[+-]\d+E[+-]\d+$/.
41 sub main'fnorm { #(string) return fnum_str
43 s/\s+//g; # strip white space
44 if (/^([+-]?)(\d*)(\.(\d*))?([Ee]([+-]?\d+))?$/ && "$2$4" ne '') {
45 &norm(($1 ? "$1$2$4" : "+$2$4"),(($4 ne '') ? $6-length($4) : $6));
51 # normalize number -- for internal use
52 sub norm { #(mantissa, exponent) return fnum_str
57 s/^([+-])0+/$1/; # strip leading zeros
58 if (length($_) == 1) {
61 $exp += length($1) if (s/(0+)$//); # strip trailing zeros
62 sprintf("%sE%+ld", $_, $exp);
68 sub main'fneg { #(fnum_str) return fnum_str
69 local($_) = &'fnorm($_[$[]);
70 vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0E+0'; # flip sign
76 sub main'fabs { #(fnum_str) return fnum_str
77 local($_) = &'fnorm($_[$[]);
83 sub main'fmul { #(fnum_str, fnum_str) return fnum_str
84 local($x,$y) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]));
85 if ($x eq 'NaN' || $y eq 'NaN') {
88 local($xm,$xe) = split('E',$x);
89 local($ym,$ye) = split('E',$y);
90 &norm(&'bmul($xm,$ym),$xe+$ye);
95 sub main'fadd { #(fnum_str, fnum_str) return fnum_str
96 local($x,$y) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]));
97 if ($x eq 'NaN' || $y eq 'NaN') {
100 local($xm,$xe) = split('E',$x);
101 local($ym,$ye) = split('E',$y);
102 ($xm,$xe,$ym,$ye) = ($ym,$ye,$xm,$xe) if ($xe < $ye);
103 &norm(&'badd($ym,$xm.('0' x ($xe-$ye))),$ye);
108 sub main'fsub { #(fnum_str, fnum_str) return fnum_str
109 &'fadd($_[$[],&'fneg($_[$[+1]));
113 # args are dividend, divisor, scale (optional)
114 # result has at most max(scale, length(dividend), length(divisor)) digits
115 sub main'fdiv #(fnum_str, fnum_str[,scale]) return fnum_str
117 local($x,$y,$scale) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]),$_[$[+2]);
118 if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0E+0') {
121 local($xm,$xe) = split('E',$x);
122 local($ym,$ye) = split('E',$y);
123 $scale = $div_scale if (!$scale);
124 $scale = length($xm)-1 if (length($xm)-1 > $scale);
125 $scale = length($ym)-1 if (length($ym)-1 > $scale);
126 $scale = $scale + length($ym) - length($xm);
127 &norm(&round(&'bdiv($xm.('0' x $scale),$ym),$ym),
132 # round int $q based on fraction $r/$base using $rnd_mode
133 sub round { #(int_str, int_str, int_str) return int_str
134 local($q,$r,$base) = @_;
135 if ($q eq 'NaN' || $r eq 'NaN') {
137 } elsif ($rnd_mode eq 'trunc') {
140 local($cmp) = &'bcmp(&'bmul($r,'+2'),$base);
143 ( $rnd_mode eq 'zero' ||
144 ($rnd_mode eq '-inf' && (substr($q,$[,1) eq '+')) ||
145 ($rnd_mode eq '+inf' && (substr($q,$[,1) eq '-')) ||
146 ($rnd_mode eq 'even' && $q =~ /[24680]$/) ||
147 ($rnd_mode eq 'odd' && $q =~ /[13579]$/) )) ) {
150 &'badd($q, ((substr($q,$[,1) eq '-') ? '-1' : '+1'));
156 # round the mantissa of $x to $scale digits
157 sub main'fround { #(fnum_str, scale) return fnum_str
158 local($x,$scale) = (&'fnorm($_[$[]),$_[$[+1]);
159 if ($x eq 'NaN' || $scale <= 0) {
162 local($xm,$xe) = split('E',$x);
163 if (length($xm)-1 <= $scale) {
166 &norm(&round(substr($xm,$[,$scale+1),
167 "+0".substr($xm,$[+$scale+1,1),"+10"),
168 $xe+length($xm)-$scale-1);
173 # round $x at the 10 to the $scale digit place
174 sub main'ffround { #(fnum_str, scale) return fnum_str
175 local($x,$scale) = (&'fnorm($_[$[]),$_[$[+1]);
179 local($xm,$xe) = split('E',$x);
183 $xe = length($xm)+$xe-$scale;
187 &norm(&round('+0',"+0".substr($xm,$[+1,1),"+10"), $scale);
189 &norm(&round(substr($xm,$[,$xe),
190 "+0".substr($xm,$[+$xe,1),"+10"), $scale);
196 # compare 2 values returns one of undef, <0, =0, >0
197 # returns undef if either or both input value are not numbers
198 sub main'fcmp #(fnum_str, fnum_str) return cond_code
200 local($x, $y) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]));
201 if ($x eq "NaN" || $y eq "NaN") {
206 ( local($xm,$xe,$ym,$ye) = split('E', $x."E$y"),
207 (($xe <=> $ye) * (substr($x,$[,1).'1')
208 || &bigint'cmp($xm,$ym))
213 # square root by Newtons method.
214 sub main'fsqrt { #(fnum_str[, scale]) return fnum_str
215 local($x, $scale) = (&'fnorm($_[$[]), $_[$[+1]);
216 if ($x eq 'NaN' || $x =~ /^-/) {
218 } elsif ($x eq '+0E+0') {
221 local($xm, $xe) = split('E',$x);
222 $scale = $div_scale if (!$scale);
223 $scale = length($xm)-1 if ($scale < length($xm)-1);
224 local($gs, $guess) = (1, sprintf("1E%+d", (length($xm)+$xe-1)/2));
225 while ($gs < 2*$scale) {
226 $guess = &'fmul(&'fadd($guess,&'fdiv($x,$guess,$gs*2)),".5");
229 &'fround($guess, $scale);