Annotation of OpenXM_contrib/gmp/mpn/generic/fib2_ui.c, Revision 1.1
1.1 ! ohara 1: /* mpn_fib2_ui -- calculate Fibonacci numbers.
! 2:
! 3: Copyright 2001, 2002 Free Software Foundation, Inc.
! 4:
! 5: This file is part of the GNU MP Library.
! 6:
! 7: The GNU MP Library is free software; you can redistribute it and/or modify
! 8: it under the terms of the GNU Lesser General Public License as published by
! 9: the Free Software Foundation; either version 2.1 of the License, or (at your
! 10: option) any later version.
! 11:
! 12: The GNU MP Library is distributed in the hope that it will be useful, but
! 13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 15: License for more details.
! 16:
! 17: You should have received a copy of the GNU Lesser General Public License
! 18: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 20: MA 02111-1307, USA. */
! 21:
! 22: #include <stdio.h>
! 23: #include "gmp.h"
! 24: #include "gmp-impl.h"
! 25: #include "longlong.h"
! 26:
! 27:
! 28: /* change this to "#define TRACE(x) x" for diagnostics */
! 29: #define TRACE(x)
! 30:
! 31:
! 32: /* The following table was generated by code at the end of this file. */
! 33:
! 34: const mp_limb_t
! 35: __gmp_fib_table[FIB_TABLE_LIMIT+2] = {
! 36:
! 37: #if GMP_NUMB_BITS >= 4
! 38: CNST_LIMB (0x1), /* -1 */
! 39: CNST_LIMB (0x0), /* 0 */
! 40: CNST_LIMB (0x1), /* 1 */
! 41: CNST_LIMB (0x1), /* 2 */
! 42: CNST_LIMB (0x2), /* 3 */
! 43: CNST_LIMB (0x3), /* 4 */
! 44: CNST_LIMB (0x5), /* 5 */
! 45: CNST_LIMB (0x8), /* 6 */
! 46: CNST_LIMB (0xD), /* 7 */
! 47: #endif
! 48: #if GMP_NUMB_BITS >= 8
! 49: CNST_LIMB (0x15), /* 8 */
! 50: CNST_LIMB (0x22), /* 9 */
! 51: CNST_LIMB (0x37), /* 10 */
! 52: CNST_LIMB (0x59), /* 11 */
! 53: CNST_LIMB (0x90), /* 12 */
! 54: CNST_LIMB (0xE9), /* 13 */
! 55: #endif
! 56: #if GMP_NUMB_BITS >= 16
! 57: CNST_LIMB (0x179), /* 14 */
! 58: CNST_LIMB (0x262), /* 15 */
! 59: CNST_LIMB (0x3DB), /* 16 */
! 60: CNST_LIMB (0x63D), /* 17 */
! 61: CNST_LIMB (0xA18), /* 18 */
! 62: CNST_LIMB (0x1055), /* 19 */
! 63: CNST_LIMB (0x1A6D), /* 20 */
! 64: CNST_LIMB (0x2AC2), /* 21 */
! 65: CNST_LIMB (0x452F), /* 22 */
! 66: CNST_LIMB (0x6FF1), /* 23 */
! 67: CNST_LIMB (0xB520), /* 24 */
! 68: #endif
! 69: #if GMP_NUMB_BITS >= 32
! 70: CNST_LIMB (0x12511), /* 25 */
! 71: CNST_LIMB (0x1DA31), /* 26 */
! 72: CNST_LIMB (0x2FF42), /* 27 */
! 73: CNST_LIMB (0x4D973), /* 28 */
! 74: CNST_LIMB (0x7D8B5), /* 29 */
! 75: CNST_LIMB (0xCB228), /* 30 */
! 76: CNST_LIMB (0x148ADD), /* 31 */
! 77: CNST_LIMB (0x213D05), /* 32 */
! 78: CNST_LIMB (0x35C7E2), /* 33 */
! 79: CNST_LIMB (0x5704E7), /* 34 */
! 80: CNST_LIMB (0x8CCCC9), /* 35 */
! 81: CNST_LIMB (0xE3D1B0), /* 36 */
! 82: CNST_LIMB (0x1709E79), /* 37 */
! 83: CNST_LIMB (0x2547029), /* 38 */
! 84: CNST_LIMB (0x3C50EA2), /* 39 */
! 85: CNST_LIMB (0x6197ECB), /* 40 */
! 86: CNST_LIMB (0x9DE8D6D), /* 41 */
! 87: CNST_LIMB (0xFF80C38), /* 42 */
! 88: CNST_LIMB (0x19D699A5), /* 43 */
! 89: CNST_LIMB (0x29CEA5DD), /* 44 */
! 90: CNST_LIMB (0x43A53F82), /* 45 */
! 91: CNST_LIMB (0x6D73E55F), /* 46 */
! 92: CNST_LIMB (0xB11924E1), /* 47 */
! 93: #endif
! 94: #if GMP_NUMB_BITS >= 64
! 95: CNST_LIMB (0x11E8D0A40), /* 48 */
! 96: CNST_LIMB (0x1CFA62F21), /* 49 */
! 97: CNST_LIMB (0x2EE333961), /* 50 */
! 98: CNST_LIMB (0x4BDD96882), /* 51 */
! 99: CNST_LIMB (0x7AC0CA1E3), /* 52 */
! 100: CNST_LIMB (0xC69E60A65), /* 53 */
! 101: CNST_LIMB (0x1415F2AC48), /* 54 */
! 102: CNST_LIMB (0x207FD8B6AD), /* 55 */
! 103: CNST_LIMB (0x3495CB62F5), /* 56 */
! 104: CNST_LIMB (0x5515A419A2), /* 57 */
! 105: CNST_LIMB (0x89AB6F7C97), /* 58 */
! 106: CNST_LIMB (0xDEC1139639), /* 59 */
! 107: CNST_LIMB (0x1686C8312D0), /* 60 */
! 108: CNST_LIMB (0x2472D96A909), /* 61 */
! 109: CNST_LIMB (0x3AF9A19BBD9), /* 62 */
! 110: CNST_LIMB (0x5F6C7B064E2), /* 63 */
! 111: CNST_LIMB (0x9A661CA20BB), /* 64 */
! 112: CNST_LIMB (0xF9D297A859D), /* 65 */
! 113: CNST_LIMB (0x19438B44A658), /* 66 */
! 114: CNST_LIMB (0x28E0B4BF2BF5), /* 67 */
! 115: CNST_LIMB (0x42244003D24D), /* 68 */
! 116: CNST_LIMB (0x6B04F4C2FE42), /* 69 */
! 117: CNST_LIMB (0xAD2934C6D08F), /* 70 */
! 118: CNST_LIMB (0x1182E2989CED1), /* 71 */
! 119: CNST_LIMB (0x1C5575E509F60), /* 72 */
! 120: CNST_LIMB (0x2DD8587DA6E31), /* 73 */
! 121: CNST_LIMB (0x4A2DCE62B0D91), /* 74 */
! 122: CNST_LIMB (0x780626E057BC2), /* 75 */
! 123: CNST_LIMB (0xC233F54308953), /* 76 */
! 124: CNST_LIMB (0x13A3A1C2360515), /* 77 */
! 125: CNST_LIMB (0x1FC6E116668E68), /* 78 */
! 126: CNST_LIMB (0x336A82D89C937D), /* 79 */
! 127: CNST_LIMB (0x533163EF0321E5), /* 80 */
! 128: CNST_LIMB (0x869BE6C79FB562), /* 81 */
! 129: CNST_LIMB (0xD9CD4AB6A2D747), /* 82 */
! 130: CNST_LIMB (0x16069317E428CA9), /* 83 */
! 131: CNST_LIMB (0x23A367C34E563F0), /* 84 */
! 132: CNST_LIMB (0x39A9FADB327F099), /* 85 */
! 133: CNST_LIMB (0x5D4D629E80D5489), /* 86 */
! 134: CNST_LIMB (0x96F75D79B354522), /* 87 */
! 135: CNST_LIMB (0xF444C01834299AB), /* 88 */
! 136: CNST_LIMB (0x18B3C1D91E77DECD), /* 89 */
! 137: CNST_LIMB (0x27F80DDAA1BA7878), /* 90 */
! 138: CNST_LIMB (0x40ABCFB3C0325745), /* 91 */
! 139: CNST_LIMB (0x68A3DD8E61ECCFBD), /* 92 */
! 140: CNST_LIMB (0xA94FAD42221F2702), /* 93 */
! 141: #endif
! 142: };
! 143:
! 144:
! 145: /* Store F[n] at fp and F[n-1] at f1p. fp and f1p should have room for
! 146: MPN_FIB2_SIZE(n) limbs.
! 147:
! 148: The return value is the actual number of limbs stored, this will be at
! 149: least 1. fp[size-1] will be non-zero, except when n==0, in which case
! 150: fp[0] is 0 and f1p[0] is 1. f1p[size-1] can be zero, since F[n-1]<F[n]
! 151: (for n>0).
! 152:
! 153: Notes:
! 154:
! 155: In F[2k+1] with k even, +2 is applied to 4*F[k]^2 just by ORing into the
! 156: low limb.
! 157:
! 158: In F[2k+1] with k odd, -2 is applied to the low limb of 4*F[k]^2 -
! 159: F[k-1]^2. This F[2k+1] is an F[4m+3] and such numbers are congruent to
! 160: 1, 2 or 5 mod 8, which means no underflow reaching it with a -2 (since
! 161: that would leave 6 or 7 mod 8).
! 162:
! 163: This property of F[4m+3] can be verified by induction on F[4m+3] =
! 164: 7*F[4m-1] - F[4m-5], that formula being a standard lucas sequence
! 165: identity U[i+j] = U[i]*V[j] - U[i-j]*Q^j.
! 166:
! 167: Enhancements:
! 168:
! 169: If there was an mpn_addlshift, it'd be possible to eliminate the yp
! 170: temporary, using xp=F[k]^2, fp=F[k-1]^2, f1p=xp+fp, fp+=4*fp, fp-=f1p,
! 171: fp+=2*(-1)^n, etc. */
! 172:
! 173: mp_size_t
! 174: mpn_fib2_ui (mp_ptr fp, mp_ptr f1p, unsigned long int n)
! 175: {
! 176: mp_ptr xp, yp;
! 177: mp_size_t size;
! 178: unsigned long nfirst, mask;
! 179: TMP_DECL (marker);
! 180:
! 181: TRACE (printf ("mpn_fib2_ui n=%lu\n", n));
! 182:
! 183: ASSERT (! MPN_OVERLAP_P (fp, MPN_FIB2_SIZE(n), f1p, MPN_FIB2_SIZE(n)));
! 184:
! 185: /* Take a starting pair from the table. */
! 186: mask = 1;
! 187: for (nfirst = n; nfirst > FIB_TABLE_LIMIT; nfirst /= 2)
! 188: mask <<= 1;
! 189: TRACE (printf ("nfirst=%lu mask=0x%lX\n", nfirst, mask));
! 190:
! 191: f1p[0] = FIB_TABLE ((int) nfirst - 1);
! 192: fp[0] = FIB_TABLE (nfirst);
! 193: size = 1;
! 194:
! 195: /* Skip to the end if the table lookup gives the final answer. */
! 196: if (mask != 1)
! 197: {
! 198: mp_size_t alloc;
! 199:
! 200: TMP_MARK (marker);
! 201: alloc = MPN_FIB2_SIZE (n);
! 202: TMP_ALLOC_LIMBS_2 (xp,alloc, yp,alloc);
! 203:
! 204: do
! 205: {
! 206: mp_limb_t c;
! 207:
! 208: /* Here fp==F[k] and f1p==F[k-1], with k being the bits of n from
! 209: n&mask upwards.
! 210:
! 211: The next bit of n is n&(mask>>1) and we'll double to the pair
! 212: fp==F[2k],f1p==F[2k-1] or fp==F[2k+1],f1p==F[2k], according as
! 213: that bit is 0 or 1 respectively. */
! 214:
! 215: TRACE (printf ("k=%lu mask=0x%lX size=%ld alloc=%ld\n",
! 216: n >> refmpn_count_trailing_zeros(mask),
! 217: mask, size, alloc);
! 218: mpn_trace ("fp ", fp, size);
! 219: mpn_trace ("f1p", f1p, size));
! 220:
! 221: /* fp normalized, f1p at most one high zero */
! 222: ASSERT (fp[size-1] != 0);
! 223: ASSERT (f1p[size-1] != 0 || f1p[size-2] != 0);
! 224:
! 225: /* f1p[size-1] might be zero, but this occurs rarely, so it's not
! 226: worth bothering checking for it */
! 227: ASSERT (alloc >= 2*size);
! 228: mpn_sqr_n (xp, fp, size);
! 229: mpn_sqr_n (yp, f1p, size);
! 230: size *= 2;
! 231:
! 232: /* Shrink if possible. Since fp was normalized there'll be at
! 233: most one high zero on xp (and if there is then there's one on
! 234: yp too). */
! 235: ASSERT (xp[size-1] != 0 || yp[size-1] == 0);
! 236: size -= (xp[size-1] == 0);
! 237: ASSERT (xp[size-1] != 0); /* only one xp high zero */
! 238:
! 239: /* Calculate F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k.
! 240: n&mask is the low bit of our implied k. */
! 241: c = mpn_lshift (fp, xp, size, 2);
! 242: fp[0] |= (n & mask ? 0 : 2); /* possible +2 */
! 243: c -= mpn_sub_n (fp, fp, yp, size);
! 244: ASSERT (n & (mask << 1) ? fp[0] != 0 && fp[0] != 1 : 1);
! 245: fp[0] -= (n & mask ? 2 : 0); /* possible -2 */
! 246: ASSERT (alloc >= size+1);
! 247: xp[size] = 0;
! 248: yp[size] = 0;
! 249: fp[size] = c;
! 250: size += (c != 0);
! 251:
! 252: /* Calculate F[2k-1] = F[k]^2 + F[k-1]^2.
! 253: F[2k-1]<F[2k+1] so no carry out of "size" limbs. */
! 254: ASSERT_NOCARRY (mpn_add_n (f1p, xp, yp, size));
! 255:
! 256: /* now n&mask is the new bit of n being considered */
! 257: mask >>= 1;
! 258:
! 259: /* Calculate F[2k] = F[2k+1] - F[2k-1], replacing the unwanted one of
! 260: F[2k+1] and F[2k-1]. */
! 261: ASSERT_NOCARRY (mpn_sub_n ((n & mask ? f1p : fp), fp, f1p, size));
! 262:
! 263: /* Can have a high zero after replacing F[2k+1] with F[2k].
! 264: f1p will have a high zero if fp does. */
! 265: ASSERT (fp[size-1] != 0 || f1p[size-1] == 0);
! 266: size -= (fp[size-1] == 0);
! 267: }
! 268: while (mask != 1);
! 269:
! 270: TMP_FREE (marker);
! 271: }
! 272:
! 273: TRACE (printf ("done size=%ld\n", size);
! 274: mpn_trace ("fp ", fp, size);
! 275: mpn_trace ("f1p", f1p, size));
! 276:
! 277: return size;
! 278: }
! 279:
! 280:
! 281:
! 282:
! 283:
! 284: /* ------------------------------------------------------------------------- */
! 285:
! 286: #if GENERATE_FIB_TABLE
! 287: /* Generate the tables of fibonacci data. This doesn't depend on the limb
! 288: size of the host, and doesn't need mpz_fib_ui working.
! 289:
! 290: The bit sizes in the table[] below will get specific setups so that a
! 291: build with GMP_NUMB_BITS equal to one of those values has as much data in
! 292: __gmp_fib_table as will fit that number of bits.
! 293:
! 294: A build with GMP_NUMB_BITS equal to some other value will effectively
! 295: fall back to the previous set of generated data. For instance if 8 and
! 296: 16 bits have been generated, but a build with 13 bits is done then
! 297: __gmp_fib_table will only contain 8 bit values, whereas it could probably
! 298: fit a few more. Everything still works, it's just that the table scheme
! 299: is not fully exploited. */
! 300:
! 301: int
! 302: main (void)
! 303: {
! 304: static struct {
! 305: int bits;
! 306: int fib_limit;
! 307: int luc_limit;
! 308: } table[] = {
! 309: { 4 },
! 310: { 8 },
! 311: { 16 },
! 312: { 32 },
! 313: { 64 },
! 314: };
! 315:
! 316: int i, t;
! 317: mpz_t f[500];
! 318: mpz_t l;
! 319:
! 320: mpz_init (l);
! 321: mpz_init_set_si (f[0], 1L); /* F[-1] */
! 322: mpz_init_set_si (f[1], 0L); /* F[0] */
! 323: for (i = 2; i < numberof(f); i++)
! 324: {
! 325: mpz_init (f[i]);
! 326: mpz_add (f[i], f[i-1], f[i-2]);
! 327: }
! 328:
! 329: for (i = 1; i < numberof (f); i++)
! 330: {
! 331: /* L[n] = F[n]+2*F[n-1] */
! 332: mpz_add (l, f[i], f[i-1]);
! 333: mpz_add (l, l, f[i-1]);
! 334:
! 335: for (t = 0; t < numberof (table); t++)
! 336: {
! 337: if (mpz_sizeinbase (f[i], 2) <= table[t].bits)
! 338: table[t].fib_limit = i-1;
! 339: if (mpz_sizeinbase (l, 2) <= table[t].bits)
! 340: table[t].luc_limit = i-1;
! 341: }
! 342: }
! 343: if (table[t].fib_limit == numberof (f) + 1)
! 344: {
! 345: printf ("Oops, need bigger f[] array\n");
! 346: abort ();
! 347: }
! 348:
! 349: for (t = numberof (table) - 1; t >= 0; t--)
! 350: {
! 351: printf ("#if GMP_NUMB_BITS >= %d\n", table[t].bits);
! 352: printf ("#define FIB_TABLE_LIMIT %d\n", table[t].fib_limit);
! 353: printf ("#define FIB_TABLE_LUCNUM_LIMIT %d\n", table[t].luc_limit);
! 354: if (t != 0)
! 355: printf ("#else\n");
! 356: }
! 357: for (t = 0; t < numberof (table); t++)
! 358: printf ("#endif /* %d */\n", table[t].bits);
! 359: printf ("\n");
! 360: printf ("\n");
! 361:
! 362: printf ("const mp_limb_t\n");
! 363: printf ("__gmp_fib_table[FIB_TABLE_LIMIT+2] = {\n");
! 364: printf ("\n");
! 365: t = 0;
! 366: i = 0;
! 367: printf ("#if GMP_NUMB_BITS >= %d\n", table[t].bits);
! 368: for (;;)
! 369: {
! 370: gmp_printf (" CNST_LIMB (0x%ZX), /* %d */\n", f[i], i-1);
! 371:
! 372: if (i-1 == table[t].fib_limit)
! 373: {
! 374: printf ("#endif\n");
! 375: do
! 376: {
! 377: t++;
! 378: if (t >= numberof (table))
! 379: goto done;
! 380: }
! 381: while (i-1 == table[t].fib_limit);
! 382: printf ("#if GMP_NUMB_BITS >= %d\n", table[t].bits);
! 383: }
! 384: i++;
! 385: }
! 386: done:
! 387: printf ("};\n");
! 388:
! 389: return 0;
! 390: }
! 391:
! 392: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to fib2_ui.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpn / generic