Exception: ghc-6.2.1: internal error: scavenge_mutable_list: strange object? 55456

[Michael Spencer Leeming]

Hi everyone,

First let me say thanks for maintaining my favourite programming language. 
I managed to crash ghci-6.2.1 for the first time tonight, with the 
attached module, and the following command:

qsFactor 1000 10000 ((nextPrime 4938203842901)*(nextPrime 4389273982791)) 4

*** Exception: ghc-6.2.1: internal error: scavenge_mutable_list: strange 
object? 55456
     Please report this as a bug to glasgow-haskell-bugs at haskell.org,
     or http://www.sourceforge.net/projects/ghc/

Cheers,

     Michael
-------------- next part --------------
module Numtheory where

import List
import Maybe
import GHC.Real
import GHC.Arr
import Bits
import Array

modA :: Integral a => a -> a -> a
modA a n
  | r < (div n 2) = r
  | otherwise = r-n
  where
    r = a `mod` n

mulPair (a,b) (c,d) = (a*c,b*d)

merge [] ys = ys
merge xs [] = xs
merge (x:xs) (y:ys)
  | x <= y = x:(merge xs (y:ys))
  | otherwise = y:(merge (x:xs) ys)

plusMod2 x y
  | even x == even y = 0
  | otherwise = 1

iSqrt n = head $ dropWhile (\x -> x*x /= n) (iterate approx n)
  where
    approx x = (x + (n `div` x)) `div` 2

-- returns the coefficients (a,b) such that a*x + b*y = gcd x y
egcd :: Integer -> Integer -> (Integer, Integer)
egcd x y
   | y == 0 = (signum x, 0)
   | otherwise = (b, a - b*(div x y))
   where
      (a, b) = egcd y (x `mod` y)

-- inverse: returns an _invertible_ k such that kx `mod` n == gcd x n 
inverse :: Integer -> Integer -> Integer
inverse x n = head results
   where
      result = fst (egcd x n)
      addNum = div n (gcd x n)
      results = [(result + k*addNum) `mod` n | k <- [0 .. ((gcd x n) -1)], 
         gcd (result + k*addNum) n == 1] 

jacobi a b
  | gcd a b > 1 = 0
  | a > b = jacobi (a `mod` b) b
  | a == 1 = 1
  | even a = let r = jacobi (div a 2) b in case (b `mod` 8) of {1 -> r; 3 -> (-r); 5 -> (-r); 7 -> r}
  | a `mod` 4 == 3 && b `mod` 4 == 3 = -k
  | odd b = k
    where k = jacobi (b `mod` a) a

modExp :: Integer -> Integer -> Integer -> Integer
modExp _ _ 0 = 1
modExp n x a 
   | even a = (result * result) `mod` n
   | otherwise = (x * result * result) `mod` n
   where result = modExp n x (div a 2)

isPrime :: Integer -> Bool
isPrime n 
   | n < 30 = n `elem` somePrimes
   | (any (==0).map (mod n)) somePrimes = False
   | otherwise = n < 960 || all isA_PSP somePrimes
   where
      somePrimes = [2,3,5,7,11,13,17,19,23,29]
      (m, odd) = twoFreeForm (n-1)
      isA_PSP a = result == 1 || any (==(n-1)) (take m (iterate sqr result))
         where 
            result = modExp n a odd
            sqr x = (x*x) `mod` n

twoFreeForm n = until (odd.snd) next (0, n)
  where next (a, s) = (a + 1, div s 2)

nextPrime :: Integer -> Integer
nextPrime n = (head.filter isPrime.enumFrom) (n+1)

primesFrom n = tail $ iterate nextPrime n
   
pollardRhoFactor :: Integer -> Integer
pollardRhoFactor n
   | even n = 2
   | otherwise = (head.filter (/= 1).map (gcd n).concat.map numbers) [5..]
   where
      q c x = ((x*x) - c) `mod` n
      numbers c = (takeWhile (/= 0).tail) (zipWith (-) 
         (iterate (q c.q c) 2) (iterate (q c) 2))

fullyRhoFactor n
  | n < 0 = (-1):(fullyRhoFactorH (-n))
  | otherwise = fullyRhoFactorH n

fullyRhoFactorH :: Integer -> [Integer]
fullyRhoFactorH n
   | n == 1 = []
   | isPrime n = [n]
   | otherwise = merge (fullyRhoFactorH a) (fullyRhoFactorH b)
   where
      a = pollardRhoFactor n
      b = div n a

repeatFromCF :: [Integer] -> [Rational]
repeatFromCF (x:xs) = tail cFracs
  where
    cFracs = infinity:(fromInteger x):(zipWith3 addMF xs (tail cFracs) cFracs)
    addMF :: Integer -> Rational -> Rational -> Rational
    addMF m p q = (numerator p + m*numerator q) % (denominator p + m*denominator q)


cfRelations n = map toPair ((repeatFromCF.cfSqrt) n)
  where
    toPair r = let a = numerator r in (a, (a*a) `modA` n)

cfSqrt :: Integer -> [Integer]
cfSqrt nI = cf (0, 1)
  where
    cf (p,q) = let k = estSurd (p,q) in k: (cf $ invertSurd (p - fromInteger k, q))
    n = fromInteger nI
    invertSurd :: (Rational, Rational) -> (Rational, Rational)
    invertSurd (p,q) = let d = p*p - q*q*n in (p/d, -q/d)
    gtZero (p,q) = signum p == signum q || p*p < q*q*n
    estSurd (p,q) = estimate (\x -> gtZero (p-fromInteger x, q))

estimate :: (Integer -> Bool) -> Integer
estimate f = foldr maybeAdd 0 vs 
  where
    vs = takeWhile f (map bit [0 ..])
    maybeAdd x v
      | f (x+v) = x+v
      | otherwise = v    

powersOf2 = iterate (2*) 1

primes = 2:[k | k <- [3,5..], all ((/=0).(mod k)) (takeWhile ((<=k).(^2)) primes)] 

getFB size integer = result
  where
    result = take size $ filter (\r -> (isFine.abs) (snd r)) (cfRelations integer) 
    z = product $ take size primes
    isFine n
      | d == n = False
      | d == 1 = True
      | otherwise = isFine d
        where d = div n (gcd n z)

bitGaussElim :: Bits a => Int -> [a] -> [a] 
bitGaussElim n vs
  | vs == [] = []
  | ones == [] = bitGaussElim (n+1) vs
  | otherwise = firstOne:reduceRest
  where
    reduceRest = bitGaussElim (n+1) ((map (xor firstOne) restOnes)++zeroes)
    (ones,zeroes) = partition ((flip testBit) n) vs
    (firstOne:restOnes) = ones

cfFactor2 size n = factorFromRelations n sPrimes fb
  where
    fb = getFB (size+1) n
    sPrimes = (-1):(take size primes)

factorFromRelations :: Integer -> [Integer] -> [(Integer, Integer)] -> Integer
factorFromRelations n primes pairs = result
  where
    l_pairs = length pairs
    l_primes = length primes
    pairArray = listArray (l_primes + 1, l_primes + l_pairs) pairs
    irs = filter (\v -> not $ any (testBit v) [1 .. l_primes]) 
      (bitGaussElim 0 $ getMatrix primes pairs) 
    rs = [[pairArray!i | i <- [(l_primes+1) .. (l_primes + l_pairs)], 
      testBit v i] | v <- irs]
    result = head [u | (x,y) <- map (foldl1 mulPair) rs, 
      u <- [gcd n (x - iSqrt y)], u /= 1, u /= n]

getMatrix primes pairs = vs 
  where
    pps = zip primes [0,1..]
    indexNumbers :: [(Integer, Int)] -> [Integer] -> Integer -> Integer
    indexNumbers list@((p,i):ps) (f:fs) v
      | f > p = indexNumbers ps (f:fs) v
      | f == p = indexNumbers list fs (complementBit v i)
    indexNumbers _ _ v = v
    vs = zipWith (indexNumbers pps) (map (fullyRhoFactor.snd) pairs) 
      (map bit [(1 + length primes) ..])

modSqrtP2 p n = (k*p+b) `mod` (p*p)
  where 
    b = modSqrt p (n `mod` p)
    k = ((n - b*b) `div` p) * ((p+1) `div` 2) * (inverse b p)

modSqrt p n
  | n >= p || n < 0 = modSqrt p (n `mod` p)
  | jacobi n p == 1 = result
    where
      (s, q) = twoFreeForm (p-1)
      w = head $ [k | k <- [5,6..], jacobi k p == -1]
      v = modExp p w q
      nI = inverse n p
      r1 = modExp p n ((q+1) `div` 2)
      sqr x = (x*x) `mod` p
      i r = length . takeWhile (/= 1) $ iterate sqr (r*r*nI)
      nextR r = (r*modExp p v (2^(s-(i r)-1))) `mod` p
      result = until (\r -> sqr r == n) nextR r1

randoms :: Integer -> [[Int]]
randoms seed = mashUp
  where
      extend list = list ++ zipWith (+) (drop 55 (extend list)) (drop 24 (extend list))
      mash x y = ((rotate x 7) `xor` y)*37 + 55
      mashUp = tail $ map (tail . extend) $ iterate (scanl1 mash) originalList
      originalList = 1:[fromInteger (mod seed p) | p <- take 55 primes]

qsFactor size m n k = error $ show relations
  where
    fbPrimes = (-1):take size [p | p <- primes, jacobi n p == 1]
    pF = estimate (\x -> (((x^k)*m)^2) < 2*n)
    aLists = choose k $ take (2*k) [p | p <- primesFrom pF, jacobi n p == 1]
    relations = take (size+5) $ concatMap (siqs fbPrimes n m) aLists

siqs fb n m aFactors = concatMap bToRs bs
  where
    fbPrimes = tail fb
    sqrtsN = map ((flip modSqrt) n) fbPrimes
    aIs = map (inverse a) fbPrimes
    a = (product aFactors)^2
    sqrts = [let r = modSqrtP2 p n in [(r, p*p), (p*p-r, p*p)] | p <- aFactors]
    bs = filter (< (div a 2)) $ map (fst . (foldl1 crtCombine)) (listProduct sqrts)
    bToC b 
      | mod (b*b - n) a == 0 = div (b*b - n) a    
    applyPoly b x = (a*x*x) + (2*b*x) + (bToC b)
    intervalM b = listArray (-m,m) (map (applyPoly b) [(-m) .. m])
    factorPairList b = [(fromInteger r, p) | (p, sqrtN, aI) <- zip3 fbPrimes sqrtsN aIs,
      r <- getSolutions b p sqrtN aI]
    getSolutions b p sqrtN aI = concat [takeWhile ((<=m).abs) (iterate (d+) v) |
      v <- [(aI*(-b + sqrtN)) `mod` p, (-aI*(b * sqrtN)) `mod` p], d <- [p,-p]]
    bToRs b = let arry = unsafeAccum removeFs (intervalM b) (factorPairList b)
      in [((a*x+b)^2, a*applyPoly b x) | x <- [(-m) .. m], arry!x == 1]  

removeFs p n
  | n `mod` p == 0 = removeFs p (div n p)
  | otherwise = n
  
choose 0 _ = [[]]
choose _ [] = []
choose n (x:xs) = (map (x:) $ choose (n-1) xs) ++ (choose n xs)

listProduct [] = [[]]
listProduct (list:lists) = [x:xs | x <- list, xs <- listProduct lists]

crtCombine (n1,m1) (n2,m2)
  | g == 1 = ((n1 + k*m1) `mod` (m1*m2), m1*m2)
    where
      g = gcd m1 m2
      k = (n2 - n1) * inverse m1 m2