home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / perl-5.003-base.tgz / perl-5.003-base.tar / fsf / perl / lib / bigrat.pl < prev    next >
Perl Script  |  1995-03-10  |  4KB  |  150 lines
  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.  
  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.  
  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;
  150.