[Haskell-cafe] Re: Importing Data.Char speeds up ghc around 70%

[Joost Behrends]

Daniel Fischer <daniel.is.fischer <at> web.de> writes:
 
> Of course, one minute after I sent my previous mail, I receive this one :(
> However, one point, it might be faster to factor out all factors p in found 
> and only then compute the intsqrt, like
> 
> found x = x{dividend = xstop, bound = intsqrt xstop, result = result x ++ 
> replicate k p}
> 	  where
> 	    p = divisor x
> 	    (xstop,k) = go (dividend x) 0
> 	    go n m
> 		| r == 0 	= go q (m+1)
> 		| otherwise = (n,m)
> 		  where
> 		    (q,r) = n `divMod` p

True - but be aware, that this will slightly slow down the computation for 
not multiple factors. And - as you recently noted - the really expensive
part are all the tried factors, which do not divide the queried number.

All this is just a first approach to the problem. When i talk of "naively
programmed", then i want to say, that number theorists might have much better 
numerical orders marching through all primes plus some more odd numbers.
I didn't search for that on the net.

The last version was some kind of resign from tries like this:

firstPrimes = [3,5,7,11,13,17]
start = last firstPrimes
pac = product firstPrimes
slen = length lsumds

lsumds = drop 1 (fst$getSummands (singleton start, start)) where
    getSummands :: (Seq Int, Int) -> (Seq Int, Int)
    getSummands r |snd r < bnd    = getSummands ((fst r)|>k, snd r + k) 
		  |otherwise      = r
        where
            bnd = 2*pac + start
            k = getNext (snd r)
            getNext n |and [(n+2)`mod`x>0 | x<-firstPrimes] = 2 
                      |otherwise                            = 2 + getNext (n+2)

smallmod :: Int -> Int -> Int
smallmod n m | n<m = n | otherwise = 0

divstep :: (DivIter,Int) -> (DivIter, Int)
divstep (x,n) | and [(fromInteger $ divisor x)<start, ximod>0] = 
                          (x {divisor = divisor x + 2}, n)
              | (fromInteger$divisor x) < start = 
                                      (x {dividend = xidiv, 
                                          bound = intsqrt(xidiv), 
                                          result = result x ++ [divisor x]}, n) 
              | ximod>0 = 
(x {divisor = divisor x + toInteger (index lsumds n)}, smallmod (n+1) slen)
              | otherwise = (x {dividend = xidiv, 
                                bound    = intsqrt(xidiv), 
                                result   = result x ++ [divisor x]}, n) 
    where 
        (xidiv, ximod) = divMod (dividend x) (divisor x)

divisions :: (DivIter, Int) -> (DivIter, Int)
divisions (y,n) | divisor y <= bound y = divisions (divstep (y,n))
                | otherwise            = (y,0)

Here the additions to divisor are taken from the sequence lsmnds (List of
SuMaNDS) - the type Seq from Data.Sequence is faster with the function index 
than Data.List with !!. getSummands is a kind of reduced sieve of 
Eratosthenes. The main improvement is the longest line:

|ximod>0 = (x {divisor = divisor x + toInteger (index lsumds n)}, 
               smallmod (n+1) slen)

I even considered converting lsmnds to ByteString and storing them - the
build of lsmnds for firstPrimes = [3,5,7,11,13,17,19,23,29] (which already
has some MB footprint) takes several minutes.  

But we have to track the number of iteration we are in. And that eats up
much more than the reduction of divisions for "failing" factors. The code works
(called slightly modificated by primefactors), but needs 5.41 minutes
for 2^61+1 :((. Also expensive might be the lookup in lsumds - the code gets
even slower with longer lists for firstPrimes.

divisions (d6$d2$d6$d4$d2$d4$d2$d4 y) is derived from

lsmnds [3,5] = [4,2,4,2,4,6,2,6].

For me the whole matter is closed for now - the 1.34 minutes are no bad result.
Amd anyway the code might represent a not too bad lower bound for efficiency of
decomposing algorithms. 

Auf Wiedersehen, Joost