Source code of Windows XP (NT5)
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

285 lines
7.8 KiB

  1. package bigint;
  2. # arbitrary size integer math package
  3. #
  4. # by Mark Biggar
  5. #
  6. # Canonical Big integer value are strings of the form
  7. # /^[+-]\d+$/ with leading zeros suppressed
  8. # Input values to these routines may be strings of the form
  9. # /^\s*[+-]?[\d\s]+$/.
  10. # Examples:
  11. # '+0' canonical zero value
  12. # ' -123 123 123' canonical value '-123123123'
  13. # '1 23 456 7890' canonical value '+1234567890'
  14. # Output values always always in canonical form
  15. #
  16. # Actual math is done in an internal format consisting of an array
  17. # whose first element is the sign (/^[+-]$/) and whose remaining
  18. # elements are base 100000 digits with the least significant digit first.
  19. # The string 'NaN' is used to represent the result when input arguments
  20. # are not numbers, as well as the result of dividing by zero
  21. #
  22. # routines provided are:
  23. #
  24. # bneg(BINT) return BINT negation
  25. # babs(BINT) return BINT absolute value
  26. # bcmp(BINT,BINT) return CODE compare numbers (undef,<0,=0,>0)
  27. # badd(BINT,BINT) return BINT addition
  28. # bsub(BINT,BINT) return BINT subtraction
  29. # bmul(BINT,BINT) return BINT multiplication
  30. # bdiv(BINT,BINT) return (BINT,BINT) division (quo,rem) just quo if scalar
  31. # bmod(BINT,BINT) return BINT modulus
  32. # bgcd(BINT,BINT) return BINT greatest common divisor
  33. # bnorm(BINT) return BINT normalization
  34. #
  35. $zero = 0;
  36. # normalize string form of number. Strip leading zeros. Strip any
  37. # white space and add a sign, if missing.
  38. # Strings that are not numbers result the value 'NaN'.
  39. sub main'bnorm { #(num_str) return num_str
  40. local($_) = @_;
  41. s/\s+//g; # strip white space
  42. if (s/^([+-]?)0*(\d+)$/$1$2/) { # test if number
  43. substr($_,$[,0) = '+' unless $1; # Add missing sign
  44. s/^-0/+0/;
  45. $_;
  46. } else {
  47. 'NaN';
  48. }
  49. }
  50. # Convert a number from string format to internal base 100000 format.
  51. # Assumes normalized value as input.
  52. sub internal { #(num_str) return int_num_array
  53. local($d) = @_;
  54. ($is,$il) = (substr($d,$[,1),length($d)-2);
  55. substr($d,$[,1) = '';
  56. ($is, reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $d)));
  57. }
  58. # Convert a number from internal base 100000 format to string format.
  59. # This routine scribbles all over input array.
  60. sub external { #(int_num_array) return num_str
  61. $es = shift;
  62. grep($_ > 9999 || ($_ = substr('0000'.$_,-5)), @_); # zero pad
  63. &'bnorm(join('', $es, reverse(@_))); # reverse concat and normalize
  64. }
  65. # Negate input value.
  66. sub main'bneg { #(num_str) return num_str
  67. local($_) = &'bnorm(@_);
  68. vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0';
  69. s/^./N/ unless /^[-+]/; # works both in ASCII and EBCDIC
  70. $_;
  71. }
  72. # Returns the absolute value of the input.
  73. sub main'babs { #(num_str) return num_str
  74. &abs(&'bnorm(@_));
  75. }
  76. sub abs { # post-normalized abs for internal use
  77. local($_) = @_;
  78. s/^-/+/;
  79. $_;
  80. }
  81. # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
  82. sub main'bcmp { #(num_str, num_str) return cond_code
  83. local($x,$y) = (&'bnorm($_[$[]),&'bnorm($_[$[+1]));
  84. if ($x eq 'NaN') {
  85. undef;
  86. } elsif ($y eq 'NaN') {
  87. undef;
  88. } else {
  89. &cmp($x,$y);
  90. }
  91. }
  92. sub cmp { # post-normalized compare for internal use
  93. local($cx, $cy) = @_;
  94. return 0 if ($cx eq $cy);
  95. local($sx, $sy) = (substr($cx, 0, 1), substr($cy, 0, 1));
  96. local($ld);
  97. if ($sx eq '+') {
  98. return 1 if ($sy eq '-' || $cy eq '+0');
  99. $ld = length($cx) - length($cy);
  100. return $ld if ($ld);
  101. return $cx cmp $cy;
  102. } else { # $sx eq '-'
  103. return -1 if ($sy eq '+');
  104. $ld = length($cy) - length($cx);
  105. return $ld if ($ld);
  106. return $cy cmp $cx;
  107. }
  108. }
  109. sub main'badd { #(num_str, num_str) return num_str
  110. local(*x, *y); ($x, $y) = (&'bnorm($_[$[]),&'bnorm($_[$[+1]));
  111. if ($x eq 'NaN') {
  112. 'NaN';
  113. } elsif ($y eq 'NaN') {
  114. 'NaN';
  115. } else {
  116. @x = &internal($x); # convert to internal form
  117. @y = &internal($y);
  118. local($sx, $sy) = (shift @x, shift @y); # get signs
  119. if ($sx eq $sy) {
  120. &external($sx, &add(*x, *y)); # if same sign add
  121. } else {
  122. ($x, $y) = (&abs($x),&abs($y)); # make abs
  123. if (&cmp($y,$x) > 0) {
  124. &external($sy, &sub(*y, *x));
  125. } else {
  126. &external($sx, &sub(*x, *y));
  127. }
  128. }
  129. }
  130. }
  131. sub main'bsub { #(num_str, num_str) return num_str
  132. &'badd($_[$[],&'bneg($_[$[+1]));
  133. }
  134. # GCD -- Euclids algorithm Knuth Vol 2 pg 296
  135. sub main'bgcd { #(num_str, num_str) return num_str
  136. local($x,$y) = (&'bnorm($_[$[]),&'bnorm($_[$[+1]));
  137. if ($x eq 'NaN' || $y eq 'NaN') {
  138. 'NaN';
  139. } else {
  140. ($x, $y) = ($y,&'bmod($x,$y)) while $y ne '+0';
  141. $x;
  142. }
  143. }
  144. # routine to add two base 1e5 numbers
  145. # stolen from Knuth Vol 2 Algorithm A pg 231
  146. # there are separate routines to add and sub as per Kunth pg 233
  147. sub add { #(int_num_array, int_num_array) return int_num_array
  148. local(*x, *y) = @_;
  149. $car = 0;
  150. for $x (@x) {
  151. last unless @y || $car;
  152. $x -= 1e5 if $car = (($x += shift(@y) + $car) >= 1e5) ? 1 : 0;
  153. }
  154. for $y (@y) {
  155. last unless $car;
  156. $y -= 1e5 if $car = (($y += $car) >= 1e5) ? 1 : 0;
  157. }
  158. (@x, @y, $car);
  159. }
  160. # subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
  161. sub sub { #(int_num_array, int_num_array) return int_num_array
  162. local(*sx, *sy) = @_;
  163. $bar = 0;
  164. for $sx (@sx) {
  165. last unless @y || $bar;
  166. $sx += 1e5 if $bar = (($sx -= shift(@sy) + $bar) < 0);
  167. }
  168. @sx;
  169. }
  170. # multiply two numbers -- stolen from Knuth Vol 2 pg 233
  171. sub main'bmul { #(num_str, num_str) return num_str
  172. local(*x, *y); ($x, $y) = (&'bnorm($_[$[]), &'bnorm($_[$[+1]));
  173. if ($x eq 'NaN') {
  174. 'NaN';
  175. } elsif ($y eq 'NaN') {
  176. 'NaN';
  177. } else {
  178. @x = &internal($x);
  179. @y = &internal($y);
  180. local($signr) = (shift @x ne shift @y) ? '-' : '+';
  181. @prod = ();
  182. for $x (@x) {
  183. ($car, $cty) = (0, $[);
  184. for $y (@y) {
  185. $prod = $x * $y + $prod[$cty] + $car;
  186. $prod[$cty++] =
  187. $prod - ($car = int($prod * 1e-5)) * 1e5;
  188. }
  189. $prod[$cty] += $car if $car;
  190. $x = shift @prod;
  191. }
  192. &external($signr, @x, @prod);
  193. }
  194. }
  195. # modulus
  196. sub main'bmod { #(num_str, num_str) return num_str
  197. (&'bdiv(@_))[$[+1];
  198. }
  199. sub main'bdiv { #(dividend: num_str, divisor: num_str) return num_str
  200. local (*x, *y); ($x, $y) = (&'bnorm($_[$[]), &'bnorm($_[$[+1]));
  201. return wantarray ? ('NaN','NaN') : 'NaN'
  202. if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0');
  203. return wantarray ? ('+0',$x) : '+0' if (&cmp(&abs($x),&abs($y)) < 0);
  204. @x = &internal($x); @y = &internal($y);
  205. $srem = $y[$[];
  206. $sr = (shift @x ne shift @y) ? '-' : '+';
  207. $car = $bar = $prd = 0;
  208. if (($dd = int(1e5/($y[$#y]+1))) != 1) {
  209. for $x (@x) {
  210. $x = $x * $dd + $car;
  211. $x -= ($car = int($x * 1e-5)) * 1e5;
  212. }
  213. push(@x, $car); $car = 0;
  214. for $y (@y) {
  215. $y = $y * $dd + $car;
  216. $y -= ($car = int($y * 1e-5)) * 1e5;
  217. }
  218. }
  219. else {
  220. push(@x, 0);
  221. }
  222. @q = (); ($v2,$v1) = @y[-2,-1];
  223. while ($#x > $#y) {
  224. ($u2,$u1,$u0) = @x[-3..-1];
  225. $q = (($u0 == $v1) ? 99999 : int(($u0*1e5+$u1)/$v1));
  226. --$q while ($v2*$q > ($u0*1e5+$u1-$q*$v1)*1e5+$u2);
  227. if ($q) {
  228. ($car, $bar) = (0,0);
  229. for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
  230. $prd = $q * $y[$y] + $car;
  231. $prd -= ($car = int($prd * 1e-5)) * 1e5;
  232. $x[$x] += 1e5 if ($bar = (($x[$x] -= $prd + $bar) < 0));
  233. }
  234. if ($x[$#x] < $car + $bar) {
  235. $car = 0; --$q;
  236. for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
  237. $x[$x] -= 1e5
  238. if ($car = (($x[$x] += $y[$y] + $car) > 1e5));
  239. }
  240. }
  241. }
  242. pop(@x); unshift(@q, $q);
  243. }
  244. if (wantarray) {
  245. @d = ();
  246. if ($dd != 1) {
  247. $car = 0;
  248. for $x (reverse @x) {
  249. $prd = $car * 1e5 + $x;
  250. $car = $prd - ($tmp = int($prd / $dd)) * $dd;
  251. unshift(@d, $tmp);
  252. }
  253. }
  254. else {
  255. @d = @x;
  256. }
  257. (&external($sr, @q), &external($srem, @d, $zero));
  258. } else {
  259. &external($sr, @q);
  260. }
  261. }
  262. 1;