contrib/perl5/lib/bigrat.pl

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