{"action":"created","issue":{"url":"https://api.github.com/repos/JuliaLang/julia/issues/11594","labels_url":"https://api.github.com/repos/JuliaLang/julia/issues/11594/labels{/name}","comments_url":"https://api.github.com/repos/JuliaLang/julia/issues/11594/comments","events_url":"https://api.github.com/repos/JuliaLang/julia/issues/11594/events","html_url":"https://github.com/JuliaLang/julia/issues/11594","id":85559992,"number":11594,"title":"improve prime sieve performance","user":{"login":"StefanKarpinski","id":153596,"avatar_url":"https://avatars.githubusercontent.com/u/153596?v=3","gravatar_id":"","url":"https://api.github.com/users/StefanKarpinski","html_url":"https://github.com/StefanKarpinski","followers_url":"https://api.github.com/users/StefanKarpinski/followers","following_url":"https://api.github.com/users/StefanKarpinski/following{/other_user}","gists_url":"https://api.github.com/users/StefanKarpinski/gists{/gist_id}","starred_url":"https://api.github.com/users/StefanKarpinski/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/StefanKarpinski/subscriptions","organizations_url":"https://api.github.com/users/StefanKarpinski/orgs","repos_url":"https://api.github.com/users/StefanKarpinski/repos","events_url":"https://api.github.com/users/StefanKarpinski/events{/privacy}","received_events_url":"https://api.github.com/users/StefanKarpinski/received_events","type":"User","site_admin":false},"labels":[{"url":"https://api.github.com/repos/JuliaLang/julia/labels/intro%20issue","name":"intro issue","color":"0052cc"},{"url":"https://api.github.com/repos/JuliaLang/julia/labels/up%20for%20grabs","name":"up for grabs","color":"0b02e1"}],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":22,"created_at":"2015-06-05T15:32:49Z","updated_at":"2015-06-13T09:14:28Z","closed_at":null,"body":"Someone emailed me directly with substantial performance improvements to our [prime sieve implementation](https://github.com/JuliaLang/julia/blob/7f55b7548d7936fe74f65a7569557c0ec998cd94/base/primes.jl#L9-L28). Unfortunately, they didn't give permission to use the code under the MIT license and haven't replied to me, so I'm going to post the non-code parts of their email so that someone else who hasn't seen their code can do an independent implementation of these improvements.\r\n\r\n> * Instead of the for loop running to `n`, we only need to go up to `isqrt(n)`\r\n> * Further, we can loop only over odd values of `p`, and set only odd multiple of `p^2` to false.\r\n> * Not finished yet – we can vectorise the last loop, with a further gain\r\n> * Finally, we can use the facts that\r\n>   * all large primes are of the form 6k±1\r\n>   * most integers are composite, so it is better to set all flags to false initially, then invert for possible primes\r\n\r\nIt would have been nice to just have a pull request fixing this (with an implied MIT license) since this individual went to the trouble of making and benchmarking these improvements, but instead here we are."},"comment":{"url":"https://api.github.com/repos/JuliaLang/julia/issues/comments/111690172","html_url":"https://github.com/JuliaLang/julia/issues/11594#issuecomment-111690172","issue_url":"https://api.github.com/repos/JuliaLang/julia/issues/11594","id":111690172,"user":{"login":"VicDrastik","id":12775830,"avatar_url":"https://avatars.githubusercontent.com/u/12775830?v=3","gravatar_id":"","url":"https://api.github.com/users/VicDrastik","html_url":"https://github.com/VicDrastik","followers_url":"https://api.github.com/users/VicDrastik/followers","following_url":"https://api.github.com/users/VicDrastik/following{/other_user}","gists_url":"https://api.github.com/users/VicDrastik/gists{/gist_id}","starred_url":"https://api.github.com/users/VicDrastik/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/VicDrastik/subscriptions","organizations_url":"https://api.github.com/users/VicDrastik/orgs","repos_url":"https://api.github.com/users/VicDrastik/repos","events_url":"https://api.github.com/users/VicDrastik/events{/privacy}","received_events_url":"https://api.github.com/users/VicDrastik/received_events","type":"User","site_admin":false},"created_at":"2015-06-13T09:14:28Z","updated_at":"2015-06-13T09:14:28Z","body":"\r\nLet's not confuse algorithms for 32-bit integers with those for 64-bit integers: Khashin's test of course works for 32-bit integers, but it is overkill to use it, and the same goes for BPSW. We can get away with using only M-R for small n i.e. both Khashin and BPSW should only be used for large n. There are 2314 MR(2) 32-bit spsp's, and any other (fixed) MR test reduces this to between 70 and 110 spsp's, normally needing yet another test to clean up (e.g. bases 2,7, and 61). My quadratic residue idea produces only 42 spsp's, but 32 of these can be quickly eliminated, and the remaining 10 can be zapped either using a table or simply 10 (fixed) trial divisions.\r\n\r\nDeinition 1.1 of Khashin is the correct one (the others are just used in the proof, not in the actual algorithm), so \r\n \r\n let c be an odd positive number (does not need to be a prime, but it always will be)\r\n let d = sqrt(c)\r\n let j = jacobisymbol(c,n)\r\n\r\n Then, if n is prime, (1+d)^n === (1+j*d) (mod n)\r\n\r\n Khashin then proves that if we choose the smallest odd c such that j==-1, then there are no Frobenius pseudoprimes less than 2^60 (assuming gcd(2.3.5.7.11.13.17,n)=1). It is also required that c<128 but this turns out to be only of theoretical interest.\r\n\r\n I have to concede your point that using the word \"provably\" is too large a claim. Also, that he does not explicitly claim this unconditional result. When I read Khashin's paper, I just automatically assumed he would always require a MR(2) test first and I immediately tested whether c<128 for base 2 strong pseudoprimes. Here is my code and the results:\r\n\r\nfunction _()\r\nf=open(\"2-SPRP-2-to-64.txt\",\"r\")\r\ng=open(\"output.txt\",\"w\")\r\ncount=zeros(Int,999)\r\nfor line in eachline(f)\r\n n=uint64(line)\r\n if n==3511^2 continue end\r\n if n==1093^2 continue end\r\n for c=3:2:999\r\n  if jacobisymbol(c,n)<1\r\n   count[c]+=1  \r\n   break\r\n  end\r\n end\r\nend\r\nfor i=1:999\r\n if count[i]>0 println(g,i,\" , \",count[i]) end\r\nend\r\nprintln(g,sum(count)+2)\r\nclose(g)\r\nclose(f)\r\nend\r\n\r\n  3 , 11515748\r\n  5 ,  7836069\r\n  7 ,  5608159\r\n 11 ,  3275113\r\n 13 ,  1727914\r\n 17 ,   926810\r\n 19 ,   487767\r\n 23 ,   254188\r\n 29 ,   130373\r\n 31 ,    65303\r\n 37 ,    33069\r\n 41 ,    16793\r\n 43 ,     8289\r\n 47 ,     4218\r\n 53 ,     2097\r\n 59 ,     1016\r\n 61 ,      548\r\n 67 ,      276\r\n 71 ,      138\r\n 73 ,       59\r\n 79 ,       32\r\n 83 ,       18\r\n 89 ,        7\r\n 97 ,        4\r\n101 ,        2\r\n103 ,        1\r\n107 ,        1\r\n\r\n31894014\r\n\r\n so, the largest c value needed is c=107. So, yes, this is not a \"proof\", because it depends on computed results rather than algebra. But, assuming Feitsma's table and my code are correct, this covers all pseudoprimes up to 2^60. Naturally, it does not cover composites up to 2^60 which are not MR(2) pseudoprimes, but it doesn't have to if you do a MR(2) test first.\r\n\r\n Several years ago, before Khashin's paper, I had a lengthy exchange of emails with Paul Underwood about his Frobenius-style test. Writing them both in a common notation, we can see the similarities\r\n\r\nUnderwood-Frobenius:\r\n for any integer p such that jacobisymbol(p^2-4,n)==-1,\r\n   (2+x)^(n+1) === 5+2*p mod(x^2-p*x+1),(mod n)\r\n\r\nKhashin-Frobenius:\r\n  for the smallest odd integer c such that jacobisymbol(c,n)==-1,\r\n   (1+x)^(n+1) === 1-c   mod(x^2-c),(mod n)\r\n \r\n  Underwood's test is not faster (if anything, it has to be slightly slower). I suggested to Paul that he should run his test on Feitsma's Fermat(2) database, and he found no counterexamples, so we should not be surprised when Khashin's version also passes.(when I have time, I will run Khashin's test on Feitsma's full data set, not just the strong pseudoprimes)\r\n\r\n Grantham's seminal paper was great, but several authors(including Khashin) have since suggested that his 1/7710 bound was far too pessimistic.\r\n\r\n The only modern test that can be compared to Khashin's is Loebenberger's version of the Crandall/Pomerance Frobenius test. Here is that algorithm in Python:\r\n\r\n\r\ndef frobenius(n):\r\n    \"\"\" returns True if n is probably prime\r\n    if n % 2 == 0 or n < 2 then   return False\r\n    #step 1\r\n    while 1:\r\n        a = randint(1, n-1)\r\n        b = randint(1, n-1)\r\n        d = a ** 2 - 4 * b\r\n        if  not (is_square(d) or gcd(2 * a  * b * d, n) != 1) then break\r\n    #step 2\r\n    m = (n - jacobi(d, n)) // 2\r\n    w = (a ** 2  * inverse(b, n) - 2) % n\r\n    u = 2\r\n    v = w \r\n    # step 3\r\n    for bit in bits_of(m):\r\n        if bit:\r\n            u = (u * v - w) % n \r\n            v = (v * v - 2) % n\r\n        else:\r\n            v = (u * v - w) % n\r\n            u = (u * u - 2) % n\r\n    # step 4\r\n    if (w * u - 2 * v) % n != 0 then return False\r\n    return (pow(b, (n - 1) // 2, n) * u) % n ==  2\r\n\r\nThe Lucas part of this takes only 2 selfridges, rather than 3 for Khashin/Underwood, but the parameters are random, and so finding pseudoprimes is not possible. \r\n"}}