lib/bigrat.pl

   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]);
  78     substr($_,0,1) = '+' unless $_ eq 'NaN';
  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;