Annotation of OpenXM_contrib/gmp/doc/multiplication, Revision 1.1
1.1 ! maekawa 1:
! 2: GMP MULTIPLICATION
! 3:
! 4:
! 5: This file describes briefly the multiplication and squaring used in GMP.
! 6: The code is likely to be hard to understand without knowing something about
! 7: the algorithms.
! 8:
! 9: GMP does NxN limb multiplications and squares using one of four algorithms,
! 10: according to the size N.
! 11:
! 12: Algorithm Sizes
! 13:
! 14: basecase < KARATSUBA_MUL_THRESHOLD
! 15: karatsuba >= KARATSUBA_MUL_THRESHOLD
! 16: toom3 >= TOOM3_MUL_THRESHOLD
! 17: fft >= FFT_MUL_THRESHOLD
! 18:
! 19: Similarly for squaring, with the SQR thresholds. Note though that the FFT
! 20: is only used if GMP is configured with --enable-fft.
! 21:
! 22: MxN multiplications of operands with different sizes are currently done by
! 23: splitting it up into various NxN pieces while above KARATSUBA_MUL_THRESHOLD.
! 24: The karatsuba and toom3 routines then operate only on equal size operands.
! 25: This is rather inefficient, and is slated for improvement in the future.
! 26:
! 27:
! 28:
! 29: BASECASE MULTIPLY
! 30:
! 31: This a straightforward rectangular set of cross-products, the same as long
! 32: multiplication done by hand and for that reason sometimes known as the
! 33: schoolbook or grammar school method.
! 34:
! 35: See Knuth (reference in gmp.texi) volume 2 section 4.3.1 algorithm M, or the
! 36: mpn/generic/mul_basecase.c code.
! 37:
! 38: Assembler implementations of mul_basecase are essentially the same as the
! 39: generic C version, but have all the usual assembler tricks and obscurities
! 40: introduced for speed.
! 41:
! 42:
! 43:
! 44: BASECASE SQUARE
! 45:
! 46: A square can be done in roughly half the time of a multiply, by using the
! 47: fact that the cross products above and below the diagonal are the same. In
! 48: practice squaring isn't 2x faster than multiplying, but it's always faster
! 49: by a decent factor.
! 50:
! 51: u0 u1 u2 u3 u4
! 52: +---+---+---+---+---+
! 53: u0 | x | | | | |
! 54: +---+---+---+---+---+
! 55: u1 | | x | | | |
! 56: +---+---+---+---+---+
! 57: u2 | | | x | | |
! 58: +---+---+---+---+---+
! 59: u3 | | | | x | |
! 60: +---+---+---+---+---+
! 61: u4 | | | | | x |
! 62: +---+---+---+---+---+
! 63:
! 64: The basic algorithm is to calculate a triangle of products below the
! 65: diagonal, double it (left shift by one bit), and add in the products on the
! 66: diagonal. This can be seen in mpn/generic/sqr_basecase.c. Again the
! 67: assembler implementations take essentially this same approach.
! 68:
! 69:
! 70:
! 71:
! 72: KARATSUBA MULTIPLY
! 73:
! 74: The Karatsuba multiplication algorithm is described in Knuth volume 2
! 75: section 4.3.3 part A, and in other texts like Geddes et al. section 4.3
! 76: (reference below). A brief description is given here.
! 77:
! 78: The Karatsuba algorithm treats its inputs x and y as each split in two parts
! 79: of equal length (or the most significant part one limb shorter if N is odd).
! 80:
! 81: high low
! 82: +----------+----------+
! 83: | x1 | x0 |
! 84: +----------+----------+
! 85:
! 86: +----------+----------+
! 87: | y1 | y0 |
! 88: +----------+----------+
! 89:
! 90: Let b be the power of 2 where the split occurs, ie. if x0 is k limbs (y0 the
! 91: same) then b=2^(k*mp_bits_per_limb). Then x=x1*b+x0 and y=y1*b+y0, and the
! 92: following holds,
! 93:
! 94: x*y = (b^2+b)*x1*y1 + b*(x1-x0)*(y1-y0) + (b+1)*x0*y0
! 95:
! 96: This formula means doing only three multiplies of (N/2)x(N/2) limbs, whereas
! 97: a basecase multiply of NxN limbs is roughly equivalent to four multiplies of
! 98: (N/2)x(N/2).
! 99:
! 100: The factors (b^2+b) etc in the formula look worse than they are,
! 101: representing simply the positions where the products must be added in.
! 102:
! 103: high low
! 104: +--------+--------+ +--------+--------+
! 105: | x1*y1 | | x0*y0 |
! 106: +--------+--------+ +--------+--------+
! 107: +--------+--------+
! 108: | x1*y1 |
! 109: +--------+--------+
! 110: +--------+--------+
! 111: | x0*y0 |
! 112: +--------+--------+
! 113: +--------+--------+
! 114: | (x1-x0)*(y1-y0) |
! 115: +--------+--------+
! 116:
! 117: The term (x1-x0)*(y1-y0) can be negative, meaning a subtraction, but the
! 118: final result is of course always positive.
! 119:
! 120: The use of three multiplies of N/2 limbs each leads to an asymptotic speed
! 121: O(N^1.585). (The exponent is log(3)/log(2).) This is a big improvement
! 122: over the basecase multiply at O(N^2) and the algorithmic advantage soon
! 123: overcomes the extra additions Karatsuba must perform.
! 124:
! 125:
! 126:
! 127:
! 128: KARATSUBA SQUARE
! 129:
! 130: A square is very similar to a multiply, but with x==y the formula reduces to
! 131: an equivalent with three squares,
! 132:
! 133: x^2 = (b^2+b)*x1^2 + b*(x1-x0)^2 + (b+1)*x0^2
! 134:
! 135: The final result is accumulated from those three squares the same way as for
! 136: the three multiplies above. The middle term (x1-x0)^2 however is now always
! 137: positive.
! 138:
! 139:
! 140:
! 141:
! 142: TOOM-COOK 3-WAY MULTIPLY
! 143:
! 144: The Karatsuba formula is part of a general approach to splitting inputs
! 145: leading to both Toom-Cook and FFT algorithms. A description of Toom-Cook
! 146: can be found in Knuth volume 2 section 4.3.3, with an example 3-way
! 147: calculation after Theorem A.
! 148:
! 149: Toom-Cook 3-way treats the operands as split into 3 pieces of equal size (or
! 150: the most significant part 1 or 2 limbs shorter than the others).
! 151:
! 152: high low
! 153: +----------+----------+----------+
! 154: | x2 | x1 | x0 |
! 155: +----------+----------+----------+
! 156:
! 157: +----------+----------+----------+
! 158: | y2 | y1 | y0 |
! 159: +----------+----------+----------+
! 160:
! 161: These parts are treated as the coefficients of two polynomials
! 162:
! 163: X(t) = x2*t^2 + x1*t + x0
! 164: Y(t) = y2*t^2 + y1*t + y0
! 165:
! 166: Again let b equal the power of 2 which is the size of the x0, x1, y0 and y1
! 167: pieces, ie. if they're k limbs each then b=2^(k*mp_bits_per_limb). With
! 168: this then x=X(b) and y=Y(b).
! 169:
! 170: Let a polynomial W(t)=X(t)*Y(t) and suppose it's coefficients are
! 171:
! 172: W(t) = w4*t^4 + w3*t^3 + w2*t^2 + w1*t + w0
! 173:
! 174: The w[i] are going to be determined, and when they are they'll give the
! 175: final result using w=W(b), since x*y=X(b)*Y(b)=W(b). The coefficients will
! 176: each be roughly b^2, so the final W(b) will be an addition like,
! 177:
! 178: high low
! 179: +-------+-------+
! 180: | w4 |
! 181: +-------+-------+
! 182: +--------+-------+
! 183: | w3 |
! 184: +--------+-------+
! 185: +--------+-------+
! 186: | w2 |
! 187: +--------+-------+
! 188: +--------+-------+
! 189: | w1 |
! 190: +--------+-------+
! 191: +-------+-------+
! 192: | w0 |
! 193: +-------+-------+
! 194: -------------------------------------------------
! 195:
! 196:
! 197: The w[i] coefficients could be formed by a simple set of cross products,
! 198: like w4=x2*x2, w3=x2*y1+x1*y2, w2=x2*y0+x1*y1+x0*y2 etc, but this would need
! 199: all nine x[i]*y[j] for i,j=0,1,2, and would be equivalent merely to a
! 200: basecase multiply. Instead the following approach is used.
! 201:
! 202: X(t) and Y(t) are evaluated and multiplied at 5 points, giving values of
! 203: W(t) at those points. The points used can be chosen in various ways, but in
! 204: GMP the following are used
! 205:
! 206: t=0 meaning x0*y0, which gives w0 immediately
! 207: t=2 meaning (4*x2+2*x1*x0)*(4*y2+2*y1+y0)
! 208: t=1 meaning (x2+x1+x0)*(y2+y1+y0)
! 209: t=1/2 meaning (x2+2*x1+4*x0)*(y2+2*y1+4*y0)
! 210: t=inf meaning x2*y2, which gives w4 immediately
! 211:
! 212: At t=1/2 the value calculated is actually 4*X(1/2)*Y(1/2), giving a value
! 213: for 16*W(1/2) (this is always an integer). At t=inf the value is actually
! 214: X(t)*Y(t)/t^2 in the limit as t approaches infinity, but it's much easier to
! 215: think of that as simply x2*y2 giving w4 immediately (much like at t=0 x0*y0
! 216: gives w0 immediately).
! 217:
! 218: Now each of the points substituted into W(t)=w4*t^4+...+w0 gives a linear
! 219: combination of the w[i] coefficients, and the value of those combinations
! 220: has just been calculated.
! 221:
! 222: W(0) = w0
! 223: 16*W(1/2) = w4 + 2*w3 + 4*w2 + 8*w1 + 16*w0
! 224: W(1) = w4 + w3 + w2 + w1 + w0
! 225: W(2) = 16*w4 + 8*w3 + 4*w2 + 2*w1 + w0
! 226: W(inf) = w4
! 227:
! 228: This is a set of five equations in five unknowns, and some elementary linear
! 229: algebra quickly isolates each w[i], by subtracting multiples of one equation
! 230: from another.
! 231:
! 232: In the code the set of five values W(0),...,W(inf) will represent those
! 233: certain linear combinations. By adding or subtracting one from another as
! 234: necessary, values which are each w[i] alone are arrived at. This involves
! 235: only a few subtractions of small multiples (some of which are powers of 2),
! 236: and so is very fast. A couple of divisions remain by powers of 2 and one
! 237: division by 3 (or by 6 rather), and that last uses the special fast
! 238: mpn_divexact_by3.
! 239:
! 240: In the code the values w4, w2 and w0 are formed in the destination, and w3
! 241: and w1 are added to them. There's an extra word at the high end of w3, w2
! 242: and w1 that are handled separately. With A=w0,B=w1,...,E=w4, the additions
! 243: are as follows.
! 244:
! 245: high low
! 246: +-------+-------+-------+-------+-------+-------+
! 247: | E | C | A |
! 248: +-------+-------+-------+-------+-------+-------+
! 249: +------+-------++------+-------+
! 250: | D || B |
! 251: +------+-------++------+-------+
! 252: -- -- --
! 253: |tD| |tC| |tB|
! 254: -- -- --
! 255: -------------------------------------------------
! 256:
! 257:
! 258: The conversion of W(t) values to the coefficients is called interpolation.
! 259: A polynomial of degree 5 like W(t) is uniquely determined by values known at
! 260: 5 different points. The points can be chosen to make the linear equations
! 261: come out with a convenient set of steps for isolating the w[i]s.
! 262:
! 263: In mpn/generic/mul_n.c the interpolate3() routine performs the
! 264: interpolation. The open-coded one-pass version may be a bit hard to
! 265: understand, the steps performed can be better seen in the USE_MORE_MPN
! 266: version.
! 267:
! 268: The use of five multiplies of N/3 limbs each leads to an asymptotic speed
! 269: O(N^1.465). (The exponent is log(5)/log(3).) This is an improvement over
! 270: Karatsuba at O(N^1.585), though Toom-Cook does more work in the evaluation
! 271: and interpolation and so it's only above a certain size that Toom-Cook
! 272: realizes its advantage.
! 273:
! 274: The formula given above for the Karatsuba algorithm has an equivalent for
! 275: Toom-Cook 3-way, involving only five multiplies, but this would be
! 276: complicated and unenlightening.
! 277:
! 278: An alternate view of Toom-Cook 3-way can be found in Zuras (reference
! 279: below). He uses a vector to represent the x and y splits and a matrix
! 280: multiplication for the evaluation and interpolation stages. The matrix
! 281: inverses are not meant to be actually used, and they have elements with
! 282: values much greater than in fact arise in the interpolation steps. The
! 283: diagram shown for the 3-way is attractive, but again it doesn't have to be
! 284: implemented that way and for example with a bit of rearrangement just one
! 285: division by 6 (not two of them) can be done.
! 286:
! 287:
! 288:
! 289:
! 290: TOOM-COOK 3-WAY SQUARE
! 291:
! 292: Toom-Cook squaring follows the same procedure as multiplication, but there's
! 293: only one X(t) and it's evaluated at 5 points, and those values squared to
! 294: give values of W(t). The interpolation is then identical, and in fact the
! 295: same interpolate3() subroutine is used for both squaring and multiplying.
! 296:
! 297:
! 298:
! 299:
! 300: FFT MULTIPLY
! 301:
! 302: At large to very large sizes a Fermat style FFT is used, meaning an FFT in a
! 303: ring of integers modulo 2^M+1. This is asymptotically fast, but overheads
! 304: mean it's only worthwhile for large products.
! 305:
! 306: Some brief notes are given here. Full explanations can be found in various
! 307: texts. Knuth section 4.3.3 part C describes the method, but using complex
! 308: numbers. In the references below Schonhage and Strassen is the original
! 309: paper, and Pollard gives some of the mathematics for a finite field as used
! 310: here.
! 311:
! 312: The FFT does its multiplication modulo 2^N+1, but by choosing
! 313: N>=bits(A)+bits(B), a full product A*B is obtained. The algorithm splits
! 314: the inputs into 2^k pieces, for a chosen k, and will recursively perform 2^k
! 315: pointwise multiplications modulo 2^M+1, where M=N/2^k. N must be a multiple
! 316: of 2^k. Those multiplications are either done by recursing into a further
! 317: FFT, or by a plain toom3 etc multiplication, whichever is optimal at the
! 318: resultant size. Note that in the current implementation M is always a
! 319: multiple of the limb size.
! 320:
! 321: The steps leading to the pointwise products are like the evaluation and
! 322: interpolation stages of the Karatsuba and Toom-Cook algorithms, but use a
! 323: convolution, which can be efficiently calculated because 2^(2N/2^k) is a
! 324: 2^k'th root of unity.
! 325:
! 326: As operand sizes gets bigger, bigger splits are used. Each time a bigger k
! 327: is used some multiplying is effectively swapped for some shifts, adds and
! 328: overheads. A table of thresholds gives the points where a k+1 FFT is faster
! 329: than a k FFT. A separate threshold gives the point where a mod 2^N+1 FFT
! 330: first becomes faster than a plain multiply of that size, and this normally
! 331: happens in the k=4 range. A further threshold gives the point where an
! 332: N/2xN/2 multiply done with an FFT mod 2^N+1 is faster than a plain multiply
! 333: of that size, and this is normally in the k=7 or k=8 range.
! 334:
! 335:
! 336:
! 337:
! 338: TOOM-COOK N-WAY (NOT USED)
! 339:
! 340: The 3-way Toom-Cook procedure described above generalizes to split into an
! 341: arbitrary number of pieces, as per Knuth volume 2 section 4.3.3 algorithm C.
! 342: Some code implementing this for GMP exists, but is not present since it has
! 343: yet to prove worthwhile. The notes here are merely for interest.
! 344:
! 345: In general a split into r+1 pieces will be made, and evaluations and
! 346: pointwise multiplications done at 2*r+1 points. So a 4-way split does 7
! 347: pointwise multiplies, 5-way does 9, etc.
! 348:
! 349: Asymptotically an r+1-way algorithm is O(n^(log(2*r+1)/log(r+1)). Only the
! 350: pointwise multiplications count towards big O() complexity, but the time
! 351: spent in the evaluate and interpolate stages grows with r and has a
! 352: significant practical impact, with the asymptotic advantage of each r
! 353: realized only at bigger and bigger sizes.
! 354:
! 355: Knuth algorithm C presents a version evaluating at points 0,1,2,...,2*r, but
! 356: exercise 4 uses -r,...,0,...,r and the latter saves some small multiplies in
! 357: the evaluate stage (or rather trades them for additions), and has a further
! 358: saving of nearly half the interpolate steps. The answer to the exercise
! 359: doesn't give full details of the interpolate, but essentially the idea is to
! 360: separate odd and even final coefficients and then perform algorithm C steps
! 361: C7 and C8 on them separately. The multipliers and divisors at step C7 then
! 362: become j^2 and 2*t*j-j*j respectively.
! 363:
! 364: In the references below, Zuras presents 3-way and 4-way Toom-Cook methods,
! 365: and compares them to small order FFTs. Hollerbach presents an N-way
! 366: algorithm from first principles. Bernstein presents the N-way algorithm in
! 367: a style that facilitates comparing its mathematics to other multiplication
! 368: algorithms.
! 369:
! 370:
! 371:
! 372:
! 373: REFERENCES
! 374:
! 375: "Algorithms for Computer Algebra", Keith O. Geddes, Stephen R. Czapor,
! 376: George Labahn, Kluwer Academic Publishers, 1992, ISBN 0-7923-9259-0.
! 377:
! 378: "Schnelle Multiplikation grosser Zahlen", by Arnold Schonhage and Volker
! 379: Strassen, Computing 7, p. 281-292, 1971.
! 380:
! 381: "The Fast Fourier Transform in a Finite Field", J.M. Pollard, Mathematics of
! 382: Computation, vol 25, num 114, April 1971.
! 383:
! 384: "On Squaring and Multiplying Large Integers", Dan Zuras, ARITH-11: IEEE
! 385: Symposium on Computer Arithmetic, 1993, pages 260 to 271. And reprinted as
! 386: "More on Multiplying and Squaring Large Integers", IEEE Transactions on
! 387: Computers, August 1994.
! 388:
! 389: "Fast Multiplication & Division of Very Large Numbers", Uwe Hollerbach, post
! 390: to sci.math.research, Jan 1996, archived at Swarthmore,
! 391: http://forum.swarthmore.edu/epigone/sci.math.research/zhouyimpzimp/x1ybdbxz5w4v@forum.swarthmore.edu
! 392:
! 393: "Multidigit Multiplication for Mathematicians", Daniel J. Bernstein,
! 394: preprint available at http://koobera.math.uic.edu/www/papers. (Every known
! 395: multiplication technique, many references.)
! 396:
! 397:
! 398:
! 399:
! 400: ----------------
! 401: Local variables:
! 402: mode: text
! 403: fill-column: 76
! 404: End:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to multiplication CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / doc