Beefy Boxes and Bandwidth Generously Provided by pair Networks
Clear questions and runnable code
get the best and fastest answer
 
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
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 musing on the Monastery: (11)
As of 2014-08-22 16:54 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    The best computer themed movie is:











    Results (161 votes), past polls