RE: [math-fun] Indoor sports: arctan relations
From: "Schroeppel, Richard" <rschroe@sandia.gov>
Speculations on finding "smooth" X^2+1:
[What was wrong with "round"? That's what Hardy & Wright use.]
That's sooooo last century!
You can generate ordinary integers in order, along with their prime factorizations, using a priority-queue-like scheme.
"Lenstra's Algorithm". By default doesn't give you the factorisations, but that part is trivial post facto.
A similar method should work for X+i, using complex primes. Unclear if it's any faster than just sieving large blocks of X+i.
The set { x^2+1 } is so sparse that sieving just along that polynomial will be fairly fast. O(R^(1/2)).
Both methods are faster than factoring N (or X^2+1) one-at-a-time.
But using Euler's factoring shortcut (the original block-GCD method), factoring is exceptionally fast. One needs to add extra prime powers, and could still miss extreme cases, but we're not after exhaustive searching are we, merely computationally efficient searching. My 1-line GP/Pari script finds all P-smooth X^2+1 numbers up to (10^7)^2 in only seconds using this not-completely-naive factoring technique. I only wen't to p=100, and found the following 10^6<X<10^7 1984933 2343692 3449051 6225244
Both methods adapt to ignoring large prime divisors, just keeping those up to some limit.
Lenstra's method is defined in those terms. Alas it's something like O(f(R)*#p) for f(R)=o(R) by default - the larger your set, the more collisions you have, and the more time you spend hopping down lists trying to find the right place to insert values. So it tends to slow down for large sets of divisors, and may take a while before the smoothness-contraint causes the unspecified o(R) to become preferable to the square+1 constraint's O(R^(1/2)). I have a 'Lenstra Engine' somewhere - if someone wants me to perform some actual hunts for smooth numbers with it, just let me know the details, and I'll get a spare machine onto it immmediately. Phil () ASCII ribbon campaign () Hopeless ribbon campaign /\ against HTML mail /\ against gratuitous bloodshed [stolen with permission from Daniel B. Cristofani] __________________________________________________ Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com
* Phil Carmody <thefatphil@yahoo.co.uk> [Apr 12. 2006 12:22]:
From: "Schroeppel, Richard" <rschroe@sandia.gov>
Speculations on finding "smooth" X^2+1:
[What was wrong with "round"? That's what Hardy & Wright use.]
That's sooooo last century!
8-))
You can generate ordinary integers in order, along with their prime factorizations, using a priority-queue-like scheme.
"Lenstra's Algorithm". By default doesn't give you the factorisations, but that part is trivial post facto.
Dijkstra's algorithm (your follow-up). Clueless me: <cite http://en.wikipedia.org/wiki/Dijkstra's_algorithm> Dijkstra's algorithm, named after its discoverer, Dutch computer scientist Edsger Dijkstra, is an algorithm that solves the single-source shortest path problem for a directed graph with nonnegative edge weights. </cite> I do not see the connection, or is there other "Dijkstra's algorithm"s?
A similar method should work for X+i, using complex primes. Unclear if it's any faster than just sieving large blocks of X+i.
The set { x^2+1 } is so sparse that sieving just along that polynomial will be fairly fast. O(R^(1/2)).
Both methods are faster than factoring N (or X^2+1) one-at-a-time.
But using Euler's factoring shortcut (the original block-GCD method), factoring is exceptionally fast.
Can you give me a reference?
One needs to add extra prime powers, and could still miss extreme cases, but we're not after exhaustive searching are we, merely computationally efficient searching.
My 1-line GP/Pari script finds all P-smooth X^2+1 numbers up to (10^7)^2 in only seconds using this not-completely-naive factoring technique. I only wen't to p=100, and found the following 10^6<X<10^7 1984933 2343692 3449051 6225244
I definitely want to learn about the method. How long does it take for p<=761 and (10^7)^2? Here are my figures: mm=10^6; \\ 10^6 == 13 sec \\ 4620 mm=10^7; \\ 10^7 == 131 sec \\ 8169 mm=10^10; \\ 10^10 == \\ 37h Share that one-liner?
Both methods adapt to ignoring large prime divisors, just keeping those up to some limit.
Lenstra's method is defined in those terms. Alas it's something like O(f(R)*#p) for f(R)=o(R) by default - the larger your set, the more collisions you have, and the more time you spend hopping down lists trying to find the right place to insert values. So it tends to slow down for large sets of divisors, and may take a while before the smoothness-contraint causes the unspecified o(R) to become preferable to the square+1 constraint's O(R^(1/2)).
I have a 'Lenstra Engine' somewhere - if someone wants me to perform some actual hunts for smooth numbers with it, just let me know the details, and I'll get a spare machine onto it immmediately.
Here! Another one I though about: For all small primes P compute the square roots of minus one mod P, mod P^2, mod P^3, etc. Do sieving as follows: For all small primes P. add log(P) at positions X where (X^2+1) divides P; again add log(P) at positions X where (X^2+1) divides P^2; etc. Let L be the log-sum at position X. Check whether log(X^2+1) equals L within precision according to prime limit. Use a fast hand-woven approx-log routine aplog() for the computation of the log(X^2+1). I just haven't found a good alog(), else this would be my next attempt. But C's float log() might just do together with: after X_j was valid (i.e. smooth) the next valid entry L_i must be strictly greater than L_j. Hmmm... this might eliminate almost all log() computation. I'll try that.
* Joerg Arndt <arndt@jjj.de> [Apr 12. 2006 13:06]:
[...]
Another one I though about: For all small primes P compute the square roots of minus one mod P, mod P^2, mod P^3, etc. Do sieving as follows: For all small primes P. add log(P) at positions X where (X^2+1) divides P; again add log(P) at positions X where (X^2+1) divides P^2; etc.
Let L be the log-sum at position X.
Check whether log(X^2+1) equals L within precision according to prime limit. Use a fast hand-woven approx-log routine aplog() for the computation of the log(X^2+1).
I just haven't found a good alog(), else this would be my next attempt. But C's float log() might just do together with: after X_j was valid (i.e. smooth) the next valid entry L_i must be strictly greater than L_j. Hmmm... this might eliminate almost all log() computation.
I'll try that.
Done. And its fast! For X<=10^10 the old code needed 37h, the new one needs 23min30sec. Updating the log as decribed reduces the number of log-computations to (exactly) the number of smooth X found. Using three priority queues, according to size of the prime power sizes. This gives fine cache locality. Now computing up to 10^12, should take 40h. Thanks a lot for the very useful ideas!
participants (2)
-
Joerg Arndt -
Phil Carmody