diff options
Diffstat (limited to 'ipl/procs/rational.icn')
| -rw-r--r-- | ipl/procs/rational.icn | 220 |
1 files changed, 220 insertions, 0 deletions
diff --git a/ipl/procs/rational.icn b/ipl/procs/rational.icn new file mode 100644 index 0000000..0f3c311 --- /dev/null +++ b/ipl/procs/rational.icn @@ -0,0 +1,220 @@ +############################################################################ +# +# File: rational.icn +# +# Subject: Procedures for arithmetic on rational numbers +# +# Author: Ralph E. Griswold +# +# Date: June 10, 2001 +# +############################################################################ +# +# This file is in the public domain. +# +############################################################################ +# +# Contributor: Gregg M. Townsend +# +############################################################################ +# +# These procedures perform arithmetic on rational numbers (fractions): +# +# addrat(r1,r2) Add rational numbers r1 and r2. +# +# divrat(r1,r2) Divide rational numbers r1 and r2. +# +# medrat(r1,r2) Form mediant of r1 and r2. +# +# mpyrat(r1,r2) Multiply rational numbers r1 and r2. +# +# negrat(r) Produce negative of rational number r. +# +# rat2real(r) Produce floating-point approximation of r +# +# rat2str(r) Convert the rational number r to its string +# representation. +# +# real2rat(v,p) Convert real to rational with precision p. +# The default precision is 1e-10. +# (Too much precision gives huge, ugly factions.) +# +# reciprat(r) Produce the reciprocal of rational number r. +# +# str2rat(s) Convert the string representation of a rational number +# (such as "3/2") to a rational number. +# +# subrat(r1,r2) Subtract rational numbers r1 and r2. +# +############################################################################ +# +# Links: numbers +# +############################################################################ + +link numbers + +record rational(numer, denom, sign) + +procedure addrat(r1, r2) #: sum of rationals + local denom, numer, div + + r1 := ratred(r1) + r2 := ratred(r2) + + denom := r1.denom * r2.denom + numer := r1.sign * r1.numer * r2.denom + + r2.sign * r2.numer * r1.denom + + if numer = 0 then return rational (0, 1, 1) + + div := gcd(numer, denom) + + return rational(abs(numer / div), abs(denom / div), numer / abs(numer)) + +end + +procedure divrat(r1, r2) #: divide rationals. + + r1 := ratred(r1) + r2 := ratred(r2) + + return mpyrat(r1, reciprat(r2)) + +end + +procedure medrat(r1, r2) #: form rational mediant + local numer, denom, div + + r1 := ratred(r1) + r2 := ratred(r2) + + numer := r1.numer + r2.numer + denom := r1.denom + r2.denom + + div := gcd(numer, denom) + + return rational(numer / div, denom / div, r1.sign * r2.sign) + +end + +procedure mpyrat(r1, r2) #: multiply rationals + local numer, denom, div + + r1 := ratred(r1) + r2 := ratred(r2) + + numer := r1.numer * r2.numer + denom := r1.denom * r2.denom + + div := gcd(numer, denom) + + return rational(numer / div, denom / div, r1.sign * r2.sign) + +end + +procedure negrat(r) #: negative of rational + + r := ratred(r) + + return rational(r.numer, r.denom, -r.sign) + +end + +procedure rat2real(r) #: floating-point approximation of rational + + r := ratred(r) + + return (real(r.numer) * r.sign) / r.denom + +end + +procedure rat2str(r) #: convert rational to string + + r := ratred(r) + + return "(" || (r.numer * r.sign) || "/" || r.denom || ")" + +end + +procedure ratred(r) #: reduce rational to lowest terms + local div + + if r.denom = 0 then runerr(501) + if abs(r.sign) ~= 1 then runerr(501) + + if r.numer = 0 then return rational(0, 1, 1) + + if r.numer < 0 then r.sign *:= -1 + if r.denom < 0 then r.sign *:= -1 + + r.numer := abs(r.numer) + r.denom := abs(r.denom) + + div := gcd(r.numer, r.denom) + + return rational(r.numer / div, r.denom / div, r.sign) + +end + +# real2rat(v, p) -- convert real to rational with precision p +# +# Originally based on a calculator algorithm posted to usenet on August 19, +# 1987, by Joseph D. Rudmin, Duke University Physics Dept. (duke!dukempd!jdr) + +$define MAXITER 40 # maximum number of iterations +$define PRECISION 1e-10 # default conversion precision + +procedure real2rat(r, p) #: convert to rational with precision p + local t, d, i, j + static x, y + initial { x := list(MAXITER); y := list(MAXITER + 2) } + + t := abs(r) + /p := PRECISION + every i := 1 to MAXITER do { + x[i] := integer(t) + y[i + 1] := 1 + y[i + 2] := 0 + every j := i to 1 by -1 do + y[j] := x[j] * y[j + 1] + y[j + 2] + if abs(y[1] / real(y[2]) - r) < p then break + d := t - integer(t) + if d < p then break + t := 1.0 / d + } + return rational(y[1], y[2], if r >= 0 then 1 else -1) + +end + +procedure reciprat(r) #: reciprocal of rational + + r := ratred(r) + + return rational(r.denom, r.numer, r.sign) + +end + +procedure str2rat(s) # convert string to rational + local div, numer, denom, sign + + s ? { + ="(" & + numer := integer(tab(upto('/'))) & + move(1) & + denom := integer(tab(upto(')'))) & + pos(-1) + } | fail + + return ratred(rational(numer, denom, 1)) + +end + +procedure subrat(r1, r2) #: difference of rationals + + r1 := ratred(r1) + r2 := ratred(r2) + + return addrat(r1, negrat(r2)) + +end |