home *** CD-ROM | disk | FTP | other *** search
/ back2roots/padua / padua.7z / padua / lang / perl4.035.V010.lzh / perl4.035 / lib / bigrat.pl < prev    next >
Encoding:
Perl Script  |  1992-08-05  |  4.1 KB  |  149 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.     if ($gcd ne '+1') { 
  59.         $num = &'bdiv($num,$gcd);
  60.         $dom = &'bdiv($dom,$gcd);
  61.     } else {
  62.         $num = &'bnorm($num);
  63.         $dom = &'bnorm($dom);
  64.     }
  65.     substr($dom,0,1) = '';
  66.     "$num/$dom";
  67.     }
  68. }
  69.  
  70. # negation
  71. sub main'rneg { #(rat_num) return rat_num
  72.     local($_) = &'rnorm($_[0]);
  73.     tr/-+/+-/ if ($_ ne '+0/1');
  74.     $_;
  75. }
  76.  
  77. # absolute value
  78. sub main'rabs { #(rat_num) return $rat_num
  79.     local($_) = &'rnorm($_[0]);
  80.     substr($_,0,1) = '+' unless $_ eq 'NaN';
  81.     $_;
  82. }
  83.  
  84. # multipication
  85. sub main'rmul { #(rat_num, rat_num) return rat_num
  86.     local($xn,$xd) = split('/',&'rnorm($_[0]));
  87.     local($yn,$yd) = split('/',&'rnorm($_[1]));
  88.     &norm(&'bmul($xn,$yn),&'bmul($xd,$yd));
  89. }
  90.  
  91. # division
  92. sub main'rdiv { #(rat_num, rat_num) return rat_num
  93.     local($xn,$xd) = split('/',&'rnorm($_[0]));
  94.     local($yn,$yd) = split('/',&'rnorm($_[1]));
  95.     &norm(&'bmul($xn,$yd),&'bmul($xd,$yn));
  96. }
  97.  
  98. # addition
  99. sub main'radd { #(rat_num, rat_num) return rat_num
  100.     local($xn,$xd) = split('/',&'rnorm($_[0]));
  101.     local($yn,$yd) = split('/',&'rnorm($_[1]));
  102.     &norm(&'badd(&'bmul($xn,$yd),&'bmul($yn,$xd)),&'bmul($xd,$yd));
  103. }
  104.  
  105. # subtraction
  106. sub main'rsub { #(rat_num, rat_num) return rat_num
  107.     local($xn,$xd) = split('/',&'rnorm($_[0]));
  108.     local($yn,$yd) = split('/',&'rnorm($_[1]));
  109.     &norm(&'bsub(&'bmul($xn,$yd),&'bmul($yn,$xd)),&'bmul($xd,$yd));
  110. }
  111.  
  112. # comparison
  113. sub main'rcmp { #(rat_num, rat_num) return cond_code
  114.     local($xn,$xd) = split('/',&'rnorm($_[0]));
  115.     local($yn,$yd) = split('/',&'rnorm($_[1]));
  116.     &bigint'cmp(&'bmul($xn,$yd),&'bmul($yn,$xd));
  117. }
  118.  
  119. # int and frac parts
  120. sub main'rmod { #(rat_num) return (rat_num,rat_num)
  121.     local($xn,$xd) = split('/',&'rnorm($_[0]));
  122.     local($i,$f) = &'bdiv($xn,$xd);
  123.     if (wantarray) {
  124.     ("$i/1", "$f/$xd");
  125.     } else {
  126.     "$i/1";
  127.     }   
  128. }
  129.  
  130. # square root by Newtons method.
  131. #   cycles specifies the number of iterations default: 5
  132. sub main'rsqrt { #(fnum_str[, cycles]) return fnum_str
  133.     local($x, $scale) = (&'rnorm($_[0]), $_[1]);
  134.     if ($x eq 'NaN') {
  135.     'NaN';
  136.     } elsif ($x =~ /^-/) {
  137.     'NaN';
  138.     } else {
  139.     local($gscale, $guess) = (0, '+1/1');
  140.     $scale = 5 if (!$scale);
  141.     while ($gscale++ < $scale) {
  142.         $guess = &'rmul(&'radd($guess,&'rdiv($x,$guess)),"+1/2");
  143.     }
  144.     "$guess";          # quotes necessary due to perl bug
  145.     }
  146. }
  147.  
  148. 1;
  149.