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