mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Miscellaneous Math

Reply
 
Thread Tools
Old 2007-07-27, 10:32   #1
rdotson
 
rdotson's Avatar
 
Jul 2005

508 Posts
Question Modification of Fermat's method

Greetings,

I believe I've discovered a way to replace one of the loops in Fermat's factorization method with a binary search. It does what was intended and seems to work, but is nothing more than a curiousity because the multiplication needed in the binary search is so much more computationally expensive than addition that the method is slower than Fermat's traditional method.

The motivation for the technique described below was a reply I received to a question on the Yahoo PrimeNumbers forum some time ago at: http://tech.groups.yahoo.com/group/p.../message/18687

The specific implementation of Fermat's method I used is described in Knuth's "Art of Computer Programming" series, in Volume 2, "Seminumerical Algorithms", Algorithm C - Fermat's method, in section 4.5.

Since I'm a programmer rather than a mathematician I would rather simply post the program itself rather than try to express it in terms of equations and botch things up. It's written in the Maple programming language.

The derivation of the predictor equation "f(n)= r-n*(y+n-1)" used in function "f2" below is another subject altogether, but is outlined in the Yahoo PrimeNumbers message mentioned above. I will be happy to provide the Maple code for the derivation if anyone is interested.

The "inner loop" that was replaced with a binary search is this code from the function "f1":
while (r > 0)
do
r:= r - y:
y:= y + 2:
od:

#================================================
# Maple Code Follows
#================================================

Code:
# Ron S. Dotson (c) July 27, 2007
restart:    # NOTE: "N" is the number to be factored
kernelopts(ASSERT=true):
#=============== Fermat Classic ================
f1:= proc(r0,x0,y0) local r, x, y, u, v, i, j;
  r:= r0:
  x:= x0:
  y:= y0:
  i:= 0:
  j:= 0:
  while (r <> 0)
    do
    while (r < 0)
      do
      r:= r + x - y:
      y:= y + 2:
      x:= x + 2:
      j:= j+1:
      od:
    while (r > 0)
      do
      r:= r - y:
      y:= y + 2:
      i:= i+1:
      od:
    od:
u:= (x-y)/2;
v:= (x+y-2)/2;
return [u,v,i,j]:
end proc:

#=============== Fermat Modification ================
f2:= proc(r0,x0,y0) local f,r,x,y,y1,rNext,rUpper,rLower,rHalfRange, rEst, yEst, u,v,i,j;
  f:= proc(n, r, y) local fn;  #p:= r-n*(y+n-1);
    fn:= r-n*(y+n-1);          #f:= unapply(p,n):
    return(fn);
  end proc:
  r:= r0:
  x:= x0:
  y:= y0:
  i:= 0:
  j:= 0:
  rLower:= 0:
  rNext:= 0:
  rUpper:= 0:
  rEst:= 0:
  while (r <> 0)
    do
    y1:= y:
    if (r < 0) then
      r:= r + x:
      x:= x + 2:
      j:= j+1:
    fi:
    #printf("% 4d: rLower= % 4d, rNext= % 6d, ", i,rLower,rNext);
    #printf("r= % 6d, y= % 6d, rUpper= % 6d, rEst= % 6d, x= % 6d\n" ,r,y,rUpper,rEst,x);
    while (r > 0)
      do
      rUpper:= trunc(r/y)+1; # first estimate is an upper bound on "r"
      rLower:= 0;
      rHalfRange:= trunc((rUpper-rLower)/2):
      rNext:= rHalfRange:
      #printf("% 4d: rLower= % 4d, rNext= % 6d, ", i,rLower,rNext);
      #printf("r= % 6d, y= % 6d, rUpper= % 6d, rEst= % 6d\n" ,r,y,rUpper,rEst);
      while (rHalfRange > 0)
        do
        i:= i+1:
        rEst:= f(rNext, r, y):
        if (rEst < 0) then
          ASSERT(rNext<rUpper, "rNext NOT < rUpper");
          rUpper:= rNext; # decrease upper bound
        else  # rEst > 0
          ASSERT(rNext>rLower, "rNext NOT > rLower");
          rLower:= rNext; # increase lower bound
        fi:
        rHalfRange:= trunc((rUpper-rLower)/2):
        rNext:= rLower + rHalfRange:
        #printf("% 4d: rLower= % 4d, rNext= % 6d, ", i,rLower,rNext);
        #printf("r= % 6d, y= % 6d, rUpper= % 6d, rEst= % 6d\n" ,r,y,rUpper,rEst);
        od:
      #printf("\n");
      if f(rLower, r, y)=0 then 
        rNext:= rLower;
      else
        rNext:= rUpper;
      fi;
      r:= f(rNext, r, y);
      y:= y1 + 2*rNext;
      #ASSERT(r=0, "r NOT = 0");
      #printf("% 4d: rLower= % 4d, rNext= % 6d, ", i,rLower,rNext);
      #printf("r= % 6d, y= % 6d, rUpper= % 6d, rEst= % 6d\n" ,r,y,rUpper,rEst);
      od:
    od:
u:= (x-y)/2;
v:= (x+y-2)/2;
return [u,v,i,j]:
end proc:

#============ Begin Program ===============
Digits:= 7:
f1Total:= 0.0:
f2Total:= 0.0:
f1Totali:=0:
f2Totali:= 0:
for i from 1 to 10
do
a:= nextprime(trunc(rand()/10^(12-Digits)));  # rand() = 12 digits
b:= nextprime(trunc(rand()/10^(12-Digits)));
N:= a*b;
printf("\n%2d Factors: a= %d, b= %d, N= %d\n", i, a, b, N);
sqrtN:= trunc(evalf[800](sqrt(N))): # floor(sqrt(N))
y:= 1:
x:= 2*sqrtN+1:   # 2 times the integer square root of N, plus 1
r:= sqrtN^2 - N: # error in using sqrtN as a factor (negative)
ASSERT(r<0,"Initialization1: r NOT < 0");
L1:= [];
L2:= [];

st:= time():
L1:= f1(r,x,y):
et:= time():

f1secs:= et-st;
f1Total:= f1Total + f1secs:
f1Totali:= f1Totali + L1[3]:
printf("    f1secs = % 5.1f secs, X-iterations= % 6d, Y-iterations= % 6d\n", f1secs,L1[4],L1[3]);

st:= time():
L2:= f2(r,x,y):
et:= time():

f2secs:= et-st;
f2Total:= f2Total + f2secs:
f2Totali:= f2Totali + L2[3]:
printf("    f2secs = % 5.1f secs, X-iterations= % 6d, Y-iterations= % 6d\n", f2secs,L2[4],L2[3]);
u:= L1[1]:
v:= L1[2]:
if (u*v <> N) then
    printf("%2d f1 FAILURE, a= %d, b= %d, N= %d\n", i,a,b,N);
fi:
x:= L2[1]:
y:= L2[2]:
if (x*y <> N) then
    printf("%2d f2 FAILURE, a= %d, b= %d, N= %d\n", i,a,b,N);
fi:
od:
printf("TOTALS: f1= % 5.0f secs (%d Y-iterations), f2=% 5.0f secs (%d Y-iterations)\n", f1Total,f1Totali, f2Total,f2Totali);
rdotson is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
Help with series convergence in Fermat method of factoring EdH Other Mathematical Topics 52 2021-01-29 21:17
B1 and B2 in P-1 method Miszka Math 13 2013-12-27 20:23
New Method Unregistered Miscellaneous Math 14 2013-05-24 10:55
Elliptic Curve Method factoring - Fermat numbers philmoore Math 131 2006-12-18 06:27
Fermat's factoring method with Gauss's inprovement Carlos Programming 0 2005-09-11 12:50

All times are UTC. The time now is 07:33.

Fri May 7 07:33:11 UTC 2021 up 29 days, 2:14, 0 users, load averages: 2.02, 2.25, 2.26

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.