Algorithms for trial division

Modular arithmetics

To find Fm=0 (mod k.2n+1), there is a simple base algorithm:
Let x=2, then it is necessary to execute a maximum of n-2 iterations:
x=x2 mod (k.2n+1)
Certainly it is possible to slightly optimize calculations, if first time to accept x=2j, where j is defined depending on n.
To implement test with the necessary efficiency, we need:

Let x2=A.2n+B and let A=q.k+r where 0<=r<k; we can determine q and r rapidly, when k is a single-precision number.
Then x2=B-q+r.2n (mod k.2n+1), so the remainder is easily obtained. (The Art Of Computer Programming, vol 2, Donald E. Knuth).

Fast algorithm for N=k.2n+1 < 95 bit

We start calculating the algorithm:
kbit=log2k, and Z=kbit+n+96
R=2Z/N (accuracy of calculation R=96bit)
Basic algorithm of calculation x=x2 mod N
Sq=x2 (6 multiplication)
A=trunc(Sq*R) (6 multiplication, high bits Sq)
B=A*N/2Z (6 multiplication)
x=Sq-B (low bits Sq)

Erastothenes approach

Certainly performance of the basic algorithm is rather slowly, but it is notmandatory to make calculations for all k. It is possible to exclude a lot of odd k by Erastosthenes method. Let p be a prime number. We know that if p divides k.2n+1 then p also divides (k+j.p).2n+1, where 0 <= j < infinity. At the start of the program, do some setting up, for example to find that 3 divides k=2(mod 3), 5 divides k=1(mod 5), 7 divides k=4(mod 7), etc.
Then use a fast implementation of the sieve of Erastosthenes to quickly mark these values as composite:
Clear all bits in an array
3 divides k=2 (mod 3), so set bits 2, 5, 8, etc.
5 divides k=1 (mod 5), so set bytes 1, 6, 11, etc
7 divides k=4 (mod 7), set set bytes 4, 11, 18, etc.
etc.

The fastest algorithm for finding value k0 , that satisfies the condition k0. 2n+1=0 (mod p), has been recently offered by Paul Jobling. He has written:
There is an approach that is independent of the size of p and has a speed that is O(log n), so the speed is fixed regardless of the size of p

p | k0 .2n+1 ==> - k0 = 2-n (mod p)

2-n = (1/2)n, and 1/2 (mod p) = (p+1)/2.

So just do a left-to-right exponentiation (mod p) of (p+1)/2 to the power n.

For example, if n=13 and p=17:
(p+1)/2 = 9, and 9^13 (mod 17) can be found by left-to-right exponentiation to be 8. So if k0 = - 8(mod 17), then 17 | k0 .2^13+1. And indeed 17 does divide (-8+17).2^13+1, (9+17).2^13+1, etc.

ECM

What is an elliptic curve? "Most" second degree equations in two variables represent conic sections, where by "most", we mean that we want to exclude certain degenerate cases, such as an equation like y^2 - x^2 = 0 which actually consists of two straight lines. "Most" third or fourth degree equations in two variables can be equipped with a structure which turns them into elliptic curves. (We can exclude bad cases by throwing out the ones where the discriminant of the equation is equal to zero.) In the good cases, we can come up with a way of "adding" points on the curve, that is combining any two points on the curve to get another point on the curve in a way that the combination is associative, there is an identity element, and each element has an inverse. Talking about elliptic curves is simplified by the fact that we can change variables in a way that transforms the equation into what is known as the Weierstrass form: y^2 = x^3 + A*x^2 + B*x + C, where A, B, and C are constants. (Some variations set either A or C to zero.) The points on this elliptic curve are the solutions of this equation along with one more point called the "point at infinity", denoted O. We can pick this point at infinity to be the identity element. Then, the way that we add two points (x1,y1), (x2,y2) on this curve is: Step 1: Find the equation of the line which passes through these two points. If these two points are the same point, take the tangent line to the curve instead.

Step 2: This straight line intersects the equation y^2 = x^3 + A*x^2 + B*x + C in three points unless it is a vertical line. Two of these points are (x1,y1) and (x2,y2). Call the third point (x3,y3). Step 3: Define the "group addition law sum" of (x1,y1) and (x2,y2) to be the point (x3,-y3). Note the minus sign in front of y3. Since y^2 = x^3 + A*x^2 + B*x + C is even in y, this point is also on the curve. If the straight line is vertical, define the sum to be the point O at infinity. This will be the case only when (x2,y2) = (x1,-y1). (We also need to define the sum of two points when one of them is the point at infinity, but since O is the identity, (x,y)+O=(x,y) and O+O=O finishes the definition.)