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