Investigate pollard's algorithm for factorising primes between 2^32 and 2^64

[noughtmare]

I just wrote a very simple naive Pollard rho algorithm for prime factorisation (not even Brent's improvement) and it turns out to beat arithmoi's factorise function for composite numbers around 2^33, the difference in performance becomes even larger if I filter out composites whose smallest prime factor is at least 2^16. Here's my implementation and benchmark:

{-# LANGUAGE BangPatterns #-}
module Main where

import           Data.Maybe
import qualified Debug.Trace
import           Gauge.Main
import           Math.NumberTheory.Primes

-- trace = Debug.Trace.trace
trace _ x = x

pollard :: Integer -> Integer -> Integer -> Integer
pollard !n !c !x0 | isJust (isPrime n) = n
                  | otherwise          = go x0 ((x0 * x0 + c) `rem` n)
 where
  go !x !y
    | g == 1 = go ((x * x + c) `rem` n)
                  (let y' = (y * y + c) `rem` n in (y' * y' + c) `rem` n)
    | otherwise = pollard' n c x0 g x y -- get the rare case out of the way
    where g = gcd (x - y) n

{-# NOINLINE pollard' #-}
pollard'
  :: Integer -> Integer -> Integer -> Integer -> Integer -> Integer -> Integer
pollard' n c x0 g x y = if g == n
  then pollard n c (x0 + 1)
  -- pollard can return multiples of prime factors, so we fall back on trial division
  else unPrime (fst (head (factorise g)))

composite :: Integer -> Maybe Integer
composite n = safeMinimum $ map (unPrime . fst) $ factorise n

safeMinimum :: Ord a => [a] -> Maybe a
safeMinimum [] = Nothing
safeMinimum xs = Just (minimum xs)

safeHead :: [a] -> Maybe a
safeHead (x : _) = Just x
safeHead _       = Nothing

main :: IO ()
main =
  let cs =
        {- this part filters composites with large prime factors 
        mapMaybe
            (\x -> composite x >>= \c -> if c < 2 ^16
              then Nothing
              else Just x
            )
           -}[2 ^ 33 .. 2 ^ 33 + 100000]
  in  do
        print (length cs)
        print -- to make sure my pollard function is correct
          (take 10 $ mapMaybe
            (\n -> if pollard n 1 2 `elem` map (unPrime . fst) (factorise n)
              then Nothing
              else Just (pollard n 1 2, map (unPrime . fst) (factorise n), n)
            )
            cs
          )
        defaultMain
          [ bench "pollard" $ nf (map (\n -> pollard n 1 2)) cs
          , bench "arithmoi" $ nf (map (head . factorise)) cs
          ]

Results without filtering composites with large prime factors:

$ cabal run fast-factorisation -- -s
Up to date
100001
[]
benchmarking pollard ... took 12.05 s, total 56 iterations
pollard                                  mean 218.3 ms  ( +- 1.277 ms  )
benchmarking arithmoi ... took 22.64 s, total 56 iterations
arithmoi                                 mean 413.3 ms  ( +- 11.14 ms  )

Results with filtering:

$ cabal run fast-factorisation -- -s
Up to date
4620
[]
pollard                                  mean 67.10 ms  ( +- 2.051 ms  )
benchmarking arithmoi ... took 19.48 s, total 56 iterations
arithmoi                                 mean 360.2 ms  ( +- 13.59 ms  )

[Bodigrim]

Nice! I was working on an overhaul of prime numbers and factorization, but unlikely to make much progress before finishing text-utf8 project.

[noughtmare]

I totally understand, in the meantime I am doing some more experiments.