https://wiki.haskell.org/api.php?action=feedcontributions&user=3Jane&feedformat=atomHaskellWiki - User contributions [en]2024-03-28T16:41:08ZUser contributionsMediaWiki 1.35.5https://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17649Talk:99 questions/31 to 412007-12-17T19:24:18Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35, given from "Another possibility is .." on.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only.<br />
<br />
b) factor*factor is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then result res ++ [z] else result res<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, snd m ++ [2] )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y |ximod > 0 = y { divisor=(divisor y) + 2 } <br />
|otherwise = y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = result y ++ [divisor y]}<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y |divisor y <= bound y = divisions (divstep y) <br />
|otherwise = y</haskell><br />
<br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement. Compiled with this version:<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 20 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17648Talk:99 questions/31 to 412007-12-17T19:23:54Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35, given from "Another possibility is .." on.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only.<br />
b) factor*factor is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then result res ++ [z] else result res<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, snd m ++ [2] )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y |ximod > 0 = y { divisor=(divisor y) + 2 } <br />
|otherwise = y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = result y ++ [divisor y]}<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y |divisor y <= bound y = divisions (divstep y) <br />
|otherwise = y</haskell><br />
<br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement. Compiled with this version:<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 20 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17647Talk:99 questions/31 to 412007-12-17T19:15:27Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35, given from "Another possibility is .." on.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then the dividend for the following divisions gets possibly smaller and - what's more - the square root of it too.<br />
<br />
b) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
c) factor*factor is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then result res ++ [z] else result res<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, snd m ++ [2] )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y |ximod > 0 = y { divisor=(divisor y) + 2 } <br />
|otherwise = y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = result y ++ [divisor y]}<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y |divisor y <= bound y = divisions (divstep y) <br />
|otherwise = y</haskell><br />
<br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement. Compiled with this version:<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 20 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17646Talk:99 questions/31 to 412007-12-17T19:14:54Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35, given from "Another possibility is .." on.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then the dividend for the following divisions gets possibly smaller and - what's more - the square root of it too.<br />
<br />
b) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
c) factor*factor is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then result res ++ [z] else result res<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, snd m ++ [2] )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y |ximod > 0 = y { divisor=(divisor y) + 2 } <br />
|otherwise = y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = result y ++ [divisor y]}<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y |divisor y <= bound y = divisions (divstep y) <br />
|otherwise = y</haskell><br />
<br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement. Compiled with this version:<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 20 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17645Talk:99 questions/31 to 412007-12-17T19:08:29Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35, given from "Another possibility is .." on.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then the dividend for the following divisions gets possibly smaller and - what's more - the square root of it too.<br />
<br />
b) The division by a current divisor is done twice in the solution, once for getting the quotient, a second time for getting the remainder. We must use divMod. This reduces the number of the expensive divisions to the half again.<br />
<br />
c) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
d) factor*factor is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
All in all these proposals should give improvement in performance in a magnitude of 1000% to the solution in the text, with code still naive and easily readable for beginners in Haskell (as me):<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then result res ++ [z] else result res<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, snd m ++ [2] )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y |ximod > 0 = y { divisor=(divisor y) + 2 } <br />
|otherwise = y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = result y ++ [divisor y]}<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y |divisor y <= bound y = divisions (divstep y) <br />
|otherwise = y</haskell><br />
<br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement. Compiled with this version:<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 20 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17644Talk:99 questions/31 to 412007-12-17T19:06:21Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35, given from "Another possibility is .." on.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then the dividend for the following divisions gets possibly smaller and - what's more - the square root of it too.<br />
<br />
b) The division by a current divisor is done twice in the solution, once for getting the quotient, a second time for getting the remainder. We must use divMod. This reduces the number of the expensive divisions to the half again.<br />
<br />
c) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
d) factor*factor is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
All in all these proposals should give improvement in performance in a magnitude of 1000% to the solution in the text, with code still naive and easily readable for beginners in Haskell (as me):<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then result res ++ [z] else result res<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, snd m ++ [2] )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y |ximod > 0 = y { divisor=(divisor y) + 2 } <br />
|otherwise = y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = result y ++ [divisor y]}<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y |divisor y <= bound y = divisions (divstep y) <br />
|otherwise = y</haskell><br />
<br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement. Compiled with this version:<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 20 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17643Talk:99 questions/31 to 412007-12-17T15:03:20Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35, given from "Another possibility is .." on.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then the dividend for the following divisions gets possibly smaller and - what's more - the square root of it too.<br />
<br />
b) The division by a current divisor is done twice in the solution, once for getting the quotient, a second time for getting the remainder. We must use divMod. This reduces the number of the expensive divisions to the half again.<br />
<br />
c) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
d) factor*factor is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
All in all these proposals should give improvement in performance in a magnitude of 1000% to the solution in the text, with code still naive and easily readable for beginners in Haskell (as me):<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then result res ++ [z] else result res<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, snd m ++ [2] )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y = if ximod>0 then <br />
y { divisor=(divisor y) + 2 } <br />
else<br />
y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = result y ++ [divisor y] }<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y = if divisor y <= bound y then divisions (divstep y) <br />
else y</haskell><br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement. Compiled with this version:<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 20 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17625Talk:99 questions/31 to 412007-12-17T00:26:49Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35, given from "Another possibility is .." on.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then the dividend for the following divisions gets possibly smaller and - what's more - the square root of it too.<br />
<br />
b) The division by a current divisor is done twice in the solution, once for getting the quotient, a second time for getting the remainder. We must use divMod. This reduces the number of the expensive divisions to the half again.<br />
<br />
c) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
d) factor*factor is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
All in all these proposals should give improvement in performance in a magnitude of 1000% to the solution in the text, with code still naive and easily readable for beginners in Haskell (as me):<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then result res ++ [z] else result res<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, snd m ++ [2] )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y = if ximod>0 then <br />
y { divisor=(divisor y) + 2 } <br />
else<br />
y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = result y ++ [divisor y] }<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y = if divisor y <= bound y then divisions (divstep y) <br />
else y</haskell><br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement. Compiled with this version:<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 30 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs.</div>3Janehttps://wiki.haskell.org/index.php?title=99_questions/31_to_41&diff=1762499 questions/31 to 412007-12-17T00:13:29Z<p>3Jane: </p>
<hr />
<div>__NOTOC__<br />
<br />
This is part of [[H-99:_Ninety-Nine_Haskell_Problems|Ninety-Nine Haskell Problems]], based on [https://prof.ti.bfh.ch/hew1/informatik3/prolog/p-99/ Ninety-Nine Prolog Problems].<br />
<br />
== Arithmetic ==<br />
<br />
== Problem 31 ==<br />
<br />
(**) Determine whether a given integer number is prime.<br />
<br />
Example:<br />
<pre><br />
* (is-prime 7)<br />
T<br />
</pre><br />
Example in Haskell:<br />
<pre><br />
P31> isPrime 7<br />
True<br />
</pre><br />
<br />
Solution:<br />
<haskell><br />
isPrime :: Integral a => a -> Bool<br />
isPrime p = p > 1 && all (\n -> p `mod` n /= 0 ) $ takeWhile (\n -> n*n <= p) [2..] <br />
</haskell><br />
<br />
Well, a natural number p is a prime number iff it is larger than 1 and no natural number n with n >= 2 and n^2 <= p is a divisor of p. That's exactly what is implemented: we take the list of all integral numbers starting with 2 as long as their square is at most p and check that for all these n there is a remainder concerning the division of p by n.<br />
<br />
== Problem 32 ==<br />
<br />
(**) Determine the greatest common divisor of two positive integer numbers.<br />
Use Euclid's algorithm.<br />
<br />
<pre><br />
Example:<br />
* (gcd 36 63)<br />
9<br />
<br />
Example in Haskell:<br />
[myGCD 36 63, myGCD (-3) (-6), myGCD (-3) 6]<br />
[9,3,3]<br />
</pre><br />
<br />
Solution:<br />
<haskell><br />
gcd' 0 y = y<br />
gcd' x y = gcd' (y `mod` x) x<br />
myGCD x y | x < 0 = myGCD (-x) y<br />
| y < 0 = myGCD x (-y)<br />
| y < x = gcd' y x<br />
| otherwise = gcd' x y<br />
</haskell><br />
<br />
The Prelude includes a gcd function, so we have to choose another name for ours. The function gcd' is a straightforward implementation of Euler's algorithm, and myGCD is just a wrapper that makes sure the arguments are positive and in increasing order.<br />
<br />
== Problem 33 ==<br />
<br />
(*) Determine whether two positive integer numbers are coprime.<br />
Two numbers are coprime if their greatest common divisor equals 1.<br />
<br />
Example:<br />
<pre><br />
* (coprime 35 64)<br />
T<br />
</pre><br />
<br />
Example in Haskell:<br />
<pre><br />
* coprime 35 64<br />
True <br />
</pre><br />
<br />
Solution:<br />
<haskell><br />
coprime a b = gcd a b == 1<br />
</haskell><br />
<br />
Here we use the prelude function for computing gcd's along with a test of the result's equality to one. <br />
<br />
== Problem 34 ==<br />
<br />
(**) Calculate Euler's totient function phi(m).<br />
Euler's so-called totient function phi(m) is defined as the number of positive integers r (1 <= r < m) that are coprime to m.<br />
Example: m = 10: r = 1,3,7,9; thus phi(m) = 4. Note the special case: phi(1) = 1.<br />
<pre><br />
Example:<br />
* (totient-phi 10)<br />
4<br />
Example in Haskell:<br />
* totient 10<br />
4<br />
</pre><br />
<br />
Solution:<br />
<haskell><br />
totient 1 = 1<br />
totient a = length $ filter (coprime a) [1..a-1]<br />
where coprime a b = gcd a b == 1<br />
</haskell><br />
<br />
We take coprime from the previous exercise and give it to filter, which applies it to each element of a list from 1 to one less than the number, returning only those that are true. length tells us how many elements are in the resulting list, and thus how many elements are coprime to n<br />
<br />
== Problem 35 ==<br />
<br />
(**) Determine the prime factors of a given positive integer.<br />
Construct a flat list containing the prime factors in ascending order.<br />
<br />
Example:<br />
<pre><br />
* (prime-factors 315)<br />
(3 3 5 7)<br />
</pre><br />
Example in Haskell:<br />
<pre><br />
> primeFactors 315<br />
[3, 3, 5, 7]<br />
</pre><br />
<br />
Solution:<br />
<haskell><br />
primeFactors :: Integer -> [Integer]<br />
primeFactors a = let (f, f1) = factorPairOf a<br />
f' = if prime f then [f] else primeFactors f<br />
f1' = if prime f1 then [f1] else primeFactors f1<br />
in f' ++ f1'<br />
where<br />
factorPairOf a = let f = head $ factors a<br />
in (f, div a f)<br />
factors a = filter (isFactor a) [2..a-1]<br />
isFactor a b = rem a b == 0<br />
prime a = (length $ factors a) == 0<br />
</haskell><br />
<br />
Kind of ugly, but it works, though it may have bugs in corner cases. This uses the factor tree method of finding prime factors of a number. factorPairOf picks a factor and takes it and the factor you multiply it by and gives them to primeFactors. primeFactors checks to make sure the factors are prime. If not it prime factorizes them. In the end a list of prime factors is returned.<br />
<br />
Another possibility is to observe that you need not ensure that potential divisors are primes, as long as you consider them in ascending order:<br />
<haskell><br />
primeFactors n = primeFactors' n 2 where<br />
primeFactors' 1 _ = []<br />
primeFactors' n factor<br />
| n `mod` factor == 0 = factor : primeFactors' (n `div` factor) factor<br />
| otherwise = primeFactors' n (factor + 1)<br />
</haskell><br />
Thus, we just loop through all possible factors and add them to the list if they divide the original number. As the primes get farther apart, though, this will do a lot of needless checks to see if composite numbers are prime factors.<br />
However we can stop as soon as the candidate factor exceeds the square root of n:<br />
<haskell><br />
primeFactors n = primeFactors' n 2 where<br />
primeFactors' n factor<br />
| factor*factor > n = [n]<br />
| n `mod` factor == 0 = factor : primeFactors' (n `div` factor) factor<br />
| otherwise = primeFactors' n (factor + 1)<br />
</haskell><br />
<br />
You can avoid the needless work by just looping through the primes:<br />
<haskell><br />
primeFactors n = factor n primes<br />
where factor n (p:ps) | p*p > n = [n]<br />
| n `mod` p /= 0 = factor n ps<br />
| otherwise = p : factor (n `div` p) (p:ps)<br />
primes = 2 : filter ((==1) . length . primeFactors) [3,5..]<br />
</haskell><br />
(see also discussion of this solution following "Discuss this page")<br />
<br />
== Problem 36 ==<br />
(**) Determine the prime factors of a given positive integer.<br />
<br />
Construct a list containing the prime factors and their multiplicity.<br />
<br />
Example:<br />
<pre><br />
* (prime-factors-mult 315)<br />
((3 2) (5 1) (7 1))<br />
</pre><br />
Example in Haskell:<br />
<pre><br />
*Main> prime_factors_mult 315<br />
[(3,2),(5,1),(7,1)]<br />
</pre><br />
<br />
Solution:<br />
<haskell><br />
prime_factors_mult n = map swap $ encode $ primeFactors n<br />
where swap (x,y) = (y,x)<br />
</haskell><br />
using <tt>primeFactors</tt> from problem 35 to generate the list of prime factors in ascending order, and <tt>encode</tt> from problem 10 to compress this list to a list of numbers paired with their multiplicity.<br />
<br />
== Problem 37 ==<br />
<br />
(**) Calculate Euler's totient function phi(m) (improved).<br />
See problem 34 for the definition of Euler's totient function. If the list of the prime factors of a number m is known in the form of problem 36 then the function phi(m) can be efficiently calculated as follows: Let ((p1 m1) (p2 m2) (p3 m3) ...) be the list of prime factors (and their multiplicities) of a given number m. Then phi(m) can be calculated with the following formula:<br />
<br />
<pre><br />
phi(m) = (p1 - 1) * p1 ** (m1 - 1) + (p2 - 1) * p2 ** (m2 - 1) + (p3 - 1) * p3 ** (m3 - 1) + ...<br />
</pre><br />
<br />
Note that a ** b stands for the b'th power of a.<br />
<i>Note</i>: Actually, the official problems show this as a sum, but it should be a product.<br />
<br />
Solution: Given prime_factors_mult from problem 36, we get<br />
<haskell><br />
totient m = product [(p - 1) * p ^ (c - 1) | (p, c) <- prime_factors_mult m]<br />
</haskell><br />
This just uses a list comprehension to build each term of the product in the formula for phi, then multiplies them all.<br />
<br />
== Problem 38 ==<br />
<br />
(*) Compare the two methods of calculating Euler's totient function.<br />
<br />
Use the solutions of problems 34 and 37 to compare the algorithms. Take the number of reductions as a measure for efficiency. Try to calculate phi(10090) as an example.<br />
<br />
(no solution required)<br />
<br />
== Problem 39 ==<br />
<br />
(*) A list of prime numbers.<br />
<br />
Given a range of integers by its lower and upper limit, construct a list of all prime numbers in that range.<br />
<br />
<pre><br />
Example in Haskell:<br />
P29> primesR 10 20<br />
[11,13,17,19]<br />
</pre><br />
<br />
Solution 1:<br />
<haskell><br />
primesR :: Integral a => a -> a -> [a]<br />
primesR a b = filter isPrime [a..b]<br />
</haskell><br />
<br />
If we are challenged to give all primes in the range between a and b we simply take all number from a up to b and filter the primes out.<br />
<br />
Solution 2:<br />
<haskell><br />
primes :: Integral a => [a]<br />
primes = let sieve (n:ns) = n:sieve [ m | m <- ns, m `mod` n /= 0 ] in sieve [2..]<br />
<br />
primesR :: Integral a => a -> a -> [a]<br />
primesR a b = takeWhile (<= b) $ dropWhile (< a) primes<br />
</haskell><br />
<br />
Another way to compute the claimed list is done by using the ''Sieve of Eratosthenes''. The <hask>primes</hask> function generates a list of all (!) prime numbers using this algorithm and <hask>primesR</hask> filter the relevant range out. [But this way is very slow and I only presented it because I wanted to show how nice the ''Sieve of Eratosthenes'' can be implemented in Haskell :)]<br />
<br />
== Problem 40 ==<br />
(**) Goldbach's conjecture.<br />
Goldbach's conjecture says that every positive even number greater than 2 is the sum of two prime numbers. Example: 28 = 5 + 23. It is one of the most famous facts in number theory that has not been proved to be correct in the general case. It has been numerically confirmed up to very large numbers (much larger than we can go with our Prolog system). Write a predicate to find the two prime numbers that sum up to a given even integer.<br />
<br />
<pre><br />
Example:<br />
* (goldbach 28)<br />
(5 23)<br />
<br />
Example in Haskell:<br />
*goldbach 28<br />
(5, 23)<br />
</pre><br />
<br />
Solution 1:<br />
<haskell><br />
goldbach a = head $<br />
filter (\e -> (isPrime $ fst e) && (isPrime $ snd e)) $<br />
map (\e -> (e, a - e)) [1,3..div a 2]<br />
where<br />
factors a = filter (isFactor a) [2..a-1]<br />
isFactor a b = rem a b == 0<br />
isPrime a = (length $ factors a) == 0<br />
</haskell><br />
<br />
Woohoo! I've solved the goldbach conjecture! Where do I collect my prize? This does the obvious thing. It makes a list of odd numbers and then uses it to make up pairs of odd numbers that sum to n. Then it looks for a pair of odd numbers where both are prime and returns it as a tuple.<br />
<br />
Solution 2:<br />
<haskell><br />
-- using the previous problem...<br />
goldbach n = head [(x,y) | x <- pr, y <- pr, x+y==n]<br />
where pr = takeWhile (< n) primes<br />
</haskell><br />
<br />
Solution 3:<br />
<haskell><br />
goldbach n = head [(x,y) | x <- takeWhile (< n) primes, let y = n-x, isPrime y]<br />
</haskell><br />
<br />
== Problem 41 ==<br />
<br />
(**) Given a range of integers by its lower and upper limit, print a list of all even numbers and their Goldbach composition.<br />
<br />
In most cases, if an even number is written as the sum of two prime numbers, one of them is very small. Very rarely, the primes are both bigger than say 50. Try to find out how many such cases there are in the range 2..3000.<br />
<br />
Example:<br />
<pre><br />
* (goldbach-list 9 20)<br />
10 = 3 + 7<br />
12 = 5 + 7<br />
14 = 3 + 11<br />
16 = 3 + 13<br />
18 = 5 + 13<br />
20 = 3 + 17<br />
* (goldbach-list 1 2000 50)<br />
992 = 73 + 919<br />
1382 = 61 + 1321<br />
1856 = 67 + 1789<br />
1928 = 61 + 1867<br />
</pre><br />
<br />
Example in Haskell:<br />
<pre><br />
*Exercises> goldbachList 9 20<br />
[(3,7),(5,7),(3,11),(3,13),(5,13),(3,17)]<br />
*Exercises> goldbachList' 4 2000 50<br />
[(73,919),(61,1321),(67,1789),(61,1867)]<br />
</pre><br />
<br />
Solution:<br />
<haskell><br />
goldbachList lb ub = map goldbach $ [even_lb,even_lb+2..ub]<br />
where even_lb = max ((lb+1) `div` 2 * 2) 4<br />
goldbachList' lb ub mv = filter (\(a,b) -> a > mv && b > mv) $<br />
goldbachList lb ub<br />
</haskell><br />
using the <tt>goldbach</tt> function from problem 40.<br />
<br />
[[Category:Tutorials]]</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17623Talk:99 questions/31 to 412007-12-17T00:09:32Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35, given from "Another possibility is .." on.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then the dividend for the following divisions gets possibly smaller and - what's more - the square root of it too.<br />
<br />
b) The division by a current divisor is done twice in the solution, once for getting the quotient, a second time for getting the remainder. We must use divMod. This reduces the number of the expensive divisions to the half again.<br />
<br />
c) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
d) factor*factor is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
All in all these proposals should give improvement in performance in a magnitude of 1000% to the solution in the text, with code still naive and easily readable for beginners in Haskell (as me):<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then result res ++ [z] else result res<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, snd m ++ [2] )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y = if ximod>0 then <br />
y { divisor=(divisor y) + 2 } <br />
else<br />
y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = result y ++ [divisor y] }<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y = if divisor y <= bound y then divisions (divstep y) <br />
else y</haskell><br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Compiled with this version for "divisions" (indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement):<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 30 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17622Talk:99 questions/31 to 412007-12-17T00:01:45Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35, given from "Another possibility is .." on.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then the dividend for the following divisions gets possibly smaller and - what's more - the square root of it too.<br />
<br />
b) The division by a current divisor is done twice in the solution, once for getting the quotient, a second time for getting the remainder. We must use divMod. This reduces the number of the expensive divisions to the half again.<br />
<br />
c) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
d) factor*factor is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
All in all these proposals should give improvement in performance in a magnitude of 1000% to the solution in the text, with code still naive and easily readable for beginners in Haskell (as me):<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then revres ++ [z] else revres<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
revres = reverse (result res)<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, 2:(snd m) )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y = if ximod>0 then <br />
y { divisor=(divisor y) + 2 } <br />
else<br />
y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = (divisor y):(result y) }<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y = if divisor y <= bound y then divisions (divstep y) <br />
else y</haskell><br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Compiled with this version for "divisions" (indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement):<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 30 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17621Talk:99 questions/31 to 412007-12-17T00:00:54Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35, given from "Another possibility is .." on.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then the dividend for the following divisions gets possibly smaller and - what's more - the square root of it too.<br />
<br />
b) The division by a current divisor is done twice in the solution, once for getting the quotient, a second time for getting the remainder. We must use divMod. This reduces the number of the expensive divisions to the half again.<br />
<br />
c) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
d) factor*factor is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
All in all these proposals should give improvement in performance in a magnitude of 1000% to the solution in the text, with code still naive and easily readable for beginners in Haskell (as me):<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then revres ++ [z] else revres<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
revres = reverse (result res)<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, 2:(snd m) )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y = if ximod>0 then <br />
y { divisor=(divisor y) + 2 } <br />
else<br />
y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = (divisor y):(result y) }<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y = if divisor y <= bound y then divisions (divstep y) <br />
else y</haskell><br />
<br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Compiled with this version for "divisions" (indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement):<br />
<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 30 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17620Talk:99 questions/31 to 412007-12-16T23:58:15Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35, given from "Another possibility is .." on.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then the dividend for the following divisions gets possibly smaller and - what's more - the square root of it too.<br />
<br />
b) The division by a current divisor is done twice in the solution, once for getting the quotient, a second time for getting the remainder. We must use divMod. This reduces the number of the expensive divisions to the half again.<br />
<br />
c) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
d) factor*factor is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
All in all these proposals should give improvement in performance in a magnitude of 1000% to the solution in the text, with code still naive and easily readable for beginners in Haskell (as me):<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then revres ++ [z] else revres<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
revres = reverse (result res)<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, 2:(snd m) )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y = if ximod>0 then <br />
y { divisor=(divisor y) + 2 } <br />
else<br />
y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = (divisor y):(result y) }<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y = if divisor y <= bound y then divisions (divstep y) <br />
else y</haskell><br />
<br />
<br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Compiled with this version for "divisions" (indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement):<br />
<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 30 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17619Talk:99 questions/31 to 412007-12-16T23:52:44Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35, given from "Another possibility is .." on.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then that divisor is done completely. Thus we must go through the divisors only once for the whole process (not start by dividing by 2 again after every "successful" division). This reduces the number of computations at least to the half.<br />
<br />
b) The division by a current divisor is done twice in the solution, once for getting the quotient, a second time for getting the remainder. We must use divMod. This reduces the number of the expensive divisions to the half again.<br />
<br />
c) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
d) n*n is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
All in all these proposals should give improvement in performance in a magnitude of 1000% to the solution in the text, with code still naive and easily readable for beginners in Haskell (as me):<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then revres ++ [z] else revres<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
revres = reverse (result res)<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, 2:(snd m) )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y = if ximod>0 then <br />
y { divisor=(divisor y) +2 } <br />
else<br />
y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = (divisor y):(result y) }<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y = if divisor y <= bound y then divisions (divstep y) <br />
else y</haskell><br />
<br />
<br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Compiled with this versions for "divisions"<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 30 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs. Indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17618Talk:99 questions/31 to 412007-12-16T23:47:01Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then that divisor is done completely. Thus we must go through the divisors only once for the whole process (not start by dividing by 2 again after every "successful" division). This reduces the number of computations at least to the half.<br />
<br />
b) The division by a current divisor is done twice in the solution, once for getting the quotient, a second time for getting the remainder. We must use divMod. This reduces the number of the expensive divisions to the half again.<br />
<br />
c) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
d) n*n is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
All in all these proposals should give improvement in performance in a magnitude of 1000% to the solution in the text, with code still naive and easily readable for beginners in Haskell (as me):<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then revres ++ [z] else revres<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
revres = reverse (result res)<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, 2:(snd m) )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y = if ximod>0 then <br />
y { divisor=(divisor y) +2 } <br />
else<br />
y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = (divisor y):(result y) }<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y = if divisor y <= bound y then divisions (divstep y) <br />
else y</haskell><br />
<br />
<br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Compiled with this versions for "divisions"<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 30 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs. Indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17617Talk:99 questions/31 to 412007-12-16T23:44:31Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then that divisor is done completely. Thus we must go through the divisors only once for the whole process (not start by dividing by 2 again after every "successful" division). This reduces the number of computations at least to the half.<br />
<br />
b) The division by a current divisor is done twice in the solution, once for getting the quotient, a second time for getting the remainder. We must use divMod. This reduces the number of the expensive divisions to the half again.<br />
<br />
c) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
d) n*n is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.<br />
<br />
All in all these proposals should give improvement in performance in a magnitude of 1000% to the solution in the text, with code still naive and easily readable for beginners in Haskell (as me):<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter { dividend :: Integer, <br />
divisor :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then revres ++ [z] else revres<br />
where <br />
res = divisions (DivIter { dividend = o1,<br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2})<br />
revres = reverse (result res)<br />
z = dividend res <br />
--indeed gets 1 in certain situations<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, 2:(snd m) )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y = if ximod>0 then <br />
y { divisor=(divisor y) +2 } <br />
else<br />
y {dividend = xidiv, <br />
bound = intsqrt xidiv, <br />
result = (divisor y):(result y) }<br />
where (xidiv, ximod) = divMod (dividend y) (divisor y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y = if divisor y <= bound y then divisions (divstep y) <br />
else y</haskell><br />
<br />
<br />
This computes <br />
primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]<br />
in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Compiled with this versions for "divisions"<br />
<br />
<haskell>divisions = do<br />
y <- get<br />
if divisor y <= bound y then do<br />
put ( divstep y )<br />
divisions<br />
else <br />
return y</haskell><br />
<br />
called from primefactors by <br />
<br />
<haskell>res = execState divisions (DivIter { dividend = o1, <br />
divisor = 3, <br />
bound = intsqrt(o1),<br />
result = o2 })</haskell><br />
<br />
it computes <br />
primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]<br />
in about 30 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs. Indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17566Talk:99 questions/31 to 412007-12-15T23:22:13Z<p>3Jane: </p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^64. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n, which exceeds 2^32. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then it is done completely. Thus we must go through the divisors only once for the whole process (not start by dividing by 2 again after every "successful" division). This reduces the number of computations at least to the half.<br />
<br />
b) The division by a current divisor is done twice in the solution, once for getting the quotient, a second time for getting the remainder. We must use divMod. This reduces the number of the expensive divisions to the half again.<br />
<br />
c) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
d)n*n is computed at each step. We should compute a bound ("intsqrt") only once for each "successful" division.<br />
<br />
All in all these proposals should give improvement in performance in a magnitude of 1000% to the solution in the text, with code still naive and easily readable for beginners in Haskell (as me):<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter {x :: Integer, <br />
divi :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then revres ++ [z] else revres<br />
where<br />
res = divisions (DivIter {x = o1, <br />
divi = 3, <br />
bound = intsqrt(o1), <br />
result = o2})<br />
revres = reverse (result res)<br />
z = x res -- yes, sometimes equal to 1<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, 2:(snd m) )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y = if ximod>0 then y { divi=(divi y) +2 } else<br />
y {x = xidiv, bound = intsqrt xidiv, result = (divi y):(result y) }<br />
where<br />
(xidiv, ximod) = divMod (x y) (divi y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y = if divi y <= bound y then divisions (divstep y) else y</haskell><br />
<br />
<br />
The reader should notice, that "twosect" is an easy model for the (still easy) function "divstep".<br />
<br />
This proposal assumes, that named fields really are syntactic sugar and will be completely factored out by compilation. Otherwise a quadruple should be used.<br />
<br />
"divisions" is separated from "divstep" here not for readability, but for possible replacement. We MUST ensure tail recursion here, or use a monadic version of "divisions" with destructive updates - but that's a more advanced topic.</div>3Janehttps://wiki.haskell.org/index.php?title=Talk:99_questions/31_to_41&diff=17564Talk:99 questions/31 to 412007-12-15T22:42:45Z<p>3Jane: Critics of Solution 35</p>
<hr />
<div>Hi,<br />
<br />
there are several problems with the solution to 35.<br />
<br />
First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^64. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n, which exceeds 2^32. And there are doubts anyway concerning the price for building such a list first.<br />
<br />
Second: <br />
<br />
a) If found a dividend, we must factor it out completely, that means divide by it, as often as possible. Then it is done completely. Thus we must go through the divisors only once for the whole process (not start by dividing by 2 again after every "successful" division). This reduces the number of computations at least to the half.<br />
<br />
b) The division by a current divisor is done twice in the solution, once for getting the quotient, a second time for getting the remainder. We must use divMod. This reduces the number of the expensive divisions to the half again.<br />
<br />
c) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.<br />
<br />
d)n*n is computed at each step. We should compute a bound ("intsqrt") only once for each "successful" division.<br />
<br />
All in all these proposals should give improvement in performance in a magnitude of 1000% to the solution in the text, with code still naive and easily readable for beginners in Haskell (as me):<br />
<br />
<br />
<haskell>data DivIter = <br />
DivIter {x :: Integer, <br />
divi :: Integer,<br />
bound :: Integer,<br />
result :: [Integer] }<br />
<br />
intsqrt m = floor (sqrt $ fromInteger m)<br />
<br />
primefactors :: Integer -> [Integer]<br />
primefactors n | n < 2 = []<br />
| even n = o2 ++ (primefactors o1)<br />
| otherwise = if z /=1 then revres ++ [z] else revres<br />
where<br />
res = divisions (DivIter {x = o1, <br />
divi = 3, <br />
bound = intsqrt(o1), <br />
result = o2})<br />
revres = reverse (result res)<br />
z = x res -- yes, sometimes equal to 1<br />
(o1, o2) = twosect (n, [])<br />
<br />
twosect :: (Integer, [Integer]) -> (Integer, [Integer])<br />
twosect m | odd (fst m) = m<br />
| even (fst m) = twosect ( div (fst m) 2, 2:(snd m) )<br />
<br />
divstep :: DivIter -> DivIter <br />
divstep y = if ximod>0 then y { divi=(divi y) +2 } else<br />
y {x = xidiv, bound = intsqrt xidiv, result = (divi y):(result y) }<br />
where<br />
(xidiv, ximod) = divMod (x y) (divi y)<br />
<br />
divisions :: DivIter -> DivIter<br />
divisions y = if divi y <= bound y then divisions (divstep y) else y</haskell><br />
<br />
<br />
The reader should notice, that "twosect" is an easy model for the (still easy) function "divstep".<br />
<br />
This proposal assumes, that named fields really are syntactic sugar and will be completely factored out by compilation. Otherwise a quadruple should be used.<br />
<br />
"divisions" is separated from "divstep" here not for readability, but for possible replacement. We MUST ensure tail recursion here, or use a monadic version of "divisions" with destructive updates - but that's a more advanced topic.</div>3Jane