-3

Complex plot of the Analytic version

I've been doing some number theory during all this Covid downtime and I think I've discovered an pretty interesting (if not novel) algorithm for detecting primality. I posted my writeup on my LinkedIn page, you don't have to signup or anything.

https://www.linkedin.com/pulse/efficient-prime-number-generation-christopher-wolfe/

I just want to know if this technique is already known, because it's quite fast (constant time) and pretty compact in size. I did a deeper dive a while back which you can read on my blog.

http://jasuto.com/ideas/primes/

It would be great to get some feedback and to verify that this is new. I have quite a bit more to release, but thought I would start light :)

Here is the demonstration code if you don't feel like reading the article. I have seen many prime number implementations, but nothing O(1) or this small...

# Copyright 2021 Christopher Wolfe (chris@jasuto.com)
# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.


# constant time primality testing, along with a fairly concise prime
# number generator.  


# bignum division handles flooring
def floor(n):
    return n


def alg_prime_q(ni):
    n = ni - 1
    return floor((n*2**(n + 1) + 2)//(2*n - (-1)**n + 3)) - floor((n*2**(n + 1) - 2)//(2*n - (-1)**n + 3))


# cos(π x) + i sin(π x)
def phasor(n):
    return ((n + 1) & 1) << 1 - 1


# this will include Fermat pseudoprimes as well
def prime_q(n):
    q = n - 1
    s = q * (1 << (q + 1))
    t = (q << 1) - phasor(q) + 3
    return floor((s + 2) // t) - floor((s - 2) // t)


# returns all primes < n
def primes(n):
    return [i for i in range(3, n, 2) if prime_q(i)]

res = primes(1000)
print(res)


print('done')
Chris
  • 1
  • 2
  • 6
    Please explain the algorithm in your question, not in a link. – kaya3 Feb 07 '21 at 16:58
  • i wrote an article about it, first link. i believe this is a new algorithm, just trying to confirm before i release some variants. – Chris Feb 07 '21 at 17:02
  • 2
    You need to include the explanation in the question, not in a link, as I said. Also this algorithm absolutely is not O(1) time. – kaya3 Feb 07 '21 at 17:04
  • to check if a number is prime, it definitely is. – Chris Feb 07 '21 at 17:06
  • the prime generation is just for reference, a constant time primality test would be novel i believe. – Chris Feb 07 '21 at 17:07
  • 2
    No, it is not. How do you expect a computer to shift an arbitrary number of bits in constant time? – kaya3 Feb 07 '21 at 17:07
  • This isn't really a normal post but here's a wikipedia article on prime algorithms if you can't find it there its probably a new way: https://en.wikipedia.org/wiki/Primality_test – Yash Feb 07 '21 at 17:09
  • Your comment in the code says it "will include Fernat pseudoprimes as well" - what exactly do you mean by "include", does it say they are prime, or does it correctly say they are not prime? – kaya3 Feb 07 '21 at 17:09
  • Unfortunately, it does not work, your generator gives non prime numbers in the result : 341, 561, 645 ! Could you explain how you did it ? Moreover you have a for loop, so I do not think it is O(1), it is rather O(n/2) I think. – Malo Feb 07 '21 at 17:13
  • you can easily remove the pseudoprimes, it just adds to the complexity, and i was mostly curious if this was a known technique. the 'errors' are actually even a subclass of the fermat primes, but it does require an additional divide... – Chris Feb 07 '21 at 17:21
  • @Malo The OP is claiming that the `prime_q` function is constant time, not the `primes` function; it is not, because arithmetic and bitwise operations on arbitrarily large integers are not constant time. – kaya3 Feb 07 '21 at 17:21
  • @Kaya3 thanks for explaining, but title and last sentence are misleading : "I have seen many prime number implementations, but nothing O(1) or this small..." – Malo Feb 07 '21 at 17:23
  • "floor" does nothing, why is it included? – Christian Sloper Feb 07 '21 at 17:24
  • 1
    *"you can easily remove the pseudoprimes"* - that's not how this works, if what you have is an algorithm that *you think* can easily be made to work, then what you actually have (until you actually make it work) is an algorithm that doesn't work. I'm voting to close this question because it is based on multiple false premises; you don't have a primality testing algorithm, and it isn't constant time. – kaya3 Feb 07 '21 at 17:24
  • @Kaya3 so we agree – Malo Feb 07 '21 at 17:25
  • 1
    I think is a variation of https://oeis.org/wiki/Chinese_hypothesis – Christian Sloper Feb 07 '21 at 17:26
  • but i agree with Kaya and Malo that it would be intersesting to see how you easily remove the pseudoprimes. – Christian Sloper Feb 07 '21 at 17:27
  • lol, didn't mean to be a bother, i just thought it might be interesting for others to see. this is just a tiny piece of some ideas i've been working on for a while. and yes, you can generate primes in a novel way with this but it needs a proper writeup. – Chris Feb 07 '21 at 17:28
  • 2
    Well, it really doesn't work. I tested it for the primes below 200000 (it gets very slow, more than linearily for large values of n). It took about 20 seconds, while the functions in https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n take a fraction of a second. Your code gives 18089 primes under this value, while there really are only 17984 - that's 95 too much... – Thierry Lathuille Feb 07 '21 at 17:28
  • 3
    As for *"I have quite a bit more to release, but thought I would start light"*, I suggest you approach your other discoveries from the starting point of "is this correct?" rather than "is this novel?". If, like in this case, you *already know* it is not correct but you think it can be "easily" made correct, then do that so-called "easy" work before asking others to tell you so. – kaya3 Feb 07 '21 at 17:29
  • It involves calculating (n-1)*2**n , so no wonder it starts to take some time – Christian Sloper Feb 07 '21 at 17:30
  • the floor is there to stop someone from algebraically simplifying the divisors. it is done within pythons division code. it might not be the case for all langs though. – Chris Feb 07 '21 at 17:30
  • @Chris, no problem, do not hesitate to show us the final version of your novel work. – Malo Feb 07 '21 at 17:30
  • Would be happy to see an extended version. – Christian Sloper Feb 07 '21 at 17:30
  • maybe a version without all the bitshifting so its easier to follow. q << 1 could be replaced by 2*q , ((n + 1) & 1) << 1 - 1 is better with "if odd then ... else .. " – Christian Sloper Feb 07 '21 at 17:33
  • there are different analytical forms, you can read a little about them on my article, this is not going to be the fastest version, especially using python's bignums... – Chris Feb 07 '21 at 17:33
  • i added the non-bit version to the code... – Chris Feb 07 '21 at 17:38
  • @ChristianSloper it's not related to the chinese-hypothesis, but the addition pseudo-primes are exactly the Sarrus numbers that are in that post. – Chris Feb 07 '21 at 17:42

1 Answers1

5

I think I can show that this is a version of the chinese hypothesis i linked in comments.

Let us simplify the version without the bit shifts.

def alg_prime_q(ni):
    n = ni - 1
    return floor((n*2**(n + 1) + 2)//(2*n - (-1)**n + 3)) - floor((n*2**(n + 1) - 2)//(2*n - (-1)**n + 3))

I am going to assume ni is odd, as this simplify the expression a lot and is also the only interesting case (determining primality of even numbers is easy).

In that case, (-1)**n = 1

floor((n*2**(n + 1) + 2)//(2*n - 1 + 3)) - floor((n*2**(n + 1) - 2)//(2*n - 1+ 3))

also 2 is a factor in both numerator and denominator so this reduces further to:

floor((n*2**n  + 1)//(n + 1)) - floor((n*2**(n ) - 1)//(n + 1))

replacing back ni for n-1 :

floor((ni-1)*2**(ni-1)  + 1)//ni - floor((ni-1)*2**(ni-1) - 1)//ni

The claim is that ni is prime if ni divides ((ni-1)*2**(ni-1) + 1) more times than ni divides ((ni-1)*2**(ni-1) - 1).

Observe that the difference of the two is 2, so for this to happen (ni-1)*2**(ni-1) + 1)%ni has to be either 1 or 0.

Let us check the two separately: Case 1,

(ni-1)*2**(ni-1) + 1 = 1 (mod ni)
ni*2**(ni-1)-*2**(ni-1) = 0 (mod ni)
2**(ni-1) = 0 (mod ni)

Since n is odd, it can never be a power of 2 so this case never happens

Case 2.

(ni-1)*2**(ni-1) + 1 = 0 (mod ni)
ni*2**(ni-1)-2**(ni-1) + 1 = 0 (mod ni)
2**(ni-1) + 1 = 0 (mod ni)
2**ni + 2 = 0 (mod ni)

Which is the chinese hypothesis.

Christian Sloper
  • 6,238
  • 2
  • 11
  • 28
  • i think you just essentially proved that the sequence contains the Poulet set, which we already knew. but it also contains all the primes as well. which the chinese conjecture does not. also i already implemented the version without the pseudoprimes using a slightly different analytical form that i go over in the article, running as a Cuda kernel no less. Let me add a picture of the plot, you can clearly see it's not a simple mersenne-ish function. it's actually strikingly close to the zeta function in certain areas. – Chris Feb 07 '21 at 20:22
  • 1
    All primes satisfy the chinese conjecture (by Fermats little theorem) – Christian Sloper Feb 07 '21 at 20:23
  • 1
    Your formula is just the chinese conjecture by extra steps. – Christian Sloper Feb 07 '21 at 20:24
  • i cant embed pictures yet, lol, so there is another link... – Chris Feb 07 '21 at 20:27
  • 1
    I think i have shown quite convincingly that your formula reduces (at least for odd numbers) directly to the chinese conjecture. It contains exactly the same numbers as the chinese conjecture (primes and Poulet set) as it is, in the end, what you are doing. – Christian Sloper Feb 07 '21 at 20:39
  • You CANNOT factor the the floors out, that is why i left them in. – Chris Feb 07 '21 at 20:41
  • try without using the integer division, which handles flooring internally like i said. – Chris Feb 07 '21 at 20:41
  • 1
    I did not factor the floors out.. – Christian Sloper Feb 07 '21 at 20:42
  • import math; def floor(n): return n # return math.floor(n) def prime_q(ni): n = ni - 1 return floor((n*2**(n + 1) + 2)/(2*n - (-1)**n + 3)) - floor((n*2**(n + 1) - 2)/(2*n - (-1)**n + 3)) for i in range(50): print(prime_q(i)) – Chris Feb 07 '21 at 20:42
  • 1
    not sure what you are on about, the only factoring i did was factoring a two from both denominator and numerator. which is perfectly fine – Christian Sloper Feb 07 '21 at 20:43
  • i see that use // divison for flooring at a point, added back floor function for clarity. changes notthign – Christian Sloper Feb 07 '21 at 20:45
  • without the floors: 2.0 2.0 0.6666666666666667 0.6666666666666665 0.40000000000000036 0.40000000000000036 0.2857142857142847 0.2857142857142847 0.22222222222222854 0.22222222222220012 0.18181818181818699 0.18181818181813014 0.15384615384618883 0.15384615384618883 0.13333333333321207 0.13333333333321207 0.1176470588252414 0.1176470588252414 0.10526315789320506 0.10526315789320506 0.09523809526581317 0.09523809526581317 0.08695652172900736 0.08695652196183801 0.0800000000745058 0.0800000000745058 0.07407407462596893 0.074074067175388 .... – Chris Feb 07 '21 at 20:47
  • 1
    I havent removed the floors man – Christian Sloper Feb 07 '21 at 20:48
  • and with: 2 2 1 1 1 1 1 1 0 0 1 1 1 1 0 0 1 1 1 1 0 0 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 0 0 1 1 1 1 0 0 1 1 0 0 – Chris Feb 07 '21 at 20:48
  • those fractions are not as clean as you are expecting – Chris Feb 07 '21 at 20:49
  • 1
    I havent removed any floors..... – Christian Sloper Feb 07 '21 at 20:50
  • remove the // and replace with / and set the count to 30 or something, overwise it will overflow – Chris Feb 07 '21 at 20:50
  • anyhow, like i said i already have this working, it simply doesnt reduce down to what you calculated – Chris Feb 07 '21 at 20:51
  • yes, it does..... show where in my reduction i am wrong. – Christian Sloper Feb 07 '21 at 20:51
  • Ok, whatever, I tried to help you. Let me know when you get this published. – Christian Sloper Feb 07 '21 at 20:52
  • yea, but he has this "n = ni -1 " , so when ni is odd, then n is even" @kaya3 – Christian Sloper Feb 08 '21 at 07:03
  • Ah, my mistake. – kaya3 Feb 08 '21 at 07:04