* Schroeppel, Richard <rschroe@sandia.gov> [Apr 12. 2006 11:05]:
Speculations on finding "smooth" X^2+1:
[What was wrong with "round"? That's what Hardy & Wright use.]
Smooth as in factoring papers (where we also have "power smooth"). Sorry if that caused confusion.
You can generate ordinary integers in order, along with their prime factorizations, using a priority-queue-like scheme. [Suppose the current integer is 6, with divisors 2 and 3. Place queue entries for "2 divides 8" and "3 divides 9". Minor adjustments for prime powers.] 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. Both methods are faster than factoring N (or X^2+1) one-at-a-time. Both methods adapt to ignoring large prime divisors, just keeping those up to some limit.
One can (an IMHO should) avoid all complex quantities. I see two basic versions of priority-queue algorithm (and the priority-queue just makes it cache friendly): 1. mark the good ones (those that must have at least one divisor 5,13,...,761 (the first 64). Problem: almost all are marked, rendering the scheme practically useless. Well, about 23% ? x=1.0;forprime(p=5,761,if(1==p%4, x*=(1-2/p))); x 0.230249577826900 2. mark the bad ones, those that, have a prime factor among the first, say, 10000 primes > 761. Problem: almost nothing is marked, bad again. With all 'big' primes p=4k+1 and 761<p<10^7 we get no practical win: ? x=1.0;forprime(p=762, 10000000,if(1==p%4, x*=(1-2/p))); x 0.415370346793571 Practically, I can factor up to a given limit (pari/gp allows that) and test whether the biggest factor is greater than my 761. This is blazing fast! A quick C program with 128-bit integers was just about 15% faster than the script! A factor of two might be feasible with hand-written (well, generated) assembler, but that is not enough.
If you choose your prime divisor limit ahead of time, you can generate all numbers < 10^12 that factor into just those primes. There's an obvious method that just starts with a list {1}, multiplies in powers of 2 to get a longer list, then multiplies that list with powers of 5, then 13, etc. This would have to be filtered by "is the number an X^2+1?", but that's cheap.
Yes, par/gp's issquare() does it extremely efficiently (as described in Cohen's book). Problem: with more than 20 primes the program runs virtually forever.
The list length could be a problem, and the numbers aren't exactly generated in order. There's a two-list variation that builds two lists from the first and second halves of the prime-set, and then generates the cross product of the two lists in numerical order (and applies the X^2+1 filter).
Will that really be practical?
I believe that one can design an algorithm that will directly generate numbers with all prime factors in a given list, in order of size, and not much memory, by keeping a list of "useful steps" along with a recent history. With list = {2,3}, starting from 1, the first steps are 2 and 3, updated to 3/2 and 4/3, then 9/8, 32/27, and so on. The two-prime case looks a lot like finding continued fraction approximants to log3/log2. I'm not sure if the number of steps required in the step-list for the more-prime case is too big.
Sounds a bit like the recursion that tries all divisors d of X^2+1 and tests whether (X+d)^2+1 is (also) smooth. The algorithm generates monsters like X=15689622322199725018 15689622322199725018^2+1= [5 2] [29 2] [41 2] [89 1] [101 2] [113 1] [197 1] [229 2] [277 2] [353 1] [569 1] [653 2] within minutes(!) However, the method misses quite a few valid X.
Another approach is to use LLL on the arctangents, directly looking for near-relations among the arctangents of small complex primes. Apparently LLL can be run on matrices of a few hundred rows.
Bad! Nullspace computations are enough: Attach to the exponents of a prime P a minus if 2*P<X, take the nullspace of the matrix obtained. E.g. the numbers [173932, 1068, 515, 239, 192] give the matrix 173932 [-5, -1, 1, -2] 1068 [ 6, 0, 1, 0] (1068^2+1==5^6*73) 515 [ 0, 1, 0, -2] 239 [ 0, -4, 0, 0] (239^2+1==2*13^4) 192 [-1, 0, 1, 1] ? K=matkerint(M) [56] [32] [-100] [-39] [-88] Thereby: +88[192] +39[239] +100[515] -32[1068] -56[173932] == 1 * Pi/4 One does have to test what multiple of Pi/4 is actually obtained, often it is just zero.
Rich
Thanks a lot for your thoughts!
-----Original Message----- [...]
The gap between what is known (well, to me) and what can be done seems huge: 1) how to prove the number of N-term relations is even finite? 2) given a set of primes smaller than some limit M, is there any known bound on the largets X so that (X^2+1) is M-smooth? 3) why do the "best" n-term relations have such a striking tendency to prefer the smallest primes? (If it wasn't for that I wouldn't be so confident that the relations up to, say, 11-terms might actually be the best possible). 4) How come that N-term relations tend to make a much smaller progress over (N-1)-term relations if N is even? (very mysterios IMHO) I'll try the sieving and/or some more tree search with added heuristics (like forbidding large powers of a single prime). Or I just let the simple search run for a long time. Or I might just see a doctor to get over the arctan addiction...