[QUOTE=mahbel;461790]I can only agree with the above. The problem will persist even if I were to provide 10^6 examples because the next example, the 10^6 +1, may not work. I realize that the method is not based on a sound theoretical basis (no derivation from first mathematical principles, no theorems, no nothing).[/QUOTE]

I'm not concerned by the number of examples -- a dozen randomly-selected examples would be plenty to convince me, a million is overkill. I'm concerned about the size of the examples. 121, 91, and 77 are really, really small. You don't need to use huge numbers, but 10 or 20 digits would be much better in terms of convincing me that the method has sufficient generality and utility.

[QUOTE=mahbel;461790]But in more than 18 months of testing, I have never encountered a number, among the small ones that were considered, that cannot be factored by a combination of the 3 combinations (squares, non-squares, and squares + non-squares). But of course, I realize that it's not enough and the numbers tried are small.[/QUOTE]

I doubt it fails for any small number. If you implemented it the way I did, though, it would take roughly time n^(3/2) to go through all sums of four squares, which is considerably worse than trial division. So even if it did work for all numbers, and not just for all small numbers, it's not clear that it could be made practical.

If you [i]do[/i] want to make it practical, we're going to need not just a faster way of generating sums of four squares, but a way to limit the number of sums we use.

[QUOTE=science_man_88;461800]doh before it's [B][U]not[/U][/B] * luck of the draw. my math has the possibility of the first semiprime this may fail for being higher than 544 million because at 8* number of divisors * 80 checks ( though I did the math with 81 including 0 multiple times) means a semiprime which has 4 divisors total that don't divide by 4 ( assuming it's not 4) could have 8*81*number of divisors = 8*81*4= 648*4 -> 23227( prime 2592) ; the nextprime that could possibly be outside of the representations then is 23251; which is only the lowest divisor of a number if it's greater than or equal to 540,609,001 potentially. and proving me wrong without math or fast code ( testing one number a second for example) you will take more than 17 years to disprove/prove it doing them in order potentially. edit: and that's probably even low with my crap math.[/QUOTE]

We need to be careful here. The Jacobi count includes sums like 0^2 + 4^2 + (-5)^2 + 6^2. Excluding the negatives takes us down to about n/2 for hard semiprimes where sigma(n) is close to n. Then the number of checks is, I think, 72 when you cut out the duplication:
[code]mahbelDoubleExtended(n)=my(a2,b2,c2,d,d2,g); for(a=sqrtup(n/4),sqrtint(n), a2=a^2; for(b=sqrtup((n-a2)/3),min(sqrtint(n-a2),a), b2=b^2; for(c=sqrtup((n-a2-b2)/2),min(sqrtint(n-a2-b2),b), c2=c^2; if(issquare(n-a2-b2-c2,&d)&&c>=d, d2=d^2; g=gcd(n,a); if(g>1 && g<n, return(g)); g=gcd(n,b); if(g>1 && g<n, return(g)); g=gcd(n,c); if(g>1 && g<n, return(g)); g=gcd(n,d); if(g>1 && g<n, return(g)); g=gcd(n,a2+b2); if(g>1 && g<n, return(g)); g=gcd(n,a2+c2); if(g>1 && g<n, return(g)); g=gcd(n,a2+d2); if(g>1 && g<n, return(g)); g=gcd(n,d); if(g>1 && g<n, return(g)); g=gcd(n,c); if(g>1 && g<n, return(g)); g=gcd(n,c+d); if(g>1 && g<n, return(g)); g=gcd(n,c+d2); if(g>1 && g<n, return(g)); g=gcd(n,c2+d); if(g>1 && g<n, return(g)); g=gcd(n,b); if(g>1 && g<n, return(g)); g=gcd(n,b+d); if(g>1 && g<n, return(g)); g=gcd(n,b+d2); if(g>1 && g<n, return(g)); g=gcd(n,b+c); if(g>1 && g<n, return(g)); g=gcd(n,b+c+d); if(g>1 && g<n, return(g)); g=gcd(n,b+c+d2); if(g>1 && g<n, return(g)); g=gcd(n,b+c2); if(g>1 && g<n, return(g)); g=gcd(n,b+c2+d); if(g>1 && g<n, return(g)); g=gcd(n,a2+b2-b); if(g>1 && g<n, return(g)); g=gcd(n,b2+d); if(g>1 && g<n, return(g)); g=gcd(n,b2+c); if(g>1 && g<n, return(g)); g=gcd(n,b2+c+d); if(g>1 && g<n, return(g)); g=gcd(n,a2+c2-c); if(g>1 && g<n, return(g)); g=gcd(n,a2+d2-d); if(g>1 && g<n, return(g)); g=gcd(n,a); if(g>1 && g<n, return(g)); g=gcd(n,a+d); if(g>1 && g<n, return(g)); g=gcd(n,a+d2); if(g>1 && g<n, return(g)); g=gcd(n,a+c); if(g>1 && g<n, return(g)); g=gcd(n,a+c+d); if(g>1 && g<n, return(g)); g=gcd(n,a+c+d2); if(g>1 && g<n, return(g)); g=gcd(n,a+c2); if(g>1 && g<n, return(g)); g=gcd(n,a+c2+d); if(g>1 && g<n, return(g)); g=gcd(n,a2-a+b2); if(g>1 && g<n, return(g)); g=gcd(n,a+b); if(g>1 && g<n, return(g)); g=gcd(n,a+b+d); if(g>1 && g<n, return(g)); g=gcd(n,a+b+d2); if(g>1 && g<n, return(g)); g=gcd(n,a+b+c); if(g>1 && g<n, return(g)); g=gcd(n,a+b+c+d); if(g>1 && g<n, return(g)); g=gcd(n,a+b+c+d2); if(g>1 && g<n, return(g)); g=gcd(n,a+b+c2); if(g>1 && g<n, return(g)); g=gcd(n,a+b+c2+d); if(g>1 && g<n, return(g)); g=gcd(n,a2-a+b2-b); if(g>1 && g<n, return(g)); g=gcd(n,a+b2); if(g>1 && g<n, return(g)); g=gcd(n,a+b2+d); if(g>1 && g<n, return(g)); g=gcd(n,a2-a+c2); if(g>1 && g<n, return(g)); g=gcd(n,a+b2+c); if(g>1 && g<n, return(g)); g=gcd(n,a+b2+c+d); if(g>1 && g<n, return(g)); g=gcd(n,a2-a+c2-c); if(g>1 && g<n, return(g)); g=gcd(n,a2-a+d2); if(g>1 && g<n, return(g)); g=gcd(n,a2-a+d2-d); if(g>1 && g<n, return(g)); g=gcd(n,a2-a); if(g>1 && g<n, return(g)); g=gcd(n,a2+d); if(g>1 && g<n, return(g)); g=gcd(n,a2+c); if(g>1 && g<n, return(g)); g=gcd(n,a2+c+d); if(g>1 && g<n, return(g)); g=gcd(n,b2+c2-c); if(g>1 && g<n, return(g)); g=gcd(n,b2+d2-d); if(g>1 && g<n, return(g)); g=gcd(n,a2+b); if(g>1 && g<n, return(g)); g=gcd(n,a2+b+d); if(g>1 && g<n, return(g)); g=gcd(n,b2-b+c2); if(g>1 && g<n, return(g)); g=gcd(n,a2+b+c); if(g>1 && g<n, return(g)); g=gcd(n,a2+b+c+d); if(g>1 && g<n, return(g)); g=gcd(n,b2-b+c2-c); if(g>1 && g<n, return(g)); g=gcd(n,b2-b+d2); if(g>1 && g<n, return(g)); g=gcd(n,b2-b+d2-d); if(g>1 && g<n, return(g)); g=gcd(n,b2-b); if(g>1 && g<n, return(g)); g=gcd(n,c2+d2-d); if(g>1 && g<n, return(g)); g=gcd(n,c2-c+d2); if(g>1 && g<n, return(g)); g=gcd(n,c2-c+d2-d); if(g>1 && g<n, return(g)); g=gcd(n,c2-c); if(g>1 && g<n, return(g)); g=gcd(n,d2-d); if(g>1 && g<n, return(g)))))); "could not factor";[/code]

So the work is about 36n which is a fair bit better than 640n, though it doesn't compare well to trial division's 2*sqrt(n)/log n. (And we still don't know that it can factor everything.)

[QUOTE=CRGreathouse;461803]We need to be careful here. The Jacobi count includes sums like 0^2 + 4^2 + (-5)^2 + 6^2. Excluding the negatives takes us down to about n/2 for hard semiprimes where sigma(n) is close to n. Then the number of checks is, I think, 72 when you cut out the duplication:
[code]mahbelDoubleExtended(n)=my(a2,b2,c2,d,d2,g); for(a=sqrtup(n/4),sqrtint(n), a2=a^2; for(b=sqrtup((n-a2)/3),min(sqrtint(n-a2),a), b2=b^2; for(c=sqrtup((n-a2-b2)/2),min(sqrtint(n-a2-b2),b), c2=c^2; if(issquare(n-a2-b2-c2,&d)&&c>=d, d2=d^2; g=gcd(n,a); if(g>1 && g<n, return(g)); g=gcd(n,b); if(g>1 && g<n, return(g)); g=gcd(n,c); if(g>1 && g<n, return(g)); g=gcd(n,d); if(g>1 && g<n, return(g)); g=gcd(n,a2+b2); if(g>1 && g<n, return(g)); g=gcd(n,a2+c2); if(g>1 && g<n, return(g)); g=gcd(n,a2+d2); if(g>1 && g<n, return(g)); g=gcd(n,d); if(g>1 && g<n, return(g)); g=gcd(n,c); if(g>1 && g<n, return(g)); g=gcd(n,c+d); if(g>1 && g<n, return(g)); g=gcd(n,c+d2); if(g>1 && g<n, return(g)); g=gcd(n,c2+d); if(g>1 && g<n, return(g)); g=gcd(n,b); if(g>1 && g<n, return(g)); g=gcd(n,b+d); if(g>1 && g<n, return(g)); g=gcd(n,b+d2); if(g>1 && g<n, return(g)); g=gcd(n,b+c); if(g>1 && g<n, return(g)); g=gcd(n,b+c+d); if(g>1 && g<n, return(g)); g=gcd(n,b+c+d2); if(g>1 && g<n, return(g)); g=gcd(n,b+c2); if(g>1 && g<n, return(g)); g=gcd(n,b+c2+d); if(g>1 && g<n, return(g)); g=gcd(n,a2+b2-b); if(g>1 && g<n, return(g)); g=gcd(n,b2+d); if(g>1 && g<n, return(g)); g=gcd(n,b2+c); if(g>1 && g<n, return(g)); g=gcd(n,b2+c+d); if(g>1 && g<n, return(g)); g=gcd(n,a2+c2-c); if(g>1 && g<n, return(g)); g=gcd(n,a2+d2-d); if(g>1 && g<n, return(g)); g=gcd(n,a); if(g>1 && g<n, return(g)); g=gcd(n,a+d); if(g>1 && g<n, return(g)); g=gcd(n,a+d2); if(g>1 && g<n, return(g)); g=gcd(n,a+c); if(g>1 && g<n, return(g)); g=gcd(n,a+c+d); if(g>1 && g<n, return(g)); g=gcd(n,a+c+d2); if(g>1 && g<n, return(g)); g=gcd(n,a+c2); if(g>1 && g<n, return(g)); g=gcd(n,a+c2+d); if(g>1 && g<n, return(g)); g=gcd(n,a2-a+b2); if(g>1 && g<n, return(g)); g=gcd(n,a+b); if(g>1 && g<n, return(g)); g=gcd(n,a+b+d); if(g>1 && g<n, return(g)); g=gcd(n,a+b+d2); if(g>1 && g<n, return(g)); g=gcd(n,a+b+c); if(g>1 && g<n, return(g)); g=gcd(n,a+b+c+d); if(g>1 && g<n, return(g)); g=gcd(n,a+b+c+d2); if(g>1 && g<n, return(g)); g=gcd(n,a+b+c2); if(g>1 && g<n, return(g)); g=gcd(n,a+b+c2+d); if(g>1 && g<n, return(g)); g=gcd(n,a2-a+b2-b); if(g>1 && g<n, return(g)); g=gcd(n,a+b2); if(g>1 && g<n, return(g)); g=gcd(n,a+b2+d); if(g>1 && g<n, return(g)); g=gcd(n,a2-a+c2); if(g>1 && g<n, return(g)); g=gcd(n,a+b2+c); if(g>1 && g<n, return(g)); g=gcd(n,a+b2+c+d); if(g>1 && g<n, return(g)); g=gcd(n,a2-a+c2-c); if(g>1 && g<n, return(g)); g=gcd(n,a2-a+d2); if(g>1 && g<n, return(g)); g=gcd(n,a2-a+d2-d); if(g>1 && g<n, return(g)); g=gcd(n,a2-a); if(g>1 && g<n, return(g)); g=gcd(n,a2+d); if(g>1 && g<n, return(g)); g=gcd(n,a2+c); if(g>1 && g<n, return(g)); g=gcd(n,a2+c+d); if(g>1 && g<n, return(g)); g=gcd(n,b2+c2-c); if(g>1 && g<n, return(g)); g=gcd(n,b2+d2-d); if(g>1 && g<n, return(g)); g=gcd(n,a2+b); if(g>1 && g<n, return(g)); g=gcd(n,a2+b+d); if(g>1 && g<n, return(g)); g=gcd(n,b2-b+c2); if(g>1 && g<n, return(g)); g=gcd(n,a2+b+c); if(g>1 && g<n, return(g)); g=gcd(n,a2+b+c+d); if(g>1 && g<n, return(g)); g=gcd(n,b2-b+c2-c); if(g>1 && g<n, return(g)); g=gcd(n,b2-b+d2); if(g>1 && g<n, return(g)); g=gcd(n,b2-b+d2-d); if(g>1 && g<n, return(g)); g=gcd(n,b2-b); if(g>1 && g<n, return(g)); g=gcd(n,c2+d2-d); if(g>1 && g<n, return(g)); g=gcd(n,c2-c+d2); if(g>1 && g<n, return(g)); g=gcd(n,c2-c+d2-d); if(g>1 && g<n, return(g)); g=gcd(n,c2-c); if(g>1 && g<n, return(g)); g=gcd(n,d2-d); if(g>1 && g<n, return(g)))))); "could not factor";[/code]

So the work is about 36n which is a fair bit better than 640n, though it doesn't compare well to trial division's 2*sqrt(n)/log n. (And we still don't know that it can factor everything.)[/QUOTE]

well if N is odd we know it's divisors are odd so we only need to check if the sum is divisible by odd values ( so if b,c and d are even for example b+c+d=even and is only possibly divisible by one of N's divisors if the sum is 2 mod 4).

[QUOTE=CRGreathouse;461803]...So the work is about 36n which is a fair bit better than 640n, though it doesn't compare well to trial division's 2*sqrt(n)/log n. (And we still don't know that it can factor everything.)[/QUOTE]
You only need to scan a few 9-digit semiprime inputs (with both factors a random 5-digit prime) and to find it fail.

[QUOTE=Batalov;461810]You only need to scan a few 9-digit semiprime inputs (with both factors a random 5-digit prime) and to find it fail.[/QUOTE]

I think I'll need to be cleverer than that.

I'm running some random tests now as you suggest, but no failures so far. But remember, the number of tests is mind-bogglingly large: we're looking for a multiple of p, which happens with frequency at least 1/sqrt(n), but we're doing 36n tests -- 36*sqrt(n) times more than the expected number needed. So if we were just picking the numbers randomly, there's essentially no way the test would fail: (1 - 1/sqrt(n))^(36n) --> exp(-36*sqrt(n)) probability of failure, which is in the neighborhood of 10^-156354 already by 10^8.

So the only way it could be expected to fail is if the numbers it picked were chosen in such a way to avoid hitting the prime factors of n. But I think this method might fit the bill. If n = a^2 + b^2 + c^2 + d^2, then a, b, c, and d are fairly small, and many of the numbers they generate will be duplicates. There may be ways of limiting their ability to combine to form the smaller prime p. As for the squares, n can be chosen such that it is not the sum of three squares so the only possible combinations involve squares and first powers. Then maybe some modular or other trickery can make it hard to hit p. Then just apply the method a lot, hoping it avoids q as well (you may have some degrees of freedom left to affect it also).

In any case this would probably take real work and for now I'm just spitballing. But it's interesting to think about.

The probability of other powers of a,b,c or d or any random number less than N (or greater than N if you allow subtraction), yielding GCD factors is high enough to justify no good reason not to try say summation of cubes (our mixed powers such as cubes and squares ) before moving to other 4 squares. Also I think there is no justification to perform exclusive combinations of squire sets. Why not mix them up?
Larger semi primes are less likely to yield, since they are less probable to have a common factor with sum of small random coprimes.

[QUOTE=a1call;461819]The probability of other powers of a,b,c or d or any random number less than N (or greater than N if you allow subtraction), yielding GCD factors is high enough to justify no good reason not to try say summation of cubes (our mixed powers such as cubes and squares ) before moving to other 4 squares. Also I think there is no justification to perform exclusive combinations of squire sets. Why not mix them up?[/QUOTE]

I don't really know why you'd use any of these methods. They seem much worse than ordinary trial division. That's why I'm interested in the geometric side of things: maybe there is some pattern that can be seen there and taken advantage of. Most probably it would only lead to the factorization of (thin) special classes of numbers, but you never know.

[QUOTE=CRGreathouse;461803]We need to be careful here. The Jacobi count includes sums like 0^2 + 4^2 + (-5)^2 + 6^2. Excluding the negatives takes us down to about n/2 for hard semiprimes where sigma(n) is close to n. Then the number of checks is, I think, 72 when you cut out the duplication:
[code]mahbelDoubleExtended(n)={my(a2,b2,c2,d,d2,g);
for(a=sqrtup(n/4),sqrtint(n), a2=a^2;
for(b=sqrtup((n-a2)/3),min(sqrtint(n-a2),a), b2=b^2;
for(c=sqrtup((n-a2-b2)/2),min(sqrtint(n-a2-b2),b), c2=c^2;
if(issquare(n-a2-b2-c2,&d)&&c>=d, d2=d^2;
[COLOR=DarkRed]g=gcd(n,a); if(g>1 && g<n, return(g));
g=gcd(n,b); if(g>1 && g<n, return(g));
g=gcd(n,c); if(g>1 && g<n, return(g));
g=gcd(n,d); if(g>1 && g<n, return(g));[/COLOR]
g=gcd(n,a2+b2); if(g>1 && g<n, return(g));
g=gcd(n,a2+c2); if(g>1 && g<n, return(g));
g=gcd(n,a2+d2); if(g>1 && g<n, return(g));
[COLOR=Red]g=gcd(n,d); if(g>1 && g<n, return(g));
g=gcd(n,c); if(g>1 && g<n, return(g));[/COLOR]
g=gcd(n,c+d); if(g>1 && g<n, return(g));
g=gcd(n,c+d2); if(g>1 && g<n, return(g));
g=gcd(n,c2+d); if(g>1 && g<n, return(g));
[COLOR=Red]g=gcd(n,b); if(g>1 && g<n, return(g));[/COLOR]
g=gcd(n,b+d); if(g>1 && g<n, return(g));
g=gcd(n,b+d2); if(g>1 && g<n, return(g));
g=gcd(n,b+c); if(g>1 && g<n, return(g));
g=gcd(n,b+c+d); if(g>1 && g<n, return(g));
g=gcd(n,b+c+d2); if(g>1 && g<n, return(g));
g=gcd(n,b+c2); if(g>1 && g<n, return(g));
g=gcd(n,b+c2+d); if(g>1 && g<n, return(g));
g=gcd(n,a2+b2-b); if(g>1 && g<n, return(g));
g=gcd(n,b2+d); if(g>1 && g<n, return(g));
g=gcd(n,b2+c); if(g>1 && g<n, return(g));
g=gcd(n,b2+c+d); if(g>1 && g<n, return(g));
g=gcd(n,a2+c2-c); if(g>1 && g<n, return(g));
g=gcd(n,a2+d2-d); if(g>1 && g<n, return(g));
[COLOR=Red]g=gcd(n,a); if(g>1 && g<n, return(g));[/COLOR]
g=gcd(n,a+d); if(g>1 && g<n, return(g));
g=gcd(n,a+d2); if(g>1 && g<n, return(g));
g=gcd(n,a+c); if(g>1 && g<n, return(g));
g=gcd(n,a+c+d); if(g>1 && g<n, return(g));
g=gcd(n,a+c+d2); if(g>1 && g<n, return(g));
g=gcd(n,a+c2); if(g>1 && g<n, return(g));
g=gcd(n,a+c2+d); if(g>1 && g<n, return(g));
g=gcd(n,a2-a+b2); if(g>1 && g<n, return(g));
g=gcd(n,a+b); if(g>1 && g<n, return(g));
g=gcd(n,a+b+d); if(g>1 && g<n, return(g));
g=gcd(n,a+b+d2); if(g>1 && g<n, return(g));
g=gcd(n,a+b+c); if(g>1 && g<n, return(g));
g=gcd(n,a+b+c+d); if(g>1 && g<n, return(g));
g=gcd(n,a+b+c+d2); if(g>1 && g<n, return(g));
g=gcd(n,a+b+c2); if(g>1 && g<n, return(g));
g=gcd(n,a+b+c2+d); if(g>1 && g<n, return(g));
g=gcd(n,a2-a+b2-b); if(g>1 && g<n, return(g));
g=gcd(n,a+b2); if(g>1 && g<n, return(g));
g=gcd(n,a+b2+d); if(g>1 && g<n, return(g));
g=gcd(n,a2-a+c2); if(g>1 && g<n, return(g));
g=gcd(n,a+b2+c); if(g>1 && g<n, return(g));
g=gcd(n,a+b2+c+d); if(g>1 && g<n, return(g));
g=gcd(n,a2-a+c2-c); if(g>1 && g<n, return(g));
g=gcd(n,a2-a+d2); if(g>1 && g<n, return(g));
g=gcd(n,a2-a+d2-d); if(g>1 && g<n, return(g));
g=gcd(n,a2-a); if(g>1 && g<n, return(g));
g=gcd(n,a2+d); if(g>1 && g<n, return(g));
g=gcd(n,a2+c); if(g>1 && g<n, return(g));
g=gcd(n,a2+c+d); if(g>1 && g<n, return(g));
g=gcd(n,b2+c2-c); if(g>1 && g<n, return(g));
g=gcd(n,b2+d2-d); if(g>1 && g<n, return(g));
g=gcd(n,a2+b); if(g>1 && g<n, return(g));
g=gcd(n,a2+b+d); if(g>1 && g<n, return(g));
g=gcd(n,b2-b+c2); if(g>1 && g<n, return(g));
g=gcd(n,a2+b+c); if(g>1 && g<n, return(g));
g=gcd(n,a2+b+c+d); if(g>1 && g<n, return(g));
g=gcd(n,b2-b+c2-c); if(g>1 && g<n, return(g));
g=gcd(n,b2-b+d2); if(g>1 && g<n, return(g));
g=gcd(n,b2-b+d2-d); if(g>1 && g<n, return(g));
g=gcd(n,b2-b); if(g>1 && g<n, return(g));
g=gcd(n,c2+d2-d); if(g>1 && g<n, return(g));
g=gcd(n,c2-c+d2); if(g>1 && g<n, return(g));
g=gcd(n,c2-c+d2-d); if(g>1 && g<n, return(g));
g=gcd(n,c2-c); if(g>1 && g<n, return(g));
g=gcd(n,d2-d); if(g>1 && g<n, return(g))))));
"could not factor";}[/code]So the work is about 36n which is a fair bit better than 640n, though it doesn't compare well to trial division's 2*sqrt(n)/log n. (And we still don't know that it can factor everything.)[/QUOTE]
There are some repeats in the gcd code... ([COLOR=DarkRed]above[/COLOR])

...and some -'s are in there already (even though, I see that they are equivalent to +'s if subtracted from n: e.g.
gcd(n,c2-c+d2-d) = gcd(n,n-(c2-c+d2-d)) = gcd(n,a2+b2+c+d))

[QUOTE=Batalov;461821]There are some repeats in the gcd code...[/QUOTE]

Thanks! Must have been from combining the existing program with the new list, or something.

[QUOTE=Batalov;461821]...and some -'s are in there already (even though, I see that they are equivalent to +'s if subtracted from n: e.g.
gcd(n,c2-c+d2-d) = gcd(n,n-(c2-c+d2-d)) = gcd(n,a2+b2+c+d))[/QUOTE]

Right, just trying to write them in the tersest form. A full expansion to subtraction would add a lot more terms.

Fixed version (thanks Serge):
[code]sqrtup(n)=n=ceil(n);if(issquare(n,&n),n,sqrtint(n)+1);
mahbelDoubleExtended(n)=my(a2,b2,c2,d,d2,g); for(a=sqrtup(n/4),sqrtint(n), a2=a^2; for(b=sqrtup((n-a2)/3),min(sqrtint(n-a2),a), b2=b^2; for(c=sqrtup((n-a2-b2)/2),min(sqrtint(n-a2-b2),b), c2=c^2; if(issquare(n-a2-b2-c2,&d)&&c>=d, d2=d^2; g=gcd(n,a); if(g>1, return(g));g=gcd(n,b); if(g>1, return(g));g=gcd(n,c); if(g>1, return(g));g=gcd(n,d); if(g>1, return(g));g=gcd(n,a2+b2); if(g>1 && g<n, return(g));g=gcd(n,a2+c2); if(g>1 && g<n, return(g));g=gcd(n,a2+d2); if(g>1 && g<n, return(g));g=gcd(n,c+d); if(g>1 && g<n, return(g));g=gcd(n,c+d2); if(g>1 && g<n, return(g));g=gcd(n,c2+d); if(g>1 && g<n, return(g));g=gcd(n,b+d); if(g>1 && g<n, return(g));g=gcd(n,b+d2); if(g>1 && g<n, return(g));g=gcd(n,b+c); if(g>1 && g<n, return(g));g=gcd(n,b+c+d); if(g>1 && g<n, return(g));g=gcd(n,b+c+d2); if(g>1 && g<n, return(g));g=gcd(n,b+c2); if(g>1 && g<n, return(g));g=gcd(n,b+c2+d); if(g>1 && g<n, return(g));g=gcd(n,a2+b2-b); if(g>1 && g<n, return(g));g=gcd(n,b2+d); if(g>1 && g<n, return(g));g=gcd(n,b2+c); if(g>1 && g<n, return(g));g=gcd(n,b2+c+d); if(g>1 && g<n, return(g));g=gcd(n,a2+c2-c); if(g>1 && g<n, return(g));g=gcd(n,a2+d2-d); if(g>1 && g<n, return(g));g=gcd(n,a+d); if(g>1 && g<n, return(g));g=gcd(n,a+d2); if(g>1 && g<n, return(g));g=gcd(n,a+c); if(g>1 && g<n, return(g));g=gcd(n,a+c+d); if(g>1 && g<n, return(g));g=gcd(n,a+c+d2); if(g>1 && g<n, return(g));g=gcd(n,a+c2); if(g>1 && g<n, return(g));g=gcd(n,a+c2+d); if(g>1 && g<n, return(g));g=gcd(n,a2-a+b2); if(g>1 && g<n, return(g));g=gcd(n,a+b); if(g>1 && g<n, return(g));g=gcd(n,a+b+d); if(g>1 && g<n, return(g));g=gcd(n,a+b+d2); if(g>1 && g<n, return(g));g=gcd(n,a+b+c); if(g>1 && g<n, return(g));g=gcd(n,a+b+c+d); if(g>1 && g<n, return(g));g=gcd(n,a+b+c+d2); if(g>1 && g<n, return(g));g=gcd(n,a+b+c2); if(g>1 && g<n, return(g));g=gcd(n,a+b+c2+d); if(g>1 && g<n, return(g));g=gcd(n,a2-a+b2-b); if(g>1 && g<n, return(g));g=gcd(n,a+b2); if(g>1 && g<n, return(g));g=gcd(n,a+b2+d); if(g>1 && g<n, return(g));g=gcd(n,a2-a+c2); if(g>1 && g<n, return(g));g=gcd(n,a+b2+c); if(g>1 && g<n, return(g));g=gcd(n,a+b2+c+d); if(g>1 && g<n, return(g));g=gcd(n,a2-a+c2-c); if(g>1 && g<n, return(g));g=gcd(n,a2-a+d2); if(g>1 && g<n, return(g));g=gcd(n,a2-a+d2-d); if(g>1 && g<n, return(g));g=gcd(n,a2-a); if(g>1 && g<n, return(g));g=gcd(n,a2+d); if(g>1 && g<n, return(g));g=gcd(n,a2+c); if(g>1 && g<n, return(g));g=gcd(n,a2+c+d); if(g>1 && g<n, return(g));g=gcd(n,b2+c2-c); if(g>1 && g<n, return(g));g=gcd(n,b2+d2-d); if(g>1 && g<n, return(g));g=gcd(n,a2+b); if(g>1 && g<n, return(g));g=gcd(n,a2+b+d); if(g>1 && g<n, return(g));g=gcd(n,b2-b+c2); if(g>1 && g<n, return(g));g=gcd(n,a2+b+c); if(g>1 && g<n, return(g));g=gcd(n,a2+b+c+d); if(g>1 && g<n, return(g));g=gcd(n,b2-b+c2-c); if(g>1 && g<n, return(g));g=gcd(n,b2-b+d2); if(g>1 && g<n, return(g));g=gcd(n,b2-b+d2-d); if(g>1 && g<n, return(g));g=gcd(n,b2-b); if(g>1 && g<n, return(g));g=gcd(n,c2+d2-d); if(g>1 && g<n, return(g));g=gcd(n,c2-c+d2); if(g>1 && g<n, return(g));g=gcd(n,c2-c+d2-d); if(g>1 && g<n, return(g));g=gcd(n,c2-c); if(g>1 && g<n, return(g));g=gcd(n,d2-d); if(g>1 && g<n, return(g)))))); "could not factor";[/code]

The count remains 72; I duplicated the cases a, b, c, and d but did not double-count them.

[QUOTE=CRGreathouse;461780]Another problem is that the examples used are small, and small numbers are weird. See Guy's papers on the Small Law of Small Numbers.[/QUOTE]
Despite talking of producing 10^6 examples, the [i]only[/i] example the OP has shown us entirely on his own, is 91; he also managed to produce a factor of 77 by changing his method, after I pointed out that his originally-described method wouldn't factor it. His repeated insistence that his "strong extended" method stumbles over a factor, doesn't make 77 get any bigger. He has declined to factor so small a challenge number as 19673, which is amenable to very simple methods.

The OP has also failed to give any reason why anyone should think his method is generally applicable.

At this point, the most interesting possibility I can see in pursuing this topic, is a faster method for producing representations of N as a sum of four squares. If there were any way to narrow the focus to those representations that are somehow likely to produce factors, that might be interesting. At this point, all I can see in that regard is maybe excluding possibilities using simple inequalities.

An amusing possible result of trying this method did occur to me: It might happen that you try so [i]many[/i] representations without success, that you obtain a proof that N is composite without finding a factor! If you find enough representations to prove that

r[sub]4[/sub](n) > 8*(N + 1),

without finding a factor, then you [i]do[/i] have a proof that N is composite!

I hasten to point out that this would probably qualify as the world's slowest compositeness test.