func-with-traceShow.hs

-- | @'sieveFactor' fs n@ finds the prime factorisation of @n@ using the 'FactorSieve' @fs@.
--   For negative @n@, a factor of @-1@ is included with multiplicity @1@.
--   After stripping any present factors @2@, the remaining cofactor @c@ (if larger
--   than @1@) is factorised with @fs@. This is most efficient of course if @c@ does not
--   exceed the bound with which @fs@ was constructed. If it does, trial division is performed
--   until either the cofactor falls below the bound or the sieve is exhausted. In the latter
--   case, the elliptic curve method is used to finish the factorisation.
sieveFactor :: FactorSieve -> Integer -> [(Integer,Int)]
sieveFactor (FS bnd sve) = check
  where
    bound = fromIntegral bnd
    check 0 = error "0 has no prime factorisation"
    check 1 = []
    check n
      | n < 0       = (-1,1) : check (-n)
      | n <= bound  = go2w (fromIntegral n)     -- avoid expensive Integer ops if possible
      | fromInteger n .&. (1 :: Int) == 1 = sieveLoop n
      | otherwise   = go2 n
    go2w n
        | n .&. 1 == 1 = intLoop ((n-3) `shiftR` 1)
        | otherwise = case shiftToOddCount n of
                        (k,m) -> (2,k) : if m == 1 then [] else intLoop ((m-3) `shiftR` 1)
    go2 n = case shiftToOddCount n of
              (k,m) -> (2,k) : if m == 1 then [] else sieveLoop m
    sieveLoop n
        | bound < n  = tdLoop n (integerSquareRoot' n) 0
        | otherwise = intLoop (fromIntegral (n `shiftR` 1)-1)
    intLoop :: Word -> [(Integer,Int)]
    intLoop !n = case unsafeAt sve (fromIntegral n) of
                  0 -> [(2*fromIntegral n+3,1)]
                  p -> let p' = fromIntegral p in countLoop p' (n `quot` p' - 1) 1
    countLoop !p !i !c
        = case unsafeAt sve (fromIntegral i) of
            0 | p-3 == 2*i -> [(fromIntegral p,c+1)]
              | otherwise  -> (fromIntegral p,c) : (2*fromIntegral i+3,1) : []
            q | fromIntegral q == p -> countLoop p (i `quot` p - 1) (c+1)
              | otherwise -> (fromIntegral p, c) : intLoop i
    lstIdx = snd (bounds sve)
    tdLoop n sr ix
        | lstIdx < ix   = curve n
        | sr < p        = trace "SR " [(n,1)]
        | pix /= 0      = tdLoop n sr (ix+1)    -- not a prime
        | otherwise     = case splitOff p n of
                            (0,_) -> tdLoop n sr (ix+1)
                            (k,m) -> (p,k) : case m of
                                               1 -> []
                                               j | j <= bound -> intLoop (fromIntegral (j `shiftR` 1) - 1)
                                                 | otherwise -> tdLoop j (integerSquareRoot' j) (ix+1)
          where
            -- traced versions of new code
--            p = traceShow (pix, ix, 2*ix+3) $ fromIntegral $ 2*ix + 3
--            pix = unsafeAt sve ix
            -- traced versions of old code
            p = traceShow (unsafeAt sve $ toPrim ix, ix, toPrim ix) $ toPrim ix
            pix = unsafeAt sve $ fromIntegral p
    curve n = trace ("CURVE " ++ show n) $ stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n

test-case.hs

import Math.NumberTheory.Primes.Factorisation
main = do
	putStrLn $ ("RETVAL "++) $ show $ sieveFactor (factorSieve 15) 75
	putStrLn $ ("RETVAL "++) $ show $ sieveFactor (factorSieve 2048) (29*73)

{-
Output from ORIGINAL.
* The 3-tuples are (pix, ix, p).
* SR from the first test case indicates that 75 was considered prime because its square root was reached.
* Because CURVE does not appear, neither test case began the (unnecessary) curve factorization step.

(0,0,7)
(0,1,11)
SR
RETVAL [(75,1)]
(0,0,7)
(5,1,11)
(0,2,13)
(0,3,17)
(0,4,19)
(7,5,23)
(0,6,29)
RETVAL [(29,1),(73,1)]


======

Output from PATCHED
* The 3-tuples are (pix, ix, p).
* Because CURVE does not appear, neither test case began the (unnecessary) curve factorization step.

(0,0,3)
(0,1,5)
RETVAL [(3,1),(5,2)]
(0,0,3)
(0,1,5)
(0,2,7)
(3,3,9)
(0,4,11)
(0,5,13)
(3,6,15)
(0,7,17)
(0,8,19)
(3,9,21)
(0,10,23)
(5,11,25)
(3,12,27)
(0,13,29)
RETVAL [(29,1),(73,1)]

-}