Prime numbers and primality testing Yahoo Group Programming Question =============================================== William F Sindelar Message 1 of 32 Jan 15, 2002 ----------------------------------------------- Hi Everybody I would welcome comments on a probably stupid question, but first I want to thank Ignacio Canestro, Jim Fougeron, Tom Hadley, and Paul Leyland for their responses to my recent post "Is the next one prime?" Guys, I truly appreciate your doing this just to satisfy my curiosity! I often thought of downloading some of the huge number crunching programs that were mentioned on this site, but I seriously doubt my ability to understand and run them successfully. Jim, your little note "caveat utilitor" further reinforces that doubt. I understand that the thousands of digits primes that are announced here are of special forms which are amenable to manipulation according to fixed rules, and the result further checked by a different method, but somehow I feel uneasy about this. Maybe I may be able to understand these programs a bit better by posing this question: Suppose I have a thousands of digits integer which does not have a 3 or a 7 as the final digit, and I wish to be 100% certain it is a square integer. Also say that I can input this integer without error into my computer. Are there programs that can do this? Are there other ways of determining squareness besides extracting the root? Thanks and regards everyone Bill Sindelar =============================================== Paul Jobling Message 2 of 32 Jan 15, 2002 ----------------------------------------------- Bill, > 7 as the final digit, and I wish to be 100% certain it is a square > integer. Also say that I can input this integer without error into my > computer. Are there programs that can do this? Are there other ways of > determining squareness besides extracting the root? You can do some quick tests that throw out numbers that are *not* squares. If you have got the number in a binary (rather than a decimal) format, look at the bottom 6 bits (ie the number mod 64) - there are only a small number of allowable values. If you can do a quick bit of long division to find the number mod 63 and mod 65 as well (63 and 65 are off the top of my head, by the way - Cohen's book has got the details), you can do a similar thing to throw out even more non-squares. Those few numbers that are left must be tested by finding an integer square root (using Newton's method), then squaring this and comparing it to the original number. Obviously a big integer package and some programming is required, though Pari is probably a lot easier if you don't want to program (and it has an issquare() test that does all of that for you). Regards, Paul. __________________________________________________ Virus checked by MessageLabs Virus Control Centre. =============================================== Paul Leyland Message 3 of 32 Jan 15, 2002 ----------------------------------------------- > computer. Are there programs that can do this? Are there other ways of > determining squareness besides extracting the root? That is an interesting question. Extracting the root is straightforward, as Paul Jobling has posted, but can we test in such a way that we know for certain whether an integer N is a perfect square, without calculating the square root? The answer is "yes", and here is one approach. Choose a set of small integers n_i, all co-prime to each other (i.e. any one has no factors in common with any other, thus 24 is coprime to 35 even though neither is prime) such that their product Pi(n_i) > N. By the Chinese Remainder Theorem, there is a unique set of r_i such that r_i = N mod n_i. If every single one of these r_i is a quadratic residue modulo n_i, N must be a perfect square. If at least one r_i is not a quadratic residue, N is not a perfect square. You made a start at this approach by noting that 3 and 7 are not quadratic residues modulo 10. For completeness, you should have included 2 and 8 but I assume you were interested only in odd numbers. Paul =============================================== Phil Carmody Message 4 of 32 Jan 15, 2002 ----------------------------------------------- On Tue, 15 January 2002, William F Sindelar wrote: > Suppose I have a thousands of digits integer which does not have a 3 or a > 7 as the final digit, and I wish to be 100% certain it is a square > integer. Also say that I can input this integer without error into my > computer. Are there programs that can do this? Are there other ways of > determining squareness besides extracting the root? It depends what you want from the 'program'. For many of us on the list the 'program' is nothing more than a library that we use in combination with whatever is our favourite programming language. I know that C and C++ programmers frequently use either the GMP (GNU Multi-Precision, www.swox.org) library or Mike Scott's Miracl library (from Shamus Software, www.shamus.ie) for these tasks, but there are many other libraries (hey, I still occasionally use Lenstra's LIP, and am proud of the fact!). Each of these libraries contain optimised issquare()-like functions, which will do what you want. For others on the list the 'program' is a fully integrated mathematical algebra system, which will contain programming primitives too. Examples include Maple, Mathematica, Pari/GP, Macsyma, and others. These too contain issquare() functions too. The particular technique each library or program uses may vary slightly, but the general principle is the same: In order: 1) throw away trivially easy cases (mod 2^n) - O(1) time 2) throw away easy cases (mod p or p^2 for small p) - O(N) time 3) use brute force. square root and resquare. - O(??) time. There are 3 choices for (1): count trailing zeroes; look at the last n bits in a lookup table; some self-indulgent bit-twiddling I discovered about a 6 months ago. This stage says 'NO' about 86% of the time. There is a lot of choice for stage (2), revolving around a) how many primes, which primes b) how you perform the modular reduction The two issues are intertwined. e.g. If you chose 3^2, 5, 7 and 13, then you can simply reduce the number mod (3^2*5*7*13) = 2^12-1. Reduction mod 2^12-1 is blazingly fast on a binary computer as it requires no multiplication or division. When the value is reduced mod 2^12-1, it's trivial to pull the mod 63 and mod 65 values out, and then look them up in a bitmask. 2^24-1 gives you even more small primes, as do 2^30-1, 2^48-1 and 2^60-1. Which number you chose will depend on your CPU architecture, and whether you want to do one of these stages or two. It's hard to do more than two of these blazingly fast stages as it's harder to isolate the primes in expressions such as the above, you end up with repeats. Once 2^n-1 type reductions are exhausted, there's always a traditional modular reduction by carefully chosen values. The values will be the largest products of small primes that you've not yet tested that you can fit into a single machine word (to make the reduction as fast as possible). Again, when you have the value reduced by the composite modulus, simply take remainders modulo the prime factors of that modulus and look up in a bitmask. Each prime or prime power gives on average just under a 50% rejection rate. (3 is poor, but 3^2 is better, 5 likewise.) So assuming that you can stick about 10 primes together, you're laughing at having handled about 99.8% of cases very quickly indeed. The final stage could in theory be very slow compared to the above as it's quite rarely needed assuming the above easy cases do reject the proportions that they are expected to do. Note, the first 2 stages are very quick to say "No", but _don't_ ever actually say "yes", you can only get the "yes" answer from the brute force stage 3. Phil Don't be fooled, CRC Press are _not_ the good guys. They've taken Wolfram's money - _don't_ give them yours. http://mathworld.wolfram.com/erics_commentary.html Find the best deals on the web at AltaVista Shopping! http://www.shopping.altavista.com =============================================== Phil Carmody Message 5 of 32 Jan 15, 2002 ----------------------------------------------- On Tue, 15 January 2002, "Paul Leyland" wrote: > > computer. Are there programs that can do this? Are there other ways of > > determining squareness besides extracting the root? > > > That is an interesting question. Extracting the root is > straightforward, as Paul Jobling has posted, but can we test in such a > way that we know for certain whether an integer N is a perfect square, > without calculating the square root? The answer is "yes", and here is > one approach. Remember that I said that the only place you can get a "yes" was from my brute force stage 3. Paul has just just prodded me to indicate that I was being sloppy and should have prefixed that with "assuming that N is larger than the product of the small primes". Because if you don't assume that... > Choose a set of small integers n_i, all co-prime to each other (i.e. any > one has no factors in common with any other, thus 24 is coprime to 35 > even though neither is prime) such that their product Pi(n_i)> N. By > the Chinese Remainder Theorem, there is a unique set of r_i such that > r_i = N mod n_i. If every single one of these r_i is a quadratic > residue modulo n_i, N must be a perfect square. If at least one r_i is > not a quadratic residue, N is not a perfect square. Spot on. Effectively this is a (deliberately wishy-washily defined) "going all the way" version of my stage 2. You can get "no" answers quite quickly, but the "yes" answers require you to go the whole way. Remembering that if N is large (and we were talking large N), then operations on all the n_i becomes ever more time consuming, and sometimes it's just easier to bite the bullet and take the sqrt hit. So in general this method isn't widely used. However, there are methods of doing bignums where Paul's method is better, methods for example where modular reductions are 'free'. Unfortunately using these bignum systems some other operations which are frequently used (div/rem type things) become much more complicated. See Knuth TAoCP vol.2 SNA, probably under "mixed radix representations", or similar. Sometimes it's worth popping into mixed-radix representation to do a whole bunch of stuff that's swiftest there, and then to pop back into a simple bignum representation for other operations. However, the conversion overhead is non-trivial. My handy hint: for version 1, K.I.S.S. Phil Don't be fooled, CRC Press are _not_ the good guys. They've taken Wolfram's money - _don't_ give them yours. http://mathworld.wolfram.com/erics_commentary.html Find the best deals on the web at AltaVista Shopping! http://www.shopping.altavista.com =============================================== Paul Leyland Message 6 of 32 Jan 15, 2002 ----------------------------------------------- Phil Carmody wrote: > Remember that I said that the only place you can get a "yes" > was from my brute force stage 3. Paul has just just prodded > me to indicate that I was being sloppy and should have > prefixed that with "assuming that N is larger than the > product of the small primes". Because if you don't assume that... Not so fast young man... It's time for me to wipe the egg off my face and issue a retraction. Marcel Martin was kind enough to send me this email: ] Let N = 79, n0 = 3, n1 = 5, n2 = 7. ] Now ] 3*5*7 = 105 and 105 > 79, ] 79 mod 3 = 1 and Jacobi(1,3) = 1, ] 79 mod 5 = 4 and Jacobi(4,5) = 1, ] 79 mod 7 = 2 and Jacobi(2,7) = 1, ] but ] 79 is not a perfect square. Oh dear. I recant. It remains an interesting question as to whether squareness can be proved without extracting a square root. I'll think more before attempting to answer it a second time. Paul =============================================== Phil Carmody Message 7 of 32 Jan 15, 2002 ----------------------------------------------- On Tue, 15 January 2002, "Paul Leyland" wrote: > Not so fast young man... Young man. At your age, at our age, young man... > It's time for me to wipe the egg off my face and issue a retraction. OK, I retract my part retraction too. > Marcel Martin was kind enough to send me this email: As always, hats off to Marcel. > ] Let N = 79, n0 = 3, n1 = 5, n2 = 7. > ] Now > ] 3*5*7 = 105 and 105 > 79, > ] 79 mod 3 = 1 and Jacobi(1,3) = 1, > ] 79 mod 5 = 4 and Jacobi(4,5) = 1, > ] 79 mod 7 = 2 and Jacobi(2,7) = 1, > ] but > ] 79 is not a perfect square. > > Oh dear. I recant. Oh good. I've have a chance to investigate and learn something. Is your glass half empty, Paul? (/me gives Paul a top-up from a spare pint of Stella...) > It remains an interesting question as to whether squareness can be > proved without extracting a square root. I'll think more before > attempting to answer it a second time. It seems clear that if Jacobi(n%b, b) = 1 for all bsqrt(n) (done without taking a square root, of course!) In Marcel's case the smallest square denier is b=2^2, for reference. Those who like to whip out their copies of Pari/GP in public may like to find some numbers which have ever increasing smallest-square-deniers. I'm at work, so you've all got a good several-hour head-start! Phil It is against US Department of Agriculture regulations to advertise or sell as "Prime Rib" any cut of meat containing a non-prime number of ribs. Find the best deals on the web at AltaVista Shopping! http://www.shopping.altavista.com =============================================== jbrennen Message 8 of 32 Jan 15, 2002 ----------------------------------------------- --- In primenumbers@y..., Phil Carmody wrote: > In other words, given non-square n, what's the size of the smallest b that is a square-denier (yes, neologism) for n. > > Is there a case where the smallest square-denier for n is >sqrt(n) (done without taking a square root, of course!) > In Marcel's case the smallest square denier is b=2^2, for reference. > Cases of non-square n with smallest square-denier > sqrt(n): n square-denier(n) --------------------- 2 3 3 4 5 3 6 4 7 4 8 3 10 4 12 5 13 5 15 4 21 8 24 7 40 7 45 7 60 9 This is likely a complete list. There are no others < 10^7. Other numbers which are "record-setters" as the smallest n with a particular square-denier: n square-denier(n) --------------------- 184 11 280 13 364 16 1456 17 3124 19 5236 23 17185 25 25249 29 49504 37 233776 41 364144 43 775369 47 3864169 53 8794864 61 So, in some sense, these numbers could be considered pseudo-squares? =============================================== Phil Carmody Message 9 of 32 Jan 15, 2002 ----------------------------------------------- On Tue, 15 January 2002, "jbrennen" wrote: > --- In primenumbers@y..., Phil Carmody wrote: > > In other words, given non-square n, what's the size of the smallest > b that is a square-denier (yes, neologism) for n. > > > > Is there a case where the smallest square-denier for n is >sqrt(n) > (done without taking a square root, of course!) > > In Marcel's case the smallest square denier is b=2^2, for reference. > > > > Cases of non-square n with smallest square-denier > sqrt(n): Hehe, I knew I'd be beaten to this. Hmmm, must buy laptop... Incidentally, what algorithm did you use? I'd guess a straight sieve would work (which would be blazingly fast)? The CRT method seems to fan out too much to be practical. What range do you think you could practically search in an hour, say? (Assuming a sieve runs at a constant rate.) > n square-denier(n) > --------------------- > 2 3 > 3 4 > 5 3 > 6 4 > 7 4 > 8 3 > 10 4 > 12 5 > 13 5 > 15 4 > 21 8 > 24 7 > 40 7 > 45 7 > 60 9 > > This is likely a complete list. There are no others < 10^7. I'd guess that's the case too. > Other numbers which are "record-setters" as the smallest n > with a particular square-denier: > > n square-denier(n) > --------------------- > 184 11 > 280 13 > 364 16 > 1456 17 > 3124 19 > 5236 23 > 17185 25 > 25249 29 > 49504 37 > 233776 41 > 364144 43 > 775369 47 > 3864169 53 > 8794864 61 > > So, in some sense, these numbers could be considered pseudo-squares? The contiguous nature of the numbers used to test them (i.e. all smallest prime powers) does make the definition a little arbitrary, but still useful. I genuinely believe that this is an interesting investigation and worth further study. It only becomes practical if the growth-rate picks up, which alas I don't believe it will. Phil Don't be fooled, CRC Press are _not_ the good guys. They've taken Wolfram's money - _don't_ give them yours. http://mathworld.wolfram.com/erics_commentary.html Find the best deals on the web at AltaVista Shopping! http://www.shopping.altavista.com =============================================== Hadley, Thomas H (Tom), ALINF Message 10 of 32 Jan 15, 2002 ----------------------------------------------- May I try to summarize what has been found so far on this topic? (for all of us who have a hard time keeping up with you blazingly fast thinkers.) A. We're trying to find out if a number, N, is a square integer, without going through a square root calculation. B. We can prove that a number is not square if, for some number b, N mod b (or N%b) is NOT a quadratic residue of b. This means that it is impossible to find a number that when squared and reduced mod b, gives N%b. For example, no number when squared and reduced mod 7 can give you a 6. Six is not a quadratic residue of 7. Thus, no number which is equivalent to 6 (mod 7) can be a square. C. The Jacobi symbol, written in ascii as Jacobi(x,y) is a convenient way to represent whether or not x is a quadratic residue of y. A value of 1 means it is a quadratic residue, a value of -1 means it is not a quadratic residue and a value of 0 means a divides b evenly(?). D. It appears that all non-squares have some smallest number such that Jacobi(b,n) = 1. I don't know if this is proven or not. E. An open question is, "How many numbers b do you need to test with Jacobi(b,N) to prove that N is square? All b < N? All b < sqrt(N)? All b < 100? F. Paul and Phil thought at first that you only need to test enough b's such that their product is bigger than N, but Marcel came up with a counterexample. G. Another open question. "Does a small subset of numbers exist that when tested can prove N is a square?" H. Are there any "Carmichael squares", numbers N for which Jacobi(b,N)=1 for all b < N, but it is not square? I. A square-denier is a number b such that Jacobi(b,N) is -1, but for all a wrote: > Hehe, I knew I'd be beaten to this. Hmmm, must buy laptop... > Incidentally, what algorithm did you use? I'd guess a straight sieve would work (which would be blazingly fast)? The CRT method seems to fan out too much to be practical. What range do you think you could practically search in an hour, say? (Assuming a sieve runs at a constant rate.) I used a crude sieve which could probably be improved in speed by an order of magnitude easily. I allocated a sieve array of length 10^7, set all entries to 1 except for the perfect squares. Then for each number N beginning with 3, find the set of quadratic non-residues modulo N. Clear to 0 all entries in the sieve array corresponding to those non-residues. The entries in the sieve array still equal to 1 are those non-squares with "square-denier" greater than N. Repeat with the next N. With suitable modifications to allow sieving on a sliding window, and some optimization of the sieving to be cache-friendly, I would guess that a run rate of at least 10^10 per hour would be possible. At some point, probably around p = 30 to 50, it probably becomes better to just look up Jacobi symbols on the few candidates which slip through the sieve to that point. Further sieving is counterproductive; remember that with each new prime, we must eliminate essentially half of the sieve array. =============================================== Phil Carmody Message 12 of 32 Jan 17, 2002 ----------------------------------------------- On Tue, 15 January 2002, "jbrennen" wrote: > With suitable modifications to allow sieving on a sliding window, > and some optimization of the sieving to be cache-friendly, I would > guess that a run rate of at least 10^10 per hour would be possible. > > At some point, probably around p = 30 to 50, it probably becomes > better to just look up Jacobi symbols on the few candidates > which slip through the sieve to that point. Further sieving is > counterproductive; remember that with each new prime, we must > eliminate essentially half of the sieve array. It appears that my generic sieve can't cope with non-prime bases in its filters. (i.e. the prime power bases). So I'm running at about 5% of my top speed currently :-( (Note to Paul - it appears that it works with -s1, i.e. no 'special' stage - any ideas.) Even at 5% speed, I can't guarantee the correctness of the following: I don't suppose you could confirm or deny, could you, please Jack? 10869664 %64 =32 32384209 %67 =27 105361344 %71 =42 173899609 %73 =31 425088976 %89 =56 2140195264 %97 =68 2805544681 %101 =11 10310263441 %103 =6 11940537961 %107 =77 38432362081 %113 =107 43395268849 %127 =39 51802119889 %131 =90 I'm looking for >131 results currently. From my experience of Euler trinomials, I think that if I can crank the sieve up to full speed results of over 200 should be possible to find. Phil Don't be fooled, CRC Press are _not_ the good guys. They've taken Wolfram's money - _don't_ give them yours. http://mathworld.wolfram.com/erics_commentary.html Find the best deals on the web at AltaVista Shopping! http://www.shopping.altavista.com =============================================== Jack Brennen Message 13 of 32 Jan 18, 2002 ----------------------------------------------- Phil Carmody wrote: > Even at 5% speed, I can't guarantee the correctness of the following: I don't suppose you could confirm or deny, could you, please Jack? > > 10869664 %64 =32 > 32384209 %67 =27 > 105361344 %71 =42 > 173899609 %73 =31 > 425088976 %89 =56 > 2140195264 %97 =68 > 2805544681 %101 =11 > 10310263441 %103 =6 > 11940537961 %107 =77 > 38432362081 %113 =107 > 43395268849 %127 =39 > 51802119889 %131 =90 > > I'm looking for >131 results currently. Well, I improved my sieve a bit, and I currently have a run rate of just over 20 million numbers a second, or about 7.5*10^10 / hour. This is on a Pentium II (350 MHz), so you can extrapolate to figure what a current model year PC might be able to do. Interestingly, my sieve was most efficient when I stopped sieving at prime power 11 (only using a sieve for 5,7,8,9,11), and then looked up precomputed values of the Jacobi() function for all prime powers >= 13. Also, the sieve "window" is small (2520 bytes) in order to fit in the P-II's level-one cache. I can confirm the last six of your results (my sieve program only considers those candidates with "square-denier" over 100), and add the next one: 299530084681 %137 =10 The sieve is still running, and is up to about 8.25*10^11, with the next "record" not found yet. =============================================== djbroadhurst Message 14 of 32 Jan 18, 2002 ----------------------------------------------- Jack Brennen wrote: > 299530084681 %137 =10 At last: the Eddington denial. > The reciprocal of the fine structure constant of > electromagnetism is close to 137. [Eddington] But not all that close... David =============================================== Phil Carmody Message 15 of 32 Jan 18, 2002 ----------------------------------------------- On Fri, 18 January 2002, Jack Brennen wrote: > Phil Carmody wrote: > > Even at 5% speed, I can't guarantee the correctness of the following: I don't suppose you could confirm or deny, could you, please Jack? > > > > 10869664 %64 =32 > > 32384209 %67 =27 > > 105361344 %71 =42 > > 173899609 %73 =31 > > 425088976 %89 =56 > > 2140195264 %97 =68 > > 2805544681 %101 =11 > > 10310263441 %103 =6 > > 11940537961 %107 =77 > > 38432362081 %113 =107 > > 43395268849 %127 =39 > > 51802119889 %131 =90 > > > > I'm looking for >131 results currently. > > Well, I improved my sieve a bit, and I currently have a run rate > of just over 20 million numbers a second, or about 7.5*10^10 / hour. > This is on a Pentium II (350 MHz), so you can extrapolate to > figure what a current model year PC might be able to do. Given the choice, I always sieve on my 533MHz alpha. > Interestingly, my sieve was most efficient when I stopped sieving > at prime power 11 (only using a sieve for 5,7,8,9,11), > and then looked up precomputed values of the Jacobi() function > for all prime powers >= 13. Also, the sieve "window" is small > (2520 bytes) in order to fit in the P-II's level-one cache. I have some crappy overheads that I need to get rid of per block so I've aimed my sieve at my 96K L2 instead (I think I use 72K of it). > I can confirm the last six of your results (my sieve program > only considers those candidates with "square-denier" over 100), > and add the next one: > > 299530084681 %137 =10 Thanks loads, I _really_ don't know why my sieve should work when cranked down to its slowest speed, but not when put into its turbo modes. It's a useful datapoint for debugging though. I ran overnight in molasses mode (dunno what speed, sorry, but it appears it's about 1% speed, not 5% speed) and got the following: Summary: 299530084681 %137 =10 (1204776760009 %139 =14) (2834563560049 %149 =59) 963619651801 %151 =3 7115492078809 %157 =80 16711966296601 %163 =105 40249953292201 %167 =73 (%169) (%173) 23150469048001 %179 =79 (paranthetical ones are 'late', or not at all) > The sieve is still running, and is up to about 8.25*10^11, with > the next "record" not found yet. This morning I removed prime powers from my sieve, only using the p^1 power, and added them into my post-processing (10^-11 of the time) stage instead. That permits me to crank it up to James Dean mode! (In theory, if I were doing a serious sieve, I probably have 2 more notches, but they would mean a block would take a couple of hours, or over a day respectively, with no mid-block restart capability). In James Dean mode, running at << 28814 phil 16 0 5736 5736 1072 R 93.6 2.2 9:22 sqrdeny.alpha >> 93.6% CPU, for 9 mins 22 seconds. I'm sieving a block of size: <<< DoSieve(79992282930560, 9223372032559808512) 98199866973121 %173 =168 99348981507241 %169 =141 #P sieved to 101409094987200 >>> at a rate of: <<< phil@megaspaz:GFN$ dc -e '101409094987200 79992282930560-60 9*20+/p' 38244307244 >>> 38 billion a second. Woh! Not even I expected that! And it's nice for that range to have filled in precisely the two values I was looking for, 169 and 173! Looks like it's gonna be a great day today! I shall run this until the evening, and then stop it, as the other 6.4% of CPU time is the more important sieve to me ;-) Phil Don't be fooled, CRC Press are _not_ the good guys. They've taken Wolfram's money - _don't_ give them yours. http://mathworld.wolfram.com/erics_commentary.html Find the best deals on the web at AltaVista Shopping! http://www.shopping.altavista.com =============================================== grmulkey Message 16 of 32 Jan 18, 2002 ----------------------------------------------- If you read msg 4790, there may be a slight misunderstanding involving the distinction between the quadratic residue concept and the Jacobi symbol used to evaluate it efficiently. Please see below. --- In primenumbers@y..., "Hadley, Thomas H (Tom), ALINF" wrote: > May I try to summarize what has been found so far on this topic? (for > all of us who have a hard time keeping up with you blazingly fast > thinkers.) > > A. We're trying to find out if a number, N, is a square integer, without > going through a square root calculation. > > B. We can prove that a number is not square if, for some number b, N mod > b (or N%b) is NOT a quadratic residue of b. This means that it is > impossible to find a number that when squared and reduced mod b, gives > N%b. For example, no number when squared and reduced mod 7 can give you > a 6. Six is not a quadratic residue of 7. Thus, no number which is > equivalent to 6 (mod 7) can be a square. > > C. The Jacobi symbol, written in ascii as Jacobi(x,y) is a convenient > way to represent whether or not x is a quadratic residue of y. A value > of 1 means it is a quadratic residue, a value of -1 means it is not a > quadratic residue and a value of 0 means a divides b evenly(?). But Jacobi(73,35)=Jacobi(3,35)=Jacobi(3,5)Jacobi(3,7)=(-1)(-1)=1, however 3 (and 73) are not quadratic residues of 35. One needs to use only prime powers for moduli, as several have in extensive calculations posted recently. > > D. It appears that all non-squares have some smallest number such that > Jacobi(b,n) = 1. I don't know if this is proven or not. Do you mean = -1 ? > > E. An open question is, "How many numbers b do you need to test with > Jacobi(b,N) to prove that N is square? All b < N? All b < sqrt(N)? > All b < 100? > > F. Paul and Phil thought at first that you only need to test enough b's > such that their product is bigger than N, but Marcel came up with a > counterexample. > > G. Another open question. "Does a small subset of numbers exist that > when tested can prove N is a square?" They would need to be prime powers. > > H. Are there any "Carmichael squares", numbers N for which Jacobi(b,N)=1 > for all b < N, but it is not square? > > I. A square-denier is a number b such that Jacobi(b,N) is -1, but for > all a quadratic residue. > > I hope I understand the topic correctly. Corrections are welcome. > > Tom Hadley Some nice questions here, enough for several PhD's, maybe. =============================================== Phil Carmody Message 17 of 32 Jan 18, 2002 ----------------------------------------------- On Fri, 18 January 2002, "grmulkey" wrote: > If you read msg 4790, there may be a slight misunderstanding > involving the distinction between the quadratic residue concept and > the Jacobi symbol used to evaluate it efficiently. Please see below. Note that 4790 had a bit of 2+2=5 logic in it. However, without that post, Jack and I would never have started hunting for the aberrant 5s! > --- In primenumbers@y..., "Hadley, Thomas H (Tom), ALINF" > wrote: > > May I try to summarize what has been found so far on this topic? > (for > > all of us who have a hard time keeping up with you blazingly fast > > thinkers.) > > > > A. We're trying to find out if a number, N, is a square integer, > without > > going through a square root calculation. > > > > B. We can prove that a number is not square if, for some number b, N > mod > > b (or N%b) is NOT a quadratic residue of b. This means that it is > > impossible to find a number that when squared and reduced mod b, > gives > > N%b. For example, no number when squared and reduced mod 7 can give > you > > a 6. Six is not a quadratic residue of 7. Thus, no number which is > > equivalent to 6 (mod 7) can be a square. > > > > C. The Jacobi symbol, written in ascii as Jacobi(x,y) is a > convenient > > way to represent whether or not x is a quadratic residue of y. A > value > > of 1 means it is a quadratic residue, a value of -1 means it is not > a > > quadratic residue and a value of 0 means a divides b evenly(?). > > > But Jacobi(73,35)=Jacobi(3,35)=Jacobi(3,5)Jacobi(3,7)=(-1)(-1)=1, > however 3 (and 73) are not quadratic residues of 35. One needs to > use only prime powers for moduli, as several have in extensive > calculations posted recently. Absolutely. Well spotted. I did read Tom's post all through, but was tainted by my own preconceptions. However, take care of the order of the parameters to the Jacobi function in comparison to the much-misliked-by-me (a,p) notation for Jacobi symbols. > > D. It appears that all non-squares have some smallest number such > that > > Jacobi(b,n) = 1. I don't know if this is proven or not. > > Do you mean = -1 ? Or does he mean a largest number? A non-square will always have a square-denier, as any prime larger than it denies its squareness. Whether is has one less than itself is unknown, but almost certainly true. > > E. An open question is, "How many numbers b do you need to test with > > Jacobi(b,N) to prove that N is square? All b < N? All b < > sqrt(N)? > > All b < 100? > > > > F. Paul and Phil thought at first that you only need to test enough > b's > > such that their product is bigger than N, but Marcel came up with a > > counterexample. > > > > G. Another open question. "Does a small subset of numbers exist > that > > when tested can prove N is a square?" > > > They would need to be prime powers. Not all numbers. Each prime power creates a roughly 50/50 filter, and all prime powers under consideration are coprime these filters can be shoved into the Chinese Remainder Theorem, and will always produce a complete set of solutions. The ratios decrease, but the sets of solutions grow without bound. e.g. my (broken but bodged) filter file for minimal square-deniers begins: 3 +2 : 0 1 5 +3 : 0 1 4 7 +4 : 0 1 2 4 11 +6 : 0 1 3 4 5 9 13 +7 : 0 1 3 4 9 10 12 17 +9 : 0 1 2 4 8 9 13 15 16 19 +10 : 0 1 4 5 6 7 9 11 16 17 ... Which means that modulo (3*5*7*11*13*17*19) there will be (2*3*4*6*7*9*10) solutions. (A working sieve would have prime powers not just primes). > > H. Are there any "Carmichael squares", numbers N for which > Jacobi(b,N)=1 > > for all b < N, but it is not square? Therefore no. It's _trivial_ to simply throw these filter criteria at the CRT, however, you have a vast number, and the chances of you finding the minimal one (and as we were looking for bounds, so minimality is intended) are next to zero. Having said that they're about as good as me and Jack trying a brute force sieving approach from zero. We have the advantage of sieves being cheaper than CRT solutions. It's just a matter of patience. Hmmm, I just noticed the following in my logs... 335953582561081 %193 =35 So I should be up to 10^15 by the evening - maybe a >200 minimal square denier is findable. (Hmmm, maybe I'll run the sieve until morning...) > > I. A square-denier is a number b such that Jacobi(b,N) is -1, but > for > > all a not a > > quadratic residue. We're only testing prime power a's though. However, that doesn't actually affect the above definition. > > I hope I understand the topic correctly. Corrections are welcome. > > > > Tom Hadley > > > Some nice questions here, enough for several PhD's, maybe. Hehe, I can fund myself if anyone wants to supervise me doing a Ph.D. on demonically fast sieving! }:-)~ Phil It is against US Department of Agriculture regulations to advertise or sell as "Prime Rib" any cut of meat containing a non-prime number of ribs. Find the best deals on the web at AltaVista Shopping! http://www.shopping.altavista.com =============================================== jbrennen Message 18 of 32 Jan 18, 2002 ----------------------------------------------- --- In primenumbers@y..., Phil Carmody wrote: > On Fri, 18 January 2002, Jack Brennen wrote: > > Well, I improved my sieve a bit, and I currently have a run rate > > of just over 20 million numbers a second, or about 7.5*10^10 / hour. > > I'm sieving a block of size: > ... > at a rate of: > ... > 38 billion a second. > I bow before the master. I should go throw my sieve program into a bottomless pit, and study at the FatPhil school of sieving, I guess. =============================================== Phil Carmody Message 19 of 32 Jan 18, 2002 ----------------------------------------------- On Fri, 18 January 2002, "jbrennen" wrote: > I should go throw my sieve program into a bottomless pit, > and study at the FatPhil school of sieving, I guess. 1) Use CRT as much as possible. (I could go to ~72 then ~130G/s by putting another prime or two in the CRT stage). 2) Saturate the pipelines for each tiny prime (i.e. do many in parallel, and maybe unroll the loop to do 2 words simultaniously). (Hint - use a machine with 31 registers ;-) ) 3) Don't waste time sieving things that are already fully zeroed. (You knew this already, but there's more to the issue than just the superficial interpretation.) 4) Balance the post-sieving precalulated rejections between bitmap-oriented techniques and sorted-list oriented techniques. My thanks as always to Paul J. for pushing me towards the bitmap technique, I was previously wastefully doing more 'tiny' primes in a fast loop, and then flipping to a sorted-list. (The function is called doPaul() in his honour!) 5) Don't waste time sieving things when you're ultra-sparse, simply come up with a constant-time yes/no answer once and for all. (Again, well known, but not really relevant for this task.) My runtime parameters were: <<< #I ZIP: 0 primes, relation = 0 mod 1 (expanded) #I SPESH: 7 primes, reduction = 37182145 -> 544320 #I SMALL@7: 7 dumb, 0 quick to 47 #I PAUL: 24 primes <=167 #I WHEEL: 0 other primes <= 0 #I FACTOR: 1.368196e-11 >>> meaning: (ZIP+SPESH) CRT reduced the search space by 544320/37182145 (SMALL) I used 7 'tiny' primes (however, I do 7 simultaniously in my loop, so only one iteration of an actual sieve), using only one version of my inner loop (there's a 'quick' version which can be faster on 2nd and subsequent runs through the array. (PAUL) The bitmap oriented post-sieving was 24 primes (WHEEL) No post-sieving was done with lists of things to reject (FACTOR) the ratio that passes through the above stages Alas the world release of the gensv code has been put back by this 'prime power base' bug. My faithful beta-testers are invited to race me to a proper solution (rather than my hack of doing prime powers later). Phil It is against US Department of Agriculture regulations to advertise or sell as "Prime Rib" any cut of meat containing a non-prime number of ribs. Find the best deals on the web at AltaVista Shopping! http://www.shopping.altavista.com =============================================== Hadley, Thomas H (Tom), ALINF Message 20 of 32 Jan 18, 2002 ----------------------------------------------- OK, I think I understand more now about the Jacobi symbol (after reading more in the Prime Pages Glossary about it). >On Fri, 18 January 2002, "grmulkey" wrote: >> If you read msg 4790, there may be a slight misunderstanding >> involving the distinction between the quadratic residue concept and >> the Jacobi symbol used to evaluate it efficiently. Please see below. > You're right. Here's what I found: The Legendre symbol, (a,p), not the Jacobi symbol, tells us whether or not a number is a quadratic residue of another number, however, the denominator (or second number in the (a,p) form) must be a prime. The Jacobi symbol, (a,n), expands the Legendre symbol to allow n to be composite, but now, if n is composite, we can't say that a result of 1 implies that a is necessarily a quadratic residue of n. Also, if and only if Legendre(a,p) == 0, then p divides a. If and only if Jacobi(a,n) == 0, then gcd(a,n) > 1. > >> > D. It appears that all non-squares have some smallest number such >> that >> > Jacobi(b,n) = 1. I don't know if this is proven or not. >> >> Do you mean = -1 ? > Yes, I should have said -1. It should have read: "It appears that all non-squares have some number b < n such that Jacobi(b,n) = -1." >> > >> > G. Another open question. "Does a small subset of numbers exist >> that when tested can prove N is a square?" >> >> >> They would need to be prime powers. > >Not all numbers. Each prime power creates a roughly 50/50 >filter, and all prime powers under consideration are coprime >these filters can be shoved into the Chinese Remainder >Theorem, and will always produce a complete set of solutions. >The ratios decrease, but the sets of solutions grow without bound. > OK, you've lost me now. Are you saying that you must use prime powers or that, for efficiency reasons, you should use prime powers? I can see that if a given n = 7(mod 15), then it cannot be a square, so couldn't 15 be a square-denier? I see now...if 15 is a square-denier, then either 3 or 5 must be a square denier, also. So no composite number can be a *minimum* square denier. Is that right? If you're just checking for squareness, why not use 105, say, instead of 3 and 5 and 7 separately? Tom Hadley =============================================== Phil Carmody Message 21 of 32 Jan 18, 2002 ----------------------------------------------- On Fri, 18 January 2002, "Hadley, Thomas H (Tom), ALINF" wrote: > >Not all numbers. Each prime power creates a roughly 50/50 > >filter, and all prime powers under consideration are coprime > >these filters can be shoved into the Chinese Remainder > >Theorem, and will always produce a complete set of solutions. > >The ratios decrease, but the sets of solutions grow without bound. > > > > OK, you've lost me now. Are you saying that you must use prime > powers or that, for efficiency reasons, you should use prime powers? Don't be fooled, any composite based tests are not Jacobi symbols. They are the 'intersection' of two Legendres. If one fails, the composite fails. > I can see that if a given n = 7(mod 15), then it cannot be a square, > so couldn't 15 be a square-denier? I see now...if 15 is a > square-denier, then either 3 or 5 must be a square denier, > also. So no composite number can be a *minimum* square denier. > Is that right? I'll cut to the toe-waggling scene: yes. > If you're just checking for squareness, why not use 105, say, > instead of 3 and 5 and 7 separately? Due to coprimality of those 3 numbers (all being prime and different helps in that regard!) the CRT guarantees a unique 1:1 correspondence between a {3,5,7} problem and a {105} problem. They are the same problem. Each odd prime rejects (p-1)/2 of the residues, which is where my 50/50 approximations come from. The CRT equivalence is why the 'three' tests modulo {63,64,65} are so common, they're actually a {2^8,3^2,5,7,13} test, i.e. 5 tests, 2 of which (the ^8 and ^2) are far stronger tests than the single power would be so the rejection rate is far better than 1:32 that 5 approximate 50/50s would provide. Phil Don't be fooled, CRC Press are _not_ the good guys. They've taken Wolfram's money - _don't_ give them yours. http://mathworld.wolfram.com/erics_commentary.html Find the best deals on the web at AltaVista Shopping! http://www.shopping.altavista.com =============================================== jfribrg@aol.com Message 22 of 32 Jan 19, 2002 ----------------------------------------------- The only reference I saw addressing this question is in Cohen's "A Course in Computational Algebraic Number Theory", which states that a number n is a square if and only if Jacobi(n,p) = 1 for all primes p not dividing n. This is completely useless from a computational standpoint. If you knew the prime factors of n, you would have no trouble determining if n is a square. At the same time, the squareness algorithm given in the above referenced book recommends calculating legendre symbols for a few well chosen moduli, and finally performing a square root if necessary. [Non-text portions of this message have been removed] =============================================== Ronald Hallam Message 23 of 32 Jan 19, 2002 ----------------------------------------------- The square detcection algorithm is in Cohen, section 1.7.2, page 39-40. Ron -----Original Message----- From: jfribrg@... To: primenumbers@yahoogroups.com Date: 19 January 2002 17:00 Subject: Re: [PrimeNumbers] Re: Programming Question >The only reference I saw addressing this question is in Cohen's "A Course in >Computational Algebraic Number Theory", which states that a number n is a >square if and only if Jacobi(n,p) = 1 for all primes p not dividing n. This >is completely useless from a computational standpoint. If you knew the prime >factors of n, you would have no trouble determining if n is a square. At >the same time, the squareness algorithm given in the above referenced book >recommends calculating legendre symbols for a few well chosen moduli, and >finally performing a square root if necessary. > > >[Non-text portions of this message have been removed] > > > >Unsubscribe by an email to: primenumbers-unsubscribe@egroups.com >The Prime Pages : http://www.primepages.org > > > >Your use of Yahoo! Groups is subject to http://docs.yahoo.com/info/terms/ > > =============================================== Phil Carmody Message 24 of 32 Jan 20, 2002 ----------------------------------------------- On Sat, 19 January 2002, jfribrg@... wrote: > The only reference I saw addressing this question is in Cohen's "A Course in > Computational Algebraic Number Theory", which states that a number n is a > square if and only if Jacobi(n,p) = 1 for all primes p not dividing n. This > is completely useless from a computational standpoint. If you knew the prime > factors of n, you would have no trouble determining if n is a square. Absolutely. I put a smallish challenge to a guy on sci.math who came up with the same observation. Finding ones which are hard to factor was hard! I finally found one which had no factor below 4e6. > At > the same time, the squareness algorithm given in the above referenced book > recommends calculating legendre symbols for a few well chosen moduli, and > finally performing a square root if necessary. Not having the book (grrr, the seller of the second hand one on Amazon won't ship internationally :-( ), I'd guess that these moduli are the small prime squares, which have the best size/rejection features. Everyone uses 2^n, some use 3^2 (as part of 63), but everyone seems to forget about 5^2. On a side note, last night I threw my sieve up another notch (74.4G/s woo woo!), and this morning it found 7841727547787169 %229 =200 So I have finally found the number which passes all Jacobi tests for Prime powers under 200 and a cherry on top! Honestly, I'll stop this sieve eventually. I think that 10^16 is probably far enough. Maybe. :-) Phil Don't be fooled, CRC Press are _not_ the good guys. They've taken Wolfram's money - _don't_ give them yours. http://mathworld.wolfram.com/erics_commentary.html Find the best deals on the web at AltaVista Shopping! http://www.shopping.altavista.com =============================================== thefatphil Message 25 of 32 Jan 20, 2002 ----------------------------------------------- --- In primenumbers@y..., Phil Carmody wrote: > So I have finally found the number which passes all Jacobi tests for Prime powers under 200 and a cherry on top! Someone beat me to the cherry, alas: http://www.ams.org/journal-getitem?pii=S0025-5718-96-00678-3 Oh dear - GIFs not cut'n'pasted... <<< Abstract: If is an odd prime, the pseudosquare is defined to be the least positive nonsquare integer such that and the Legendre symbol for all odd primes . In this paper we first discuss the connection between pseudosquares and primality testing. We then describe a new numerical sieving device which was used to extend the table of known pseudosquares up to . We also present several numerical results concerning the growth rate of the pseudosquares, results which so far confirm that , an inequality that must hold under the extended Riemann Hypothesis. >>> The crazy thing about these guys is that they decided to not use prime powers, only primes, but then throw in a p==1 mod 8. Now when I did sums at university, 8 was considered to be a prime power, not a prime. However, they appear to have had a 500G/s sieve... Ouch! Almost certainly custom hardware though, and factors of a thousand are possible when you've got custom hardware. (c.f. the DES cracking chip, for instance.) Phil =============================================== Paul Jobling Message 26 of 32 Jan 21, 2002 ----------------------------------------------- > Not having the book (grrr, the seller of the second hand one > on Amazon won't ship internationally :-( ), I'd guess that > these moduli are the small prime squares, which have the best > size/rejection features. Everyone uses 2^n, some use 3^2 (as > part of 63), but everyone seems to forget about 5^2. The values that Cohen uses are 64, 63, 65, and 11. 99% of non-squares fail ones of these tests. After this, an integer square root must be calculated. Cohen suggests that the tests are performed in the order above. Regards, Paul. __________________________________________________ Virus checked by MessageLabs Virus Control Centre. =============================================== Phil Carmody Message 27 of 32 Jan 21, 2002 ----------------------------------------------- On Mon, 21 January 2002, "Paul Jobling" wrote: > > Not having the book (grrr, the seller of the second hand one > > on Amazon won't ship internationally :-( ), I'd guess that > > these moduli are the small prime squares, which have the best > > size/rejection features. Everyone uses 2^n, some use 3^2 (as > > part of 63), but everyone seems to forget about 5^2. > > The values that Cohen uses are 64, 63, 65, and 11. 99% of non-squares fail > ones of these tests. After this, an integer square root must be calculated. > Cohen suggests that the tests are performed in the order above. Does he say why? It's obvious really: 64 is 'instant' But why not do 64 bits simultaniously, rather than just 6? - see my 4-clock-tick no-lookup-table pseudocoded algorithm some time last year. Several % better than using just base 64. Evaluate residues modulo 63 and 65 _simultaniously_, it's no slower than doing just one reduction for a numbers of more than a few words in length. (If you can't see how, try doing mod 9 and mod 11 simultaniously in base 10) In handwaving terms, 63 contains 3 rejection criteria (3, 9, and 7), but 65 has only 2 (5 and 13), so use the most effective one first. However, exact ratios should be used, and the different rejection ratios are not the same. When worked out exactly, 63 does indeed reject more than 65. 11? pah, if you're going to throw an extra single prime in then you should be made to justify it. For _large_ numbers, you want many small prime tests. For small numbers the square root isn't hugely expensive. Therefore the number of extra primes you use after 63,64,65 should _depend on the size of the number being tested_. To arbitrarily chose to include 1 is a bit, erm, arbitrary. My guess is that he needed it for the mystical magical decadigital 99% rejection ratio, providing some kind of 10-fingered feel-good factor. But why 11???? Both 17 and 31 can be done faster on most typical architectures, and both of them have a better rejection ratio. Phil It is against US Department of Agriculture regulations to advertise or sell as "Prime Rib" any cut of meat containing a non-prime number of ribs. Find the best deals on the web at AltaVista Shopping! http://www.shopping.altavista.com =============================================== Paul Jobling Message 28 of 32 Jan 21, 2002 ----------------------------------------------- > > The values that Cohen uses are 64, 63, 65, and 11. 99% of > non-squares fail > > ones of these tests. After this, an integer square root > must be calculated. > > Cohen suggests that the tests are performed in the order above. > > Does he say why? Because the rejection ratios are greatest in that order, as you say. Cohen suggests that it might be better to do it all in one go (save for the mod 64, obviously) , and as 6K is nothing these days then working mod 45760 would be just a long division and a lookup into a 5720 byte bitmap. And IMHO I do feel that 11 was included just to get it below 1% - without it it is about 1.2% (AFAIR). Regards, Paul. __________________________________________________ Virus checked by MessageLabs Virus Control Centre. =============================================== Phil Carmody Message 29 of 32 Jan 21, 2002 ----------------------------------------------- On Mon, 21 January 2002, Marcel Martin wrote: > [*] As a matter of fact, we can do considerably better than a division > of a big integer N by a small integer p. > > Since p is odd, we can use an operation that is much faster than the > division. > 1) Compute the inverse U of p modulo B (ordinary, B = 2^32). This can > be done in, say, 50-60 CPU cycles using a lookup table of inverses > mod 2^8. > 2) 'Divide' N using p and its inverse U. This gives us a carry C such > that 0 <= C < p. > > At this point, we have N = Q*p - C*B^i (i is equal to Size(N) or > Size(N)-1). BTW, if C = 0, we know that p divides N. > > Now, Jacobi(N, p) = Jacobi(Q*p - C*B^i, p) = Jacobi(-C*B^i, p) > but > Jacobi(-C*B^i, p) = Jacobi(-C, p) * Jacobi(B^i, p) > and > Jacobi(B^i, p) = 1 (since B is a square and gcd(B,p) = 1). > Finally, it comes > Jacobi(N, p) = Jacobi(-C, p) > > Dividing the bigint N by p requires Size(N) single divisions. > 'Dividing' it using U and p requires 2*Size(N) single multiplications. > Practically, this is 3.5 times faster on AMD processors. > > I put 'dividing' between '' because this is not really a division. > The 'quotient' is equal to the one obtained with the standard division > if and only if N is a multiple of p. Yeah, but wouldn't you rather 'divide' it by 16777215 = 0xFFFFFF? Take 3 words (p,q,r) from memory and partition the bytes as follows into 4 24-bit quantities: ppppqqqqrrrr wwwxxxyyyzzz Now accumulate 4 32-bit sums of www xxx yyy and zzz Wrap carry bytes every 256*3 words (=3K bytes) This 'division' takes about 5 shifts and a few masks per 3 words, plus a a few operations at the end to merge the 4 sums. Et voila, residues mod 63, 65 and 17 (and 241 if that takes your fancy). Phil Don't be fooled, CRC Press are _not_ the good guys. They've taken Wolfram's money - _don't_ give them yours. http://mathworld.wolfram.com/erics_commentary.html Find the best deals on the web at AltaVista Shopping! http://www.shopping.altavista.com =============================================== Phil Carmody Message 30 of 32 Jan 21, 2002 ----------------------------------------------- On Mon, 21 January 2002, Hagen Reddmann wrote: > Or You can use this with 2^96-1. We partitionate at 3 x 32Bit Words and can > use it then with all products of > 3,3,5,7,13,17,97,193,241,257,673,65537,22253377 > This give us 1222 possible modulis. A 96Bit Digit have the advantage that we > use ONLY additions, no shifts or masking. I saw someone else (I think it was the GMP guys) recommend the idea for using 4 separate registers rahter than one massive 96-bit value requiring carries between the three words, and they claimed to have a lightning fast x86 machine code inner loop in mind. They're the GMP programmers, so when they say it's fast, I trust them! In truth they probably run at a similar speed, as it's probably memory-bound. It is also possible to extract those moduli from my scheme too if you treat the carry bytes differently. However, on the whole yours is simpler. Note however, that even with all those moduli, it's still not necessarily the best, as it might be possible to do two smaller values in parallel. I think that x=30 and x=48 or something like that make a good pair. (but shifts are required!) > If You interesst I have wrote some code to produce all modulis to any X of > the form 2^x -1 and all corresponending primefactors. I took it up to 64 in the past. Please feel free to post your results up to 96 to the list. It will hopefully jog my memory what the magic pair of powers was! Cheers for your input, Phil Don't be fooled, CRC Press are _not_ the good guys. They've taken Wolfram's money - _don't_ give them yours. http://mathworld.wolfram.com/erics_commentary.html Find the best deals on the web at AltaVista Shopping! http://www.shopping.altavista.com =============================================== Phil Carmody Message 31 of 32 Jan 21, 2002 ----------------------------------------------- On Mon, 21 January 2002, Marcel Martin wrote: > Checking with your method takes T1 (running time). > Checking with a mod operator (not necessarily Cohen's one) takes T2. > Computing a square root takes T3. > We always have T1 < T2. In fact, since the complexity is the same, > we can compute a constant c = T2/T1 > > 1) With your method > > Running time = T1 + 0.0075 T3 > > 2) With a modulo operator of which the proba is p (with the one I use > in Primo, p = 0.00028. No typo, there are 3 zeros) > > Running time = T2 + p T3 = c T1 + p T3 > > Your method is better if and only if (c - 1) T1 > (0.0075 - p) T3. > This is true for a given size of numbers. Obviously T3 increases much > more quickly than T1 when the number size increases. As I said in http://groups.yahoo.com/group/primenumbers/message/4855 when I first lambasted Henri's arbitrary choice of 11, <<< Therefore the number of extra primes you use after 63,64,65 should _depend on the size of the number being tested_. >>> There's no one correct solution. > Clearly, your test is always better than Cohen's one (0.0075 < p). > But, a priori, it is not sure whether it is better than the one I use > in Primo. In that case p is very small and c is close to 1. I will > try to see... :-) There will be better ones, certainly, 2^24-1 was only the first number out of the hat. You must also remember that my programming environment of choice is a Dec Alpha, where there isn't even a divide instruction, and the multiply is not pipelined (on the integer side, that is). So I always favour bit-twiddles over multiply and divide. (But no, I wouldn't use the 32-bit oriented scheme on the alpha, obviously.) With pipelined division via multiplication, one can probably extract two finely-crafted sets of primes. I'd recommed sticking 5^2 in the mix certainly, something that only a few 2^n+/-1 give you (2^10+1 => 2^(10n-1) too). 7^2 is better than any single prime, and you need 2^21-1 to get that. So I'd guess that (3.5.7.11.13)^2.17 would be a good start, but that's a guess. Whether (11.13)^2 is better than 3 other primes I can't tell you off the top of my head. Maybe bumping up the power of 3 to 3^4 at the expense of 11 or 13 improves things, too? And this is all working on the assumption that you're doing a 32-bit modulus. If I know you - you've probably been able to crank it up to 52 or even 64 bits at no loss of speed by judicious use of the FPU. The optimal value for each size could in theory be calculated, I'm sure, just by brute force (no more than about a hundred thousand combinations for the 32-bit problem), but alas I'm not in a programming mood at the moment (come on, it's 8am, who's in a programming mood at 8am?). Anyone want to brute force this for 31, 32, 51, 52, 63, 64 bits? (working on the assumtion that division via multiplication sacrifices 1 bit of precision). Phil Don't be fooled, CRC Press are _not_ the good guys. They've taken Wolfram's money - _don't_ give them yours. http://mathworld.wolfram.com/erics_commentary.html Find the best deals on the web at AltaVista Shopping! http://www.shopping.altavista.com =============================================== Paul Jobling Message 32 of 32 Jan 22, 2002 ----------------------------------------------- Marcel Martin asked: > >Cohen suggests that it might be better to do it all in one go (save > >for the mod 64, obviously), and as 6K is nothing these days then > >working mod 45760 would be just a long division and a lookup into a > >5720 byte bitmap. > > Why 45760? This was a typo I guess. 45045 was what I meant. Regards, Paul. __________________________________________________ Virus checked by MessageLabs Virus Control Centre. =============================================== Cached by Georg Fischer at Nov 14 2019 12:47 with clean_yahoo.pl V1.4