Problem 127

http://projecteuler.net/index.php?section=problems&id=127

いまいちぱっとしないコードである。

import Number
import Data.List
import Data.Array.IArray

upper = 110000

rad = listArray (1,upper).map (product.nub.factors) $ [1..upper]::Array Integer Integer

hamming [1] = [1]
hamming ns = hs
    where hs = (product ns):foldl1 merge [map (n*) hs|n<-ns\\[1]]

gather low lim (x:xs) | lim' <= 1 && low <  1 = [[]]
                      | lim' <= 1 && low >= 1 = []
                      | otherwise = map (x:) (gather low' lim' xs) ++ gather low lim xs
    where low' = low / fromIntegral x
          lim' = lim / fromIntegral x

factorAB = [(factorA,factorB) |
            factorA <- ([1]:).gather 1 upper' $ primes,
            let radA =fromIntegral. product $factorA,
            factorB <- gather radA (upper' /radA)$ primes\\factorA]
    where upper' = fromIntegral upper

findC (fa,fb) = [a+b |
               a <-takeWhile(<upper) $ hamming fa,
               b <-takeWhile(<upper-a).dropWhile(<2*radAB-a)$ hamming fb,
               radAB * (rad!) (a+b) < a+b]
    where radAB = product fa * product fb

main = print.sum.map (sum.findC)$factorAB

a,bをそのまま動かしただけだと、探索範囲が広い（広そう）なので、
a,bの素因数で考えている。

まぁ、速くなったが・・・そんなに速くないぞ

別の方法、シンプル。

import Number
import Data.Array.IArray
import Data.List
import Data.Ord

upper = 110000

rad = listArray (1,upper). map (product.nub.factors) $ [1..upper]::Array Integer Integer

abcHit = [c | 
          (c,radC) <- assocs rad,
          (a,radA) <- takeWhile(\(a,ra)-> ra*radC*2<c) sortRad,
          let b = c -a,
          a < b,
          gcd radC radA == 1,
          radA*(rad!b)*radC < c]
    where sortRad = sortBy(comparing snd). assocs $ rad

main = print.sum$ abcHit

速さは大して変わらなかった。