summaryrefslogtreecommitdiff
path: root/lib/bigrat.pl
blob: 2d3738f805b9ac7b596a7fa1e83d87ef11bc1080 (plain)
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
150
151
152
153
154
155
package bigrat;
require "bigint.pl";
#
# This library is no longer being maintained, and is included for backward
# compatibility with Perl 4 programs which may require it.
#
# In particular, this should not be used as an example of modern Perl
# programming techniques.
#
# 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;