Beefy Boxes and Bandwidth Generously Provided by pair Networks
Syntactic Confectionery Delight


by gumby (Scribe)
on Jun 11, 2002 at 18:19 UTC ( #173606=note: print w/replies, xml ) Need Help??

in reply to Vampire Numbers

A more efficient algorithm (of Gauss and later Kraitchik) would be to find solutions to the congruence x^2=y^2 (mod p) e.g.

38^2=3^2 (mod 1435)

It is trivial to then check if the factor (38+3 in this case) is a permutation of n/2 digits from the original number.

Replies are listed 'Best First'.
Re: Factors
by kvale (Monsignor) on Jun 12, 2002 at 21:01 UTC
    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.


      The Quadratic Sieve has running time


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


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


      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.

      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.


      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% 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.

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

Log In?

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://173606]
[marto]: interesting

How do I use this? | Other CB clients
Other Users?
Others having an uproarious good time at the Monastery: (5)
As of 2018-04-20 05:25 GMT
Find Nodes?
    Voting Booth?