]]> 2001 American Mathematical Society Mathematics of Computation Math. Comp. mcom 71 238 20010803 July 23, 1999 August 8, 2000 August 3, 2001 863 863-871 9 S0025-5718-01-01341-2 10.1090/S0025-5718-01-01341-2 0025-5718 1088-6842 English 11Y40 11Y60 11Y16 11B68 11R42 11R29 94A60 11R18 2000 http://www.ams.org/jourcgi/jour-getitem?pii=S0025-5718-01-01341-2 Definitions Let k 0 be a totally real number field, and let p be an odd prime. Let k 1 k 0( p) , where p n will denote a primitive p n -th root of unity. Let k 1 k 0 , and let . Let p e be the largest power of p such that p e k 0( p) . definition Let k 0 be the zeta function for k 0 . We say that p is k 0 -regular if p is relatively prime to k 0 (1-2m) for all integers m such that 2 2m - 2 and also p is relatively prime to p e k 0 (1 - ) . The number of such zeta-values that are divisible by p will be the index of k 0 -irregularity of p . definition According to a well-known theorem of Kummer, p divides the order of the class group of ( p) if and only if p divides the numerator of a Bernoulli number B m for some even m such that 2 m p-3 . Such primes are called irregular; the others are called regular. In the setting we have described above, the author proved in his thesis ( dissert , see also Holden99 ), building on work of Greenberg and Kudo, that under a certain technical condition Kummer's criterion can be extended to give information about whether p divides the class group of k 0( p) . To be exact, let k 1 denote the maximal real subfield of k 1 , which is equal to k 0( p p -1 ) . Let h(k 1) denote the class number of k 1 and h (k 1) denote the class number of k 1 . It is known that h (k 1) h(k 1) ; we let the relative class number h - (k 1) be the quotient. theorem Greenberg, Holden GH Assume that no prime of the field k 1 lying over p splits in k 1 . Then p divides h - (k 1) if and only if p is not k 0 -regular. theorem As an application, we note that one common way of constructing public-key cryptographic systems is to utilize the problem of finding a discrete logarithm in some abelian group. In order to make sure that the discrete logarithm problem is computationally hard, one needs to know something about the structure of the group involved, e.g., that it is divisible by a large prime. Theorem GH shows that if p is a large k 0 -irregular prime and the conditions of the theorem are met, then the class group of k 0( p) may be especially suitable for cryptography. (One should see BP97 for more on the use of class groups in cryptography.) For the case we consider, k 0 will be a real quadratic field () , with D a positive fundamental discriminant. For such a k 0 , we will say that primes are D -regular or have given index of D -irregularity, and we will let the zeta function k 0 be also denoted by D . (More generally, we may refer to the concept as quadratic irregularity''.) In this case will be equal to p-1 unless D p , in which case (p-1) 2 . Also, e is always equal to 1 when p does not divide the order of k 0 over , which is true in this case since p is odd. For the condition in Theorem GH that no prime of the field k 1 lying over p splits in k 1 to be satisfied, it is sufficient that p should not divide D , and we should also note that since p does not divide the degree of k 0 () over , a theorem of Leopoldt shows that p divides h(k 1) if and only if p divides h - (k 1) . In general, we will consider three cost models for the time of multiplication: first using naive multiplication ( O(tt') ), second using Schonhage-Strassen fast multiplication or a similar method ( O(t (t') O(1) )), and third using a model where multiplication (or addition) takes constant time regardless of the size of the factors ( O(1) ). We do not expect constant time multiplication to occur asymptotically in the real world, but it can provide useful estimates in situations where the size of the numbers involved is small compared to the word size of the actual computer in question. (In these running time bounds, t is the number of bits in the larger multiplicand and t' the number of bits in the smaller.) First formula Siegel's first formula to compute D(1-2m) for m 1 an integer is analogous to the formula (1-2m) -B 2m (2m). Using elementary methods, Siegel showed that similarly equation eq1 D(1-2m) B 2m 4m 2 D 2m-1 j 1 D (j) B 2m (j D). equation Here (j) , the Kronecker symbol, and B 2m (j D) indicates the 2m -th Bernoulli polynomial evaluated at the fraction j D . The Bernoulli polynomial B r(x) can be computed from the Bernoulli numbers as equation B r(x) s 0 rB r-s x s. equation It is not difficult to estimate the sizes of the numbers involved. We will assume throughout that B m , 1 m M , are precomputed over the common denominator of the final result, and stored in this fashion each has size O(m (m D)) bits for a total table size of O(M 2(M D)) bits. (The precomputation does not in fact add to the asymptotic running time.) The rational numbers B 2m (j D) can also be stored in O(m(m D)) bits, as can the total. See Holden98 for more details. A first attempt at an algorithm based on ( eq1 ) might compute B 0(), , B M() naively from the formula. The time taken for this would be dominated by the powerings. For a b some rational number, the total time with naive multiplication would be O(M 4 (M a b) O(1) ). Using fast multiplication instead of naive multiplication would improve this to O(M 3 (M a b) O(1) ), while with constant time multiplication we need only time O(M 2M) regardless of a and b . However we can do better than this, using a cross between Horner's method of evaluating polynomials and an algorithm used by Brent to calculate Bernoulli numbers, as was previously discussed by the author in Holden98 . This method gives a total time of O(M 3 (M a b) O(1) ) using either constant or fast multiplication, and time O(M 2) using constant time multiplication. Using either of these algorithms to compute D(1-2m) , 2 2m M , isthen relatively straightforward. Note that the Kronecker symbol (j) can becomputed in time O( 2 D) . The slower version of the algorithm has time O(M 4D (M D) O(1) ) with naive multiplication, O(M 3D (M D) O(1) ) with fast multiplication, and O(M 2D (M D) O(1) ) with constant time multiplication. The faster version runs in time O(M 3D (D M) O(1) ) with either naive or fast multiplication (the O(1) factor is different, of course) and again in time O(M 2D (M D) O(1) ) with constant time multiplication. Second formula Siegel's second formula is, as I said, derived from the theory of modular forms. In general, for k 0 a totally real number field as above, it says that k 0 (1-2m) - 2 n c 2mn -1 l 1 r c 2mn,l s l k 0 (2m), where n k 0: , c 2mn c 2mn,0 and c 2mn,l are rational integers depending only on 2mn and l (given by explicit formulas which we will discuss), r cases mn 6 if 2mn 2 modulo 12 , mn 6 1 otherwise , cases and s l k 0 is a sum over norms of ideals in the ring of integers of k 0 , namely, s l k 0 (2m) () -1 , 0, () l 2m-1 (()), where 2m-1 () N() 2m-1 is a generalization of the usual sum of powers function and is the different of k 0 . In the quadratic case this all becomes much easier: equation eq2 D(1-2m) - 4 c 4m -1 l 1 r c 4m,l s l D(2m), r m 3 1, equation s l D(2m) () -1 , 0, () l 2m-1 (()), and s l D(2m) can also be expressed in terms of a purely arithmetic function e 2m-1 (n) , as follows: s l k 0 (2m) jl D(j) j 2m-1 e 2m-1 ((l j) 2D) and e 2m-1 (n) x 2n 4 x 2m-1 ( n-x 2 4) where 2m-1 (n) dn d 2m-1 is the usual sum-of-powers function. (See Siegel68 , Zagier , Cohen75 and Cohen76 for more detailed descriptions of these formulas, and for their derivations.) The coefficients c 4m,l are most easily expressed as the coefficients of a certain power series, and can be computed as needed without adding to the asymptotic running time. We will give explicit formulas for the power series and discuss its computation in Section coefficients . It is not hard to prove that in this form c 4m,l is of size O(m) . The running time for calculating the function e 2m-1 (n) is complicated by the need for factoring; we use here an estimate based on the elliptic curve factoring method (we would expect something very similar with any of the other standard subexponential methods) to get an expected running time involving the function (x) e x x . Given this, we get an expected running time to compute e 2m-1 (n) of O((n) 1 o(1) (2m-1) 2 2 n) using naive multiplication. If we now applied ( eq2 ) as it was originally stated to compute all D(1-2m) , 2 2m M , we would get a running time of O(M 3 (M) O(1) (D) O(1) M M 5 M (M D) O(1) ), again using naive multiplication. However, it is more efficient to rearrange the terms of the formula as follows: eqnarray D(1-2m) -4 c 4m -1 l 1 r c 4m,l s l k 0 (2m) -4 c 4m -1 l 1 r c 4m,l jl D(j) j 2m-1 e 2m-1 ((l j) 2D) eq3 -4 c 4m -1 k 1 r ( j 1 r k D(j) j 2m-1 c 4m,jk ) e 2m-1 (k 2D). eqnarray This rearrangement of the formula requires fewer calls to compute e 2m-1 by a factor of m . Using this version of the formula, the time necessary to compute all D(1-2m) , 2 2m M , using naive multiplication is O(M 3 (M) O(1) (D) O(1) M 5 (M D) O(1) ). This is much worse than the best algorithm based on ( eq1 ) in terms of M , but it is better in terms of D . Also, except for one final division by c 4m -1 , all of the arithmetic in this formula deals only with rational integers, unlike the previous formulas. Note that the first term comes from the factoring process, while the second term comes from multiplications. It should be noted that the asymptotic running time of this algorithm is greatly improved by using Schonhage-Strassen fast multiplication or constant time multiplication, in which cases the second term becomes smaller than the first and the running time becomes O(M 3 (M) O(1) (D) O(1) ). This is still worse than using ( eq1 ) in terms of M , but only by a subexponential factor. It should also be noted that ( eq2 ) and ( eq3 ) also present opportunities for time savings when computing zeta-values for multiple D in the same range of M , at a sacrifice of memory space. The controlling factor in the speed of the algorithm is the number of times that 2m-1 (n) must be calculated. Note that in computing all d(1-2m) , 5 d D , there can only be O(m 2D) different values of n . However, following the algorithm strictly, we would ordinarily make O(m 2D 3 2 ) calls to the subroutine that calculates this function. Thus if we compute all d(1-2m) , 5 d D , storing values of 2m-1 (n) as we compute them, and then repeat this process for each m in the range 2 2m M , the running time should be O(D (D) O(1) ) in terms of D , rather than O(D 3 2 (D) O(1) ) as one would obtain following the algorithm strictly. This compares very favorably with the time of O(D 2) in terms of D which holds for algorithms using ( eq1 ). Since the exponent 2m-1 used in the 2m-1 (n) function changes as m does, we can dispose of the table when we change m . The table that we need to keep requires at most O(M 3D(M D)) bits of storage, which could be a significant barrier. More efficient storage of the important information may be valuable here; we will discuss this somewhat more in Section conclusion . Computing the numbers c 4m,l coefficients The integers c 4m,l are defined as follows. Let G k 1- 2k B k n 1 k-1 (n) , q n be the (normalized) Eisenstein series of order k for k 6, 10 , and 14 . (For the general c 2mn,l one also needs k 0, 4 , and 8 .) Let aligned t q n 1 (1-q n) 24 q ( n 0 (-1) n (2n 1) , q n(n 1) 2 ) 8 aligned be the discriminant series. Let r m 3 1 as before, and let T 4m G 12r-4m 2 -r n -r c 4m, -n , q n. Then c 4m c 4m, 0 , and the other c 4m,l for 1 l r can also be read off as coefficients of T 4m . Luckily, the expression 12r-4m 2 only takes on the values 6 , 10 , and 14 . (In the more general case we can define T 2mn similarly; the expression 12r-2mn 2 can take on the values 0 , 4 , and 8 in addition to those above.) The best algorithm known to the author for calculating these coefficients goes roughly as follows. At the start of the computations for D(1-2m) , 2 2m M , calculate G 6 , G 10 , G 14 , and -1 with the maximum number of coefficients necessary (about M 12 ). Instead of trying to compute all of the needed series -r at once, we calculate it as a running product which only needs to be updated when r changes. Then, whenever m changes, we multiply truncated versions G 12r-4m 2 and -r (with about m 6 coefficients each) to find the required coefficients of T 4m . The series -1 can also be expressed as aligned t -1 q n 1 (1-q n) -24 q ( n 1 1 (1-q n) ) 24 q ( n 0 p(n) ,q n) 24 aligned where p(n) takes on integer values and is well-known as the partition function from additive number theory. Hardy and Ramanujan proved an asymptotic expression for p(n) which shows that p(n) is of order . (See, for example, Chapter 14 of Apostol .) From this it is easy to show that the coefficients of -r take at most O(M 1.5 ) bits of storage each, as do the coefficients of G 12r-4m 2 and T 4m . Thus the storage required for all of the computation of c 4m, l necessary for a fixed m is O(M 2.5 ) . The resulting table, which is of course independent of d , can be used to compute d(1-2m) for any range of d and then disposed of when m is changed. As far as the time for this algorithm is concerned, computing the tables necessary for all d(1-2m) , 2 2m M , 5 d D , takes time O(M 5) with naive multiplication of integers and also of polynomials. By using FFT methods to multiply polynomials the time can be reduced to O(M 3.5 M M) with fast multiplication of integers and O(M 2 M) with constant time multiplication of integers. (The use of FFT methods here was suggested by A.O.L. Atkin and Will Galway.) This is within our previously established time bounds in the naive and constant time cases; however in the fast multiplication case it could add to the total asymptotic time in terms of M , which for the previous parts of the algorithm was established as O(M 3 (M) O(1) ) in terms of M . On the other hand it should be noted that these calculations only need to be done once per value of m no matter how many values of d one is examining. Also the constant involved in the O(M 3.5 M M) seems to be quite good in practice compared to that in the O(M 3 (M) O(1) ) . Summary and future work conclusion Table algstable presents the various algorithms, for comparison. We present the asymptotic order of the running time to compute D(1-2m) , 2 2m M , for a given D , using the three methods of multiplication discussed earlier. table Comparison of algorithms for D(1-2m) , 2 2m M algstable tabular lll Equation used Multiplication Time order ( eq1 ) naive M 4D (M D) O(1) ( eq1 ) fast M 3D (M D) O(1) ( eq1 ) constant M 2 D (M D) O(1) ( eq1 ) from Holden98 naive M 3D (D M) O(1) ( eq1 ) from Holden98 fast M 3D (D M) O(1) ( eq1 ) from Holden98 constant M 2 D (M D) O(1) ( eq2 ) naive M 3 (M) O(1) (D) O(1) M M 5 (M D) O(1) ( eq2 ) fast M 3 (M) O(1) (D) O(1) M ( eq2 ) constant M 3 (M) O(1) (D) O(1) M ( eq3 ) naive M 3 (M) O(1) (D) O(1) M 5 (M D) O(1) ( eq3 ) fast M 3 (M) O(1) (D) O(1) ( eq3 ) constant M 3 (M) O(1) (D) O(1) tabular table The factor of M in the times for algorithms based on ( eq2 ) has been included to emphasize that these algorithms are slower than those based on ( eq3 ), even though the factor of M could be absorbed into that of (M) O(1) . We also provide, in Tables timetable and timetable2 , tables of actual timings for some of the algorithms, using naive multiplication. The times were measured on a Sun SPARC Ultra-1 computer using the GP-Pari interpreted language. (See PARImanual .) The data computed by these programs will be analyzed in a future paper. Table timetable measures the time to compute D(1-2m) , 2 2m M , for a given D , in hours, minutes, and seconds. The number in parentheses indicates the size of the stack used, in bytes. These numbers are only a very rough guide to the actual amount of memory used. table Calculating D(1-2m) , 2 2m M timetable tabular r r rr rr D M time ( eq1 ) from Holden98 (stacksize) time ( eq3 ) (stacksize) 5 100 3.155 (10M) .838 (10M) 101 100 55.561 (10M) 3.758 (10M) 501 100 4:50.282 (10M) 12.480 (10M) 1001 100 10:52.411 (10M) 20.670 (10M) 5001 100 48:27.107 (12M) 1:43.548 (10M) 5 500 2:05.670 (4M) 8:11.603 (4M) 101 500 42:38.612 (4M) 1:18:26.615 (4M) 501 500 3:48:18.615 (8M) 4:45:20.438 (4M) 1001 500 7:49:56.048 (12M) 7:49:00.783 (4M) 5 1000 10:05.903 (4M) 3:55:46.908 (4M) 5 2000 1:14:01.992 (16M) 118:17:46.020 (4M) tabular table Table timetable2 measures the time to compute d(1-2m) , 2 2m M , 5 d D . We use the algorithm based on ( eq3 ), both with and without keeping a table of 2m-1 (n) as described earlier. The units and stack size numbers should be interpreted as in Table timetable . table Calculating d(1-2m) , 2 2m M , 5 d D timetable2 tabular r r rr rr D M time ( eq3 ) (stacksize) time ( eq3 ) with table (stacksize) 100 100 1:19.638 (4M) 1:27.983 (4M) 500 100 18:34.377 (4M) 13:57.166 (8M) 1000 100 1:03:50.464 (4M) 37:52.552 (12M) 5000 100 26:39:46.109 (4M) 7:10:44.870 (64M) tabular table As one can see, memory usage for algorithms based on ( eq3 ) with a table of values 2m-1 (n) goes up quite quickly, and even so a lot of redundant work is being done. For one thing, the same numbers n will have to be factored repeatedly for different values of 2m-1 , but they will lead to different values of 2m-1 (n) so the actual factorization would need to be stored and not just a function value. This would be even more memory intensive. Carl Pomerance has suggested some approaches to this, including using a cache rather than a complete table, and storing only the largest prime factor rather than a complete factorization. The question of a cache leads naturally to the question of which numbers will appear as integer values of (k D-x 2) 4 as k and D vary, and how often. Henri Cohen, in Cohen76 , gives some variations on Siegel's formula which could cut down on the number of times 2m-1 (n) needs to be computed. Algorithms based on these might provide a linear speedup over the algorithms presented here, but the asymptotic behavior would probably be the same. The anonymous reviewer has suggested that some of the arithmetic could possibly be speeded up by the use of modular techniques and the Chinese remainder theorem. This certainly deserves more consideration. Another prospect for future work is the analysis of first-hit'' versions of these algorithms, namely determining how long we should expect to search before finding, say, the first D -irregular prime larger than a certain bound for D in a given range. This might be particularly useful for cryptographic applications, in which we would be explicitly looking for a class group (or a small number of them) with a hard discrete logarithm problem. This will be explained in more detail in mathcomp2 . Acknowledgments The author would like to thank A.O.L. Atkin, Will Galway, and the members of the NMBRTHY electronic mailing list for their very helpful suggestions, and Robert Harley for generating some of the precomputed tables used in the computations reported in this paper. He would also like to especially thank Johannes Buchmann, Henri Cohen, and Carl Pomerance for suggestions and encouragement, and Gary Walsh for his encouraging remarks on an early version of this paper. Finally, he would like to thank the anonymous reviewer for encouragement and several helpful suggestions.