# My Weblog

## Elliptic curve prime factorisation using Montgomery curve

This is almost similar to previous example except using Montgomery curve rather than Lenstra curve. Programmingpraxis used Montgomery curve to find the prime factorisation. I personally feel that Lenstra algorithm worked fine when factors were small how ever Montgomery return the number as prime. Run this program against a small input (digits 15) and chances are that you will get small digit as itself prime and large input as prime factor. This code uses some bitwise operations , upper bound of number base 2 and Montgomery Ladder algorithm.

```--y^2 = x^3+ax^2+x mod n
import Random
import System.Random
import Data.Bits
data Elliptic = Cone Integer Integer Integer deriving (Eq,Show)
data Point = Conp Integer Integer  deriving (Eq,Show)

addPoints (Conp 0 0) p2 _ _ = p2
addPoints p1 (Conp 0 0) _ _ = p1
addPoints (Conp x_1 z_1) (Conp x_2 z_2) (Conp x_0 z_0) n =
let
x_3= mod ((((x_1-z_1)*(x_2+z_2)+(x_1+z_1)*(x_2-z_2))^2)*z_0) n
z_3= mod ((((x_1-z_1)*(x_2+z_2)-(x_1+z_1)*(x_2-z_2))^2)*x_0) n
in (Conp x_3 z_3)

doublePoint::Point->Elliptic->Point
doublePoint (Conp x z) (Cone an ad n) = Conp (mod x' n) (mod y' n)
where
where
t=(x+z)^2-(x-z)^2

ilognew::Integer->Integer->Int
ilognew b n = f (div ubound 2)  ubound where
ubound = until (\e -> b ^ e > n) (* 2) 1
f lo hi | lo>=hi = div (lo+hi) 2
| 2 ^ mid <= n = f (mid+1) hi
| otherwise = f lo mid
where mid = div (lo + hi) 2

montMul::Integer->Elliptic->Point->Integer->Point
montMul _ _ _ 0= (Conp 0 0)
montMul n eCurve pOint k = tmpFun  (Conp 0 0)  pOint (ilognew 2 k  -1) where
tmpFun r0 r1  0 = if (.&.) k ( shiftL (1::Integer) 0 ) /=0 then  addPoints r0 r1 pOint n
else  doublePoint r0 eCurve
tmpFun r0 r1  cnt |(.&.) k ( shiftL (1::Integer) cnt ) /=0  =
tmpFun (addPoints r0 r1 pOint n) (doublePoint r1 eCurve)  (cnt-1)
|otherwise =
tmpFun (doublePoint r0 eCurve) (addPoints r0 r1 pOint n)  (cnt-1)

montgomeryCurve::Integer-> IO [Integer]
montgomeryCurve  n = do
s<-getStdRandom (randomR (6,n-1))
let u=mod (s*s-5) n
let v=mod (4*s) n
let an=mod (((v-u)^3)*(3*u+v)) n
return [an,ad,mod (u*u*u) n,mod (v*v*v) n]

prime = 2:3:5:filter isPrime [6*k+r | k<-[1..],r<-[1,5]]
isPrime n = all (notDiv n ) \$ takeWhile (\p->p*p<=n) prime where
notDiv n p = mod n p /=0

allp::Integer->Integer->[Integer]
allp p_1 p_2 = takeWhile (<=p_2) \$ dropWhile (<p_1) prime

--http://programmingpraxis.com/2010/04/27/modern-elliptic-curve-factorization-part-2/
ellipticFac::Integer->Integer->Elliptic->Point->Integer
ellipticFac  b_1 b_2 (Cone an ad n)  (Conp u v) =
let
c= foldl lcm 1 [1..b_1]
(Conp x z)= montMul n (Cone an ad n) (Conp u v) c
g=gcd n z
in case (and [1<g, g<n]) of
True -> g
_ ->
let
(Conp x' z')= foldl ( montMul n (Cone an ad n) ) (Conp x z) (allp b_1 b_2)
g'=gcd (mod (g*z') n) n
in case (and [1<g',g'<n]) of
True -> g'
_ -> n

factor::Integer->Integer->IO [Integer]
factor 1 _ = return []
factor n cnt | cnt>=5 = return [n]
| otherwise = do
let p=ellipticFac 10000 20000 (Cone an ad n) (Conp u v)
case (and [1<p,p<n]) of
False -> factor n (cnt+1)
_ -> do
ps<-factor (div n p) (0)
ls<-factor p (0)
return (ls++ps)

main = do
line<-getLine
lst<-factor n 0
print lst
return ()
```

## Elliptic Curve prime factorisation

Lenstra algorithm explains about Lenstra elliptic curve prime factorisation and wiki page about Elliptic curve is good enough to start about elliptic curves.
Here is algorithm which i borrowed from Primality Testing and Integer Factorization in Public-Key Cryptography . Its very good book specially for those who wants to implement the different algorithm about prime numbers. This implementation is highly influenced by Basic Number Theory in Haskell. I borrowed multiEC from this library. Excellent library of number theory in Haskell.

```[1] (Choose an Elliptic Curve) Choose a random pair (E, P ),
where E is an   elliptic curve y^2 = x^3 +ax+b over Z/nZ, and
P (x, y) ∈ E(Z/nZ) is a point on E. That is, choose a, x, y ∈ Z/nZ at random, and set  b ← y^2 − x^3 − ax. If gcd(4a^3 + 27b^2 , n) = 1, then E is not an  elliptic curve, start all over and choose another pair
(E, P ).
[2] (Choose an Integer k) Just as in the “p−1” method, select a positive integer k that is divisible by many prime powers, for example, k = lcm(1, 2, . . . , B) or k = B! for a suitable bound B; the larger B is the more likely the method  will succeed in producing a factor, but the longer the method will take to work.
[3] (Calculate kP ) Calculate the point kP ∈ E(Z/nZ). We use the following formula to compute
P3 (x3 , y3 ) = P1 (x1 , y1 ) + P2 (x2 , y2 ) modulo n:
(x3 , y3 ) = (s^2 − x1 − x2) mod n, s(x1 − x3 ) − y1 mod n
where
s=3x1^2+a/2y1 mod n     if P1=P2
s=y1-y2/x1-x2 mod n

The computation of kP mod n can be done in O(log k) doublings and
[4] (Computing GCD) If kP ≡ OE (mod n), then compute d = gcd(m2 , n), where m2=2y1 or x1-x2 else go to Step [1] to make a new choice for “a” or even for a new pair  (E, P ).
[5] (Factor Found?) If 1 < d < n, then d is a nontrivial factor of n, output d  and go to Step [7].
[6] (Start Over?) If d is not a nontrivial factor of n and if you still wish to try more elliptic curves, then go to Step [1] to start all over again, else go to Step [7].
[7] (Exit) Terminate the algorithm.
```
```--y^2= x^3+ax+b mod n
import Random
data Elliptic = Conelliptic Integer Integer Integer deriving (Eq,Show)
data Point = Conpoint Integer Integer | Identity deriving (Eq,Show)

addPoints _ Identity p_2 = Left p_2
addPoints _ p_1 Identity = Left p_1
addPoints (Conelliptic a b n) (Conpoint x_p y_p) (Conpoint x_q y_q)
| x_p/=x_q = let
[u,v,d]= extended_gcd (x_p-x_q) n
s=mod ((y_p-y_q)*u) n
x_r=mod (s*s-x_p-x_q) n
y_r=mod (-y_p-s*(x_r-x_p)) n
in case ((Conpoint x_r y_r),d) of
(_,1)->Left (Conpoint x_r y_r)
(_,d')->Right d'
| otherwise = if mod (y_p+y_q) n ==0 then Left Identity
else  let
[u,v,d]= extended_gcd (2*y_p) n
s=mod ((3*x_p*x_p+a)*u) n
x_r=mod (s*s-2*x_p) n
y_r=mod (-y_p-s*(x_r-x_p)) n
in case ((Conpoint x_r y_r),d) of
(_,1)->Left (Conpoint x_r y_r)
(_,d')->Right d'

extended_gcd::Integer->Integer->[Integer]
extended_gcd a b=
let
[u,v,d]=exgcd a b
where
exgcd a b | mod a b == 0 = [0,1,b]
| otherwise = [y,x-y*(div a b),z]
where
[x,y,z] = exgcd b (mod a b)
in if d<0 then [-u,-v,-d] else [u,v,d]

multiEC::Elliptic->Point->Integer->Either Point Integer
multiEC _ _ 0 = Left Identity
multiEC ecurve point k | k>0 = fun Identity point k
where
fun p _ 0 = Left p
fun p q n = let
p'=if odd n then addPoints ecurve p q else Left p
in case (p',q') of
(Left p'',Left q'')-> fun p'' q'' (div n 2)
(Right p'',_)-> Right p''
(_,Right q'')->Right q''

dscrmntEC a b = 4 * a * a * a + 27 * b * b

randomCurve ::Integer->  IO [Integer]
randomCurve n= do
a<-getStdRandom (randomR (1,n))
u<-getStdRandom (randomR (1,n))
v<-getStdRandom (randomR (1,n))
let b=mod (v*v-u*u*u-a*u) n
return [a,b,n,u,v]

factor::Integer->Integer->Integer->IO [Integer]
factor 1 _ _ = return []
factor n m cnt | cnt>=5 = return [n]
| otherwise = do
[a,b,n,u,v]<-randomCurve n
let p=multiEC (Conelliptic a b n) (Conpoint u v ) m
case p of
Left p' -> factor n m (cnt+1)
Right p'-> do
ps<-factor (div n p') m 0 --start again factoring
ls<-factor p' m 0 --start again factoring
return (ls++ps)

main = do
line<-getLine
m<-factor n (foldl lcm 1 [1..10000]) 0
print m
return ()
```

Here is my python implementation which is faster than Haskell. Kind of surprising for me. When i started learning haskell then some where it was given that because of type checking haskell is faster than python.

```#!/usr/local/bin/python
# -*- coding: utf-8 -*-
import math
import random
#y^2=x^3+ax+b mod n

prime=[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271 ]

# ax+by=gcd(a,b). This function returns [gcd(a,b),x,y]. Source Wikipedia
def extended_gcd(a,b):
x,y,lastx,lasty=0,1,1,0
while b!=0:
q=a/b
a,b=b,a%b
x,lastx=(lastx-q*x,x)
y,lasty=(lasty-q*y,y)
if a<0:
return (-a,-lastx,-lasty)
else:
return (a,lastx,lasty)

# pick first a point P=(u,v) with random non-zero coordinates u,v        #(mod N), then pick a random non-zero A (mod N),
# then take B = u^2 - v^3 - Ax (mod N).
# http://en.wikipedia.org/wiki/Lenstra_elliptic_curve_factorization

def randomCurve(N):
A=random.randrange(N)
u=random.randrange(N)
v=random.randrange(N)
B=(v*v-u*u*u-A*u)%N
return [(A,B,N),(u,v)]

# Given the curve y^2 = x^3 + ax + b over the field K
#(whose  characteristic we assume to be neither 2 nor 3),
#and points
# P = (xP, yP) and Q = (xQ, yQ) on the curve, assume first that
# xP != xQ. Let the slope of the line s = (yP − yQ)/(xP − xQ);
#since K  is a field, s is well-defined. Then we can define
# R = P + Q = (xR, − yR) by
# s=(xP-xQ)/(yP-yQ) Mod N
#	xR=s^2-xP-xQ	Mod N
#	yR=yP+s(xR-xP)	Mod N
#  If xP = xQ, then there are two options: if yP = −yQ, including
# the case where yP = yQ = 0, then the sum is defined as
#  0[Identity]. 	 thus, the inverse of each point
# on the curve is
# found by reflecting it across the x-axis. If yP = yQ != 0,
# then R = P + P = 2P =(xR, −yR) is given by
#	s=3xP^2+a/(2yP)	Mod N
#	xR=s^2-2xP	Mod N
#	yR=yP+s(xR-xP)	Mod N
#	http://en.wikipedia.org/wiki/Elliptic_curve#The_group_law

if p_1=="Identity": return [p_2,1]
if p_2=="Identity": return [p_1,1]
a,b,n=E
(x_1,y_1)=p_1
(x_2,y_2)=p_2
x_1%=n
y_1%=n
x_2%=n
y_2%=n
if x_1 != x_2 :
d,u,v=extended_gcd(x_1-x_2,n)
s=((y_1-y_2)*u)%n
x_3=(s*s-x_1-x_2)%n
y_3=(-y_1-s*(x_3-x_1))%n
else:
if (y_1+y_2)%n==0:return ["Identity",1]
else:
d,u,v=extended_gcd(2*y_1,n)
s=((3*x_1*x_1+a)*u)%n
x_3=(s*s-2*x_1)%n
y_3=(-y_1-s*(x_3-x_1))%n

return [(x_3,y_3),d]

# http://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication
#	Q=0 [Identity element]
#	while m:
#		if (m is odd) Q+=P
#		P+=P
#		m/=2
#	return Q')

def mulPoint(E,P,m):
Ret="Identity"
d=1
while m!=0:
if d!=1 : return [Ret,d]
if d!=1 : return [Ret,d]
m>>=1
return [Ret,d]

def ellipticFactor(N,m,times=5):
for i in xrange(times):
E,P=randomCurve(N);
Q,d=mulPoint(E,P,m)
if d!=1 : return d
return N

if __name__=="__main__":
n=input()
for p in prime:#preprocessing
while n%p==0:
print p
n/=p
m=int(math.factorial(2000))
while n!=1:
k=ellipticFactor(n,m)
n/=k
print k
```

It could be possible that these codes may have indentation problem due to formatting and both languages are indentation sensitive.

January 17, 2011 Posted by | Programming | 6 Comments

## SPOJ FIBOSUM

SPOJ Fibonacci Sum is about adding fibonacci number in given range [n,m]. The Wikipedia article says “The sum of the first n Fibonacci numbers is the (n + 2)nd Fibonacci number minus 1″ [second identiy]. Rest is simple.

```#include<cstdio>
#include<iostream>
#include<cstring>
using namespace std;
typedef unsigned long long ll;
ll  MOD=ll( 1000000007);
ll Fibo(int n)
{
ll  fib[2][2]={{1,1},{1,0}},ret[2][2]={{1,0},{0,1}},tmp[2][2]={{0,0},{0,0}};
while(n)
{
if(n&1)
{
memset(tmp,0,sizeof tmp);
for(int i=0;i<2;i++) for(int j=0;j<2;j++) for(int k=0;k<2;k++)
tmp[i][j]=(tmp[i][j]+ret[i][k]*fib[k][j])%MOD;
for(int i=0;i<2;i++) for(int j=0;j<2;j++) ret[i][j]=tmp[i][j];
}
memset(tmp,0,sizeof tmp);
for(int i=0;i<2;i++) for(int j=0;j<2;j++) for(int k=0;k<2;k++)
tmp[i][j]=(tmp[i][j]+fib[i][k]*fib[k][j])%MOD;
for(int i=0;i<2;i++) for(int j=0;j<2;j++) fib[i][j]=tmp[i][j];
n/=2;

}
return (ret[0][1])%MOD;
}

int main()
{
int t,m,n;
for(scanf("%d",&t);t>0;t--)
{
scanf("%d%d",&n,&m);
cout<<(Fibo(m+2)-Fibo(n+1)+MOD)%MOD<<endl;;
}
}
```

January 9, 2011 Posted by | Programming | 3 Comments