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

Re: Factors

by kvale (Monsignor)
on Jun 12, 2002 at 21:01 UTC ( #174006=note: print w/replies, xml ) Need Help??

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.



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


    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.

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.


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

Log In?

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://174006]
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others contemplating the Monastery: (6)
As of 2017-09-21 18:04 GMT
Find Nodes?
    Voting Booth?
    During the recent solar eclipse, I:

    Results (251 votes). Check out past polls.