我必须划分两个非常大的整数(7000+数字是常态)。第一个总是小于第二个:划分返回NaN
> let a = prod fp'
> let b = prod fq'
> a/b
NaN
> length . show $ a
7309
> length . show $ b
7309
> logBase 10 a -- Ehm what? Whatever...
Infinity
我确定有一种方法可以分割这些整数。实际的代码是有点长:
module Main
where
import Data.List (transpose, sort, groupBy)
import Math.Sieve.Factor (factor, sieve)
thueMorse :: [Int]
thueMorse = 0 : interleave (map (1-) thueMorse) (tail thueMorse)
where interleave (x:xs) ys = x : interleave ys xs
cc :: [[Integer]]
cc = zipWith (\n b -> [[2*n+1, 2*n+2], [2*n+2,2*n+1]]!!b) [0..] thueMorse
mapGroup f = map f . groupBy eq . sort
where eq a b = fst a == fst b
takeOut xs = takeOut' xs [] []
where takeOut' [] fs qs = [fs,qs]
takeOut' ([b,e]:xs) fs qs
| e > 0 = takeOut' xs ([b,e]:fs) qs
| e < 0 = takeOut' xs fs ([b,-e]:qs)
| otherwise = takeOut' xs fs qs
sub x y = takeOut . mapGroup collapse $ x ++ y'
where y' = [(b, -e) | (b,e) <- y]
collapse l = [fst $ head l, sum [e | (b,e) <- l]]
factor' :: (Integral a) => Int -> a -> [(a, a)]
factor' l = factor (sieve l)
factorize :: (Integral a) => Int -> [a] -> [(a, a)]
factorize l ns = mapGroup getsum $ concatMap (factor' l) ns
where getsum l = (fst $ head l, sum . map snd $ l)
approx :: Int -> Double
approx n = prod fp'/prod fq' -- <-- The problem is here
where limit = 2^(n-1)
[1:ps,qs] = transpose $ take limit cc
[fp,fq] = map (factorize (2*limit)) [ps,qs]
[fp',fq'] = sub fp fqb
prod l = fromIntegral . product $ [b^e | [b,e] <- l]
功能运行是approx
。从1
到10
是没有问题的,但approx 11
结果是NaN
开始...
> map approx [1..13]
[0.5,0.6666666666666666,0.7,0.7061728395061728,0.7070239390108867,0.7071021244800849,0.7071066220582782,0.707106777975181,0.7071067811490775,0.7071067811862988,NaN,NaN,NaN]
如果你写了'approx n = fromRational $ prod fp'/ prod fq'',那么你得到一个很好的Double值。当然,通过'Rational'比直接将整数转换为'Double'需要更多的时间,但是它解决了你的'overflow〜> NaN'问题。 – 2014-09-11 14:03:28
@DanielFischer:谢谢!这样可行! – rubik 2014-09-11 15:46:12