Long Precision Arithmetic

1
Long Precision Arithmetic
Dr. Peter-Michael Seidel Computer Science and
Engineering Southern Methodist University Dallas,
TX, 75275
2
Basic Data Structures (1)

• Data Structure
• Single digits fixed to 32 bits (1 Digit
  Word)
• Long Precision Data Structure
• numbers capable of fixed digits 4096 bits
  (multiplication of two 2048 bit numbers)
  -gt 128 digits,
  fractions ?
• access to units of Words(32bit) DWord(64bit)
• taking care of machine precision and available
  data types
• memory allocation has large impact on running
  times
• Operations
• basic arithmetic operations (Long Precision and
  Modular Arithmetic)
• Constructor, Destructor
• GetDigit() to access an arbitrary digit of a
  number

3
Variable Precision Data Structure

• Digit
• Implemented based on Long Precision Digit (32
  or 16 bits)
• Long Precision Number
• Number Components
• magnitude Arraym_iMaxLen of Word
• sign Bool
• actual length m_MSD pos. Integer in 0
  m_iMaxLen

4
Operations

• Constructors (Destructor)
• Four different types possible (including copying
  existing numbers)
• Import/Export Functions
• To convert numbers and interface
• Addition
• Use Assembler Macro to add 2 Digits Carry
• Addition only necessary for digits up to
  minA.MSD, B.MSD
• Propagation of Carries for digits up to
  maxA.MSD, B.MSD
• Adjustment of actual result length addition
  maxA.MSD, B.MSD, maxA.MSD, B.MSD1
  subtraction 0 maxA.MSD, B.MSD
• Multiplication
• Consider recursive methods that reduce problem to
  addition, e.g. Karatsubas Algorithm
• Comparing fixed solutions at 16 bits and
  at 32 bits
• Considering School Method for lowest level(s)
  (fewer additions)

5
FP Expansions

• Represent long precision as sum of FP numbers
• aligned significands
• Range
• Single precision (exponent range
  -126,127)
• Could represent numbers up to 2128 with
  resolution of 2-139 (266 bits)
• Exact representation requires up to 12 FP numbers
• Double precision (exponent range
  -1022,1023)
• Could represent numbers up to 21024 with
  resolution of 2-1074 (2096 bits)
• Exact representation requires up to 40 FP numbers
• Extension Could be considered as digits in
  larger system

6
Exact Addition of 2 FP numbers

• Find C,D, so that
• AB CD
• C,D non-overlapping
• D lt C2-p
• Computation in 3 steps (assume A gt B)
• C RNE(AB)
• E RNE(C-A) (representable)
• D RNE(B-E) (representable)
• Achieves most compact representation for FP
  Expansions
• Sort, process from right to left (effort can be
  reduced to 2 adds/FP number)

A
B
C
D
7
Multiplication

• Of two n-digit numbers

8
How to multiply 2 n-digit numbers.

X





2





9
Multiplication of 2 n-bit numbers

• X
• Y
• X a 2n/2 b Y c 2n/2 d
• XY ac 2n (adbc) 2n/2 bd

a
b
c
d
10
Mult of 2 n-bit(digit) numbers

• X
• Y
• XY ac 2n (adbc) 2n/2 bd

a
b
c
d
MULT(X,Y) If X Y 1 then RETURN XY
Break X into ab and Y into cd RETURN
MULT(a,c) 2n (MULT(a,d) MULT(b,c)) 2n/2
MULT(b,d)
11
Time required by MULT

• T(n) 4T(n/2) 3TAdd(n)
• T(1) 1
• Form T(n) an2 bn c for some a,b,c
• Calculate T(1) 1 ? a b c 1
• T(n) 4 T(n/2) 3n
• ? an2 bn c 4 a(n/2)2 b(n/2) c 3n
• an2 2bn 4c 3n
• ? (b3)n 3c 0
• Therefore b-3 c0 a4
• gt T(n) 4n2 - 3n

12
MULT computation
MULT(X,Y) If X Y 1 then RETURN XY
Break X into ab and Y into cd RETURN
MULT(a,c) 2n (MULT(a,d) MULT(b,c)) 2n/2
MULT(b,d)

• MULT is recursively called 4 times.

13
Karatsuba-Ofman Multiplication

• Input a,b,c,d Output ac, adbc, bd
• X1 a b
• X2 c d
• X3 X1 X2 ac ad bc bd
• X4 ac
• X5 bd
• X6 X3 X4 - X5 ad bc
• T(n) 3T(n/2) 3TAdd(n)
• T(1) 1

14
Karatsuba-Ofman Mult 1962
MULT(X,Y) If X Y 1 then RETURN XY
Break X into ab and Y into cd e MULT(a,c)
and f MULT(b,d) RETURN e2n (MULT(ab, cd)
e - f) 2n/2 f

• T(n) 3 T(n/2) 3n
• Actually T(n) 2 T(n/2) T(n/2 1) 3n

15
T(n)
3n

T(n/2)
T(n/2)
T(n/2)
16
T(n)
3n

T(n/2)
T(n/2)
17
T(n)
3n

18
3n
3n/2 3n/2 3n/2


Level i is the sum of 3i copies of 3n/2i

. . . . . . . . . . . . . . . . . . . . . . . . . .
11111111111111111111111111111111111111111111
0
1
2
3n/4 3n/4 3n/4 3n/4 3n/4 3n/4 3n/4
3n/4 3n/4
i
19
3n
9/2n
27/4n
3(3/2)in
3n(13/2(3/2)2 . . . (3/2) ) -2n(3/2)
20
The Geometric Sum
21
Toom-Cook Multiplication

• Recursive reduction of n bit/digit multiplication
• to 5 multiplications of size n/3 some
  additions
• gt T(n) O(nlog5/log3) O(n1.465)

22
Multiplication Algorithms
conventional n2
Karatsuba n1.58
Toom - Cook n1.465
FFT (Brigham) n logn2
(Schönhage Strassen) n logn loglogn
23
Practical Comparisons

• Implementation in GMP / NTL (GNU Multiple
  Precision)

24
REFERENCES

• Borodin, A. and Munro, I. The Computational
  Complexity of Algebraic and Numeric Problems. New
  York American Elsevier, 1975.
• Borwein, J. M. Borwein, P. B. and Bailey, D. H.
  "Ramanujan, Modular Equations, and Approximations
  to Pi, or How to Compute One Billion Digits of
  Pi." Amer. Math. Monthly 96, 201-219, 1989.
• Brigham, E. O. The Fast Fourier Transform.
  Englewood Cliffs, NJ Prentice-Hall, 1974.
• Brigham, E. O. Fast Fourier Transform and
  Applications. Englewood Cliffs, NJ
  Prentice-Hall, 1988.
• Cook, S. A. On the Minimum Computation Time of
  Functions. Ph.D. Thesis. Cambridge, MA Harvard
  University, pp. 51-77, 1966.
• Hollerbach, U. "Fast Multiplication Division of
  Very Large Numbers." sci.math.research posting,
  Jan. 23, 1996.
• Karatsuba, A. and Ofman, Yu. "Multiplication of
  Many-Digital Numbers by Automatic Computers."
  Doklady Akad. Nauk SSSR 145, 293-294, 1962.
  Translation in Physics-Doklady 7, 595-596, 1963.
• Knuth, D. E. The Art of Computer Programming,
  Vol. 2 Seminumerical Algorithms, 3rd ed.
  Reading, MA Addison-Wesley, pp. 278-286, 1998.
• Schönhage, A. and Strassen, V. "Schnelle
  Multiplikation Grosser Zahlen." Computing 7,
  281-292, 1971.
• Toom, A. L. "The Complexity of a Scheme of
  Functional Elements Simulating the Multiplication
  of Integers." Dokl. Akad. Nauk SSSR 150, 496-498,
  1963. English translation in Soviet Mathematics
  3, 714-716, 1963.
• Zuras, D. "More on Squaring and Multiplying Large
  Integers." IEEE Trans. Comput. 43, 899-908, 1994.