Re: [Axiom-developer] Re: factor bug

[Waldek Hebisch]

William Sit wrote:
> Tim:
> 
> Could you please re-verify the NAG version?
> 
> According to your report, NAG version already has the 
> problem for factoring 119646463, which is a prime (so in 
> this case, 'factor' is correct and 'prime?' is wrong, 
> according to your report for NAG). Also, my Windows 
> version reports prime? 119646463 => true, so there is no 
> error here.
> 
> Note that in my earlier report, the number is 119643463, 
> which is not prime and is 10111 * 11833.
> 

Francois Maltey wrote:
> I test :
>
> [#factors (117661597+0*i) for i in 1..1000]   -- there are 1 and 2 factors
> reduce (+, [#factors (117661597+0*i) for i in 1..1000])     -- around 1650

The exact solution is 2000 of corse.

AFAICS the problem is that the map x -> (x^2+5) mod p has cycle
of the same length for both factors of 119643463 and 117661597. For
many examples just running Pollard method more times is sufficient.
However there are some nasty examples like 93281*97499 = 9094804219
or 90149*93229*97387 = 818489150670827 where probability of getting
right answer (using x^2+5 map) seem to be lower then 0.03.
Heuristically, using different map for succesive trials should solve
the problem, but at least at first glance it looks hard to prove
that changing map will work for _all_ numbers that we want to factor.

Below I implemented a different approach: I we detect that all factors
give cycle of the same length we look at the _first_ place when cycle
begin.  It is easy to prove that with probability at least 0.5 
cycles for different factors begin in different places. If this is the
case we can still find a proper factor.

In practise, it seems that that probability of success is much higher
than 0.5, but there is still non-negligable probability of failure, so
I also decided to repeat Pollard method up to 20 times.  I have not good
justification for choosing 20 -- for all bad examples that I was able to
produce 20 repetitions bring probability of failure way below hardware
error rate.

We probably should use better factorization method -- Pollard rho needs
a lot of time to find 17-digit factor and will probably never finish
if factors have more than 20 digits.  For comparison, I tried ecm-6.1.1
program and it found 37-digit factor of 69-digit number in about 2000
seconds, which is probably faster then Axiom Pollard for 17-digit 
factors.  But even simplest ECM is more compilcated then Pollard and
AFAICS ECM becomes really competitive only when one adds a bunch of
improvements which significantly increase complexity.  Other fast
methods look more complicated then ECM...

Anyway, the following patch fixed for me factorisation problems:

--- wh-sandbox/src/algebra/intfact.spad.pamphlet        2006-11-17 
19:30:20.000000000 +0100
+++ wh-sandbox2/src/algebra/intfact.spad.pamphlet       2006-12-20 
02:05:36.284242184 +0100
@@ -401,27 +401,51 @@
        r:I := 1
        q:I := 1
        G:I := 1
+       l:I
+       k:I
        until G > 1 repeat
           x := y
+          ys := y
           for i in 1..convert(r)@Integer repeat
              y := (y*y+5::I) rem n
              q := (q*abs(x-y)) rem n
-             k:I := 0
+          k := 0::I
+          G := gcd(q,n)
           until (k>=r) or (G>1) repeat
              ys := y
              for i in 1..convert(min(m,r-k))@Integer repeat
                 y := (y*y+5::I) rem n
-                q := q*abs(x-y) rem n
+                q := (q*abs(x-y)) rem n
              G := gcd(q,n)
              k := k+m
+          k := k + r
           r := 2*r
        if G=n then
+          l := k - m
+          G := 1::I
           until G>1 repeat
              ys := (ys*ys+5::I) rem n
              G := gcd(abs(x-ys),n)
+             l := l + 1
+          if G = n then
+             y := x0
+             x := x0
+             for i in 1..convert(l)@Integer repeat
+               y := (y*y+5::I) rem n
+             G := gcd(abs(x-y), n)
+             until G>1 repeat
+               y := (y*y+5::I) rem n
+               x := (x*x+5::I) rem n
+               G := gcd(abs(x-y), n)
        G=n => "failed"
        G
 
+    PollardSmallFactor20(n:I):Union(I,"failed") ==
+       for i in 1..20 repeat
+          r := PollardSmallFactor n
+          r case I => return r
+       r
+
     BasicSieve(r, lim) ==
        l:List(I) :=
           [1::I,2::I,2::I,4::I,2::I,4::I,2::I,4::I,6::I,2::I,6::I]
@@ -470,7 +494,7 @@
           (y:=perfectSqrt (x**2-n)) case I =>
                 insert_!(x+y,a,c)
                 insert_!(x-y,a,c)
-          (d := PollardSmallFactor n) case I =>
+          (d := PollardSmallFactor20 n) case I =>
              for m in 0.. while zero?(n rem d) repeat n := n quo d
              insert_!(d, a, m * c)
              if n > 1 then insert_!(n, a, c)


-- 
                              Waldek Hebisch
address@hidden