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