Annotation of OpenXM_contrib/gmp/PROJECTS, Revision 1.1
1.1 ! maekawa 1: IDEAS ABOUT THINGS TO WORK ON
! 2:
! 3: * mpq_cmp: Maybe the most sensible thing to do would be to multiply the, say,
! 4: 4 most significant limbs of each operand and compare them. If that is not
! 5: sufficient, do the same for 8 limbs, etc.
! 6:
! 7: * Write mpi, the Multiple Precision Interval Arithmetic layer.
! 8:
! 9: * Write `mpX_eval' that take lambda-like expressions and a list of operands.
! 10:
! 11: * As a general rule, recognize special operand values in mpz and mpf, and
! 12: use shortcuts for speed. Examples: Recognize (small or all) 2^n in
! 13: multiplication and division. Recognize small bases in mpz_pow_ui.
! 14:
! 15: * Implement lazy allocation? mpz->d == 0 would mean no allocation made yet.
! 16:
! 17: * Maybe store one-limb numbers according to Per Bothner's idea:
! 18: struct {
! 19: mp_ptr d;
! 20: union {
! 21: mp_limb val; /* if (d == NULL). */
! 22: mp_size size; /* Length of data array, if (d != NULL). */
! 23: } u;
! 24: };
! 25: Problem: We can't normalize to that format unless we free the space
! 26: pointed to by d, and therefore small values will not be stored in a
! 27: canonical way.
! 28:
! 29: * Document complexity of all functions.
! 30:
! 31: * Add predicate functions mpz_fits_signedlong_p, mpz_fits_unsignedlong_p,
! 32: mpz_fits_signedint_p, etc.
! 33:
! 34: mpz_floor (mpz, mpq), mpz_trunc (mpz, mpq), mpz_round (mpz, mpq).
! 35:
! 36: * Better random number generators. There should be fast (like mpz_random),
! 37: very good (mpz_veryrandom), and special purpose (like mpz_random2). Sizes
! 38: in *bits*, not in limbs.
! 39:
! 40: * It'd be possible to have an interface "s = add(a,b)" with automatic GC.
! 41: If the mpz_xinit routine remembers the address of the variable we could
! 42: walk-and-mark the list of remembered variables, and free the space
! 43: occupied by the remembered variables that didn't get marked. Fairly
! 44: standard.
! 45:
! 46: * Improve speed for non-gcc compilers by defining umul_ppmm, udiv_qrnnd,
! 47: etc, to call __umul_ppmm, __udiv_qrnnd. A typical definition for
! 48: umul_ppmm would be
! 49: #define umul_ppmm(ph,pl,m0,m1) \
! 50: {unsigned long __ph; (pl) = __umul_ppmm (&__ph, (m0), (m1)); (ph) = __ph;}
! 51: In order to maintain just one version of longlong.h (gmp and gcc), this
! 52: has to be done outside of longlong.h.
! 53:
! 54: Bennet Yee at CMU proposes:
! 55: * mpz_{put,get}_raw for memory oriented I/O like other *_raw functions.
! 56: * A function mpfatal that is called for exceptions. Let the user override
! 57: a default definition.
! 58:
! 59: * Make all computation mpz_* functions return a signed int indicating if the
! 60: result was zero, positive, or negative?
! 61:
! 62: * Implement mpz_cmpabs, mpz_xor, mpz_to_double, mpz_to_si, mpz_lcm, mpz_dpb,
! 63: mpz_ldb, various bit string operations. Also mpz_@_si for most @??
! 64:
! 65: * Add macros for looping efficiently over a number's limbs:
! 66: MPZ_LOOP_OVER_LIMBS_INCREASING(num,limb)
! 67: { user code manipulating limb}
! 68: MPZ_LOOP_OVER_LIMBS_DECREASING(num,limb)
! 69: { user code manipulating limb}
! 70:
! 71: Brian Beuning proposes:
! 72: 1. An array of small primes
! 73: 3. A function to factor a mpz_t. [How do we return the factors? Maybe
! 74: we just return one arbitrary factor? In the latter case, we have to
! 75: use a data structure that records the state of the factoring routine.]
! 76: 4. A routine to look for "small" divisors of an mpz_t
! 77: 5. A 'multiply mod n' routine based on Montgomery's algorithm.
! 78:
! 79: Dough Lea proposes:
! 80: 1. A way to find out if an integer fits into a signed int, and if so, a
! 81: way to convert it out.
! 82: 2. Similarly for double precision float conversion.
! 83: 3. A function to convert the ratio of two integers to a double. This
! 84: can be useful for mixed mode operations with integers, rationals, and
! 85: doubles.
! 86:
! 87: Elliptic curve method description in the Chapter `Algorithms in Number
! 88: Theory' in the Handbook of Theoretical Computer Science, Elsevier,
! 89: Amsterdam, 1990. Also in Carl Pomerance's lecture notes on Cryptology and
! 90: Computational Number Theory, 1990.
! 91:
! 92: * Harald Kirsh suggests:
! 93: mpq_set_str (MP_RAT *r, char *numerator, char *denominator).
! 94:
! 95: * New function: mpq_get_ifstr (int_str, frac_str, base,
! 96: precision_in_som_way, rational_number). Convert RATIONAL_NUMBER to a
! 97: string in BASE and put the integer part in INT_STR and the fraction part
! 98: in FRAC_STR. (This function would do a division of the numerator and the
! 99: denominator.)
! 100:
! 101: * Should mpz_powm* handle negative exponents?
! 102:
! 103: * udiv_qrnnd: If the denominator is normalized, the n0 argument has very
! 104: little effect on the quotient. Maybe we can assume it is 0, and
! 105: compensate at a later stage?
! 106:
! 107: * Better sqrt: First calculate the reciprocal square root, then multiply by
! 108: the operand to get the square root. The reciprocal square root can be
! 109: obtained through Newton-Raphson without division. To compute sqrt(A), the
! 110: iteration is,
! 111:
! 112: 2
! 113: x = x (3 - A x )/2.
! 114: i+1 i i
! 115:
! 116: The final result can be computed without division using,
! 117:
! 118: sqrt(A) = A x .
! 119: n
! 120:
! 121: * Newton-Raphson using multiplication: We get twice as many correct digits
! 122: in each iteration. So if we square x(k) as part of the iteration, the
! 123: result will have the leading digits in common with the entire result from
! 124: iteration k-1. A _mpn_mul_lowpart could help us take advantage of this.
! 125:
! 126: * Peter Montgomery: If 0 <= a, b < p < 2^31 and I want a modular product
! 127: a*b modulo p and the long long type is unavailable, then I can write
! 128:
! 129: typedef signed long slong;
! 130: typedef unsigned long ulong;
! 131: slong a, b, p, quot, rem;
! 132:
! 133: quot = (slong) (0.5 + (double)a * (double)b / (double)p);
! 134: rem = (slong)((ulong)a * (ulong)b - (ulong)p * (ulong)quot);
! 135: if (rem < 0} {rem += p; quot--;}
! 136:
! 137: * Speed modulo arithmetic, using Montgomery's method or my pre-inversion
! 138: method. In either case, special arithmetic calls would be needed,
! 139: mpz_mmmul, mpz_mmadd, mpz_mmsub, plus some kind of initialization
! 140: functions. Better yet: Write a new mpr layer.
! 141:
! 142: * mpz_powm* should not use division to reduce the result in the loop, but
! 143: instead pre-compute the reciprocal of the MOD argument and do reduced_val
! 144: = val-val*reciprocal(MOD)*MOD, or use Montgomery's method.
! 145:
! 146: * mpz_mod_2expplussi -- to reduce a bignum modulo (2**n)+s
! 147:
! 148: * It would be a quite important feature never to allocate more memory than
! 149: really necessary for a result. Sometimes we can achieve this cheaply, by
! 150: deferring reallocation until the result size is known.
! 151:
! 152: * New macro in longlong.h: shift_rhl that extracts a word by shifting two
! 153: words as a unit. (Supported by i386, i860, HP-PA, POWER, 29k.) Useful
! 154: for shifting multiple precision numbers.
! 155:
! 156: * The installation procedure should make a test run of multiplication to
! 157: decide the threshold values for algorithm switching between the available
! 158: methods.
! 159:
! 160: * Fast output conversion of x to base B:
! 161: 1. Find n, such that (B^n > x).
! 162: 2. Set y to (x*2^m)/(B^n), where m large enough to make 2^n ~~ B^n
! 163: 3. Multiply the low half of y by B^(n/2), and recursively convert the
! 164: result. Truncate the low half of y and convert that recursively.
! 165: Complexity: O(M(n)log(n))+O(D(n))!
! 166:
! 167: * Improve division using Newton-Raphson. Check out "Newton Iteration and
! 168: Integer Division" by Stephen Tate in "Synthesis of Parallel Algorithms",
! 169: Morgan Kaufmann, 1993 ("beware of some errors"...)
! 170:
! 171: * Improve implementation of Karatsuba's algorithm. For most operand sizes,
! 172: we can reduce the number of operations by splitting differently.
! 173:
! 174: * Faster multiplication: The best approach is to first implement Toom-Cook.
! 175: People report that it beats Karatsuba's algorithm already at about 100
! 176: limbs. FFT would probably never beat a well-written Toom-Cook (not even for
! 177: millions of bits).
! 178:
! 179: FFT:
! 180: {
! 181: * Multiplication could be done with Montgomery's method combined with
! 182: the "three primes" method described in Lipson. Maybe this would be
! 183: faster than to Nussbaumer's method with 3 (simple) moduli?
! 184:
! 185: * Maybe the modular tricks below are not needed: We are using very
! 186: special numbers, Fermat numbers with a small base and a large exponent,
! 187: and maybe it's possible to just subtract and add?
! 188:
! 189: * Modify Nussbaumer's convolution algorithm, to use 3 words for each
! 190: coefficient, calculating in 3 relatively prime moduli (e.g.
! 191: 0xffffffff, 0x100000000, and 0x7fff on a 32-bit computer). Both all
! 192: operations and CRR would be very fast with such numbers.
! 193:
! 194: * Optimize the Schoenhage-Stassen multiplication algorithm. Take advantage
! 195: of the real valued input to save half of the operations and half of the
! 196: memory. Use recursive FFT with large base cases, since recursive FFT has
! 197: better memory locality. A normal FFT get 100% cache misses for large
! 198: enough operands.
! 199:
! 200: * In the 3-prime convolution method, it might sometimes be a win to use 2,
! 201: 3, or 5 primes. Imagine that using 3 primes would require a transform
! 202: length of 2^n. But 2 primes might still sometimes give us correct
! 203: results with that same transform length, or 5 primes might allow us to
! 204: decrease the transform size to 2^(n-1).
! 205:
! 206: To optimize floating-point based complex FFT we have to think of:
! 207:
! 208: 1. The normal implementation accesses all input exactly once for each of
! 209: the log(n) passes. This means that we will get 0% cache hit when n >
! 210: our cache. Remedy: Reorganize computation to compute partial passes,
! 211: maybe similar to a standard recursive FFT implementation. Use a large
! 212: `base case' to make any extra overhead of this organization negligible.
! 213:
! 214: 2. Use base-4, base-8 and base-16 FFT instead of just radix-2. This can
! 215: reduce the number of operations by 2x.
! 216:
! 217: 3. Inputs are real-valued. According to Knuth's "Seminumerical
! 218: Algorithms", exercise 4.6.4-14, we can save half the memory and half
! 219: the operations if we take advantage of that.
! 220:
! 221: 4. Maybe make it possible to write the innermost loop in assembly, since
! 222: that could win us another 2x speedup. (If we write our FFT to avoid
! 223: cache-miss (see #1 above) it might be logical to write the `base case'
! 224: in assembly.)
! 225:
! 226: 5. Avoid multiplication by 1, i, -1, -i. Similarly, optimize
! 227: multiplication by (+-\/2 +- i\/2).
! 228:
! 229: 6. Put as many bits as possible in each double (but don't waste time if
! 230: that doesn't make the transform size become smaller).
! 231:
! 232: 7. For n > some large number, we will get accuracy problems because of the
! 233: limited precision of our floating point arithmetic. This can easily be
! 234: solved by using the Karatsuba trick a few times until our operands
! 235: become small enough.
! 236:
! 237: 8. Precompute the roots-of-unity and store them in a vector.
! 238: }
! 239:
! 240: * When a division result is going to be just one limb, (i.e. nsize-dsize is
! 241: small) normalization could be done in the division loop.
! 242:
! 243: * Never allocate temporary space for a source param that overlaps with a
! 244: destination param needing reallocation. Instead malloc a new block for
! 245: the destination (and free the source before returning to the caller).
! 246:
! 247: * Parallel addition. Since each processors have to tell it is ready to the
! 248: next processor, we can use simplified synchronization, and actually write
! 249: it in C: For each processor (apart from the least significant):
! 250:
! 251: while (*svar != my_number)
! 252: ;
! 253: *svar = my_number + 1;
! 254:
! 255: The least significant processor does this:
! 256:
! 257: *svar = my_number + 1; /* i.e., *svar = 1 */
! 258:
! 259: Before starting the addition, one processor has to store 0 in *svar.
! 260:
! 261: Other things to think about for parallel addition: To avoid false
! 262: (cache-line) sharing, allocate blocks on cache-line boundaries.
! 263:
! 264:
! 265: Local Variables:
! 266: mode: text
! 267: fill-column: 77
! 268: fill-prefix: " "
! 269: version-control: never
! 270: End:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to PROJECTS CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp