LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 510.84 I-£6r no. 171-187 cop. a. Report No. 175 n>\0- 17 3 April 22, 196; THE NUMBER OF LATTICE POINTS IN A K -DIMENSIONAL HYPERSPHERE by W. C. Mitchell JW ■■•• -V 11... Digitized by the Internet Archive in 2013 http://archive.org/details/numberoflatticep175mitc Report No. 175 THE NUMBER OF LATTICE POINTS IN A K-DIMENSIONAL HYPERSPHERE by W. C. Mitchell April 22, 1965 Department of Computer Science University of Illinois Urbana, Illinois ACKNOWLEDGMENTS The author would like to express his appreciation to Stephen Czuchlewski of Manhattan College, J. R. Ehrman of the University of Illinois, and T. M. Apostol of the California Institute of Technology for assistance in this work. This work was supported by National Science Foundation Grant Gl6489 during the summer of I96U and was continued with some financial support for auxiliary computing time from the Computing Center of Miami University and from the Department of Mathematics of the California Institute of Technology. 1. Introduction One of the most interesting problems of analytic number theory- involves the difference between the number of lattice points in a K-dimensional hyper sphere and the "volume" of the hyper sphere. Define the set L ( x ) as follows: L k (x) = j {J V J 2 , . . . ,J k ) I E j\ S x f (i) v. i=l J where the J. are integers. Let A, (x) be the number of distinct points in L, (x). Thus \.(x) is the number of lattice points in a K-dimensional hyper- 1/2 sphere of radius x ' . Define V, (x) to be the "volume" of a K-dimensional K. hypersphere of radius x ' . k/2 k/2 [k/2] k/2 r( | +1) I- where [z] is the integer part of z and k! ! = k(k - 2) (k - *+)... (l or 2). The problem of primary interest is to find the Greatest Lower Bound 0. of the set of values 6 for which k P k (x) ^ A^x) - V k (x) = 0(x S ) (3) Walfisz [1] gives the following general results: P k (x) = 0(x (k - l} / 2 ) p k ( x ) = fi ( x k /2-l) {k) P i+ (x) = 0(x log 2 x) k > 5 P k (x) = 0(x k / 2 " 1 ) -1- Thus for k > h k 2 The value of k which has received the greatest attention is k = 2, the number of lattice points in a circle. Wilton [2] gives an account of the early work in this problem. Since that time there have been several results given establishing new values of 6 for which P p (x) = 0(x ). One of the most recent is Hua's proof [3] that P p (x) = 0(x ' ). It has also been shown by Hardy (see [2]) that P p (x) = fi(x ' ). It is a common conjucture that P 2 (x) = 0(x 1 / i++€ ), G >0. With the advent of high-speed computers it has become possible to evaluate P^(x) for "large" x in order to see if the calculated results are consistent with theoretical results or if it is reasonable to make any new conjectures concerning 6 . There have been three previous papers published on this subject. Fraser and Gotlieb [h] calculated isolated values of P p (x) and 1/2 P (x) for x ' < 2000 on an IBM 65O. However, their conclusions differ with the present paper for P (x). Harry Mitchell [5] calculated P p (x) for 1/2 1/2 x ' < 200,000 on an IBM 7090 but the results for x' > 3000 are incorrect, as pointed out by Keller and Swenson [6]. Keller and Swenson determined P p (x) 1/2 for many values of x ' < 260,000 on an IBM 7090 and their method of inter- pretation leads them to suggest that 6 may be as small as .3. This seems unlikely from the method of interpretation used in the present paper. The present calculations on an IBM 709^ include P,(x) for ri k = 2, 3, k, 5, 6 and all integer x < 250,000 (x 1 / 2 < 500); isolated values of P 2 (x) for x 1 / 2 < 10,000,000; and a few values of P_(x) for x 1 / 2 < 9000. The results of this work show that the calculated values follow the theoretical limits quite closely. The results for k = 2 fail to indicate that e p is less than Hua's bound of 13/^0. For k = 3 the most reasonable conclusion is 5 < e < .6, -2- Efficient algorithms for various combinations of k and x are presented in section 2. Section 3 is composed of computing methods and section k contains conclusions. 2. Counting Algorithms The most efficient method of evaluating A, (x) depends on the range of k and x and upon whether isolated values or a large number of consecutive values are desired. The following formula is similar to one given by Walfisz [1]. A 1 (0) = 1 (5) A 1 (x) = A 1 (x - 1) + 26 ( x ) , c/ \ Jl for x a perfect square where o(x) =i_ ,, . [0 otherwise [nTx] A R (x) = A^x) + 2 I A k _ 1 ( x " i ) i=l Formula (5) provides the basic method of calculating A (x). For large values of x some terms of the summation may be larger than the fixed- point single-word capacity of the computer (2 - 1 on the IBM 709^). This difficulty can often be remedied by defining R, (x) to be the number of points K. (j , J , . . . ,J ) such that k p Z J^ = x (6) i=l X It is evident that ■3- r -\(x) - A v( x - l) for x an integer V x ) = j (7) ^ otherwise [x] A-([x]) = Z R (i) * 1=0 k Also R 1 (o) = l (8) R 1 (x) = 26 (x) [>/x] p V X) = R k-l (x) + 2 E R k-l (x ~ i } i=l Note the similarities of (5) and (8). By changing initial values the same procedure may be used for either A,(x) or R (x). The next formula makes use of the symmetries involves in the set L k (x) resulting from permutations and negatives of ordered k -tuples. Define L k (x) ={(J r J 2 >...,J k ) 6 \(*)\0< J x < J 2 < ... < J k | (9) and let MCJ^, J^, . . . , J fc ) be the number of distinct permutations and negations of J^, J^, ..., J . Then we have M(J r J 2 , . . . , J k ) = 2 k - n (°) -^JJ (10) k n n(p)! p=0 where n(p) is the number of i for which J. = p. Thus if Y = (j J J ) 1 ~m 1' 2' ' ' ' ' m V x) E 1=2 M(Y, ) (11) -k- By rearranging (9) and (ll) it follows that [x] \(x) = 1 + E E ^Ir-i'V (12) Jk=1 4_ lG L. .(x-J 2 ) >k-l k-l x k-1— k Now if J < J then M(j, , J , . . . , J ) = 2kM(j, , J , . . . , J .). Similarly, if k-i k-i+1 : ... = J k , then M(J x ,J 2 ,...,J k ) = 2 ( . )M( J ± , Jg, . . . , J R _.) . Thus we have [X] ^ , V \(x) - 1 + Z Z 2 ± (") E M(Y k _.) (13) J * =1 1=1 2 k -i £L i-i< x - iJ k 2 ) J k-i Z < It is convenient to note that S k (x,oo) = A^x) (16) S k (oo,j) = (2J - l) k Formula (15) used with (l6) is the basic method for taking advantage of symmetries among the points of L,(x). However (15) can be efficiently modified still farther by noting that for mJ < Z, S (Z,J) = S (00, j). Thus (15) becomes J-l m . N m . s m (z,j) = 1 + E E 2 1 ( m )s .(z - ij 2 ,j ) + E E 2 x ( m )s .(»,j ) (it) J =N + 1 i=l X m " 1 m m J =1 1=1 V ^ m m m where N = min( [s/z/m '], J - l). Now applying (l6) and changing an order of summation we get m J-l . N m s a (z,j) - 1 + E E ^(^.(z - i^j ) + E E (- ) 2 1 (2j m - ir (18) i=l J =N+1 x m x m m j =1 i=1 1 m m m By the binomial theorem and the fact that Z < -* S (Z,j) = we have -6- m J-l . Q N S (Z,J) = 1+ 2 Z 2 1 ( m )S .(Z - iJ 2 ,J ) + Z (2J + l) m - (2 J - l) m i=u mwL V m " 1 m m j =i m m m MIN. m i = 1 + Z Z 2 X ( m )S .(Z - i«lV ) + (2N + 1)" - 1 (19) . . _ „ _ l m-i nr m 1=1 J =M+1 m where MIN. = min( [n/z/1*], J - l). However, when (19) is used recursively to compute A (x) as defined in (16), the following relation holds: N s ([-v/z/m' ],J - l) = [v/z/m' ]. Clearly the relation holds for m - k. Now assume that it holds for id m m < k. Then Jz/m' /x72r:r + 8 Z [v/x - J 2 '] (21) J=[n/x/2]+1 This formula was known to Gauss. Formula (20) is of course ideally suited to programming for an algorithmic compiler; however, it was thought that a machine-oriented language would be more efficient. Consequently it was coded in SCATRE for the 709U and used for isolated values of A (x). -7- There is an extension of (2l) which should be very valuable for computing isolated values of Ao(x). It is possible, for certain x, to compute A (x + u) for all u subject to |u| < 2.8 x ' in nearly the same time necessary to compute A (x) alone. For x = 10 , the largest argument used in this work, this would have made available over 17,000 results in about twice the time required for the single value. Needless to say this is a significant improvement. Unfortunately this method was not known when the computations were done for this paper, but the method has been used for small x. Define the following: 2 2 x = r + 2r = (r + l) - 1, r an integer (22) U < 2 s/2r = 2.8: 1/k -U + 2 < u < U and W = [six - J 2 ], where J is used in the context of (2l). We shall establish the following theorem: "W + 1, u > (W + l) 2 - (x - J 2 ) [six + u - J 2 ' ] = ^x - J 2 ' - 1 > n/x - r 2 ' - 1 = \/2r - 1 -8- Thus -2W < u < 2W + 2 and because u and W are integers -2W + 1 < u < 2W + 1 And thus x + u - J and -x-J 2 -2W+l>x-J 2 -2 n/x - J 2 ' + 1 = (n/x - J 2 " - l) 2 > (W - l) 2 and so W - 1 < six + u - J d ' < W + 2 Thus W - 1 < [sjx + u - J 2 ' ] < W + 1 and the theorem is proved. It is clear when each of the three possible cases holds. The true value of this method lies in computing A~(x + u) for all defined u simultaneously. Define Q(0) = q(u) =0, -U+2 u ( , [>/(x+u)/2'] Q(u) = Q(0) + Z q(v) + [\/u - 1 ] - Z [n/x + u - J 2 '] v=l J=[n/x/2 ']+l Q(-u) = Q(0) + Z q(-v) + Z [n/x - u - J^ ] v=l J=[n/(x-u)/2']+1 And then, for all u, A 2 (x + u) = 1 + k[slx + u'] + 4[n/(x + u)/^] 2 + 8Q(u) (23) 3. Computer Methods and Numerical Results The time-consuming part of computing A (x) by any of the methods mentioned in this paper is evaluating [s/x]. However from the logical order of summation it happens that successive arguments are often close together. Also, most square root routines are floating-point whereas exact fixed-point results are necessary for this program. Thus it may he possible to improve upon the efficiency of the square root routine by using fixed-point operations and by using the previous result as a first approximation for the current argument. Using this and the identity (N + l) 2 = N 2 + 2N + 1 -10- A, (x) may be calculated on a binary machine with no multiplications and no 1/2 numbers except the sum larger than x ' . This procedure was developed independent by Keller and Swenson [6] and the present author. It is particularly useful for employing (23). Keller and Swenson present the necessary algorithm and a basic derivation of the process. For the current paper two methods were used for determining A, (x). For calculating isolated values (20) and the above method were used. Special square root routines were used throughout. The time to compute \.( x ) was on the order of T(k,w) = w k x (k-l) / 2 where 1 - l//k w, , = r- 1 w. k+1 k k It should be noted that w, decreases as k increases. k When a large quantity of consecutive values was desired (5) and (8) were used. An IBM 1301 Disc File was available for additional storage. A particularly desirable aspect of this Disc File is the opportunity of using one portion of core storage for computation while moving data between the Disc File and another portion of core. For the present problem this effectively created one million words of core storage. The time required to compute A (x) for 2 < k < K and all integer x < X is T(K,X) = w(k - l)x 3 / 2 • (23) Using this method 3-1/2 hours were required for T (6,250000). Formula (8) was used with the assumption that IV. (x) < 2 . This assumption was violated near Rg (Uo,000). -11- Integer arithmetic was used exclusively for A, (x) in all programs and it is expected that all values are correct. Complete agreement was noted for all values published in [6]. Similar agreement existed with [k] except for A (l800 ), the largest argument published "by [*+]. This value was calculated twice for this paper, each calculation requiring three minutes. k. Conclusions Q The problem of establishing that P, (x) = 0(x ) is equivalent to finding X a sequence <(X.,Y. )>. , such that 11 1=1 |P k (x) I < Y. + 0(x y ) for x < X. +1 (2k) and lim — r < +°° i-*» X. 1 Since L-^(x) is composed only of k-tuples of integers, \(x) is piecewise constant over [n,n + 1], where n is an integer. Thus for x € [n,n + 1] P k (x) = (P k (x) - P k (n)) + P k (n) = C(x k / 2 - n k / 2 ) + P k (n) = P k (n) + 0(x k / 2 - l) (25) However by (h) , 6 > k/2 - 1, the sequence of "extreme" points (N , I P (N ) I ) such that |R (n) | < |p. (N. ) | for n < N. , satisfies the first 1 xi 'k 1 — ' k 1 ' l+l requirement of (2k) for all 6 of interest. This sequence is uniquely determined, given an initial element, and is equivalent to the set of (N, |p ,(n)|) for which |P k (n)| < |P k (N)| for n < N. Thus N. is the first integer for which ' P k^ N i^' K l P k^ N i+l^* A table of the first 50 extreme points for k = 2, 3 is given in Table 1. -12- TABLE 1. FIRST 50 EXTREME POINTS FOR k = 2, 3 X. 1 A 2 (X.) p 2 ( Xi ) X. 1 A 3 (X.) P 3 ( Xi ) 1 5 2 1 7 3 2 9 3 2 19 7 5 21 5 5 57 10 10 37 6 6 81 19 20 69 6 11+ 251 32 24 69 6 21 437 34 26 89 7 29 691 37 1+1 137 8 30 739 51 53 177 10 5^ 17^3 81 130 1+21 13 90 3695 119 lk9 1+81 13 13^ 6619 122 205 657 13 155 8217 134 23U 7I+9 li+ 17U 9771 157 287 885 17 230 14771 160 340 1085 17 23l+ 15155 l6l 410 1305 17 251 16831 174 1+25 1353 18 270 18805 221 480 1I+89 19 342 26745 252 586 1861 20 37^ 30551 254 81+0 2617 22 l+6l 41755 294 850 2693 23 1+9^ 46297 305 986 3125 27 550 5^339 309 1680 521+9 29 666 72359 364 181+3 5761 29 750 86407 371 2260 7129 29 810 96969 405 2591 8109 31 990 131059 580 3023 9I+65 32 1890 344859 682 3021+ 9I+65 35 2070 395231 734 34oo 10717 36 21+86 519963 756 3539 121+01 37 2757 607141 763 3960 121+01 1+0 2966 677397 776 5182 16237 43 3150 741509 959 5183 16237 1+6 3566 893019 1028 7920 2U833 48 3630 917217 1105 9796 30725 50 4554 1288415 1120 11233 35237 53 1+829 14-06811 1170 1I+883 1+6701 55 5670 1789599 1205 15119 I+7I+1+I 57 5750 1827927 1550 15120 l+71+l+l 60 8154 3085785 1570 19593 61U93 60 8382 3216051 1576 21600 67797 61 8774 3444439 1851 21603 67805 63 8910 3524869 1930 2160I+ 67805 66 10350 4412643 2028 22177 69605 66 10710 4645127 24o4 28559 89653 68 1573^ 8269399 2411 28560 89653 71 15750 8282167 2565 31679 99I+I+9 74 16302 8721339 2675 31680 99I+I+9 77 17550 9741669 2895 38015 119349 79 23310 14910309 2905 38016 119349 82 2389I+ 15474065 2940 38017 119349 85 21+174 15746999 3133 -13- 1/2 The number of extreme points for x < 250,000 (x ' < 500 ) is Number of Extreme Points 2 3 76 80 1+ 5 6 171 i+36 k7k (x < U0,000) This reflects the increasing value of 6 and perhaps more "regularity" for the K. higher values of k. The problem of showing P,(N.) k 1 lim — < +c i-w> N 1 is equivalent to showing lim (log |P (N.)| - 6 log N.) < +00 (26) . K 1 1 1-900 Graphically this corresponds to finding a straight line with slope 6 which majorizes the points (log N., log |p (N.)|). The sequence of extreme 1 K. 1 points for x < 250,000 is shown in Fig. 1. Only a sample of the points for k = 5> 6 are shown. The lines drawn represent the minimum slopes which parallel the extreme points. In addition for k = 2, 3, ^ the theoretical minimums for 6 (see (h)) are shown. If 6 is estimated from these points the results are e 2 = 0.325 = 13A0 (27) e 3 =0.60 = 3/5 e^ = 1.06 e 5 = 1.52 e 6 = 2.01 -ik- Q FIGURE 1. EXTREME POINTS FOR x < 250,000 AND k = 2, 3, b, 5, 6 -15- The accuracy of visual estimation limits this method to a precision of at most + .01. For instance in Fig. 1 for k = 2 a line with slope 13/^+0 would be indistinguishable from one of slope l/3. From these results it seems likely that will approach (k/2 - l) for large k, as expected from (h) . Indeed the 1/2 extreme points for k = h above log x = h (x ' = 100) are majorized by the line with slope 1.00. That the extreme points for k = k, 5, 6 eventually approach slope (k/2 - l) suggests the possibility that 9 and 9 are also smaller than the above estimates. With this in mind P p (x) was calculated from some 250 values of x 1 ' 2 < 10,000,000 and P_(x) was determined for about 20 values of x ' < 9000. Some of these values are given in Table 2. From these values approximate extreme points were chosen in the same manner as the true extreme points. The approximate extreme points are not a subset of the true extreme points. The 1/2 approximate extreme points for x ' < 10,000,000 are shown with the true extreme 1/2 points for x ' < 500 in Fig. 2. If only the approximate extreme points are considered one is led to the conclusion (as was done in [k]) that "0 = l/h is not inconsistent with observed results." It is also apparent that 6 500. Thus it is not reasonable to conjecture from these results that Q is appreciably less than Hua's bound of 13/^0. A logical conjecture based upon these results is p > .30. The approximate extreme points for k = 3 suggest that . 5 < 9^ < .6, but there are too few points from which to extrapolate with assurance. For instance half of the points qualify as approximate extreme points. The time required to calculate more values of P^(x) would be prohibitive. -16. TABLE 2 1/2' x ' A 2 (x) l? 2 (x)J 1/2 x ' A 3 (x) |P 3 (x)| 1000000 311+159261+9625 3965 1000 1+1887811+37 8768 1500000 70685831+6591+5 1+632 1200 7238202017 271+57 2000000 12566370610285 1+071+ 11+00 lll+9l+026l89 1*033 2500000 19631+951+076697 8239 1600 17157266213 lQk66 3000000 28271+33387381+1 81+67 1800 21+1+2898221+9 1+2225 3500000 381+81+509999277 7198 2000 33510290993 3061+5 1+000000 502651+821+51357 6080 2500 65I+I+9818205 2871+5 1+500000 63617251226505 8688 3000 113097275709 59820 5000000 78539816333093 6652 3500 17959^325^65 51+565 5500000 950331777621+29 8662 1+000 2680821+71+393 98713 6000000 113097335520185 901+8 1+500 3817031+53381 5^030 6500000 13273228960621+1 7928 5000 523598707861 67737 7000000 15393801+0012805 13095 5500 696909887157 83161+ 7500000 176711+586751+1+01 10025 6000 901+7785253^5 158889 8000000 201061929820913 8831+ 65OO 115031+61+27953 82036 8500000 226980069212125 9738 7000 11+3675^8853 91389 9000000 251+1+69001+93081+5 9928 7500 17671^5772565 95079 9500000 283528736973257 13222 8000 21I+I+660I+22929 161922 9600000 28952917891+1+573 10262 85OO 25721+1+0705977 78537 9700000 2955921+52772029 ^235 9000 3O5362785I+38I 201+908 9800000 3017185 58U38929 11835 9900000 307907^95961+805 13531 10000000 31^159265350589 8390 -17- Q^3 CD O ^^ * ' Y^ 2 =l3 /40 • • 2 = 1/4^ • • 2 ^ \ i > ( 3 l( D 12 14 LOG x FIGURE 2. DISTRIBUTION OF TRUE AND APPROXIMATE EXTREME POINTS FOR k = 2 BIBLIOGRAPHY 1. A. Walfisz. Gitterpunkte in Mehrdimensionalen Kuglen , Monografie Matematyczne, v. 33> Panstwowe Wydawnictwo Naukowe, Warsaw, 1957. (Reviewed in Mathematical Reviews , v. 20, n. 6, 1959> p. 8326. ) 2. J. R. Wilton. "The Lattice Points of a Circle: an Historical Account of the Problem." Messenger of Mathematics , v. 58, 1929 j p. 67-80. 3. L. K. Hua. "The Lattice-points in a Circle." Quart. J. Math., Oxford Ser . , v. 13, 19^2, p. 18-29. h. W. Fraser and C. C. Gotlieb. "A Calculation of the Number of Lattice Points in the Circle and Sphere." Math. Comp . , v. l6, 1962, p. 282-290. 5. H. L. Mitchell, III. Numerical Experiments on the Number of Lattice Points in the Circle . Appl. Math, and Stat. Labs., Stanford University, Stanford, California, 1961. 6. H. B. Keller and J. R. Swenson. "Experiments on the Lattice Problem of Gauss." Math. Comp . , v. 17, 1963, p. 223-230. -19- T^EBS.TYOF.LL.NO.S-URBANA 510.»4ll6«nD.C0M^« lUiac l» « processor ot visub ■PI 30112088398208