2014-09-11 106 views
2

我必须划分两个非常大的整数(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。从110是没有问题的,但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] 
+4

如果你写了'approx n = fromRational $ prod fp'/ prod fq'',那么你得到一个很好的Double值。当然,通过'Rational'比直接将整数转换为'Double'需要更多的时间,但是它解决了你的'overflow〜> NaN'问题。 – 2014-09-11 14:03:28

+0

@DanielFischer:谢谢!这样可行! – rubik 2014-09-11 15:46:12

回答

5

ab不是整数。

Prelude> :t (/) 
(/) :: Fractional a => a -> a -> a 
Prelude> :i Fractional 
class Num a => Fractional a where 
    (/) :: a -> a -> a 
    recip :: a -> a 
    fromRational :: Rational -> a 
     -- Defined in ‘GHC.Real’ 
instance Fractional Float -- Defined in ‘GHC.Float’ 
instance Fractional Double -- Defined in ‘GHC.Float’ 

所以你正在运行到的问题是,浮点除法可能最终产生NaNInfinity

你可以使用不同类型的数字,不会失去精度,或以某种方式缩放结果。

+0

也许这是指数溢出,不是精确损失? – tmyklebu 2014-09-11 12:58:01

+0

嗯,你可能是对的,但是如何工作? 1 /无穷大应该是0,对吧? – Sarah 2014-09-11 13:11:43

+3

'a'和'b'对于'Double'来说都太大了,所以你得到的是Infinity/Infinity,它是NaN。 – 2014-09-11 13:57:21