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)*(xy) = 0 (mod p) and
so (x+y) or (xy) may be factors of p.
In practice, one starts with an x that has the smallest
square larger than p, then computes x^2p 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 ddigit 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
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 100digits 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.  [reply] 
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 * xy = 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);
}
 [reply] [d/l] 

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 (xy) 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  [reply] [d/l] [select] 

If you want to find a square root of a (mod p) use the ShanksTonelli algorithm.
# Calculate the least nonnegative 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;
}
# ShanksTonelli 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 nonresidue mod p
$s = $q;
$e = 0;
while ($s % 2 == 0) {
$s = $s/2;
$e++;
}
# p1 = 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 nonresidue mod p
$g = mpow($n, $s, $p);
$r = $e;
OUT: {
do {
$y = $b;
$m = 0;
IN: while ($m <= $r1) {
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;
}
}
 [reply] [d/l] 

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**(($p1)/2) % $p == 1;
$s = $p1;
$e = 0;
while ($s % 2 == 0) {
$s = $s / 2;
$e++;
}
for ($n = 2; $n >= 2; $n++) {
last if $n**(($p1)/2) % $p == $p1;
}
$x = $a**(($s+1)/2) % $p;
$b = $a**$s % $p;
$g = $n**$s % $p;
$r = $e;
$t = $b;
while (1) {
for ($m = 0; $m <= $r1; $m++) {
last if $t % $p == 1;
$t = $t**2;
}
return $x if $m == 0;
$x = $x * $g**(2*($r$m1)) % $p;
$b = $b * $g**(2*($r$m1)) % $p;
$g = $g**(2*($r$m)) % $p;
$r = $m;
}
}
 [reply] [d/l] 

