Annotation of OpenXM_contrib/gmp/gmp.info-6, Revision 1.1.1.1
1.1 ohara 1: This is gmp.info, produced by makeinfo version 4.2 from gmp.texi.
2:
3: This manual describes how to install and use the GNU multiple precision
4: arithmetic library, version 4.1.2.
5:
6: Copyright 1991, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000,
7: 2001, 2002 Free Software Foundation, Inc.
8:
9: Permission is granted to copy, distribute and/or modify this
10: document under the terms of the GNU Free Documentation License, Version
11: 1.1 or any later version published by the Free Software Foundation;
12: with no Invariant Sections, with the Front-Cover Texts being "A GNU
13: Manual", and with the Back-Cover Texts being "You have freedom to copy
14: and modify this GNU Manual, like GNU software". A copy of the license
15: is included in *Note GNU Free Documentation License::.
16: INFO-DIR-SECTION GNU libraries
17: START-INFO-DIR-ENTRY
18: * gmp: (gmp). GNU Multiple Precision Arithmetic Library.
19: END-INFO-DIR-ENTRY
20:
21:
22: File: gmp.info, Node: FFT Multiplication, Next: Other Multiplication, Prev: Toom-Cook 3-Way Multiplication, Up: Multiplication Algorithms
23:
24: FFT Multiplication
25: ------------------
26:
27: At large to very large sizes a Fermat style FFT multiplication is
28: used, following Scho"nhage and Strassen (*note References::).
29: Descriptions of FFTs in various forms can be found in many textbooks,
30: for instance Knuth section 4.3.3 part C or Lipson chapter IX. A brief
31: description of the form used in GMP is given here.
32:
33: The multiplication done is x*y mod 2^N+1, for a given N. A full
34: product x*y is obtained by choosing N>=bits(x)+bits(y) and padding x
35: and y with high zero limbs. The modular product is the native form for
36: the algorithm, so padding to get a full product is unavoidable.
37:
38: The algorithm follows a split, evaluate, pointwise multiply,
39: interpolate and combine similar to that described above for Karatsuba
40: and Toom-3. A k parameter controls the split, with an FFT-k splitting
41: into 2^k pieces of M=N/2^k bits each. N must be a multiple of
42: (2^k)*mp_bits_per_limb so the split falls on limb boundaries, avoiding
43: bit shifts in the split and combine stages.
44:
45: The evaluations, pointwise multiplications, and interpolation, are
46: all done modulo 2^N'+1 where N' is 2M+k+3 rounded up to a multiple of
47: 2^k and of `mp_bits_per_limb'. The results of interpolation will be
48: the following negacyclic convolution of the input pieces, and the
49: choice of N' ensures these sums aren't truncated.
50:
51: ---
52: \ b
53: w[n] = / (-1) * x[i] * y[j]
54: ---
55: i+j==b*2^k+n
56: b=0,1
57:
58: The points used for the evaluation are g^i for i=0 to 2^k-1 where
59: g=2^(2N'/2^k). g is a 2^k'th root of unity mod 2^N'+1, which produces
60: necessary cancellations at the interpolation stage, and it's also a
61: power of 2 so the fast fourier transforms used for the evaluation and
62: interpolation do only shifts, adds and negations.
63:
64: The pointwise multiplications are done modulo 2^N'+1 and either
65: recurse into a further FFT or use a plain multiplication (Toom-3,
66: Karatsuba or basecase), whichever is optimal at the size N'. The
67: interpolation is an inverse fast fourier transform. The resulting set
68: of sums of x[i]*y[j] are added at appropriate offsets to give the final
69: result.
70:
71: Squaring is the same, but x is the only input so it's one transform
72: at the evaluate stage and the pointwise multiplies are squares. The
73: interpolation is the same.
74:
75: For a mod 2^N+1 product, an FFT-k is an O(N^(k/(k-1))) algorithm,
76: the exponent representing 2^k recursed modular multiplies each
77: 1/2^(k-1) the size of the original. Each successive k is an asymptotic
78: improvement, but overheads mean each is only faster at bigger and
79: bigger sizes. In the code, `MUL_FFT_TABLE' and `SQR_FFT_TABLE' are the
80: thresholds where each k is used. Each new k effectively swaps some
81: multiplying for some shifts, adds and overheads.
82:
83: A mod 2^N+1 product can be formed with a normal NxN->2N bit multiply
84: plus a subtraction, so an FFT and Toom-3 etc can be compared directly.
85: A k=4 FFT at O(N^1.333) can be expected to be the first faster than
86: Toom-3 at O(N^1.465). In practice this is what's found, with
87: `MUL_FFT_MODF_THRESHOLD' and `SQR_FFT_MODF_THRESHOLD' being between 300
88: and 1000 limbs, depending on the CPU. So far it's been found that only
89: very large FFTs recurse into pointwise multiplies above these sizes.
90:
91: When an FFT is to give a full product, the change of N to 2N doesn't
92: alter the theoretical complexity for a given k, but for the purposes of
93: considering where an FFT might be first used it can be assumed that the
94: FFT is recursing into a normal multiply and that on that basis it's
95: doing 2^k recursed multiplies each 1/2^(k-2) the size of the inputs,
96: making it O(N^(k/(k-2))). This would mean k=7 at O(N^1.4) would be the
97: first FFT faster than Toom-3. In practice `MUL_FFT_THRESHOLD' and
98: `SQR_FFT_THRESHOLD' have been found to be in the k=8 range, somewhere
99: between 3000 and 10000 limbs.
100:
101: The way N is split into 2^k pieces and then 2M+k+3 is rounded up to
102: a multiple of 2^k and `mp_bits_per_limb' means that when
103: 2^k>=mp_bits_per_limb the effective N is a multiple of 2^(2k-1) bits.
104: The +k+3 means some values of N just under such a multiple will be
105: rounded to the next. The complexity calculations above assume that a
106: favourable size is used, meaning one which isn't padded through
107: rounding, and it's also assumed that the extra +k+3 bits are negligible
108: at typical FFT sizes.
109:
110: The practical effect of the 2^(2k-1) constraint is to introduce a
111: step-effect into measured speeds. For example k=8 will round N up to a
112: multiple of 32768 bits, so for a 32-bit limb there'll be 512 limb
113: groups of sizes for which `mpn_mul_n' runs at the same speed. Or for
114: k=9 groups of 2048 limbs, k=10 groups of 8192 limbs, etc. In practice
115: it's been found each k is used at quite small multiples of its size
116: constraint and so the step effect is quite noticeable in a time versus
117: size graph.
118:
119: The threshold determinations currently measure at the mid-points of
120: size steps, but this is sub-optimal since at the start of a new step it
121: can happen that it's better to go back to the previous k for a while.
122: Something more sophisticated for `MUL_FFT_TABLE' and `SQR_FFT_TABLE'
123: will be needed.
124:
125:
126: File: gmp.info, Node: Other Multiplication, Prev: FFT Multiplication, Up: Multiplication Algorithms
127:
128: Other Multiplication
129: --------------------
130:
131: The 3-way Toom-Cook algorithm described above (*note Toom-Cook 3-Way
132: Multiplication::) generalizes to split into an arbitrary number of
133: pieces, as per Knuth section 4.3.3 algorithm C. This is not currently
134: used, though it's possible a Toom-4 might fit in between Toom-3 and the
135: FFTs. The notes here are merely for interest.
136:
137: In general a split into r+1 pieces is made, and evaluations and
138: pointwise multiplications done at 2*r+1 points. A 4-way split does 7
139: pointwise multiplies, 5-way does 9, etc. Asymptotically an (r+1)-way
140: algorithm is O(N^(log(2*r+1)/log(r+1))). Only the pointwise
141: multiplications count towards big-O complexity, but the time spent in
142: the evaluate and interpolate stages grows with r and has a significant
143: practical impact, with the asymptotic advantage of each r realized only
144: at bigger and bigger sizes. The overheads grow as O(N*r), whereas in
145: an r=2^k FFT they grow only as O(N*log(r)).
146:
147: Knuth algorithm C evaluates at points 0,1,2,...,2*r, but exercise 4
148: uses -r,...,0,...,r and the latter saves some small multiplies in the
149: evaluate stage (or rather trades them for additions), and has a further
150: saving of nearly half the interpolate steps. The idea is to separate
151: odd and even final coefficients and then perform algorithm C steps C7
152: and C8 on them separately. The divisors at step C7 become j^2 and the
153: multipliers at C8 become 2*t*j-j^2.
154:
155: Splitting odd and even parts through positive and negative points
156: can be thought of as using -1 as a square root of unity. If a 4th root
157: of unity was available then a further split and speedup would be
158: possible, but no such root exists for plain integers. Going to complex
159: integers with i=sqrt(-1) doesn't help, essentially because in cartesian
160: form it takes three real multiplies to do a complex multiply. The
161: existence of 2^k'th roots of unity in a suitable ring or field lets the
162: fast fourier transform keep splitting and get to O(N*log(r)).
163:
164: Floating point FFTs use complex numbers approximating Nth roots of
165: unity. Some processors have special support for such FFTs. But these
166: are not used in GMP since it's very difficult to guarantee an exact
167: result (to some number of bits). An occasional difference of 1 in the
168: last bit might not matter to a typical signal processing algorithm, but
169: is of course of vital importance to GMP.
170:
171:
172: File: gmp.info, Node: Division Algorithms, Next: Greatest Common Divisor Algorithms, Prev: Multiplication Algorithms, Up: Algorithms
173:
174: Division Algorithms
175: ===================
176:
177: * Menu:
178:
179: * Single Limb Division::
180: * Basecase Division::
181: * Divide and Conquer Division::
182: * Exact Division::
183: * Exact Remainder::
184: * Small Quotient Division::
185:
186:
187: File: gmp.info, Node: Single Limb Division, Next: Basecase Division, Prev: Division Algorithms, Up: Division Algorithms
188:
189: Single Limb Division
190: --------------------
191:
192: Nx1 division is implemented using repeated 2x1 divisions from high
193: to low, either with a hardware divide instruction or a multiplication by
194: inverse, whichever is best on a given CPU.
195:
196: The multiply by inverse follows section 8 of "Division by Invariant
197: Integers using Multiplication" by Granlund and Montgomery (*note
198: References::) and is implemented as `udiv_qrnnd_preinv' in
199: `gmp-impl.h'. The idea is to have a fixed-point approximation to 1/d
200: (see `invert_limb') and then multiply by the high limb (plus one bit)
201: of the dividend to get a quotient q. With d normalized (high bit set),
202: q is no more than 1 too small. Subtracting q*d from the dividend gives
203: a remainder, and reveals whether q or q-1 is correct.
204:
205: The result is a division done with two multiplications and four or
206: five arithmetic operations. On CPUs with low latency multipliers this
207: can be much faster than a hardware divide, though the cost of
208: calculating the inverse at the start may mean it's only better on
209: inputs bigger than say 4 or 5 limbs.
210:
211: When a divisor must be normalized, either for the generic C
212: `__udiv_qrnnd_c' or the multiply by inverse, the division performed is
213: actually a*2^k by d*2^k where a is the dividend and k is the power
214: necessary to have the high bit of d*2^k set. The bit shifts for the
215: dividend are usually accomplished "on the fly" meaning by extracting
216: the appropriate bits at each step. Done this way the quotient limbs
217: come out aligned ready to store. When only the remainder is wanted, an
218: alternative is to take the dividend limbs unshifted and calculate r = a
219: mod d*2^k followed by an extra final step r*2^k mod d*2^k. This can
220: help on CPUs with poor bit shifts or few registers.
221:
222: The multiply by inverse can be done two limbs at a time. The
223: calculation is basically the same, but the inverse is two limbs and the
224: divisor treated as if padded with a low zero limb. This means more
225: work, since the inverse will need a 2x2 multiply, but the four 1x1s to
226: do that are independent and can therefore be done partly or wholly in
227: parallel. Likewise for a 2x1 calculating q*d. The net effect is to
228: process two limbs with roughly the same two multiplies worth of latency
229: that one limb at a time gives. This extends to 3 or 4 limbs at a time,
230: though the extra work to apply the inverse will almost certainly soon
231: reach the limits of multiplier throughput.
232:
233: A similar approach in reverse can be taken to process just half a
234: limb at a time if the divisor is only a half limb. In this case the
235: 1x1 multiply for the inverse effectively becomes two (1/2)x1 for each
236: limb, which can be a saving on CPUs with a fast half limb multiply, or
237: in fact if the only multiply is a half limb, and especially if it's not
238: pipelined.
239:
240:
241: File: gmp.info, Node: Basecase Division, Next: Divide and Conquer Division, Prev: Single Limb Division, Up: Division Algorithms
242:
243: Basecase Division
244: -----------------
245:
246: Basecase NxM division is like long division done by hand, but in base
247: 2^mp_bits_per_limb. See Knuth section 4.3.1 algorithm D, and
248: `mpn/generic/sb_divrem_mn.c'.
249:
250: Briefly stated, while the dividend remains larger than the divisor,
251: a high quotient limb is formed and the Nx1 product q*d subtracted at
252: the top end of the dividend. With a normalized divisor (most
253: significant bit set), each quotient limb can be formed with a 2x1
254: division and a 1x1 multiplication plus some subtractions. The 2x1
255: division is by the high limb of the divisor and is done either with a
256: hardware divide or a multiply by inverse (the same as in *Note Single
257: Limb Division::) whichever is faster. Such a quotient is sometimes one
258: too big, requiring an addback of the divisor, but that happens rarely.
259:
260: With Q=N-M being the number of quotient limbs, this is an O(Q*M)
261: algorithm and will run at a speed similar to a basecase QxM
262: multiplication, differing in fact only in the extra multiply and divide
263: for each of the Q quotient limbs.
264:
265:
266: File: gmp.info, Node: Divide and Conquer Division, Next: Exact Division, Prev: Basecase Division, Up: Division Algorithms
267:
268: Divide and Conquer Division
269: ---------------------------
270:
271: For divisors larger than `DIV_DC_THRESHOLD', division is done by
272: dividing. Or to be precise by a recursive divide and conquer algorithm
273: based on work by Moenck and Borodin, Jebelean, and Burnikel and Ziegler
274: (*note References::).
275:
276: The algorithm consists essentially of recognising that a 2NxN
277: division can be done with the basecase division algorithm (*note
278: Basecase Division::), but using N/2 limbs as a base, not just a single
279: limb. This way the multiplications that arise are (N/2)x(N/2) and can
280: take advantage of Karatsuba and higher multiplication algorithms (*note
281: Multiplication Algorithms::). The "digits" of the quotient are formed
282: by recursive Nx(N/2) divisions.
283:
284: If the (N/2)x(N/2) multiplies are done with a basecase multiplication
285: then the work is about the same as a basecase division, but with more
286: function call overheads and with some subtractions separated from the
287: multiplies. These overheads mean that it's only when N/2 is above
288: `MUL_KARATSUBA_THRESHOLD' that divide and conquer is of use.
289:
290: `DIV_DC_THRESHOLD' is based on the divisor size N, so it will be
291: somewhere above twice `MUL_KARATSUBA_THRESHOLD', but how much above
292: depends on the CPU. An optimized `mpn_mul_basecase' can lower
293: `DIV_DC_THRESHOLD' a little by offering a ready-made advantage over
294: repeated `mpn_submul_1' calls.
295:
296: Divide and conquer is asymptotically O(M(N)*log(N)) where M(N) is
297: the time for an NxN multiplication done with FFTs. The actual time is
298: a sum over multiplications of the recursed sizes, as can be seen near
299: the end of section 2.2 of Burnikel and Ziegler. For example, within
300: the Toom-3 range, divide and conquer is 2.63*M(N). With higher
301: algorithms the M(N) term improves and the multiplier tends to log(N).
302: In practice, at moderate to large sizes, a 2NxN division is about 2 to
303: 4 times slower than an NxN multiplication.
304:
305: Newton's method used for division is asymptotically O(M(N)) and
306: should therefore be superior to divide and conquer, but it's believed
307: this would only be for large to very large N.
308:
309:
310: File: gmp.info, Node: Exact Division, Next: Exact Remainder, Prev: Divide and Conquer Division, Up: Division Algorithms
311:
312: Exact Division
313: --------------
314:
315: A so-called exact division is when the dividend is known to be an
316: exact multiple of the divisor. Jebelean's exact division algorithm
317: uses this knowledge to make some significant optimizations (*note
318: References::).
319:
320: The idea can be illustrated in decimal for example with 368154
321: divided by 543. Because the low digit of the dividend is 4, the low
322: digit of the quotient must be 8. This is arrived at from 4*7 mod 10,
323: using the fact 7 is the modular inverse of 3 (the low digit of the
324: divisor), since 3*7 == 1 mod 10. So 8*543=4344 can be subtracted from
325: the dividend leaving 363810. Notice the low digit has become zero.
326:
327: The procedure is repeated at the second digit, with the next
328: quotient digit 7 (7 == 1*7 mod 10), subtracting 7*543=3801, leaving
329: 325800. And finally at the third digit with quotient digit 6 (8*7 mod
330: 10), subtracting 6*543=3258 leaving 0. So the quotient is 678.
331:
332: Notice however that the multiplies and subtractions don't need to
333: extend past the low three digits of the dividend, since that's enough
334: to determine the three quotient digits. For the last quotient digit no
335: subtraction is needed at all. On a 2NxN division like this one, only
336: about half the work of a normal basecase division is necessary.
337:
338: For an NxM exact division producing Q=N-M quotient limbs, the saving
339: over a normal basecase division is in two parts. Firstly, each of the
340: Q quotient limbs needs only one multiply, not a 2x1 divide and
341: multiply. Secondly, the crossproducts are reduced when Q>M to
342: Q*M-M*(M+1)/2, or when Q<=M to Q*(Q-1)/2. Notice the savings are
343: complementary. If Q is big then many divisions are saved, or if Q is
344: small then the crossproducts reduce to a small number.
345:
346: The modular inverse used is calculated efficiently by
347: `modlimb_invert' in `gmp-impl.h'. This does four multiplies for a
348: 32-bit limb, or six for a 64-bit limb. `tune/modlinv.c' has some
349: alternate implementations that might suit processors better at bit
350: twiddling than multiplying.
351:
352: The sub-quadratic exact division described by Jebelean in "Exact
353: Division with Karatsuba Complexity" is not currently implemented. It
354: uses a rearrangement similar to the divide and conquer for normal
355: division (*note Divide and Conquer Division::), but operating from low
356: to high. A further possibility not currently implemented is
357: "Bidirectional Exact Integer Division" by Krandick and Jebelean which
358: forms quotient limbs from both the high and low ends of the dividend,
359: and can halve once more the number of crossproducts needed in a 2NxN
360: division.
361:
362: A special case exact division by 3 exists in `mpn_divexact_by3',
363: supporting Toom-3 multiplication and `mpq' canonicalizations. It forms
364: quotient digits with a multiply by the modular inverse of 3 (which is
365: `0xAA..AAB') and uses two comparisons to determine a borrow for the next
366: limb. The multiplications don't need to be on the dependent chain, as
367: long as the effect of the borrows is applied. Only a few optimized
368: assembler implementations currently exist.
369:
370:
371: File: gmp.info, Node: Exact Remainder, Next: Small Quotient Division, Prev: Exact Division, Up: Division Algorithms
372:
373: Exact Remainder
374: ---------------
375:
376: If the exact division algorithm is done with a full subtraction at
377: each stage and the dividend isn't a multiple of the divisor, then low
378: zero limbs are produced but with a remainder in the high limbs. For
379: dividend a, divisor d, quotient q, and b = 2^mp_bits_per_limb, then this
380: remainder r is of the form
381:
382: a = q*d + r*b^n
383:
384: n represents the number of zero limbs produced by the subtractions,
385: that being the number of limbs produced for q. r will be in the range
386: 0<=r<d and can be viewed as a remainder, but one shifted up by a factor
387: of b^n.
388:
389: Carrying out full subtractions at each stage means the same number
390: of cross products must be done as a normal division, but there's still
391: some single limb divisions saved. When d is a single limb some
392: simplifications arise, providing good speedups on a number of
393: processors.
394:
395: `mpn_bdivmod', `mpn_divexact_by3', `mpn_modexact_1_odd' and the
396: `redc' function in `mpz_powm' differ subtly in how they return r,
397: leading to some negations in the above formula, but all are essentially
398: the same.
399:
400: Clearly r is zero when a is a multiple of d, and this leads to
401: divisibility or congruence tests which are potentially more efficient
402: than a normal division.
403:
404: The factor of b^n on r can be ignored in a GCD when d is odd, hence
405: the use of `mpn_bdivmod' in `mpn_gcd', and the use of
406: `mpn_modexact_1_odd' by `mpn_gcd_1' and `mpz_kronecker_ui' etc (*note
407: Greatest Common Divisor Algorithms::).
408:
409: Montgomery's REDC method for modular multiplications uses operands
410: of the form of x*b^-n and y*b^-n and on calculating (x*b^-n)*(y*b^-n)
411: uses the factor of b^n in the exact remainder to reach a product in the
412: same form (x*y)*b^-n (*note Modular Powering Algorithm::).
413:
414: Notice that r generally gives no useful information about the
415: ordinary remainder a mod d since b^n mod d could be anything. If
416: however b^n == 1 mod d, then r is the negative of the ordinary
417: remainder. This occurs whenever d is a factor of b^n-1, as for example
418: with 3 in `mpn_divexact_by3'. Other such factors include 5, 17 and
419: 257, but no particular use has been found for this.
420:
421:
422: File: gmp.info, Node: Small Quotient Division, Prev: Exact Remainder, Up: Division Algorithms
423:
424: Small Quotient Division
425: -----------------------
426:
427: An NxM division where the number of quotient limbs Q=N-M is small
428: can be optimized somewhat.
429:
430: An ordinary basecase division normalizes the divisor by shifting it
431: to make the high bit set, shifting the dividend accordingly, and
432: shifting the remainder back down at the end of the calculation. This
433: is wasteful if only a few quotient limbs are to be formed. Instead a
434: division of just the top 2*Q limbs of the dividend by the top Q limbs
435: of the divisor can be used to form a trial quotient. This requires
436: only those limbs normalized, not the whole of the divisor and dividend.
437:
438: A multiply and subtract then applies the trial quotient to the M-Q
439: unused limbs of the divisor and N-Q dividend limbs (which includes Q
440: limbs remaining from the trial quotient division). The starting trial
441: quotient can be 1 or 2 too big, but all cases of 2 too big and most
442: cases of 1 too big are detected by first comparing the most significant
443: limbs that will arise from the subtraction. An addback is done if the
444: quotient still turns out to be 1 too big.
445:
446: This whole procedure is essentially the same as one step of the
447: basecase algorithm done in a Q limb base, though with the trial
448: quotient test done only with the high limbs, not an entire Q limb
449: "digit" product. The correctness of this weaker test can be
450: established by following the argument of Knuth section 4.3.1 exercise
451: 20 but with the v2*q>b*r+u2 condition appropriately relaxed.
452:
453:
454: File: gmp.info, Node: Greatest Common Divisor Algorithms, Next: Powering Algorithms, Prev: Division Algorithms, Up: Algorithms
455:
456: Greatest Common Divisor
457: =======================
458:
459: * Menu:
460:
461: * Binary GCD::
462: * Accelerated GCD::
463: * Extended GCD::
464: * Jacobi Symbol::
465:
466:
467: File: gmp.info, Node: Binary GCD, Next: Accelerated GCD, Prev: Greatest Common Divisor Algorithms, Up: Greatest Common Divisor Algorithms
468:
469: Binary GCD
470: ----------
471:
472: At small sizes GMP uses an O(N^2) binary style GCD. This is
473: described in many textbooks, for example Knuth section 4.5.2 algorithm
474: B. It simply consists of successively reducing operands a and b using
475: gcd(a,b) = gcd(min(a,b),abs(a-b)), and also that if a and b are first
476: made odd then abs(a-b) is even and factors of two can be discarded.
477:
478: Variants like letting a-b become negative and doing a different next
479: step are of interest only as far as they suit particular CPUs, since on
480: small operands it's machine dependent factors that determine
481: performance.
482:
483: The Euclidean GCD algorithm, as per Knuth algorithms E and A,
484: reduces using a mod b but this has so far been found to be slower
485: everywhere. One reason the binary method does well is that the implied
486: quotient at each step is usually small, so often only one or two
487: subtractions are needed to get the same effect as a division.
488: Quotients 1, 2 and 3 for example occur 67.7% of the time, see Knuth
489: section 4.5.3 Theorem E.
490:
491: When the implied quotient is large, meaning b is much smaller than
492: a, then a division is worthwhile. This is the basis for the initial a
493: mod b reductions in `mpn_gcd' and `mpn_gcd_1' (the latter for both Nx1
494: and 1x1 cases). But after that initial reduction, big quotients occur
495: too rarely to make it worth checking for them.
496:
497:
498: File: gmp.info, Node: Accelerated GCD, Next: Extended GCD, Prev: Binary GCD, Up: Greatest Common Divisor Algorithms
499:
500: Accelerated GCD
501: ---------------
502:
503: For sizes above `GCD_ACCEL_THRESHOLD', GMP uses the Accelerated GCD
504: algorithm described independently by Weber and Jebelean (the latter as
505: the "Generalized Binary" algorithm), *note References::. This
506: algorithm is still O(N^2), but is much faster than the binary algorithm
507: since it does fewer multi-precision operations. It consists of
508: alternating the k-ary reduction by Sorenson, and a "dmod" exact
509: remainder reduction.
510:
511: For operands u and v the k-ary reduction replaces u with n*v-d*u
512: where n and d are single limb values chosen to give two trailing zero
513: limbs on that value, which can be stripped. n and d are calculated
514: using an algorithm similar to half of a two limb GCD (see `find_a' in
515: `mpn/generic/gcd.c').
516:
517: When u and v differ in size by more than a certain number of bits, a
518: dmod is performed to zero out bits at the low end of the larger. It
519: consists of an exact remainder style division applied to an appropriate
520: number of bits (*note Exact Division::, and *note Exact Remainder::).
521: This is faster than a k-ary reduction but useful only when the operands
522: differ in size. There's a dmod after each k-ary reduction, and if the
523: dmod leaves the operands still differing in size then it's repeated.
524:
525: The k-ary reduction step can introduce spurious factors into the GCD
526: calculated, and these are eliminated at the end by taking GCDs with the
527: original inputs gcd(u,gcd(v,g)) using the binary algorithm. Since g is
528: almost always small this takes very little time.
529:
530: At small sizes the algorithm needs a good implementation of
531: `find_a'. At larger sizes it's dominated by `mpn_addmul_1' applying n
532: and d.
533:
534:
535: File: gmp.info, Node: Extended GCD, Next: Jacobi Symbol, Prev: Accelerated GCD, Up: Greatest Common Divisor Algorithms
536:
537: Extended GCD
538: ------------
539:
540: The extended GCD calculates gcd(a,b) and also cofactors x and y
541: satisfying a*x+b*y=gcd(a,b). Lehmer's multi-step improvement of the
542: extended Euclidean algorithm is used. See Knuth section 4.5.2
543: algorithm L, and `mpn/generic/gcdext.c'. This is an O(N^2) algorithm.
544:
545: The multipliers at each step are found using single limb
546: calculations for sizes up to `GCDEXT_THRESHOLD', or double limb
547: calculations above that. The single limb code is faster but doesn't
548: produce full-limb multipliers, hence not making full use of the
549: `mpn_addmul_1' calls.
550:
551: When a CPU has a data-dependent multiplier, meaning one which is
552: faster on operands with fewer bits, the extra work in the double-limb
553: calculation might only save some looping overheads, leading to a large
554: `GCDEXT_THRESHOLD'.
555:
556: Currently the single limb calculation doesn't optimize for the small
557: quotients that often occur, and this can lead to unusually low values of
558: `GCDEXT_THRESHOLD', depending on the CPU.
559:
560: An analysis of double-limb calculations can be found in "A
561: Double-Digit Lehmer-Euclid Algorithm" by Jebelean (*note References::).
562: The code in GMP was developed independently.
563:
564: It should be noted that when a double limb calculation is used, it's
565: used for the whole of that GCD, it doesn't fall back to single limb
566: part way through. This is because as the algorithm proceeds, the
567: inputs a and b are reduced, but the cofactors x and y grow, so the
568: multipliers at each step are applied to a roughly constant total number
569: of limbs.
570:
571:
572: File: gmp.info, Node: Jacobi Symbol, Prev: Extended GCD, Up: Greatest Common Divisor Algorithms
573:
574: Jacobi Symbol
575: -------------
576:
577: `mpz_jacobi' and `mpz_kronecker' are currently implemented with a
578: simple binary algorithm similar to that described for the GCDs (*note
579: Binary GCD::). They're not very fast when both inputs are large.
580: Lehmer's multi-step improvement or a binary based multi-step algorithm
581: is likely to be better.
582:
583: When one operand fits a single limb, and that includes
584: `mpz_kronecker_ui' and friends, an initial reduction is done with
585: either `mpn_mod_1' or `mpn_modexact_1_odd', followed by the binary
586: algorithm on a single limb. The binary algorithm is well suited to a
587: single limb, and the whole calculation in this case is quite efficient.
588:
589: In all the routines sign changes for the result are accumulated
590: using some bit twiddling, avoiding table lookups or conditional jumps.
591:
592:
593: File: gmp.info, Node: Powering Algorithms, Next: Root Extraction Algorithms, Prev: Greatest Common Divisor Algorithms, Up: Algorithms
594:
595: Powering Algorithms
596: ===================
597:
598: * Menu:
599:
600: * Normal Powering Algorithm::
601: * Modular Powering Algorithm::
602:
603:
604: File: gmp.info, Node: Normal Powering Algorithm, Next: Modular Powering Algorithm, Prev: Powering Algorithms, Up: Powering Algorithms
605:
606: Normal Powering
607: ---------------
608:
609: Normal `mpz' or `mpf' powering uses a simple binary algorithm,
610: successively squaring and then multiplying by the base when a 1 bit is
611: seen in the exponent, as per Knuth section 4.6.3. The "left to right"
612: variant described there is used rather than algorithm A, since it's
613: just as easy and can be done with somewhat less temporary memory.
614:
615:
616: File: gmp.info, Node: Modular Powering Algorithm, Prev: Normal Powering Algorithm, Up: Powering Algorithms
617:
618: Modular Powering
619: ----------------
620:
621: Modular powering is implemented using a 2^k-ary sliding window
622: algorithm, as per "Handbook of Applied Cryptography" algorithm 14.85
623: (*note References::). k is chosen according to the size of the
624: exponent. Larger exponents use larger values of k, the choice being
625: made to minimize the average number of multiplications that must
626: supplement the squaring.
627:
628: The modular multiplies and squares use either a simple division or
629: the REDC method by Montgomery (*note References::). REDC is a little
630: faster, essentially saving N single limb divisions in a fashion similar
631: to an exact remainder (*note Exact Remainder::). The current REDC has
632: some limitations. It's only O(N^2) so above `POWM_THRESHOLD' division
633: becomes faster and is used. It doesn't attempt to detect small bases,
634: but rather always uses a REDC form, which is usually a full size
635: operand. And lastly it's only applied to odd moduli.
636:
637:
638: File: gmp.info, Node: Root Extraction Algorithms, Next: Radix Conversion Algorithms, Prev: Powering Algorithms, Up: Algorithms
639:
640: Root Extraction Algorithms
641: ==========================
642:
643: * Menu:
644:
645: * Square Root Algorithm::
646: * Nth Root Algorithm::
647: * Perfect Square Algorithm::
648: * Perfect Power Algorithm::
649:
650:
651: File: gmp.info, Node: Square Root Algorithm, Next: Nth Root Algorithm, Prev: Root Extraction Algorithms, Up: Root Extraction Algorithms
652:
653: Square Root
654: -----------
655:
656: Square roots are taken using the "Karatsuba Square Root" algorithm
657: by Paul Zimmermann (*note References::). This is expressed in a divide
658: and conquer form, but as noted in the paper it can also be viewed as a
659: discrete variant of Newton's method.
660:
661: In the Karatsuba multiplication range this is an O(1.5*M(N/2))
662: algorithm, where M(n) is the time to multiply two numbers of n limbs.
663: In the FFT multiplication range this grows to a bound of O(6*M(N/2)).
664: In practice a factor of about 1.5 to 1.8 is found in the Karatsuba and
665: Toom-3 ranges, growing to 2 or 3 in the FFT range.
666:
667: The algorithm does all its calculations in integers and the resulting
668: `mpn_sqrtrem' is used for both `mpz_sqrt' and `mpf_sqrt'. The extended
669: precision given by `mpf_sqrt_ui' is obtained by padding with zero limbs.
670:
671:
672: File: gmp.info, Node: Nth Root Algorithm, Next: Perfect Square Algorithm, Prev: Square Root Algorithm, Up: Root Extraction Algorithms
673:
674: Nth Root
675: --------
676:
677: Integer Nth roots are taken using Newton's method with the following
678: iteration, where A is the input and n is the root to be taken.
679:
680: 1 A
681: a[i+1] = - * ( --------- + (n-1)*a[i] )
682: n a[i]^(n-1)
683:
684: The initial approximation a[1] is generated bitwise by successively
685: powering a trial root with or without new 1 bits, aiming to be just
686: above the true root. The iteration converges quadratically when
687: started from a good approximation. When n is large more initial bits
688: are needed to get good convergence. The current implementation is not
689: particularly well optimized.
690:
691:
692: File: gmp.info, Node: Perfect Square Algorithm, Next: Perfect Power Algorithm, Prev: Nth Root Algorithm, Up: Root Extraction Algorithms
693:
694: Perfect Square
695: --------------
696:
697: `mpz_perfect_square_p' is able to quickly exclude most non-squares by
698: checking whether the input is a quadratic residue modulo some small
699: integers.
700:
701: The first test is modulo 256 which means simply examining the least
702: significant byte. Only 44 different values occur as the low byte of a
703: square, so 82.8% of non-squares can be immediately excluded. Similar
704: tests modulo primes from 3 to 29 exclude 99.5% of those remaining, or
705: if a limb is 64 bits then primes up to 53 are used, excluding 99.99%.
706: A single Nx1 remainder using `PP' from `gmp-impl.h' quickly gives all
707: these remainders.
708:
709: A square root must still be taken for any value that passes the
710: residue tests, to verify it's really a square and not one of the 0.086%
711: (or 0.000156% for 64 bits) non-squares that get through. *Note Square
712: Root Algorithm::.
713:
714:
715: File: gmp.info, Node: Perfect Power Algorithm, Prev: Perfect Square Algorithm, Up: Root Extraction Algorithms
716:
717: Perfect Power
718: -------------
719:
720: Detecting perfect powers is required by some factorization
721: algorithms. Currently `mpz_perfect_power_p' is implemented using
722: repeated Nth root extractions, though naturally only prime roots need
723: to be considered. (*Note Nth Root Algorithm::.)
724:
725: If a prime divisor p with multiplicity e can be found, then only
726: roots which are divisors of e need to be considered, much reducing the
727: work necessary. To this end divisibility by a set of small primes is
728: checked.
729:
730:
731: File: gmp.info, Node: Radix Conversion Algorithms, Next: Other Algorithms, Prev: Root Extraction Algorithms, Up: Algorithms
732:
733: Radix Conversion
734: ================
735:
736: Radix conversions are less important than other algorithms. A
737: program dominated by conversions should probably use a different data
738: representation.
739:
740: * Menu:
741:
742: * Binary to Radix::
743: * Radix to Binary::
744:
745:
746: File: gmp.info, Node: Binary to Radix, Next: Radix to Binary, Prev: Radix Conversion Algorithms, Up: Radix Conversion Algorithms
747:
748: Binary to Radix
749: ---------------
750:
751: Conversions from binary to a power-of-2 radix use a simple and fast
752: O(N) bit extraction algorithm.
753:
754: Conversions from binary to other radices use one of two algorithms.
755: Sizes below `GET_STR_PRECOMPUTE_THRESHOLD' use a basic O(N^2) method.
756: Repeated divisions by b^n are made, where b is the radix and n is the
757: biggest power that fits in a limb. But instead of simply using the
758: remainder r from such divisions, an extra divide step is done to give a
759: fractional limb representing r/b^n. The digits of r can then be
760: extracted using multiplications by b rather than divisions. Special
761: case code is provided for decimal, allowing multiplications by 10 to
762: optimize to shifts and adds.
763:
764: Above `GET_STR_PRECOMPUTE_THRESHOLD' a sub-quadratic algorithm is
765: used. For an input t, powers b^(n*2^i) of the radix are calculated,
766: until a power between t and sqrt(t) is reached. t is then divided by
767: that largest power, giving a quotient which is the digits above that
768: power, and a remainder which is those below. These two parts are in
769: turn divided by the second highest power, and so on recursively. When
770: a piece has been divided down to less than `GET_STR_DC_THRESHOLD'
771: limbs, the basecase algorithm described above is used.
772:
773: The advantage of this algorithm is that big divisions can make use
774: of the sub-quadratic divide and conquer division (*note Divide and
775: Conquer Division::), and big divisions tend to have less overheads than
776: lots of separate single limb divisions anyway. But in any case the
777: cost of calculating the powers b^(n*2^i) must first be overcome.
778:
779: `GET_STR_PRECOMPUTE_THRESHOLD' and `GET_STR_DC_THRESHOLD' represent
780: the same basic thing, the point where it becomes worth doing a big
781: division to cut the input in half. `GET_STR_PRECOMPUTE_THRESHOLD'
782: includes the cost of calculating the radix power required, whereas
783: `GET_STR_DC_THRESHOLD' assumes that's already available, which is the
784: case when recursing.
785:
786: Since the base case produces digits from least to most significant
787: but they want to be stored from most to least, it's necessary to
788: calculate in advance how many digits there will be, or at least be sure
789: not to underestimate that. For GMP the number of input bits is
790: multiplied by `chars_per_bit_exactly' from `mp_bases', rounding up.
791: The result is either correct or one too big.
792:
793: Examining some of the high bits of the input could increase the
794: chance of getting the exact number of digits, but an exact result every
795: time would not be practical, since in general the difference between
796: numbers 100... and 99... is only in the last few bits and the work to
797: identify 99... might well be almost as much as a full conversion.
798:
799: `mpf_get_str' doesn't currently use the algorithm described here, it
800: multiplies or divides by a power of b to move the radix point to the
801: just above the highest non-zero digit (or at worst one above that
802: location), then multiplies by b^n to bring out digits. This is O(N^2)
803: and is certainly not optimal.
804:
805: The r/b^n scheme described above for using multiplications to bring
806: out digits might be useful for more than a single limb. Some brief
807: experiments with it on the base case when recursing didn't give a
808: noticable improvement, but perhaps that was only due to the
809: implementation. Something similar would work for the sub-quadratic
810: divisions too, though there would be the cost of calculating a bigger
811: radix power.
812:
813: Another possible improvement for the sub-quadratic part would be to
814: arrange for radix powers that balanced the sizes of quotient and
815: remainder produced, ie. the highest power would be an b^(n*k)
816: approximately equal to sqrt(t), not restricted to a 2^i factor. That
817: ought to smooth out a graph of times against sizes, but may or may not
818: be a net speedup.
819:
820:
821: File: gmp.info, Node: Radix to Binary, Prev: Binary to Radix, Up: Radix Conversion Algorithms
822:
823: Radix to Binary
824: ---------------
825:
826: Conversions from a power-of-2 radix into binary use a simple and fast
827: O(N) bitwise concatenation algorithm.
828:
829: Conversions from other radices use one of two algorithms. Sizes
830: below `SET_STR_THRESHOLD' use a basic O(N^2) method. Groups of n
831: digits are converted to limbs, where n is the biggest power of the base
832: b which will fit in a limb, then those groups are accumulated into the
833: result by multiplying by b^n and adding. This saves multi-precision
834: operations, as per Knuth section 4.4 part E (*note References::). Some
835: special case code is provided for decimal, giving the compiler a chance
836: to optimize multiplications by 10.
837:
838: Above `SET_STR_THRESHOLD' a sub-quadratic algorithm is used. First
839: groups of n digits are converted into limbs. Then adjacent limbs are
840: combined into limb pairs with x*b^n+y, where x and y are the limbs.
841: Adjacent limb pairs are combined into quads similarly with x*b^(2n)+y.
842: This continues until a single block remains, that being the result.
843:
844: The advantage of this method is that the multiplications for each x
845: are big blocks, allowing Karatsuba and higher algorithms to be used.
846: But the cost of calculating the powers b^(n*2^i) must be overcome.
847: `SET_STR_THRESHOLD' usually ends up quite big, around 5000 digits, and
848: on some processors much bigger still.
849:
850: `SET_STR_THRESHOLD' is based on the input digits (and tuned for
851: decimal), though it might be better based on a limb count, so as to be
852: independent of the base. But that sort of count isn't used by the base
853: case and so would need some sort of initial calculation or estimate.
854:
855: The main reason `SET_STR_THRESHOLD' is so much bigger than the
856: corresponding `GET_STR_PRECOMPUTE_THRESHOLD' is that `mpn_mul_1' is
857: much faster than `mpn_divrem_1' (often by a factor of 10, or more).
858:
859:
860: File: gmp.info, Node: Other Algorithms, Next: Assembler Coding, Prev: Radix Conversion Algorithms, Up: Algorithms
861:
862: Other Algorithms
863: ================
864:
865: * Menu:
866:
867: * Factorial Algorithm::
868: * Binomial Coefficients Algorithm::
869: * Fibonacci Numbers Algorithm::
870: * Lucas Numbers Algorithm::
871:
872:
873: File: gmp.info, Node: Factorial Algorithm, Next: Binomial Coefficients Algorithm, Prev: Other Algorithms, Up: Other Algorithms
874:
875: Factorial
876: ---------
877:
878: Factorials n! are calculated by a simple product from 1 to n, but
879: arranged into certain sub-products.
880:
881: First as many factors as fit in a limb are accumulated, then two of
882: those multiplied to give a 2-limb product. When two 2-limb products
883: are ready they're multiplied to a 4-limb product, and when two 4-limbs
884: are ready they're multiplied to an 8-limb product, etc. A stack of
885: outstanding products is built up, with two of the same size multiplied
886: together when ready.
887:
888: Arranging for multiplications to have operands the same (or nearly
889: the same) size means the Karatsuba and higher multiplication algorithms
890: can be used. And even on sizes below the Karatsuba threshold an NxN
891: multiply will give a basecase multiply more to work on.
892:
893: An obvious improvement not currently implemented would be to strip
894: factors of 2 from the products and apply them at the end with a bit
895: shift. Another possibility would be to determine the prime
896: factorization of the result (which can be done easily), and use a
897: powering method, at each stage squaring then multiplying in those
898: primes with a 1 in their exponent at that point. The advantage would
899: be some multiplies turned into squares.
900:
901:
902: File: gmp.info, Node: Binomial Coefficients Algorithm, Next: Fibonacci Numbers Algorithm, Prev: Factorial Algorithm, Up: Other Algorithms
903:
904: Binomial Coefficients
905: ---------------------
906:
907: Binomial coefficients C(n,k) are calculated by first arranging k <=
908: n/2 using C(n,k) = C(n,n-k) if necessary, and then evaluating the
909: following product simply from i=2 to i=k.
910:
911: k (n-k+i)
912: C(n,k) = (n-k+1) * prod -------
913: i=2 i
914:
915: It's easy to show that each denominator i will divide the product so
916: far, so the exact division algorithm is used (*note Exact Division::).
917:
918: The numerators n-k+i and denominators i are first accumulated into
919: as many fit a limb, to save multi-precision operations, though for
920: `mpz_bin_ui' this applies only to the divisors, since n is an `mpz_t'
921: and n-k+i in general won't fit in a limb at all.
922:
923: An obvious improvement would be to strip factors of 2 from each
924: multiplier and divisor and count them separately, to be applied with a
925: bit shift at the end. Factors of 3 and perhaps 5 could even be handled
926: similarly. Another possibility, if n is not too big, would be to
927: determine the prime factorization of the result based on the factorials
928: involved, and power up those primes appropriately. This would help
929: most when k is near n/2.
930:
931:
932: File: gmp.info, Node: Fibonacci Numbers Algorithm, Next: Lucas Numbers Algorithm, Prev: Binomial Coefficients Algorithm, Up: Other Algorithms
933:
934: Fibonacci Numbers
935: -----------------
936:
937: The Fibonacci functions `mpz_fib_ui' and `mpz_fib2_ui' are designed
938: for calculating isolated F[n] or F[n],F[n-1] values efficiently.
939:
940: For small n, a table of single limb values in `__gmp_fib_table' is
941: used. On a 32-bit limb this goes up to F[47], or on a 64-bit limb up
942: to F[93]. For convenience the table starts at F[-1].
943:
944: Beyond the table, values are generated with a binary powering
945: algorithm, calculating a pair F[n] and F[n-1] working from high to low
946: across the bits of n. The formulas used are
947:
948: F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k
949: F[2k-1] = F[k]^2 + F[k-1]^2
950:
951: F[2k] = F[2k+1] - F[2k-1]
952:
953: At each step, k is the high b bits of n. If the next bit of n is 0
954: then F[2k],F[2k-1] is used, or if it's a 1 then F[2k+1],F[2k] is used,
955: and the process repeated until all bits of n are incorporated. Notice
956: these formulas require just two squares per bit of n.
957:
958: It'd be possible to handle the first few n above the single limb
959: table with simple additions, using the defining Fibonacci recurrence
960: F[k+1]=F[k]+F[k-1], but this is not done since it usually turns out to
961: be faster for only about 10 or 20 values of n, and including a block of
962: code for just those doesn't seem worthwhile. If they really mattered
963: it'd be better to extend the data table.
964:
965: Using a table avoids lots of calculations on small numbers, and
966: makes small n go fast. A bigger table would make more small n go fast,
967: it's just a question of balancing size against desired speed. For GMP
968: the code is kept compact, with the emphasis primarily on a good
969: powering algorithm.
970:
971: `mpz_fib2_ui' returns both F[n] and F[n-1], but `mpz_fib_ui' is only
972: interested in F[n]. In this case the last step of the algorithm can
973: become one multiply instead of two squares. One of the following two
974: formulas is used, according as n is odd or even.
975:
976: F[2k] = F[k]*(F[k]+2F[k-1])
977:
978: F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k
979:
980: F[2k+1] here is the same as above, just rearranged to be a multiply.
981: For interest, the 2*(-1)^k term both here and above can be applied
982: just to the low limb of the calculation, without a carry or borrow into
983: further limbs, which saves some code size. See comments with
984: `mpz_fib_ui' and the internal `mpn_fib2_ui' for how this is done.
985:
986:
987: File: gmp.info, Node: Lucas Numbers Algorithm, Prev: Fibonacci Numbers Algorithm, Up: Other Algorithms
988:
989: Lucas Numbers
990: -------------
991:
992: `mpz_lucnum2_ui' derives a pair of Lucas numbers from a pair of
993: Fibonacci numbers with the following simple formulas.
994:
995: L[k] = F[k] + 2*F[k-1]
996: L[k-1] = 2*F[k] - F[k-1]
997:
998: `mpz_lucnum_ui' is only interested in L[n], and some work can be
999: saved. Trailing zero bits on n can be handled with a single square
1000: each.
1001:
1002: L[2k] = L[k]^2 - 2*(-1)^k
1003:
1004: And the lowest 1 bit can be handled with one multiply of a pair of
1005: Fibonacci numbers, similar to what `mpz_fib_ui' does.
1006:
1007: L[2k+1] = 5*F[k-1]*(2*F[k]+F[k-1]) - 4*(-1)^k
1008:
1009:
1010: File: gmp.info, Node: Assembler Coding, Prev: Other Algorithms, Up: Algorithms
1011:
1012: Assembler Coding
1013: ================
1014:
1015: The assembler subroutines in GMP are the most significant source of
1016: speed at small to moderate sizes. At larger sizes algorithm selection
1017: becomes more important, but of course speedups in low level routines
1018: will still speed up everything proportionally.
1019:
1020: Carry handling and widening multiplies that are important for GMP
1021: can't be easily expressed in C. GCC `asm' blocks help a lot and are
1022: provided in `longlong.h', but hand coding low level routines invariably
1023: offers a speedup over generic C by a factor of anything from 2 to 10.
1024:
1025: * Menu:
1026:
1027: * Assembler Code Organisation::
1028: * Assembler Basics::
1029: * Assembler Carry Propagation::
1030: * Assembler Cache Handling::
1031: * Assembler Floating Point::
1032: * Assembler SIMD Instructions::
1033: * Assembler Software Pipelining::
1034: * Assembler Loop Unrolling::
1035:
1036:
1037: File: gmp.info, Node: Assembler Code Organisation, Next: Assembler Basics, Prev: Assembler Coding, Up: Assembler Coding
1038:
1039: Code Organisation
1040: -----------------
1041:
1042: The various `mpn' subdirectories contain machine-dependent code,
1043: written in C or assembler. The `mpn/generic' subdirectory contains
1044: default code, used when there's no machine-specific version of a
1045: particular file.
1046:
1047: Each `mpn' subdirectory is for an ISA family. Generally 32-bit and
1048: 64-bit variants in a family cannot share code and will have separate
1049: directories. Within a family further subdirectories may exist for CPU
1050: variants.
1051:
1052:
1053: File: gmp.info, Node: Assembler Basics, Next: Assembler Carry Propagation, Prev: Assembler Code Organisation, Up: Assembler Coding
1054:
1055: Assembler Basics
1056: ----------------
1057:
1058: `mpn_addmul_1' and `mpn_submul_1' are the most important routines
1059: for overall GMP performance. All multiplications and divisions come
1060: down to repeated calls to these. `mpn_add_n', `mpn_sub_n',
1061: `mpn_lshift' and `mpn_rshift' are next most important.
1062:
1063: On some CPUs assembler versions of the internal functions
1064: `mpn_mul_basecase' and `mpn_sqr_basecase' give significant speedups,
1065: mainly through avoiding function call overheads. They can also
1066: potentially make better use of a wide superscalar processor.
1067:
1068: The restrictions on overlaps between sources and destinations (*note
1069: Low-level Functions::) are designed to facilitate a variety of
1070: implementations. For example, knowing `mpn_add_n' won't have partly
1071: overlapping sources and destination means reading can be done far ahead
1072: of writing on superscalar processors, and loops can be vectorized on a
1073: vector processor, depending on the carry handling.
1074:
1075:
1076: File: gmp.info, Node: Assembler Carry Propagation, Next: Assembler Cache Handling, Prev: Assembler Basics, Up: Assembler Coding
1077:
1078: Carry Propagation
1079: -----------------
1080:
1081: The problem that presents most challenges in GMP is propagating
1082: carries from one limb to the next. In functions like `mpn_addmul_1' and
1083: `mpn_add_n', carries are the only dependencies between limb operations.
1084:
1085: On processors with carry flags, a straightforward CISC style `adc' is
1086: generally best. AMD K6 `mpn_addmul_1' however is an example of an
1087: unusual set of circumstances where a branch works out better.
1088:
1089: On RISC processors generally an add and compare for overflow is
1090: used. This sort of thing can be seen in `mpn/generic/aors_n.c'. Some
1091: carry propagation schemes require 4 instructions, meaning at least 4
1092: cycles per limb, but other schemes may use just 1 or 2. On wide
1093: superscalar processors performance may be completely determined by the
1094: number of dependent instructions between carry-in and carry-out for
1095: each limb.
1096:
1097: On vector processors good use can be made of the fact that a carry
1098: bit only very rarely propagates more than one limb. When adding a
1099: single bit to a limb, there's only a carry out if that limb was
1100: `0xFF...FF' which on random data will be only 1 in 2^mp_bits_per_limb.
1101: `mpn/cray/add_n.c' is an example of this, it adds all limbs in
1102: parallel, adds one set of carry bits in parallel and then only rarely
1103: needs to fall through to a loop propagating further carries.
1104:
1105: On the x86s, GCC (as of version 2.95.2) doesn't generate
1106: particularly good code for the RISC style idioms that are necessary to
1107: handle carry bits in C. Often conditional jumps are generated where
1108: `adc' or `sbb' forms would be better. And so unfortunately almost any
1109: loop involving carry bits needs to be coded in assembler for best
1110: results.
1111:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gmp.info-6 CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp