http://www.perlmonks.org?node_id=174006

in reply to Factors
in thread Vampire Numbers

I think what is being suggested here is that instead of creating all possible pairs of numbers (x,y) from the digits of a number p and testing if x*y == p, it is more efficient to factor the number p first and and then test all possible factorizations to see if they are just a permuation of the digits of p.

Finidng a quadratic residue x^2=y^2 (mod p) is useful, because then x^2 - y^2 = (x+y)*(x-y) = 0 (mod p) and so (x+y) or (x-y) may be factors of p.

In practice, one starts with an x that has the smallest square larger than p, then computes x^2-p to see if it is a perfect square. If not, increment x and repeat.

This is pretty fast relative to testing by division from 2 to sqrt(p), but there exist even faster methods based on the quadratic residue. The Quadratic Sieve, invented by Pomerance in 1981 is a direct derivative and the Number Field Sieve is a related method.

Getting back to Vampire numbers, what are the relative efficiencies for a number p with 2d digits? For the first method, a naive approach generates all permutations of 2*d digits and splits each into two d-digit numbers to test. Thus each 2d digit number requires about (2*d)! == (2*log(p))! tests, or (2*log p)^(2*log p) by Stirling's approximation. For the second method, a quadratic sieve has running time of order exp(sqrt(log(p)*log(log(p)))) to find a single factor of p. So finding a single factor using QS is more efficient that testing all factors using the first method. But I don't know the running time for finding all the factors using the QS and I don't know the typical multiplier for the scaling result of the QS method.

Anyone?

-Mark

Replies are listed 'Best First'.
Re: Re: Factors
by gumby (Scribe) on Jun 12, 2002 at 21:45 UTC
The Quadratic Sieve has running time

O(e^(sqrt(1.125ln(n)ln(ln(n)))))

But with certain improvements (such as using Wiedemann's Gaussian Elimination algorithm over GF(2)), it's running time is given by

O(e^(sqrt(ln(n)ln(ln(n)))))

The Number Field Sieve is better for numbers over 100-digits and has running time

O(e^(1.9223((ln(n))^1/3(ln(ln(n))))^2/3))

Excuse my notation (but TeX would be overkill).

Note: O(stuff) means the magnitude of the running time is < A * stuff for some constant A and all values of n.

Re: Re: Factors
by YuckFoo (Abbot) on Jun 12, 2002 at 22:30 UTC
It seemed like a good idea so I cooked it up according to your description. Correct me if I'm wrong but it seems quadratic residues won't account for all factors.

Program correctly factors 24 to 4*6 and 2*12 but misses 3*8. Program finds no factors for 54.

YuckFoo

Update: With a little thought one can see the program will only find factor pairs that are both even or both odd. x+y * x-y = p.

#!/usr/bin/perl use strict; my (\$num) = @ARGV; my \$x = int(sqrt(\$num)) + 1; my \$y = sqrt(\$x * \$x - \$num); while (\$x - \$y >= 2) { if (\$y == sprintf("%d", \$y)) { print "Ok \$num = ", \$x-\$y, " * ", \$y+\$x, "\n"; } else { print "No x = \$x, y = \$y\n"; } \$x++; \$y = sqrt(\$x * \$x - \$num); }
Good observation! It looks like Fermat's method (the method I mentioned above) doesn't necessarily find all, or even any of the factors. It is just a heuristic. On the other hand, if you have an even number, you may extract factors of 2 until you get an odd number, in which case (x+y) and (x-y) must both be necessarily odd.

The generalization that gumby mentioned, x^2 = y^2 (mod p) may have more sucess, but I don't know if it is exhaustive. Here is code that implements it:
my \$p = shift; my %residues; my %factors; foreach my \$x (2..\$p) { my \$mod = (\$x*\$x) % \$p; if (exists \$residues{\$mod}) { \$factors{\$x-\$residues{\$mod}}++ if \$p % (\$x-\$residues{\$mod}) == 0 +; \$factors{\$x+\$residues{\$mod}}++ if \$p % (\$x+\$residues{\$mod}) == 0 +; } else { \$residues{\$mod} = \$x; } } foreach (sort {\$a <=> \$b} keys %factors) { print "\$p = \$_*", \$p/\$_, "\n"; }
For 54, this yields
1021% factor.pl 54 54 = 2*27 54 = 6*9 54 = 18*3 54 = 54*1
But this code is unsatisfying because it takes longer than the simple trial by division. Thinking from a divide and conquer point of view, it is probably faster to use these methods to find a factorization, then recursively find factorizations of the factors untill all are prime.

-Mark
If you want to find a square root of a (mod p) use the Shanks-Tonelli algorithm.
# Calculate the least non-negative remainder when an integer a # is divided by a positive integer b. sub mod { my (\$a, \$b) = @_; my \$c = \$a % \$b; return \$c if (\$a >= 0); return 0 if (\$c == 0); return (\$c + \$b); } # Calculate a^b (mod c), where a,b and c are integers and a,b>=0, c>=1 + sub mpow { my (\$a, \$b, \$c) = @_; my (\$x, \$y, \$z) = (\$a, \$b, 1); while (\$y > 0) { while (\$y % 2 == 0) { \$y = \$y / 2; \$x = \$x**2 % \$c; } \$y--; \$z = mod(\$z * \$x, \$c); } return \$z; } # Shanks-Tonelli algorithm to calculate y^2 = a (mod p) for p an odd p +rime sub tonelli { my (\$a, \$p) = @_; my (\$b, \$e, \$g, \$h, \$i, \$m, \$n, \$q, \$r, \$s, \$t, \$x, \$y, \$z); \$q = \$p - 1; \$t = mpow(\$a, \$q/2, \$p); return 0 if (\$t = \$q); # a is a quadratic non-residue mod p \$s = \$q; \$e = 0; while (\$s % 2 == 0) { \$s = \$s/2; \$e++; } # p-1 = s * 2^e; \$x = mpow(\$a, ++\$s/2, \$p); \$b = mpow(\$a, \$s, \$p); return \$x if (\$b == 1); \$n = 2; RES: while (\$n >= 2) { \$t = mpow(\$n, \$q/2, \$p); last RES if (\$t == \$q); } continue { \$n++; } # n is the least quadratic non-residue mod p \$g = mpow(\$n, \$s, \$p); \$r = \$e; OUT: { do { \$y = \$b; \$m = 0; IN: while (\$m <= \$r-1) { last IN if (\$y == 1); \$y = \$y**2 % \$p; } continue { \$m++; } return \$x if (\$m == 0); \$z = \$r - \$m - 1; \$h = \$g; for (\$i = 1; \$i < \$z; \$i++) {\$h = \$h**2 % \$p;} \$x = (\$x * \$h) % \$p; \$b = (\$b * \$h**2) % \$p; \$g = \$h**2 % \$p; \$r = \$m; } while 1; } }
Here's a more 'reference' style implementation. Basically, this is the algorithm you'll find in the paper 'Square Roots, From 1; 24, 51, 10 To Dan Shanks' by Ezra Brown, but at the cost of using Math::BigInt.
#!/usr/bin/perl -w use Math::BigInt ':constant'; use strict; sub tonelli { my (\$s, \$e, \$n, \$x, \$b, \$g, \$r, \$m, \$t); my (\$a, \$p) = @_; die "\$a has no square roots (mod \$p)" if \$a**((\$p-1)/2) % \$p == -1; \$s = \$p-1; \$e = 0; while (\$s % 2 == 0) { \$s = \$s / 2; \$e++; } for (\$n = 2; \$n >= 2; \$n++) { last if \$n**((\$p-1)/2) % \$p == \$p-1; } \$x = \$a**((\$s+1)/2) % \$p; \$b = \$a**\$s % \$p; \$g = \$n**\$s % \$p; \$r = \$e; \$t = \$b; while (1) { for (\$m = 0; \$m <= \$r-1; \$m++) { last if \$t % \$p == 1; \$t = \$t**2; } return \$x if \$m == 0; \$x = \$x * \$g**(2*(\$r-\$m-1)) % \$p; \$b = \$b * \$g**(2*(\$r-\$m-1)) % \$p; \$g = \$g**(2*(\$r-\$m)) % \$p; \$r = \$m; } }