RE: Combining UTF-16 output with :crlf is awkward
[perl.git] / lib / bigfloat.pl
1 package bigfloat;
2 require "bigint.pl";
3 #
4 # This library is no longer being maintained, and is included for backward
5 # compatibility with Perl 4 programs which may require it.
6 #
7 # In particular, this should not be used as an example of modern Perl
8 # programming techniques.
9 #
10 # Suggested alternative: Math::BigFloat
11 #
12 # Arbitrary length float math package
13 #
14 # by Mark Biggar
15 #
16 # number format
17 #   canonical strings have the form /[+-]\d+E[+-]\d+/
18 #   Input values can have embedded whitespace
19 # Error returns
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)) 
24 #   digits by default.
25 # Also used for default sqrt scale
26
27 $div_scale = 40;
28
29 # Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
30
31 $rnd_mode = 'even';
32
33 #   bigfloat routines
34 #
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
46 \f
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
51     local($_) = @_;
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));
57     } else {
58         'NaN';
59     }
60 }
61
62 # normalize number -- for internal use
63 sub norm { #(mantissa, exponent) return fnum_str
64     local($_, $exp) = @_;
65     if ($_ eq 'NaN') {
66         'NaN';
67     } else {
68         s/^([+-])0+/$1/;                        # strip leading zeros
69         if (length($_) == 1) {
70             '+0E+0';
71         } else {
72             $exp += length($1) if (s/(0+)$//);  # strip trailing zeros
73             sprintf("%sE%+ld", $_, $exp);
74         }
75     }
76 }
77
78 # negation
79 sub main'fneg { #(fnum_str) return fnum_str
80     local($_) = &'fnorm($_[$[]);
81     vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0E+0'; # flip sign
82     if ( ord("\t") == 9 ) { # ascii
83         s/^H/N/;
84     }
85     else { # ebcdic character set
86         s/\373/N/;
87     }
88     $_;
89 }
90
91 # absolute value
92 sub main'fabs { #(fnum_str) return fnum_str
93     local($_) = &'fnorm($_[$[]);
94     s/^-/+/;                                   # mash sign
95     $_;
96 }
97
98 # multiplication
99 sub main'fmul { #(fnum_str, fnum_str) return fnum_str
100     local($x,$y) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]));
101     if ($x eq 'NaN' || $y eq 'NaN') {
102         'NaN';
103     } else {
104         local($xm,$xe) = split('E',$x);
105         local($ym,$ye) = split('E',$y);
106         &norm(&'bmul($xm,$ym),$xe+$ye);
107     }
108 }
109 \f
110 # addition
111 sub main'fadd { #(fnum_str, fnum_str) return fnum_str
112     local($x,$y) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]));
113     if ($x eq 'NaN' || $y eq 'NaN') {
114         'NaN';
115     } else {
116         local($xm,$xe) = split('E',$x);
117         local($ym,$ye) = split('E',$y);
118         ($xm,$xe,$ym,$ye) = ($ym,$ye,$xm,$xe) if ($xe < $ye);
119         &norm(&'badd($ym,$xm.('0' x ($xe-$ye))),$ye);
120     }
121 }
122
123 # subtraction
124 sub main'fsub { #(fnum_str, fnum_str) return fnum_str
125     &'fadd($_[$[],&'fneg($_[$[+1]));    
126 }
127
128 # division
129 #   args are dividend, divisor, scale (optional)
130 #   result has at most max(scale, length(dividend), length(divisor)) digits
131 sub main'fdiv #(fnum_str, fnum_str[,scale]) return fnum_str
132 {
133     local($x,$y,$scale) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]),$_[$[+2]);
134     if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0E+0') {
135         'NaN';
136     } else {
137         local($xm,$xe) = split('E',$x);
138         local($ym,$ye) = split('E',$y);
139         $scale = $div_scale if (!$scale);
140         $scale = length($xm)-1 if (length($xm)-1 > $scale);
141         $scale = length($ym)-1 if (length($ym)-1 > $scale);
142         $scale = $scale + length($ym) - length($xm);
143         &norm(&round(&'bdiv($xm.('0' x $scale),$ym),&'babs($ym)),
144             $xe-$ye-$scale);
145     }
146 }
147 \f
148 # round int $q based on fraction $r/$base using $rnd_mode
149 sub round { #(int_str, int_str, int_str) return int_str
150     local($q,$r,$base) = @_;
151     if ($q eq 'NaN' || $r eq 'NaN') {
152         'NaN';
153     } elsif ($rnd_mode eq 'trunc') {
154         $q;                         # just truncate
155     } else {
156         local($cmp) = &'bcmp(&'bmul($r,'+2'),$base);
157         if ( $cmp < 0 ||
158                  ($cmp == 0 &&
159                   ( $rnd_mode eq 'zero'                             ||
160                    ($rnd_mode eq '-inf' && (substr($q,$[,1) eq '+')) ||
161                    ($rnd_mode eq '+inf' && (substr($q,$[,1) eq '-')) ||
162                    ($rnd_mode eq 'even' && $q =~ /[24680]$/)        ||
163                    ($rnd_mode eq 'odd'  && $q =~ /[13579]$/)        )) ) {
164             $q;                     # round down
165         } else {
166             &'badd($q, ((substr($q,$[,1) eq '-') ? '-1' : '+1'));
167                                     # round up
168         }
169     }
170 }
171
172 # round the mantissa of $x to $scale digits
173 sub main'fround { #(fnum_str, scale) return fnum_str
174     local($x,$scale) = (&'fnorm($_[$[]),$_[$[+1]);
175     if ($x eq 'NaN' || $scale <= 0) {
176         $x;
177     } else {
178         local($xm,$xe) = split('E',$x);
179         if (length($xm)-1 <= $scale) {
180             $x;
181         } else {
182             &norm(&round(substr($xm,$[,$scale+1),
183                          "+0".substr($xm,$[+$scale+1,1),"+10"),
184                   $xe+length($xm)-$scale-1);
185         }
186     }
187 }
188 \f
189 # round $x at the 10 to the $scale digit place
190 sub main'ffround { #(fnum_str, scale) return fnum_str
191     local($x,$scale) = (&'fnorm($_[$[]),$_[$[+1]);
192     if ($x eq 'NaN') {
193         'NaN';
194     } else {
195         local($xm,$xe) = split('E',$x);
196         if ($xe >= $scale) {
197             $x;
198         } else {
199             $xe = length($xm)+$xe-$scale;
200             if ($xe < 1) {
201                 '+0E+0';
202             } elsif ($xe == 1) {
203                 # The first substr preserves the sign, which means that
204                 # we'll pass a non-normalized "-0" to &round when rounding
205                 # -0.006 (for example), purely so that &round won't lose
206                 # the sign.
207                 &norm(&round(substr($xm,$[,1).'0',
208                       "+0".substr($xm,$[+1,1),"+10"), $scale);
209             } else {
210                 &norm(&round(substr($xm,$[,$xe),
211                       "+0".substr($xm,$[+$xe,1),"+10"), $scale);
212             }
213         }
214     }
215 }
216     
217 # compare 2 values returns one of undef, <0, =0, >0
218 #   returns undef if either or both input value are not numbers
219 sub main'fcmp #(fnum_str, fnum_str) return cond_code
220 {
221     local($x, $y) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]));
222     if ($x eq "NaN" || $y eq "NaN") {
223         undef;
224     } else {
225         ord($y) <=> ord($x)
226         ||
227         (  local($xm,$xe,$ym,$ye) = split('E', $x."E$y"),
228              (($xe <=> $ye) * (substr($x,$[,1).'1')
229              || &bigint'cmp($xm,$ym))
230         );
231     }
232 }
233 \f
234 # square root by Newtons method.
235 sub main'fsqrt { #(fnum_str[, scale]) return fnum_str
236     local($x, $scale) = (&'fnorm($_[$[]), $_[$[+1]);
237     if ($x eq 'NaN' || $x =~ /^-/) {
238         'NaN';
239     } elsif ($x eq '+0E+0') {
240         '+0E+0';
241     } else {
242         local($xm, $xe) = split('E',$x);
243         $scale = $div_scale if (!$scale);
244         $scale = length($xm)-1 if ($scale < length($xm)-1);
245         local($gs, $guess) = (1, sprintf("1E%+d", (length($xm)+$xe-1)/2));
246         while ($gs < 2*$scale) {
247             $guess = &'fmul(&'fadd($guess,&'fdiv($x,$guess,$gs*2)),".5");
248             $gs *= 2;
249         }
250         &'fround($guess, $scale);
251     }
252 }
253
254 1;