Beefy Boxes and Bandwidth Generously Provided by pair Networks
The stupid question is the question not asked
 
PerlMonks  

Re: Re: Factors

by YuckFoo (Abbot)
on Jun 12, 2002 at 22:30 UTC ( #174047=note: print w/ replies, xml ) Need Help??


in reply to Re: Factors
in thread Vampire Numbers

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


Comment on Re: Re: Factors
Download Code
Replies are listed 'Best First'.
Re: Factors
by kvale (Monsignor) on Jun 13, 2002 at 04:51 UTC
    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
Re: Re: Re: Factors
by gumby (Scribe) on Jun 13, 2002 at 15:43 UTC
    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?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://174047]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others perusing the Monastery: (13)
As of 2015-07-31 12:44 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    The top three priorities of my open tasks are (in descending order of likelihood to be worked on) ...









    Results (277 votes), past polls