Annotation of OpenXM/src/kan96xx/gmp-2.0.2-ssh-2/mpn/generic/sqrtrem.c, Revision 1.1.1.1
1.1 takayama 1: /* mpn_sqrtrem (root_ptr, rem_ptr, op_ptr, op_size)
2:
3: Write the square root of {OP_PTR, OP_SIZE} at ROOT_PTR.
4: Write the remainder at REM_PTR, if REM_PTR != NULL.
5: Return the size of the remainder.
6: (The size of the root is always half of the size of the operand.)
7:
8: OP_PTR and ROOT_PTR may not point to the same object.
9: OP_PTR and REM_PTR may point to the same object.
10:
11: If REM_PTR is NULL, only the root is computed and the return value of
12: the function is 0 if OP is a perfect square, and *any* non-zero number
13: otherwise.
14:
15: Copyright (C) 1993, 1994, 1996 Free Software Foundation, Inc.
16:
17: This file is part of the GNU MP Library.
18:
19: The GNU MP Library is free software; you can redistribute it and/or modify
20: it under the terms of the GNU Library General Public License as published by
21: the Free Software Foundation; either version 2 of the License, or (at your
22: option) any later version.
23:
24: The GNU MP Library is distributed in the hope that it will be useful, but
25: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
26: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
27: License for more details.
28:
29: You should have received a copy of the GNU Library General Public License
30: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
31: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
32: MA 02111-1307, USA. */
33:
34: /* This code is just correct if "unsigned char" has at least 8 bits. It
35: doesn't help to use CHAR_BIT from limits.h, as the real problem is
36: the static arrays. */
37:
38: #include "gmp.h"
39: #include "gmp-impl.h"
40: #include "longlong.h"
41:
42: /* Square root algorithm:
43:
44: 1. Shift OP (the input) to the left an even number of bits s.t. there
45: are an even number of words and either (or both) of the most
46: significant bits are set. This way, sqrt(OP) has exactly half as
47: many words as OP, and has its most significant bit set.
48:
49: 2. Get a 9-bit approximation to sqrt(OP) using the pre-computed tables.
50: This approximation is used for the first single-precision
51: iterations of Newton's method, yielding a full-word approximation
52: to sqrt(OP).
53:
54: 3. Perform multiple-precision Newton iteration until we have the
55: exact result. Only about half of the input operand is used in
56: this calculation, as the square root is perfectly determinable
57: from just the higher half of a number. */
58:
59: /* Define this macro for IEEE P854 machines with a fast sqrt instruction. */
60: #if defined __GNUC__ && ! defined __SOFT_FLOAT
61:
62: #if defined __sparc__
63: #define SQRT(a) \
64: ({ \
65: double __sqrt_res; \
66: asm ("fsqrtd %1,%0" : "=f" (__sqrt_res) : "f" (a)); \
67: __sqrt_res; \
68: })
69: #endif
70:
71: #if defined __HAVE_68881__
72: #define SQRT(a) \
73: ({ \
74: double __sqrt_res; \
75: asm ("fsqrtx %1,%0" : "=f" (__sqrt_res) : "f" (a)); \
76: __sqrt_res; \
77: })
78: #endif
79:
80: #if defined __hppa
81: #define SQRT(a) \
82: ({ \
83: double __sqrt_res; \
84: asm ("fsqrt,dbl %1,%0" : "=fx" (__sqrt_res) : "fx" (a)); \
85: __sqrt_res; \
86: })
87: #endif
88:
89: #if defined _ARCH_PWR2
90: #define SQRT(a) \
91: ({ \
92: double __sqrt_res; \
93: asm ("fsqrt %0,%1" : "=f" (__sqrt_res) : "f" (a)); \
94: __sqrt_res; \
95: })
96: #endif
97:
98: #endif
99:
100: #ifndef SQRT
101:
102: /* Tables for initial approximation of the square root. These are
103: indexed with bits 1-8 of the operand for which the square root is
104: calculated, where bit 0 is the most significant non-zero bit. I.e.
105: the most significant one-bit is not used, since that per definition
106: is one. Likewise, the tables don't return the highest bit of the
107: result. That bit must be inserted by or:ing the returned value with
108: 0x100. This way, we get a 9-bit approximation from 8-bit tables! */
109:
110: /* Table to be used for operands with an even total number of bits.
111: (Exactly as in the decimal system there are similarities between the
112: square root of numbers with the same initial digits and an even
113: difference in the total number of digits. Consider the square root
114: of 1, 10, 100, 1000, ...) */
115: static unsigned char even_approx_tab[256] =
116: {
117: 0x6a, 0x6a, 0x6b, 0x6c, 0x6c, 0x6d, 0x6e, 0x6e,
118: 0x6f, 0x70, 0x71, 0x71, 0x72, 0x73, 0x73, 0x74,
119: 0x75, 0x75, 0x76, 0x77, 0x77, 0x78, 0x79, 0x79,
120: 0x7a, 0x7b, 0x7b, 0x7c, 0x7d, 0x7d, 0x7e, 0x7f,
121: 0x80, 0x80, 0x81, 0x81, 0x82, 0x83, 0x83, 0x84,
122: 0x85, 0x85, 0x86, 0x87, 0x87, 0x88, 0x89, 0x89,
123: 0x8a, 0x8b, 0x8b, 0x8c, 0x8d, 0x8d, 0x8e, 0x8f,
124: 0x8f, 0x90, 0x90, 0x91, 0x92, 0x92, 0x93, 0x94,
125: 0x94, 0x95, 0x96, 0x96, 0x97, 0x97, 0x98, 0x99,
126: 0x99, 0x9a, 0x9b, 0x9b, 0x9c, 0x9c, 0x9d, 0x9e,
127: 0x9e, 0x9f, 0xa0, 0xa0, 0xa1, 0xa1, 0xa2, 0xa3,
128: 0xa3, 0xa4, 0xa4, 0xa5, 0xa6, 0xa6, 0xa7, 0xa7,
129: 0xa8, 0xa9, 0xa9, 0xaa, 0xaa, 0xab, 0xac, 0xac,
130: 0xad, 0xad, 0xae, 0xaf, 0xaf, 0xb0, 0xb0, 0xb1,
131: 0xb2, 0xb2, 0xb3, 0xb3, 0xb4, 0xb5, 0xb5, 0xb6,
132: 0xb6, 0xb7, 0xb7, 0xb8, 0xb9, 0xb9, 0xba, 0xba,
133: 0xbb, 0xbb, 0xbc, 0xbd, 0xbd, 0xbe, 0xbe, 0xbf,
134: 0xc0, 0xc0, 0xc1, 0xc1, 0xc2, 0xc2, 0xc3, 0xc3,
135: 0xc4, 0xc5, 0xc5, 0xc6, 0xc6, 0xc7, 0xc7, 0xc8,
136: 0xc9, 0xc9, 0xca, 0xca, 0xcb, 0xcb, 0xcc, 0xcc,
137: 0xcd, 0xce, 0xce, 0xcf, 0xcf, 0xd0, 0xd0, 0xd1,
138: 0xd1, 0xd2, 0xd3, 0xd3, 0xd4, 0xd4, 0xd5, 0xd5,
139: 0xd6, 0xd6, 0xd7, 0xd7, 0xd8, 0xd9, 0xd9, 0xda,
140: 0xda, 0xdb, 0xdb, 0xdc, 0xdc, 0xdd, 0xdd, 0xde,
141: 0xde, 0xdf, 0xe0, 0xe0, 0xe1, 0xe1, 0xe2, 0xe2,
142: 0xe3, 0xe3, 0xe4, 0xe4, 0xe5, 0xe5, 0xe6, 0xe6,
143: 0xe7, 0xe7, 0xe8, 0xe8, 0xe9, 0xea, 0xea, 0xeb,
144: 0xeb, 0xec, 0xec, 0xed, 0xed, 0xee, 0xee, 0xef,
145: 0xef, 0xf0, 0xf0, 0xf1, 0xf1, 0xf2, 0xf2, 0xf3,
146: 0xf3, 0xf4, 0xf4, 0xf5, 0xf5, 0xf6, 0xf6, 0xf7,
147: 0xf7, 0xf8, 0xf8, 0xf9, 0xf9, 0xfa, 0xfa, 0xfb,
148: 0xfb, 0xfc, 0xfc, 0xfd, 0xfd, 0xfe, 0xfe, 0xff,
149: };
150:
151: /* Table to be used for operands with an odd total number of bits.
152: (Further comments before previous table.) */
153: static unsigned char odd_approx_tab[256] =
154: {
155: 0x00, 0x00, 0x00, 0x01, 0x01, 0x02, 0x02, 0x03,
156: 0x03, 0x04, 0x04, 0x05, 0x05, 0x06, 0x06, 0x07,
157: 0x07, 0x08, 0x08, 0x09, 0x09, 0x0a, 0x0a, 0x0b,
158: 0x0b, 0x0c, 0x0c, 0x0d, 0x0d, 0x0e, 0x0e, 0x0f,
159: 0x0f, 0x10, 0x10, 0x10, 0x11, 0x11, 0x12, 0x12,
160: 0x13, 0x13, 0x14, 0x14, 0x15, 0x15, 0x16, 0x16,
161: 0x16, 0x17, 0x17, 0x18, 0x18, 0x19, 0x19, 0x1a,
162: 0x1a, 0x1b, 0x1b, 0x1b, 0x1c, 0x1c, 0x1d, 0x1d,
163: 0x1e, 0x1e, 0x1f, 0x1f, 0x20, 0x20, 0x20, 0x21,
164: 0x21, 0x22, 0x22, 0x23, 0x23, 0x23, 0x24, 0x24,
165: 0x25, 0x25, 0x26, 0x26, 0x27, 0x27, 0x27, 0x28,
166: 0x28, 0x29, 0x29, 0x2a, 0x2a, 0x2a, 0x2b, 0x2b,
167: 0x2c, 0x2c, 0x2d, 0x2d, 0x2d, 0x2e, 0x2e, 0x2f,
168: 0x2f, 0x30, 0x30, 0x30, 0x31, 0x31, 0x32, 0x32,
169: 0x32, 0x33, 0x33, 0x34, 0x34, 0x35, 0x35, 0x35,
170: 0x36, 0x36, 0x37, 0x37, 0x37, 0x38, 0x38, 0x39,
171: 0x39, 0x39, 0x3a, 0x3a, 0x3b, 0x3b, 0x3b, 0x3c,
172: 0x3c, 0x3d, 0x3d, 0x3d, 0x3e, 0x3e, 0x3f, 0x3f,
173: 0x40, 0x40, 0x40, 0x41, 0x41, 0x41, 0x42, 0x42,
174: 0x43, 0x43, 0x43, 0x44, 0x44, 0x45, 0x45, 0x45,
175: 0x46, 0x46, 0x47, 0x47, 0x47, 0x48, 0x48, 0x49,
176: 0x49, 0x49, 0x4a, 0x4a, 0x4b, 0x4b, 0x4b, 0x4c,
177: 0x4c, 0x4c, 0x4d, 0x4d, 0x4e, 0x4e, 0x4e, 0x4f,
178: 0x4f, 0x50, 0x50, 0x50, 0x51, 0x51, 0x51, 0x52,
179: 0x52, 0x53, 0x53, 0x53, 0x54, 0x54, 0x54, 0x55,
180: 0x55, 0x56, 0x56, 0x56, 0x57, 0x57, 0x57, 0x58,
181: 0x58, 0x59, 0x59, 0x59, 0x5a, 0x5a, 0x5a, 0x5b,
182: 0x5b, 0x5b, 0x5c, 0x5c, 0x5d, 0x5d, 0x5d, 0x5e,
183: 0x5e, 0x5e, 0x5f, 0x5f, 0x60, 0x60, 0x60, 0x61,
184: 0x61, 0x61, 0x62, 0x62, 0x62, 0x63, 0x63, 0x63,
185: 0x64, 0x64, 0x65, 0x65, 0x65, 0x66, 0x66, 0x66,
186: 0x67, 0x67, 0x67, 0x68, 0x68, 0x68, 0x69, 0x69,
187: };
188: #endif
189:
190:
191: mp_size_t
192: #if __STDC__
193: mpn_sqrtrem (mp_ptr root_ptr, mp_ptr rem_ptr, mp_srcptr op_ptr, mp_size_t op_size)
194: #else
195: mpn_sqrtrem (root_ptr, rem_ptr, op_ptr, op_size)
196: mp_ptr root_ptr;
197: mp_ptr rem_ptr;
198: mp_srcptr op_ptr;
199: mp_size_t op_size;
200: #endif
201: {
202: /* R (root result) */
203: mp_ptr rp; /* Pointer to least significant word */
204: mp_size_t rsize; /* The size in words */
205:
206: /* T (OP shifted to the left a.k.a. normalized) */
207: mp_ptr tp; /* Pointer to least significant word */
208: mp_size_t tsize; /* The size in words */
209: mp_ptr t_end_ptr; /* Pointer right beyond most sign. word */
210: mp_limb_t t_high0, t_high1; /* The two most significant words */
211:
212: /* TT (temporary for numerator/remainder) */
213: mp_ptr ttp; /* Pointer to least significant word */
214:
215: /* X (temporary for quotient in main loop) */
216: mp_ptr xp; /* Pointer to least significant word */
217: mp_size_t xsize; /* The size in words */
218:
219: unsigned cnt;
220: mp_limb_t initial_approx; /* Initially made approximation */
221: mp_size_t tsizes[BITS_PER_MP_LIMB]; /* Successive calculation precisions */
222: mp_size_t tmp;
223: mp_size_t i;
224:
225: mp_limb_t cy_limb;
226: TMP_DECL (marker);
227:
228: /* If OP is zero, both results are zero. */
229: if (op_size == 0)
230: return 0;
231:
232: count_leading_zeros (cnt, op_ptr[op_size - 1]);
233: tsize = op_size;
234: if ((tsize & 1) != 0)
235: {
236: cnt += BITS_PER_MP_LIMB;
237: tsize++;
238: }
239:
240: rsize = tsize / 2;
241: rp = root_ptr;
242:
243: TMP_MARK (marker);
244:
245: /* Shift OP an even number of bits into T, such that either the most or
246: the second most significant bit is set, and such that the number of
247: words in T becomes even. This way, the number of words in R=sqrt(OP)
248: is exactly half as many as in OP, and the most significant bit of R
249: is set.
250:
251: Also, the initial approximation is simplified by this up-shifted OP.
252:
253: Finally, the Newtonian iteration which is the main part of this
254: program performs division by R. The fast division routine expects
255: the divisor to be "normalized" in exactly the sense of having the
256: most significant bit set. */
257:
258: tp = (mp_ptr) TMP_ALLOC (tsize * BYTES_PER_MP_LIMB);
259:
260: if ((cnt & ~1) % BITS_PER_MP_LIMB != 0)
261: t_high0 = mpn_lshift (tp + cnt / BITS_PER_MP_LIMB, op_ptr, op_size,
262: (cnt & ~1) % BITS_PER_MP_LIMB);
263: else
264: MPN_COPY (tp + cnt / BITS_PER_MP_LIMB, op_ptr, op_size);
265:
266: if (cnt >= BITS_PER_MP_LIMB)
267: tp[0] = 0;
268:
269: t_high0 = tp[tsize - 1];
270: t_high1 = tp[tsize - 2]; /* Never stray. TSIZE is >= 2. */
271:
272: /* Is there a fast sqrt instruction defined for this machine? */
273: #ifdef SQRT
274: {
275: initial_approx = SQRT (t_high0 * 2.0
276: * ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1))
277: + t_high1);
278: /* If t_high0,,t_high1 is big, the result in INITIAL_APPROX might have
279: become incorrect due to overflow in the conversion from double to
280: mp_limb_t above. It will typically be zero in that case, but might be
281: a small number on some machines. The most significant bit of
282: INITIAL_APPROX should be set, so that bit is a good overflow
283: indication. */
284: if ((mp_limb_signed_t) initial_approx >= 0)
285: initial_approx = ~(mp_limb_t)0;
286: }
287: #else
288: /* Get a 9 bit approximation from the tables. The tables expect to
289: be indexed with the 8 high bits right below the highest bit.
290: Also, the highest result bit is not returned by the tables, and
291: must be or:ed into the result. The scheme gives 9 bits of start
292: approximation with just 256-entry 8 bit tables. */
293:
294: if ((cnt & 1) == 0)
295: {
296: /* The most sign bit of t_high0 is set. */
297: initial_approx = t_high0 >> (BITS_PER_MP_LIMB - 8 - 1);
298: initial_approx &= 0xff;
299: initial_approx = even_approx_tab[initial_approx];
300: }
301: else
302: {
303: /* The most significant bit of T_HIGH0 is unset,
304: the second most significant is set. */
305: initial_approx = t_high0 >> (BITS_PER_MP_LIMB - 8 - 2);
306: initial_approx &= 0xff;
307: initial_approx = odd_approx_tab[initial_approx];
308: }
309: initial_approx |= 0x100;
310: initial_approx <<= BITS_PER_MP_LIMB - 8 - 1;
311:
312: /* Perform small precision Newtonian iterations to get a full word
313: approximation. For small operands, these iteration will make the
314: entire job. */
315: if (t_high0 == ~(mp_limb_t)0)
316: initial_approx = t_high0;
317: else
318: {
319: mp_limb_t quot;
320:
321: if (t_high0 >= initial_approx)
322: initial_approx = t_high0 + 1;
323:
324: /* First get about 18 bits with pure C arithmetics. */
325: quot = t_high0 / (initial_approx >> BITS_PER_MP_LIMB/2) << BITS_PER_MP_LIMB/2;
326: initial_approx = (initial_approx + quot) / 2;
327: initial_approx |= (mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1);
328:
329: /* Now get a full word by one (or for > 36 bit machines) several
330: iterations. */
331: for (i = 16; i < BITS_PER_MP_LIMB; i <<= 1)
332: {
333: mp_limb_t ignored_remainder;
334:
335: udiv_qrnnd (quot, ignored_remainder,
336: t_high0, t_high1, initial_approx);
337: initial_approx = (initial_approx + quot) / 2;
338: initial_approx |= (mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1);
339: }
340: }
341: #endif
342:
343: rp[0] = initial_approx;
344: rsize = 1;
345:
346: #ifdef DEBUG
347: printf ("\n\nT = ");
348: mpn_dump (tp, tsize);
349: #endif
350:
351: if (tsize > 2)
352: {
353: /* Determine the successive precisions to use in the iteration. We
354: minimize the precisions, beginning with the highest (i.e. last
355: iteration) to the lowest (i.e. first iteration). */
356:
357: xp = (mp_ptr) TMP_ALLOC (tsize * BYTES_PER_MP_LIMB);
358: ttp = (mp_ptr) TMP_ALLOC (tsize * BYTES_PER_MP_LIMB);
359:
360: t_end_ptr = tp + tsize;
361:
362: tmp = tsize / 2;
363: for (i = 0;; i++)
364: {
365: tsize = (tmp + 1) / 2;
366: if (tmp == tsize)
367: break;
368: tsizes[i] = tsize + tmp;
369: tmp = tsize;
370: }
371:
372: /* Main Newton iteration loop. For big arguments, most of the
373: time is spent here. */
374:
375: /* It is possible to do a great optimization here. The successive
376: divisors in the mpn_divmod call below has more and more leading
377: words equal to its predecessor. Therefore the beginning of
378: each division will repeat the same work as did the last
379: division. If we could guarantee that the leading words of two
380: consecutive divisors are the same (i.e. in this case, a later
381: divisor has just more digits at the end) it would be a simple
382: matter of just using the old remainder of the last division in
383: a subsequent division, to take care of this optimization. This
384: idea would surely make a difference even for small arguments. */
385:
386: /* Loop invariants:
387:
388: R <= shiftdown_to_same_size(floor(sqrt(OP))) < R + 1.
389: X - 1 < shiftdown_to_same_size(floor(sqrt(OP))) <= X.
390: R <= shiftdown_to_same_size(X). */
391:
392: while (--i >= 0)
393: {
394: mp_limb_t cy;
395: #ifdef DEBUG
396: mp_limb_t old_least_sign_r = rp[0];
397: mp_size_t old_rsize = rsize;
398:
399: printf ("R = ");
400: mpn_dump (rp, rsize);
401: #endif
402: tsize = tsizes[i];
403:
404: /* Need to copy the numerator into temporary space, as
405: mpn_divmod overwrites its numerator argument with the
406: remainder (which we currently ignore). */
407: MPN_COPY (ttp, t_end_ptr - tsize, tsize);
408: cy = mpn_divmod (xp, ttp, tsize, rp, rsize);
409: xsize = tsize - rsize;
410:
411: #ifdef DEBUG
412: printf ("X =%d ", cy);
413: mpn_dump (xp, xsize);
414: #endif
415:
416: /* Add X and R with the most significant limbs aligned,
417: temporarily ignoring at least one limb at the low end of X. */
418: tmp = xsize - rsize;
419: cy += mpn_add_n (xp + tmp, rp, xp + tmp, rsize);
420:
421: /* If T begins with more than 2 x BITS_PER_MP_LIMB of ones, we get
422: intermediate roots that'd need an extra bit. We don't want to
423: handle that since it would make the subsequent divisor
424: non-normalized, so round such roots down to be only ones in the
425: current precision. */
426: if (cy == 2)
427: {
428: mp_size_t j;
429: for (j = xsize; j >= 0; j--)
430: xp[j] = ~(mp_limb_t)0;
431: }
432:
433: /* Divide X by 2 and put the result in R. This is the new
434: approximation. Shift in the carry from the addition. */
435: mpn_rshift (rp, xp, xsize, 1);
436: rp[xsize - 1] |= ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1));
437: rsize = xsize;
438: #ifdef DEBUG
439: if (old_least_sign_r != rp[rsize - old_rsize])
440: printf (">>>>>>>> %d: %0*lX, %0*lX <<<<<<<<\n",
441: i, 2 * BYTES_PER_MP_LIMB, old_least_sign_r,
442: 2 * BYTES_PER_MP_LIMB, rp[rsize - old_rsize]);
443: #endif
444: }
445: }
446:
447: #ifdef DEBUG
448: printf ("(final) R = ");
449: mpn_dump (rp, rsize);
450: #endif
451:
452: /* We computed the square root of OP * 2**(2*floor(cnt/2)).
453: This has resulted in R being 2**floor(cnt/2) to large.
454: Shift it down here to fix that. */
455: if (cnt / 2 != 0)
456: {
457: mpn_rshift (rp, rp, rsize, cnt/2);
458: rsize -= rp[rsize - 1] == 0;
459: }
460:
461: /* Calculate the remainder. */
462: mpn_mul_n (tp, rp, rp, rsize);
463: tsize = rsize + rsize;
464: tsize -= tp[tsize - 1] == 0;
465: if (op_size < tsize
466: || (op_size == tsize && mpn_cmp (op_ptr, tp, op_size) < 0))
467: {
468: /* R is too large. Decrement it. */
469:
470: /* These operations can't overflow. */
471: cy_limb = mpn_sub_n (tp, tp, rp, rsize);
472: cy_limb += mpn_sub_n (tp, tp, rp, rsize);
473: mpn_sub_1 (tp + rsize, tp + rsize, tsize - rsize, cy_limb);
474: mpn_add_1 (tp, tp, tsize, (mp_limb_t) 1);
475:
476: mpn_sub_1 (rp, rp, rsize, (mp_limb_t) 1);
477:
478: #ifdef DEBUG
479: printf ("(adjusted) R = ");
480: mpn_dump (rp, rsize);
481: #endif
482: }
483:
484: if (rem_ptr != NULL)
485: {
486: cy_limb = mpn_sub (rem_ptr, op_ptr, op_size, tp, tsize);
487: MPN_NORMALIZE (rem_ptr, op_size);
488: TMP_FREE (marker);
489: return op_size;
490: }
491: else
492: {
493: int res;
494: res = op_size != tsize || mpn_cmp (op_ptr, tp, op_size);
495: TMP_FREE (marker);
496: return res;
497: }
498: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to sqrtrem.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / kan96xx / gmp-2.0.2-ssh-2 / mpn / generic