CODEWITHSUNDEEP
haskellarbitrary-precision
Mick Francis
Mick Francis

Reputation: 139

'Infinite' precision reciprocal in Haskell

I need to generate an infinite Haskell list containing all bits (or words) from the fractional part of the reciprocal of an Integer, in MSB-first order. Is there a straightforward way to do this from standard libraries, or do I need to implement a Newton iteration function or similar?

I considered using CReal, but can't find a way of extracting bits/words.

Upvotes: 2

Views: 457

Answers (3)

Mick Francis
Mick Francis

Reputation: 139

After a bit of tinkering, I managed to do this without resorting to Rational / CReal:

recipBits :: Integer -> [Bool]
recipBits n = dblAndMod2 2
    where dblAndMod2 :: Integer -> [Bool]
          dblAndMod2 !r = let bit = r >= n
                              r'  = 2 * (if bit then r - n else r)
                          in  bit : dblAndMod2 r'

Using explicit recursion gains a moderate speedup over using iterate. I'm not too worried about where it goes into a loop, as I'm using it with such large numbers that I'll never reach the loop. Thanks again for the help!

Upvotes: 2

Daniel Wagner
Daniel Wagner

Reputation: 152997

In my other answer, I show how to do this with CReal to answer your question of whether that can be done. But I don't actually think that's a good idea, because it's invoking more power than is needed. Just to give a flavor of what I have in mind for a "more principled" approach, here's what I'm thinking:

bits :: Rational -> [Int]
bits n = whole : bits (2*frac) where
    (whole, frac) = properFraction n

In action:

> take 65 . bits $ 1/3
[0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1]

You'll notice that this is significantly faster than the CReal approach in the other answer. It should also have better memory performance, even if you switch the other answer from CReal to Rational, because it keeps the fraction upper-bounded by 1 (whereas the other solution adds ~1 bit per iteration). It can be made faster still by noticing when it starts to loop. Here's a function that returns the loop in an easily-observable way:

import Data.Set (Set)
import qualified Data.Set as S

data BitsRep
    = Loop Rational [Int]
    | Lollipop [Int] [Int]
    deriving (Eq, Ord, Read, Show)

-- always returns a Lollipop when given an empty set
bitsRaw :: Set Rational -> Rational -> BitsRep
bitsRaw s n = case S.member n s of
    True -> Loop n []
    False -> case bitsRaw (S.insert n s) (2*frac) of
        Lollipop prefix loop -> Lollipop (whole:prefix) loop
        Loop n' loop -> (if n == n' then Lollipop [] else Loop n') (whole:loop)
        where
        (whole, frac) = properFraction n

If you really want it as an infinite list, a short wrapper will do, and significantly reduce the computation needed once the loop point is reached:

bits :: Rational -> [Int]
bits n = prefix ++ cycle loop where
    Lollipop prefix loop = bitsRaw S.empty n

Upvotes: 2

Daniel Wagner
Daniel Wagner

Reputation: 152997

There are more principled ways, but CReal can certainly do it.

> bit n = (`mod` 2) . floor . (*2**n)
> [bit i (pi :: CReal) | i <- [-1..10]]
[1,1,0,0,1,0,0,1,0,0,0,0]
> [bit i (5/8 :: CReal) | i <- [1..10]]
[0,1,0,1,0,0,0,0,0,0,0]

Swap out 2 everywhere for whatever base you like. For an infinite list it's probably cheaper to iterate the multiplication, so:

> bits = map ((`mod` 2) . floor) . iterate (2*)
> take 10 (bits (5/8 :: CReal))
[0,1,0,1,0,0,0,0,0,0]

Again you can swap out 2 for any base you like.

Upvotes: 2

Related Questions