extra code in pp_concat, Take 2
[p5sagit/p5-mst-13.2.git] / lib / bigrat.pl
1 package bigrat;
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 # Arbitrary size rational math package
11 #
12 # by Mark Biggar
13 #
14 # Input values to these routines consist of strings of the form 
15 #   m|^\s*[+-]?[\d\s]+(/[\d\s]+)?$|.
16 # Examples:
17 #   "+0/1"                          canonical zero value
18 #   "3"                             canonical value "+3/1"
19 #   "   -123/123 123"               canonical value "-1/1001"
20 #   "123 456/7890"                  canonical value "+20576/1315"
21 # Output values always include a sign and no leading zeros or
22 #   white space.
23 # This package makes use of the bigint package.
24 # The string 'NaN' is used to represent the result when input arguments 
25 #   that are not numbers, as well as the result of dividing by zero and
26 #       the sqrt of a negative number.
27 # Extreamly naive algorthims are used.
28 #
29 # Routines provided are:
30 #
31 #   rneg(RAT) return RAT                negation
32 #   rabs(RAT) return RAT                absolute value
33 #   rcmp(RAT,RAT) return CODE           compare numbers (undef,<0,=0,>0)
34 #   radd(RAT,RAT) return RAT            addition
35 #   rsub(RAT,RAT) return RAT            subtraction
36 #   rmul(RAT,RAT) return RAT            multiplication
37 #   rdiv(RAT,RAT) return RAT            division
38 #   rmod(RAT) return (RAT,RAT)          integer and fractional parts
39 #   rnorm(RAT) return RAT               normalization
40 #   rsqrt(RAT, cycles) return RAT       square root
41 \f
42 # Convert a number to the canonical string form m|^[+-]\d+/\d+|.
43 sub main'rnorm { #(string) return rat_num
44     local($_) = @_;
45     s/\s+//g;
46     if (m#^([+-]?\d+)(/(\d*[1-9]0*))?$#) {
47         &norm($1, $3 ? $3 : '+1');
48     } else {
49         'NaN';
50     }
51 }
52
53 # Normalize by reducing to lowest terms
54 sub norm { #(bint, bint) return rat_num
55     local($num,$dom) = @_;
56     if ($num eq 'NaN') {
57         'NaN';
58     } elsif ($dom eq 'NaN') {
59         'NaN';
60     } elsif ($dom =~ /^[+-]?0+$/) {
61         'NaN';
62     } else {
63         local($gcd) = &'bgcd($num,$dom);
64         $gcd =~ s/^-/+/;
65         if ($gcd ne '+1') { 
66             $num = &'bdiv($num,$gcd);
67             $dom = &'bdiv($dom,$gcd);
68         } else {
69             $num = &'bnorm($num);
70             $dom = &'bnorm($dom);
71         }
72         substr($dom,$[,1) = '';
73         "$num/$dom";
74     }
75 }
76
77 # negation
78 sub main'rneg { #(rat_num) return rat_num
79     local($_) = &'rnorm(@_);
80     tr/-+/+-/ if ($_ ne '+0/1');
81     $_;
82 }
83
84 # absolute value
85 sub main'rabs { #(rat_num) return $rat_num
86     local($_) = &'rnorm(@_);
87     substr($_,$[,1) = '+' unless $_ eq 'NaN';
88     $_;
89 }
90
91 # multipication
92 sub main'rmul { #(rat_num, rat_num) return rat_num
93     local($xn,$xd) = split('/',&'rnorm($_[$[]));
94     local($yn,$yd) = split('/',&'rnorm($_[$[+1]));
95     &norm(&'bmul($xn,$yn),&'bmul($xd,$yd));
96 }
97
98 # division
99 sub main'rdiv { #(rat_num, rat_num) return rat_num
100     local($xn,$xd) = split('/',&'rnorm($_[$[]));
101     local($yn,$yd) = split('/',&'rnorm($_[$[+1]));
102     &norm(&'bmul($xn,$yd),&'bmul($xd,$yn));
103 }
104 \f
105 # addition
106 sub main'radd { #(rat_num, rat_num) return rat_num
107     local($xn,$xd) = split('/',&'rnorm($_[$[]));
108     local($yn,$yd) = split('/',&'rnorm($_[$[+1]));
109     &norm(&'badd(&'bmul($xn,$yd),&'bmul($yn,$xd)),&'bmul($xd,$yd));
110 }
111
112 # subtraction
113 sub main'rsub { #(rat_num, rat_num) return rat_num
114     local($xn,$xd) = split('/',&'rnorm($_[$[]));
115     local($yn,$yd) = split('/',&'rnorm($_[$[+1]));
116     &norm(&'bsub(&'bmul($xn,$yd),&'bmul($yn,$xd)),&'bmul($xd,$yd));
117 }
118
119 # comparison
120 sub main'rcmp { #(rat_num, rat_num) return cond_code
121     local($xn,$xd) = split('/',&'rnorm($_[$[]));
122     local($yn,$yd) = split('/',&'rnorm($_[$[+1]));
123     &bigint'cmp(&'bmul($xn,$yd),&'bmul($yn,$xd));
124 }
125
126 # int and frac parts
127 sub main'rmod { #(rat_num) return (rat_num,rat_num)
128     local($xn,$xd) = split('/',&'rnorm(@_));
129     local($i,$f) = &'bdiv($xn,$xd);
130     if (wantarray) {
131         ("$i/1", "$f/$xd");
132     } else {
133         "$i/1";
134     }   
135 }
136
137 # square root by Newtons method.
138 #   cycles specifies the number of iterations default: 5
139 sub main'rsqrt { #(fnum_str[, cycles]) return fnum_str
140     local($x, $scale) = (&'rnorm($_[$[]), $_[$[+1]);
141     if ($x eq 'NaN') {
142         'NaN';
143     } elsif ($x =~ /^-/) {
144         'NaN';
145     } else {
146         local($gscale, $guess) = (0, '+1/1');
147         $scale = 5 if (!$scale);
148         while ($gscale++ < $scale) {
149             $guess = &'rmul(&'radd($guess,&'rdiv($x,$guess)),"+1/2");
150         }
151         "$guess";          # quotes necessary due to perl bug
152     }
153 }
154
155 1;