-- Copyright (c) 2006-2011, David Amos. All rights reserved.

{-# LANGUAGE BangPatterns #-}

-- |A module for finding prime factors.
module Math.NumberTheory.Factor (module Math.NumberTheory.Prime,
                                 pfactors, ppfactors, pfactorsTo, ppfactorsTo) where

import Control.Arrow (second, (&&&))
import Data.Either (lefts)
import Data.List as L
import Math.Core.Utils (multisetSumAsc)
import Math.NumberTheory.Prime

-- |List the prime factors of n (with multiplicity). For example:
-- >>> pfactors 60
-- [2,2,3,5]
--
-- This says that 60 = 2 * 2 * 3 * 5
-- 
-- The algorithm uses trial division to find small factors,
-- followed if necessary by the elliptic curve method to find larger factors.
-- The running time increases with the size of the second largest prime factor of n.
-- It can find 10-digit prime factors in seconds, but can struggle with 20-digit prime factors.
pfactors :: Integer -> [Integer]
pfactors :: Integer -> [Integer]
pfactors n :: Integer
n | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> 0 = Integer -> [Integer] -> [Integer]
pfactors' Integer
n ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$ (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< 10000) [Integer]
primes
           | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< 0 = -1 Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: Integer -> [Integer] -> [Integer]
pfactors' (-Integer
n) ((Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< 10000) [Integer]
primes)
    where pfactors' :: Integer -> [Integer] -> [Integer]
pfactors' n :: Integer
n (d :: Integer
d:ds :: [Integer]
ds) | Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 1 = []
                             | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
dInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
d = [Integer
n]
                             | Integer
r Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 0 = Integer
d Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: Integer -> [Integer] -> [Integer]
pfactors' Integer
q (Integer
dInteger -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:[Integer]
ds)
                             | Bool
otherwise = Integer -> [Integer] -> [Integer]
pfactors' Integer
n [Integer]
ds
                             where (q :: Integer
q,r :: Integer
r) = Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
quotRem Integer
n Integer
d
          pfactors' n :: Integer
n [] = Integer -> [Integer]
pfactors'' Integer
n
          pfactors'' :: Integer -> [Integer]
pfactors'' n :: Integer
n = if Integer -> Bool
forall a. (Integral a, Random a) => a -> Bool
isMillerRabinPrime Integer
n then [Integer
n]
                         else let d :: Integer
d = Integer -> Integer
findFactorParallelECM Integer
n -- findFactorECM n
                              in [Integer] -> [Integer] -> [Integer]
forall a. Ord a => [a] -> [a] -> [a]
multisetSumAsc (Integer -> [Integer]
pfactors'' Integer
d) (Integer -> [Integer]
pfactors'' (Integer
n Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Integer
d))

-- |List the prime power factors of n. For example:
-- >>> ppfactors 60
-- [(2,2),(3,1),(5,1)]
--
-- This says that 60 = 2^2 * 3^1 * 5^1
ppfactors :: Integer -> [(Integer,Int)]
ppfactors :: Integer -> [(Integer, Int)]
ppfactors = ([Integer] -> (Integer, Int)) -> [[Integer]] -> [(Integer, Int)]
forall a b. (a -> b) -> [a] -> [b]
map ([Integer] -> Integer
forall a. [a] -> a
head ([Integer] -> Integer)
-> ([Integer] -> Int) -> [Integer] -> (Integer, Int)
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& [Integer] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) ([[Integer]] -> [(Integer, Int)])
-> (Integer -> [[Integer]]) -> Integer -> [(Integer, Int)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Integer] -> [[Integer]]
forall a. Eq a => [a] -> [[a]]
L.group ([Integer] -> [[Integer]])
-> (Integer -> [Integer]) -> Integer -> [[Integer]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> [Integer]
pfactors
-- ppfactors = map (\xs -> (head xs, length xs)) . L.group . pfactors

-- |Find the prime factors of all numbers up to n. Thus @pfactorsTo n@ is equivalent to @[(m, pfactors m) | m <- [1..n]]@,
-- except that the results are not returned in order. For example:
-- >>> pfactorsTo 10
-- [(8,[2,2,2]),(4,[2,2]),(6,[3,2]),(10,[5,2]),(2,[2]),(9,[3,3]),(3,[3]),(5,[5]),(7,[7]),(1,[])]
--
-- @pfactorsTo n@ is significantly faster than @map pfactors [1..n]@ for larger n.
pfactorsTo :: Integer -> [(Integer, [Integer])]
pfactorsTo n :: Integer
n = (Integer, [Integer]) -> [Integer] -> [(Integer, [Integer])]
pfactorsTo' (1,[]) [Integer]
primes where
    pfactorsTo' :: (Integer, [Integer]) -> [Integer] -> [(Integer, [Integer])]
pfactorsTo' (!Integer
m,![Integer]
qs) ps :: [Integer]
ps@(ph :: Integer
ph:pt :: [Integer]
pt) | Integer
m' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
n = [(Integer
m,[Integer]
qs)]
                                    | Bool
otherwise = (Integer, [Integer]) -> [Integer] -> [(Integer, [Integer])]
pfactorsTo' (Integer
m',Integer
phInteger -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:[Integer]
qs) [Integer]
ps [(Integer, [Integer])]
-> [(Integer, [Integer])] -> [(Integer, [Integer])]
forall a. [a] -> [a] -> [a]
++ (Integer, [Integer]) -> [Integer] -> [(Integer, [Integer])]
pfactorsTo' (Integer
m,[Integer]
qs) [Integer]
pt
        where m' :: Integer
m' = Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
ph
-- We avoid a reverse call, because it does make a noticeable difference to the speed.

-- |Find the prime power factors of all numbers up to n. Thus @ppfactorsTo n@ is equivalent to @[(m, ppfactors m) | m <- [1..n]]@,
-- except that the results are not returned in order. For example:
-- >>> ppfactorsTo 10
-- [(8,[(2,3)]),(4,[(2,2)]),(6,[(3,1),(2,1)]),(10,[(5,1),(2,1)]),(2,[(2,1)]),(9,[(3,2)]),(3,[(3,1)]),(5,[(5,1)]),(7,[(7,1)]),(1,[])]
--
-- @ppfactorsTo n@ is significantly faster than @map ppfactors [1..n]@ for larger n.
ppfactorsTo :: Integer -> [(Integer, [(Integer, Int)])]
ppfactorsTo = ((Integer, [Integer]) -> (Integer, [(Integer, Int)]))
-> [(Integer, [Integer])] -> [(Integer, [(Integer, Int)])]
forall a b. (a -> b) -> [a] -> [b]
map (([Integer] -> [(Integer, Int)])
-> (Integer, [Integer]) -> (Integer, [(Integer, Int)])
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (d, b) (d, c)
second (([Integer] -> (Integer, Int)) -> [[Integer]] -> [(Integer, Int)]
forall a b. (a -> b) -> [a] -> [b]
map ([Integer] -> Integer
forall a. [a] -> a
head ([Integer] -> Integer)
-> ([Integer] -> Int) -> [Integer] -> (Integer, Int)
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& [Integer] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) ([[Integer]] -> [(Integer, Int)])
-> ([Integer] -> [[Integer]]) -> [Integer] -> [(Integer, Int)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Integer] -> [[Integer]]
forall a. Eq a => [a] -> [[a]]
L.group)) ([(Integer, [Integer])] -> [(Integer, [(Integer, Int)])])
-> (Integer -> [(Integer, [Integer])])
-> Integer
-> [(Integer, [(Integer, Int)])]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> [(Integer, [Integer])]
pfactorsTo


-- Cohen, A Course in Computational Algebraic Number Theory, p488

-- return (u,v,d) s.t ua+vb = d, with d = gcd a b
extendedEuclid :: b -> b -> (b, b, b)
extendedEuclid a :: b
a b :: b
b
    | b
b b -> b -> Bool
forall a. Eq a => a -> a -> Bool
== 0 = (1,0,b
a)
    | Bool
otherwise = let (q :: b
q,r :: b
r) = b
a b -> b -> (b, b)
forall a. Integral a => a -> a -> (a, a)
`quotRem` b
b        -- a == qb + r
                      (s :: b
s,t :: b
t,d :: b
d) = b -> b -> (b, b, b)
extendedEuclid b
b b
r -- s*b+t*r == d
                  in (b
t,b
sb -> b -> b
forall a. Num a => a -> a -> a
-b
qb -> b -> b
forall a. Num a => a -> a -> a
*b
t,b
d)                   -- s*b+t*(a-q*b) == d

-- ELLIPTIC CURVE ARITHMETIC

data EllipticCurve = EC Integer Integer Integer deriving (EllipticCurve -> EllipticCurve -> Bool
(EllipticCurve -> EllipticCurve -> Bool)
-> (EllipticCurve -> EllipticCurve -> Bool) -> Eq EllipticCurve
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: EllipticCurve -> EllipticCurve -> Bool
$c/= :: EllipticCurve -> EllipticCurve -> Bool
== :: EllipticCurve -> EllipticCurve -> Bool
$c== :: EllipticCurve -> EllipticCurve -> Bool
Eq, Int -> EllipticCurve -> ShowS
[EllipticCurve] -> ShowS
EllipticCurve -> String
(Int -> EllipticCurve -> ShowS)
-> (EllipticCurve -> String)
-> ([EllipticCurve] -> ShowS)
-> Show EllipticCurve
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [EllipticCurve] -> ShowS
$cshowList :: [EllipticCurve] -> ShowS
show :: EllipticCurve -> String
$cshow :: EllipticCurve -> String
showsPrec :: Int -> EllipticCurve -> ShowS
$cshowsPrec :: Int -> EllipticCurve -> ShowS
Show)
-- EC p a b represents the curve y^2 == x^3+ax+b over Fp

data EllipticCurvePt = Inf | P Integer Integer deriving (EllipticCurvePt -> EllipticCurvePt -> Bool
(EllipticCurvePt -> EllipticCurvePt -> Bool)
-> (EllipticCurvePt -> EllipticCurvePt -> Bool)
-> Eq EllipticCurvePt
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: EllipticCurvePt -> EllipticCurvePt -> Bool
$c/= :: EllipticCurvePt -> EllipticCurvePt -> Bool
== :: EllipticCurvePt -> EllipticCurvePt -> Bool
$c== :: EllipticCurvePt -> EllipticCurvePt -> Bool
Eq, Int -> EllipticCurvePt -> ShowS
[EllipticCurvePt] -> ShowS
EllipticCurvePt -> String
(Int -> EllipticCurvePt -> ShowS)
-> (EllipticCurvePt -> String)
-> ([EllipticCurvePt] -> ShowS)
-> Show EllipticCurvePt
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [EllipticCurvePt] -> ShowS
$cshowList :: [EllipticCurvePt] -> ShowS
show :: EllipticCurvePt -> String
$cshow :: EllipticCurvePt -> String
showsPrec :: Int -> EllipticCurvePt -> ShowS
$cshowsPrec :: Int -> EllipticCurvePt -> ShowS
Show)
-- P x y

isEltEC :: EllipticCurve -> EllipticCurvePt -> Bool
isEltEC _ Inf = Bool
True
isEltEC (EC n :: Integer
n a :: Integer
a b :: Integer
b) (P x :: Integer
x y :: Integer
y) = (Integer
yInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
y Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
xInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
xInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
x Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
x Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
b) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 0


-- Koblitz p34

-- Addition in an elliptic curve, with bailout if the arithmetic fails (giving a potential factor of n)
ecAdd :: EllipticCurve
-> EllipticCurvePt
-> EllipticCurvePt
-> Either Integer EllipticCurvePt
ecAdd _ Inf pt :: EllipticCurvePt
pt = EllipticCurvePt -> Either Integer EllipticCurvePt
forall a b. b -> Either a b
Right EllipticCurvePt
pt
ecAdd _ pt :: EllipticCurvePt
pt Inf = EllipticCurvePt -> Either Integer EllipticCurvePt
forall a b. b -> Either a b
Right EllipticCurvePt
pt
ecAdd (EC n :: Integer
n a :: Integer
a b :: Integer
b) (P x1 :: Integer
x1 y1 :: Integer
y1) (P x2 :: Integer
x2 y2 :: Integer
y2)
    | Integer
x1 Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
x2 = let (_,v :: Integer
v,d :: Integer
d) = Integer -> Integer -> (Integer, Integer, Integer)
forall b. Integral b => b -> b -> (b, b, b)
extendedEuclid Integer
n ((Integer
x1Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
x2) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n)  -- we're expecting d == 1, v == 1/(x1-x2) (mod n)
                     m :: Integer
m = (Integer
y1Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
y2) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
v Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n
                     x3 :: Integer
x3 = (Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
x1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
x2) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n
                     y3 :: Integer
y3 = (-Integer
y1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer
x1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
x3)) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n
                  in if Integer
d Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 1 then EllipticCurvePt -> Either Integer EllipticCurvePt
forall a b. b -> Either a b
Right (Integer -> Integer -> EllipticCurvePt
P Integer
x3 Integer
y3) else Integer -> Either Integer EllipticCurvePt
forall a b. a -> Either a b
Left Integer
d
    | Integer
x1 Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
x2 = if (Integer
y1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
y2) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 0  -- includes the case y1 == y2 == 0
                 then EllipticCurvePt -> Either Integer EllipticCurvePt
forall a b. b -> Either a b
Right EllipticCurvePt
Inf
                 else let (_,v :: Integer
v,d :: Integer
d) = Integer -> Integer -> (Integer, Integer, Integer)
forall b. Integral b => b -> b -> (b, b, b)
extendedEuclid Integer
n ((2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
y1) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n)  -- we're expecting d == 1, v == 1/(2*y1) (mod n)
                          m :: Integer
m = (3Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
x1Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
x1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
a) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
v Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n
                          x3 :: Integer
x3 = (Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- 2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
x1) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n
                          y3 :: Integer
y3 = (-Integer
y1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer
x1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
x3)) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n
                      in if Integer
d Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 1 then EllipticCurvePt -> Either Integer EllipticCurvePt
forall a b. b -> Either a b
Right (Integer -> Integer -> EllipticCurvePt
P Integer
x3 Integer
y3) else Integer -> Either Integer EllipticCurvePt
forall a b. a -> Either a b
Left Integer
d
-- Note that b isn't actually used anywhere

-- Note, only the final `mod` n calls when calculating x3, y3 are necessary
-- and the code is faster if the others are removed

-- Scalar multiplication in an elliptic curve
ecSmult :: EllipticCurve
-> t -> EllipticCurvePt -> Either Integer EllipticCurvePt
ecSmult _ 0 _ = EllipticCurvePt -> Either Integer EllipticCurvePt
forall a b. b -> Either a b
Right EllipticCurvePt
Inf
ecSmult ec :: EllipticCurve
ec k :: t
k pt :: EllipticCurvePt
pt | t
k t -> t -> Bool
forall a. Ord a => a -> a -> Bool
> 0 = t
-> EllipticCurvePt
-> EllipticCurvePt
-> Either Integer EllipticCurvePt
forall t.
Integral t =>
t
-> EllipticCurvePt
-> EllipticCurvePt
-> Either Integer EllipticCurvePt
ecSmult' t
k EllipticCurvePt
pt EllipticCurvePt
Inf
    where -- ecSmult' k q p = k * q + p
          ecSmult' :: t
-> EllipticCurvePt
-> EllipticCurvePt
-> Either Integer EllipticCurvePt
ecSmult' 0 _ p :: EllipticCurvePt
p = EllipticCurvePt -> Either Integer EllipticCurvePt
forall a b. b -> Either a b
Right EllipticCurvePt
p
          ecSmult' i :: t
i q :: EllipticCurvePt
q p :: EllipticCurvePt
p = let p' :: Either Integer EllipticCurvePt
p' = if t -> Bool
forall a. Integral a => a -> Bool
odd t
i then EllipticCurve
-> EllipticCurvePt
-> EllipticCurvePt
-> Either Integer EllipticCurvePt
ecAdd EllipticCurve
ec EllipticCurvePt
p EllipticCurvePt
q else EllipticCurvePt -> Either Integer EllipticCurvePt
forall a b. b -> Either a b
Right EllipticCurvePt
p
                               q' :: Either Integer EllipticCurvePt
q' = EllipticCurve
-> EllipticCurvePt
-> EllipticCurvePt
-> Either Integer EllipticCurvePt
ecAdd EllipticCurve
ec EllipticCurvePt
q EllipticCurvePt
q
                           in case (Either Integer EllipticCurvePt
p',Either Integer EllipticCurvePt
q') of
                              (Right p'' :: EllipticCurvePt
p'', Right q'' :: EllipticCurvePt
q'') -> t
-> EllipticCurvePt
-> EllipticCurvePt
-> Either Integer EllipticCurvePt
ecSmult' (t
i t -> t -> t
forall a. Integral a => a -> a -> a
`div` 2) EllipticCurvePt
q'' EllipticCurvePt
p''
                              (Left _, _) -> Either Integer EllipticCurvePt
p'
                              (_, Left _) -> Either Integer EllipticCurvePt
q'

-- ELLIPTIC CURVE FACTORISATION

-- We choose an elliptic curve E over Zn, and a point P on the curve
-- We then try to calculate kP, for suitably chosen k
-- What we are hoping is that at some stage we will fail because we can't invert an element in Zn
-- This will lead to finding a non-trivial factor of n

discriminantEC :: a -> a -> a
discriminantEC a :: a
a b :: a
b = 4 a -> a -> a
forall a. Num a => a -> a -> a
* a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
a a -> a -> a
forall a. Num a => a -> a -> a
+ 27 a -> a -> a
forall a. Num a => a -> a -> a
* a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
b

-- perform a sequence of scalar multiplications in the elliptic curve, hoping for a bailout
ecTrial :: EllipticCurve
-> [t] -> EllipticCurvePt -> Either Integer EllipticCurvePt
ecTrial ec :: EllipticCurve
ec@(EC n :: Integer
n a :: Integer
a b :: Integer
b) ms :: [t]
ms pt :: EllipticCurvePt
pt
    | Integer
d Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 1 = [t] -> EllipticCurvePt -> Either Integer EllipticCurvePt
forall t.
Integral t =>
[t] -> EllipticCurvePt -> Either Integer EllipticCurvePt
ecTrial' [t]
ms EllipticCurvePt
pt
    | Bool
otherwise = Integer -> Either Integer EllipticCurvePt
forall a b. a -> Either a b
Left Integer
d
    where d :: Integer
d = Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
gcd Integer
n (Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
discriminantEC Integer
a Integer
b)
          ecTrial' :: [t] -> EllipticCurvePt -> Either Integer EllipticCurvePt
ecTrial' [] pt :: EllipticCurvePt
pt = EllipticCurvePt -> Either Integer EllipticCurvePt
forall a b. b -> Either a b
Right EllipticCurvePt
pt
          ecTrial' (m :: t
m:ms :: [t]
ms) pt :: EllipticCurvePt
pt = case EllipticCurve
-> t -> EllipticCurvePt -> Either Integer EllipticCurvePt
forall t.
Integral t =>
EllipticCurve
-> t -> EllipticCurvePt -> Either Integer EllipticCurvePt
ecSmult EllipticCurve
ec t
m EllipticCurvePt
pt of
                               Right pt' :: EllipticCurvePt
pt' -> [t] -> EllipticCurvePt -> Either Integer EllipticCurvePt
ecTrial' [t]
ms EllipticCurvePt
pt'
                               Left d :: Integer
d -> Integer -> Either Integer EllipticCurvePt
forall a b. a -> Either a b
Left Integer
d
-- In effect, we're calculating ecSmult ec (product ms) pt, but an m at a time

l :: a -> a
l n :: a
n = a -> a
forall a. Floating a => a -> a
exp (a -> a
forall a. Floating a => a -> a
sqrt (a -> a
forall a. Floating a => a -> a
log a
n a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
log (a -> a
forall a. Floating a => a -> a
log a
n)))
-- L(n) is some sort of measure of the average smoothness of numbers up to n
-- # [x <= n | x is L(n)^a-smooth] = n L(n)^(-1/2a+o(1))  -- Cohen p 482

-- q is the largest prime we're looking for - normally sqrt n
-- the b figure here is from Cohen p488
multipliers :: a -> [Integer]
multipliers q :: a
q = [Integer
p' | Integer
p <- (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
b) [Integer]
primes, let p' :: Integer
p' = [Integer] -> Integer
forall a. [a] -> a
last ((Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
b) (Integer -> [Integer]
forall a. Num a => a -> [a]
powers Integer
p))]
    where b :: Integer
b = a -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round ((a -> a
forall a. Floating a => a -> a
l a
q) a -> a -> a
forall a. Floating a => a -> a -> a
** (1a -> a -> a
forall a. Fractional a => a -> a -> a
/a -> a
forall a. Floating a => a -> a
sqrt 2))
          powers :: a -> [a]
powers x :: a
x = (a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a -> a -> a
forall a. Num a => a -> a -> a
*a
x) a
x

findFactorECM :: Integer -> Integer
findFactorECM n :: Integer
n | Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
gcd Integer
n 6 Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 1 =
    let ms :: [Integer]
ms = Double -> [Integer]
forall a. (RealFrac a, Floating a) => a -> [Integer]
multipliers (Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Integer -> Double
forall a. Num a => Integer -> a
fromInteger Integer
n)
    in [Integer] -> Integer
forall a. [a] -> a
head ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
filter ( (Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= 0) (Integer -> Bool) -> (Integer -> Integer) -> Integer -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n) ) ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$
       [Either Integer EllipticCurvePt] -> [Integer]
forall a b. [Either a b] -> [a]
lefts [EllipticCurve
-> [Integer] -> EllipticCurvePt -> Either Integer EllipticCurvePt
forall t.
Integral t =>
EllipticCurve
-> [t] -> EllipticCurvePt -> Either Integer EllipticCurvePt
ecTrial (Integer -> Integer -> Integer -> EllipticCurve
EC Integer
n Integer
a 1) [Integer]
ms (Integer -> Integer -> EllipticCurvePt
P 0 1) | Integer
a <- [1..] ]
    -- the filter is because d might be a multiple of n,
    -- for example if the problem was that the discriminant was divisible by n

-- TESTING MULTIPLE CURVES IN PARALLEL

-- Cohen p489
-- find inverse of as mod n in parallel, or a non-trivial factor of n
parallelInverse :: a -> [a] -> Either a [a]
parallelInverse n :: a
n as :: [a]
as = if a
d a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== 1 then [a] -> Either a [a]
forall a b. b -> Either a b
Right [a]
bs else a -> Either a [a]
forall a b. a -> Either a b
Left (a -> Either a [a]) -> a -> Either a [a]
forall a b. (a -> b) -> a -> b
$ [a] -> a
forall a. [a] -> a
head [a
d' | a
a <- [a]
as, let d' :: a
d' = a -> a -> a
forall a. Integral a => a -> a -> a
gcd a
a a
n, a
d' a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= 1]
    where c :: a
c:cs :: [a]
cs = [a] -> [a]
forall a. [a] -> [a]
reverse ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\x :: a
x y :: a
y -> a
xa -> a -> a
forall a. Num a => a -> a -> a
*a
y a -> a -> a
forall a. Integral a => a -> a -> a
`mod` a
n) 1 [a]
as
          ds :: [a]
ds = (a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\x :: a
x y :: a
y -> a
xa -> a -> a
forall a. Num a => a -> a -> a
*a
y a -> a -> a
forall a. Integral a => a -> a -> a
`mod` a
n) 1 ([a] -> [a]
forall a. [a] -> [a]
reverse [a]
as)
          (u :: a
u,_,d :: a
d) = a -> a -> (a, a, a)
forall b. Integral b => b -> b -> (b, b, b)
extendedEuclid a
c a
n
          bs :: [a]
bs = [a] -> [a]
forall a. [a] -> [a]
reverse [ a
ua -> a -> a
forall a. Num a => a -> a -> a
*a
nota a -> a -> a
forall a. Integral a => a -> a -> a
`mod` a
n | a
nota <- (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
cs [a]
ds]
-- let m = length as
-- then the above code requires O(m) mod calls - in fact 3m-3 calls (?)

parallelEcAdd :: Integer
-> [EllipticCurve]
-> [EllipticCurvePt]
-> [EllipticCurvePt]
-> Either Integer [EllipticCurvePt]
parallelEcAdd n :: Integer
n ecs :: [EllipticCurve]
ecs ps1 :: [EllipticCurvePt]
ps1 ps2 :: [EllipticCurvePt]
ps2 =
    case Integer -> [Integer] -> Either Integer [Integer]
forall a. Integral a => a -> [a] -> Either a [a]
parallelInverse Integer
n ((EllipticCurvePt -> EllipticCurvePt -> Integer)
-> [EllipticCurvePt] -> [EllipticCurvePt] -> [Integer]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith EllipticCurvePt -> EllipticCurvePt -> Integer
f [EllipticCurvePt]
ps1 [EllipticCurvePt]
ps2) of
    Right invs :: [Integer]
invs -> [EllipticCurvePt] -> Either Integer [EllipticCurvePt]
forall a b. b -> Either a b
Right [EllipticCurve
-> EllipticCurvePt -> EllipticCurvePt -> Integer -> EllipticCurvePt
g EllipticCurve
ec EllipticCurvePt
p1 EllipticCurvePt
p2 Integer
inv | (ec :: EllipticCurve
ec,p1 :: EllipticCurvePt
p1,p2 :: EllipticCurvePt
p2,inv :: Integer
inv) <- [EllipticCurve]
-> [EllipticCurvePt]
-> [EllipticCurvePt]
-> [Integer]
-> [(EllipticCurve, EllipticCurvePt, EllipticCurvePt, Integer)]
forall a b c d. [a] -> [b] -> [c] -> [d] -> [(a, b, c, d)]
L.zip4 [EllipticCurve]
ecs [EllipticCurvePt]
ps1 [EllipticCurvePt]
ps2 [Integer]
invs]
    Left d :: Integer
d -> Integer -> Either Integer [EllipticCurvePt]
forall a b. a -> Either a b
Left Integer
d
    where f :: EllipticCurvePt -> EllipticCurvePt -> Integer
f Inf pt :: EllipticCurvePt
pt = 1
          f pt :: EllipticCurvePt
pt Inf = 1
          f (P x1 :: Integer
x1 y1 :: Integer
y1) (P x2 :: Integer
x2 y2 :: Integer
y2) | Integer
x1 Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
x2 = Integer
x1Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
x2 -- slightly faster not to `mod` n here
                                | Integer
x1 Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
x2 = 2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
y1 -- slightly faster not to `mod` n here
          -- inverses = parallelInverse n $ zipWith f ps1 ps2
          g :: EllipticCurve
-> EllipticCurvePt -> EllipticCurvePt -> Integer -> EllipticCurvePt
g _ Inf pt :: EllipticCurvePt
pt _ = EllipticCurvePt
pt
          g _ pt :: EllipticCurvePt
pt Inf _ = EllipticCurvePt
pt
          g (EC n :: Integer
n a :: Integer
a b :: Integer
b) (P x1 :: Integer
x1 y1 :: Integer
y1) (P x2 :: Integer
x2 y2 :: Integer
y2) inv :: Integer
inv
              | Integer
x1 Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
x2 = let m :: Integer
m = (Integer
y1Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
y2) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
inv -- slightly faster not to `mod` n here
                               x3 :: Integer
x3 = (Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
x1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
x2) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n
                               y3 :: Integer
y3 = (-Integer
y1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer
x1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
x3)) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n
                           in Integer -> Integer -> EllipticCurvePt
P Integer
x3 Integer
y3
              | Integer
x1 Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
x2 = if (Integer
y1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
y2) Integer -> [Integer] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [0,Integer
n] -- `mod` n == 0  -- includes the case y1 == y2 == 0
                           then EllipticCurvePt
Inf
                           else let m :: Integer
m = (3Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
x1Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
x1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
a) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
inv -- slightly faster not to `mod` n here
                                    x3 :: Integer
x3 = (Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- 2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
x1) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n
                                    y3 :: Integer
y3 = (-Integer
y1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer
x1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
x3)) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n
                                 in Integer -> Integer -> EllipticCurvePt
P Integer
x3 Integer
y3

parallelEcSmult :: Integer
-> [EllipticCurve]
-> t
-> [EllipticCurvePt]
-> Either Integer [EllipticCurvePt]
parallelEcSmult _ _ 0 pts :: [EllipticCurvePt]
pts = [EllipticCurvePt] -> Either Integer [EllipticCurvePt]
forall a b. b -> Either a b
Right ([EllipticCurvePt] -> Either Integer [EllipticCurvePt])
-> [EllipticCurvePt] -> Either Integer [EllipticCurvePt]
forall a b. (a -> b) -> a -> b
$ (EllipticCurvePt -> EllipticCurvePt)
-> [EllipticCurvePt] -> [EllipticCurvePt]
forall a b. (a -> b) -> [a] -> [b]
map (EllipticCurvePt -> EllipticCurvePt -> EllipticCurvePt
forall a b. a -> b -> a
const EllipticCurvePt
Inf) [EllipticCurvePt]
pts
parallelEcSmult n :: Integer
n ecs :: [EllipticCurve]
ecs k :: t
k pts :: [EllipticCurvePt]
pts | t
k t -> t -> Bool
forall a. Ord a => a -> a -> Bool
> 0 = t
-> [EllipticCurvePt]
-> [EllipticCurvePt]
-> Either Integer [EllipticCurvePt]
forall t.
Integral t =>
t
-> [EllipticCurvePt]
-> [EllipticCurvePt]
-> Either Integer [EllipticCurvePt]
ecSmult' t
k [EllipticCurvePt]
pts ((EllipticCurvePt -> EllipticCurvePt)
-> [EllipticCurvePt] -> [EllipticCurvePt]
forall a b. (a -> b) -> [a] -> [b]
map (EllipticCurvePt -> EllipticCurvePt -> EllipticCurvePt
forall a b. a -> b -> a
const EllipticCurvePt
Inf) [EllipticCurvePt]
pts)
    where -- ecSmult' k qs ps = k * qs + ps
          ecSmult' :: t
-> [EllipticCurvePt]
-> [EllipticCurvePt]
-> Either Integer [EllipticCurvePt]
ecSmult' 0 _ ps :: [EllipticCurvePt]
ps = [EllipticCurvePt] -> Either Integer [EllipticCurvePt]
forall a b. b -> Either a b
Right [EllipticCurvePt]
ps
          ecSmult' k :: t
k qs :: [EllipticCurvePt]
qs ps :: [EllipticCurvePt]
ps = let ps' :: Either Integer [EllipticCurvePt]
ps' = if t -> Bool
forall a. Integral a => a -> Bool
odd t
k then Integer
-> [EllipticCurve]
-> [EllipticCurvePt]
-> [EllipticCurvePt]
-> Either Integer [EllipticCurvePt]
parallelEcAdd Integer
n [EllipticCurve]
ecs [EllipticCurvePt]
ps [EllipticCurvePt]
qs else [EllipticCurvePt] -> Either Integer [EllipticCurvePt]
forall a b. b -> Either a b
Right [EllipticCurvePt]
ps
                                 qs' :: Either Integer [EllipticCurvePt]
qs' = Integer
-> [EllipticCurve]
-> [EllipticCurvePt]
-> [EllipticCurvePt]
-> Either Integer [EllipticCurvePt]
parallelEcAdd Integer
n [EllipticCurve]
ecs [EllipticCurvePt]
qs [EllipticCurvePt]
qs
                             in case (Either Integer [EllipticCurvePt]
ps',Either Integer [EllipticCurvePt]
qs') of
                                (Right ps'' :: [EllipticCurvePt]
ps'', Right qs'' :: [EllipticCurvePt]
qs'') -> t
-> [EllipticCurvePt]
-> [EllipticCurvePt]
-> Either Integer [EllipticCurvePt]
ecSmult' (t
k t -> t -> t
forall a. Integral a => a -> a -> a
`div` 2) [EllipticCurvePt]
qs'' [EllipticCurvePt]
ps''
                                (Left _, _) -> Either Integer [EllipticCurvePt]
ps'
                                (_, Left _) -> Either Integer [EllipticCurvePt]
qs'

parallelEcTrial :: Integer
-> [EllipticCurve]
-> [t]
-> [EllipticCurvePt]
-> Either Integer [EllipticCurvePt]
parallelEcTrial n :: Integer
n ecs :: [EllipticCurve]
ecs ms :: [t]
ms pts :: [EllipticCurvePt]
pts
    | (Integer -> Bool) -> [Integer] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==1) [Integer]
ds = [t] -> [EllipticCurvePt] -> Either Integer [EllipticCurvePt]
forall t.
Integral t =>
[t] -> [EllipticCurvePt] -> Either Integer [EllipticCurvePt]
ecTrial' [t]
ms [EllipticCurvePt]
pts
    | Bool
otherwise = Integer -> Either Integer [EllipticCurvePt]
forall a b. a -> Either a b
Left (Integer -> Either Integer [EllipticCurvePt])
-> Integer -> Either Integer [EllipticCurvePt]
forall a b. (a -> b) -> a -> b
$ [Integer] -> Integer
forall a. [a] -> a
head ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
filter (Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/=1) [Integer]
ds
    where ds :: [Integer]
ds = [Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
gcd Integer
n (Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
discriminantEC Integer
a Integer
b) | EC n :: Integer
n a :: Integer
a b :: Integer
b <- [EllipticCurve]
ecs]
          ecTrial' :: [t] -> [EllipticCurvePt] -> Either Integer [EllipticCurvePt]
ecTrial' [] pts :: [EllipticCurvePt]
pts = [EllipticCurvePt] -> Either Integer [EllipticCurvePt]
forall a b. b -> Either a b
Right [EllipticCurvePt]
pts
          ecTrial' (m :: t
m:ms :: [t]
ms) pts :: [EllipticCurvePt]
pts = case Integer
-> [EllipticCurve]
-> t
-> [EllipticCurvePt]
-> Either Integer [EllipticCurvePt]
forall t.
Integral t =>
Integer
-> [EllipticCurve]
-> t
-> [EllipticCurvePt]
-> Either Integer [EllipticCurvePt]
parallelEcSmult Integer
n [EllipticCurve]
ecs t
m [EllipticCurvePt]
pts of
                                Right pts' :: [EllipticCurvePt]
pts' -> [t] -> [EllipticCurvePt] -> Either Integer [EllipticCurvePt]
ecTrial' [t]
ms [EllipticCurvePt]
pts'
                                Left d :: Integer
d -> Integer -> Either Integer [EllipticCurvePt]
forall a b. a -> Either a b
Left Integer
d

findFactorParallelECM :: Integer -> Integer
findFactorParallelECM n :: Integer
n | Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
gcd Integer
n 6 Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 1 =
    let ms :: [Integer]
ms = Double -> [Integer]
forall a. (RealFrac a, Floating a) => a -> [Integer]
multipliers (Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Integer -> Double
forall a. Num a => Integer -> a
fromInteger Integer
n)
    in [Integer] -> Integer
forall a. [a] -> a
head ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ (Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
filter ( (Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= 0) (Integer -> Bool) -> (Integer -> Integer) -> Integer -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
n) ) ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$
       [Either Integer [EllipticCurvePt]] -> [Integer]
forall a b. [Either a b] -> [a]
lefts [Integer
-> [EllipticCurve]
-> [Integer]
-> [EllipticCurvePt]
-> Either Integer [EllipticCurvePt]
forall t.
Integral t =>
Integer
-> [EllipticCurve]
-> [t]
-> [EllipticCurvePt]
-> Either Integer [EllipticCurvePt]
parallelEcTrial Integer
n [Integer -> Integer -> Integer -> EllipticCurve
EC Integer
n (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
i) 1 | Integer
i <- [1..100]] [Integer]
ms (Int -> EllipticCurvePt -> [EllipticCurvePt]
forall a. Int -> a -> [a]
replicate 100 (Integer -> Integer -> EllipticCurvePt
P 0 1)) | Integer
a <- [0,100..] ]
-- 100 at a time is chosen heuristically.