Posts Tagged ‘factorization’

Programming Praxis – Modern Elliptic Curve Factorization, Part 1

April 23, 2010

Today marks the 1-year anniversary of this blog. What better way to celebrate than to do another programming exercise?

In today’s Programming Praxis exercise we have to implement an algorithm for curve factorization. The Scheme solution is 30 lines. Let’s see if we can do any better.

Since the formulas for calculating the x and z coordinates of the addition result are nearly identical, we can factor it out into a function.

add :: Integral a => (a, a) -> (a, a) -> (a, a) -> a -> (a, a)
add (x1, z1) (x2, z2) (x0, z0) n = (f (+) z0, f (-) x0) where
    f g m = (g ((x1-z1)*(x2+z2)) ((x1+z1)*(x2-z2)))^2 * m `mod` n

For doubling, we could factor out the first three terms of the two multiplications, but it would barely save any space.

double :: Integral a => (a, a) -> a -> a -> a -> (a, a)
double (x,z) an ad n = (mod x' n, mod z' n) where
    x' = 4 * ad * (x-z)^2 * (x+z)^2
    z' = (4 * ad * (x-z)^2 + t * an) * t
    t = (x+z)^2 - (x-z)^2

Multiplication can be done through some simple recursion.

multiply :: Integral a => Int -> (a, a) -> a -> a -> a -> (a, a)
multiply 0 _ _  _  _ = (0,0)
multiply 1 p _  _  _ = p
multiply k p an ad n = let half = multiply (div k 2) p an ad n in
    if odd k then add half (multiply (div k 2 + 1) p an ad n) p n
             else double half an ad n

Calculating the curve is just executable math.

curve12 :: Integral a => a -> a -> (a, a, (a, a))
curve12 n s = ((v-u)^3 * (3 * u + v) `mod` n,
               4 * u^3 * v `mod` n,
               (u^3 `mod` n, v^3 `mod` n)) where
              u = s * s - 5 `mod` n
              v = 4 * s `mod` n

And, as usual, a test to see if everything is working correctly.

main :: IO ()
main = do let (n, s) = (5569379, 4921134)
          let (an, ad, p) = curve12 n s
          let p2 = double p an ad n
          let p3 = add p2 p p n
          let p4 = double p2 an ad n
          let p6 = double p3 an ad n
          let p7 = add p4 p3 p n
          let p13 = add p7 p6 p n
          print $ p2 == (3539269, 4477480)
          print $ p3 == (2517168, 4993956)
          print $ p4 == (4683984, 2280932)
          print $ p6 == (3440206, 1480034)
          print $ p7 == (2544042, 2445346)
          print $ p13 == (577485, 4434222)
          print $ multiply 13 p an ad n == p13

That leaves us with 16 lines, just over half the size of the Scheme solution, which is better than I anticipated since math problems tend to be hard to achieve major savings on.

Programming Praxis – Extending Pollard’s P-1 Factorization Algorithm

March 19, 2010

In today’s Programming Praxis exercise we need to write an improved version of a factorization algorithm. I was on vacation when the original exercise was posted, so let’s see what we can do with it.

As usual, some imports:

import Data.Bits
import Data.List

We need the same expm function we have used in several previous exercises. Alternatively we could use the expmod function from Codec.Encryption.RSA.NumberTheory, but it’s a lot slower than this version.

expm :: Integer -> Integer -> Integer -> Integer
expm b e m = foldl' (\r (b', _) -> mod (r * b') m) 1 .
             filter (flip testBit 0 . snd) .
             zip (iterate (flip mod m . (^ 2)) b) .
             takeWhile (> 0) $ iterate (`shiftR` 1) e

The scheme solution has some duplication in it: the start of the pollard1 and pollard2 functions is nearly identical. Since programmers hate repeating themselves, let’s factor that out into a separate function.

pollard :: (Integer -> t) -> (Integer -> t) -> Integer -> Integer -> t
pollard found notFound n b1 = f 2 2 where
    f a i | i < b1         = f (expm a i n) (i + 1)
          | 1 < d && d < n = found d
          | otherwise      = notFound a
          where d = gcd (a - 1) n

pollard1 then becomes very simple: if we don’t find anything we stop, otherwise we return the result.

pollard1 :: Integer -> Integer -> Maybe Integer
pollard1 = pollard Just (const Nothing)

pollard2 is a bit more involved, because we now have an extra step if we don’t find anything. The structure of this is very similar to the pollard function, but there are enough differences that it’s not worth the bother of abstracting it.

pollard2 :: Integer -> Integer -> Integer -> Maybe (String, Integer)
pollard2 n b1 b2 = pollard (Just . (,) "stage1") (f b1) n b1 where
    f j a | j == b2        = Nothing
          | 1 < d && d < n = Just ("stage2", d)
          | otherwise      = f (j + 1) a
          where d = gcd (expm a j n - 1) n

And of course the test to see if everything’s working correctly:

main :: IO ()
main = do print $ pollard1 15770708441 150
          print $ pollard1 15770708441 180
          print $ pollard2 15770708441 150 180

Programming Praxis – Monte Carlo factorization

June 19, 2009

In today’s Programming Praxis problem we have to implement John Pollard’s factorization algorithm. Our target is 16 lines (I’m not counting the code for the primality test, since we did that already).

First, we’re going to need to reuse the code from the Miller-Rabin primality test exercise, since we need to determine whether or not the number is prime:

import Control.Arrow
import Data.Bits
import Data.List
import System.Random

isPrime :: Integer -> StdGen -> Bool
isPrime n g =
    let (s, d) = (length *** head) . span even $ iterate (`div` 2) (n-1)
        xs = map (expm n d) . take 50 $ randomRs (2, n - 2) g
    in all (\x -> elem x [1, n - 1] ||
                  any (== n-1) (take s $ iterate (expm n 2) x)) xs

expm :: Integer -> Integer -> Integer -> Integer
expm m e b = foldl' (\r (b', _) -> mod (r * b') m) 1 .
             filter (flip testBit 0 . snd) .
             zip (iterate (flip mod m . (^ 2)) b) $
             takeWhile (> 0) $ iterate (`shiftR` 1) e

There’s not much more to the factor function than simply writing down the algorithm. It’s just math with Haskell syntax, really. The only thing you have to take care of is not to calculate gcd (x-y) n when x and y are still 2, since that will give you an incorrect result.

factor :: Integer -> Integer -> Integer
factor c n = factor' 2 2 1 where
    f x = mod (x * x + c) n
    factor' x y 1 = factor' x' y' (gcd (x' - y') n) where
                        (x', y') = (f x, f $ f y)
    factor' _ _ d = if d == n then factor (c + 1) n else d

And factors does little more than recursively calling factor, while filtering out factors of two and making sure not to call factor on a prime (since the algorithm is not designed for that).

factors :: Integer -> StdGen -> [Integer]
factors n g = sort $ fs n where
    fs x | even x      = 2 : fs (div x 2)
         | isPrime x g = [x]
         | otherwise   = f : fs (div x f) where f = factor 1 x

And to test:

main :: IO ()
main = print . factors (2^98 - 1) =<< getStdGen

The end result is 9 lines of code, which was to be expected, given that it’s just a question of writing math in your language’s syntax.

Programming Praxis – Fermat’s Method

May 19, 2009

Today’s Programming Praxis problem brings us yet another factorization problem, this time using Fermat’s method. Originally I tried using the algorithm from the assignment, but that resulted in pretty much a direct translation of his Scheme solution, since I couldn’t find a way to express it that was both more elegant and not slower. Since that would be boring, we’re going to use a slightly different approach to Fermat’s method that is slightly shorter and faster.

First we need a way to take the square root of integers:

isqrt :: Integer -> Integer
isqrt = ceiling . sqrt . fromIntegral

What this algorithm does is check whether x2 – n is square for any xs between sqrt(n) and n. If so, we return the factors of n (which are further factorized if needed).

factor :: Integer -> [Integer] -> [Integer]
factor _ []     = []
factor n (x:xs) | y*y /= x*x - n = factor n xs
                | x - y == 1     = [x + y]
                | otherwise      = concatMap factors [x - y, x + y]
                where y = isqrt (x*x - n)

Like the Scheme solution, we divide by 2 until we have an odd number to work with, since Fermat’s method only works on odd numbers.

factors :: Integer -> [Integer]
factors n | even n    = 2 : factors (div n 2)
          | otherwise = factor n [isqrt n..n]

And, as usual, we test our code:

main :: IO ()
main = print $ factors 600851475143

Nine lines of code. Not bad.

Programming Praxis – Wheel Factorization

May 8, 2009

In today’s Programming Praxis problem we’re supposed to factor numbers both the naive way and using Wheel factorization. Let’s go.

First, let’s define a function that finds factors using trial division:

tdFactors :: Integer -> [Integer] -> [Integer]
tdFactors _ [] = []
tdFactors r (x:xs) | x * x > r    = [r]
                   | mod r x == 0 = x : tdFactors (div r x) (x:xs)
                   | otherwise    = tdFactors r xs

Using that, the function for naive factorization is trivial:

factorNaive :: Integer -> [Integer]
factorNaive n = tdFactors n [2..]

Now for the wheel factorization. First we define a function that tells us which columns of our wheel to use given the circumference.

spokes :: Integer -> [Bool]
spokes c = map ((== 1) . gcd c) [1..c]

To construct our wheel, we keep repeating this pattern and use only the numbers from the correct columns (and correct the start of the wheel).

wheel :: [Integer] -> [Integer]
wheel ps = (ps ++) . drop 1 . map snd . filter fst $
           zip (cycle . spokes $ product ps) [1..]

The function to do wheel factorization is practically identical to the naive one:

factorWheel :: Integer -> [Integer]
factorWheel n = tdFactors n $ wheel [2,3,5,7]

And of course we need to test if everything works:

main :: IO ()
main = do print $ factorNaive 600851475143
          print $ factorWheel 600851475143

And we’re done. Pretty simple.