home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #1 / NN_1993_1.iso / spool / comp / programm / 3420 < prev    next >
Encoding:
Internet Message Format  |  1993-01-06  |  4.4 KB
  1. Path: sparky!uunet!math.fu-berlin.de!mailgzrz.TU-Berlin.DE!news.netmbx.de!Germany.EU.net!urmel.informatik.rwth-aachen.de!martin
  2. From: martin@math.rwth-aachen.de (  Martin Schoenert)
  3. Newsgroups: comp.programming
  4. Subject: Re: Arbitrary precision with large bases
  5. Date: 6 Jan 93 09:10:24 GMT
  6. Organization: Rechnerbetrieb Informatik - RWTH Aachen
  7. Lines: 104
  8. Message-ID: <martin.726311424@bert>
  9. References: <1icne6INNpjb@uwm.edu>
  10. NNTP-Posting-Host: bert.math.rwth-aachen.de
  11.  
  12. markh@csd4.csd.uwm.edu (Mark) writes:
  13.  
  14.        A while back, I came up with an efficient arbitrary precision BCD
  15.     arithmetic package that uses base 1,000,000,000 (which happens to be
  16.     close to 2^30).
  17.  
  18.        I'm looking for other packages for comparison that use large bases,
  19.     either BCD or binary.  The non-trivial aspects to such as package are
  20.     the following:
  21.  
  22.        (1) An algorithm to carry out A, B -> AB div BASE, AB mod BASE
  23.        (2) A division algorithm that minimizes the amount of backtracking.
  24.            The grade school algorithm is a backtracking algorithm, and only
  25.            in base 2 is backtracking avoided.
  26.        (3) A square root algorithm, like the one taught in school, that likewise
  27.            minimizes backtracking.
  28.  
  29. I have looked at quite a few large integer arithmetic packages and have
  30. written two myself.  Mark Riordan has collected several, they are
  31. available for anonymous 'ftp' on 'cl-next2.cl.msu.edu' (35.8.4.22).
  32.  
  33. Most packages use a base that is a power of 2.  The most common choices
  34. are 2^16 and 2^32.  I assume that 2^16 does not count as a large base
  35. because with base 2^16 computing AB div BASE and AB mod BASE is not a
  36. problem at all.
  37.  
  38. Base 2^32 has the following advantages.
  39.  
  40. a)  the larger the base the fewer digits any given integer has in this
  41.     base.  Since (for example) the performance of the multiplication
  42.     is proportional to the product of the number of digits of the
  43.     operands it is always a win to use larger bases (though the win
  44.     with base 2^32 over base 10^9 is only about 12.5%).
  45.  
  46. b)  many modern processors have a 32bit x 32bit -> 64bit instruction.
  47.     This instruction cannot directly be used from within C, so one needs
  48.     a little bit of assembler.  However if one does that one gets
  49.     AB div BASE and AB mod BASE with one instruction.
  50.  
  51. Of course base 2^32 has the disadvantage that converting an integer to
  52. human readable format (i.e., base 10) is more work than in base 10^9.
  53.  
  54. With respect to aspect (2).  Knuth shows that if you normalize the
  55. divisor such that its leading digit is larger than or equal to BASE/2 and
  56. compute an appoximation of a digit of the quotient by dividing the
  57. leading three digits of the remainder by the leading two digits of the
  58. divisor the approximation is at most too large by 1 *independent* of the
  59. base.  It is easy to detect if the digit is too large and make the
  60. appropriate corrections.
  61.  
  62. With respect to aspect (3).  The  usual method to compute square roots is
  63. Newton's.
  64.  
  65. Here is another method, which I think is just a careful implementation of
  66. what you refer to as the scool method.  This method is considerably
  67. faster than Newton's, in principle it is as fast as squaring the result,
  68. in practice it is about a factor of 2 slower in my arithmetic.
  69.  
  70. Let BB be a power of BASE and let
  71.  
  72.     A  =  A_1 BB^2 + A_0
  73.  
  74. such that BB^2/4 <= A_1.  Compute the upper half of the root recursively
  75.  
  76.     X_1  =  Floor( Root( A_1 ) );
  77.     R  =  A - (X_1 BB)^2;
  78.  
  79. now compute the approximation for the rest of the root
  80.  
  81.     X_0  =  R / (2 X_1 BB);
  82.     S  =  R - 2 X_0 X_1 BB;
  83.  
  84. this approximation may be off by one, correct if neccessary
  85.  
  86.     if S < X_0^2  then
  87.         X_0  =  X_0 - 1;
  88.         S  =  S + 2 X_1 BB;  /* = R - 2 X_0 X_1 BB again */
  89.     fi;
  90.  
  91. then with
  92.  
  93.     X  =  X_1 BB + X_0;
  94.     T  =  S - X_0^2;
  95.  
  96. we have
  97.  
  98.     X  =  Floor( Root( A ) );
  99.     T  =  A - X^2;
  100.  
  101. You can either choose BB to be BASE, in which case each iteration
  102. computes one more digit of the root, or you can choose BB such that A_1
  103. and A_0 have the same number of digits, in which case each iteration
  104. doubles the number of digits of the root.  The second is faster for very
  105. large integers (~1000 decimal digits) if one uses Karatsuba's trick for
  106. the computation of X_0^2.
  107.  
  108. This algorithm is again independent of the value of BASE.
  109.  
  110. Hope this helps,  Martin.
  111.  
  112. -- .- .-. - .. -.  .-.. --- ...- . ...  .- -. -. .. -.- .-
  113. Martin Sch"onert,   Martin.Schoenert@Math.RWTH-Aachen.DE,  +49 241 804551
  114. Lehrstuhl D f"ur Mathematik, Templergraben 64, RWTH, D 51 Aachen, Germany
  115.  
  116.