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::Point->Point->Point->Integer->Point
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 
					   x'= mod (4*ad*(x+z)^2*(x-z)^2) n
					   y'= mod ((4*ad*(x-z)^2+t*an)*t) n 
					    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
			let ad=mod (4*u^3*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 
				[an,ad,u,v] <-montgomeryCurve n
				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 
		let n=read line ::Integer
		lst<-factor n 0
		print lst
		return ()