Re: Factors

[kvale]

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

[gumby]

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.

[YuckFoo]

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);
   }

[kvale]

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

[gumby]

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;
   }
}

[gumby]

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;
   }
}