gistfile1.hs

import Data.Array

data Op = DN | UP deriving (Show, Eq, Ord)

paso :: (Integral a) => (a -> a) -> a -> a
paso hazpar x = if odd x then div (hazpar x) 2 else div x 2

serieInfinita :: (Integral a) => (a -> a) -> a -> [Op]
serieInfinita f x = (if odd x then UP else DN) : serieInfinita f siguiente where
    siguiente = paso f x

serie 1 = [1]
serie x = (if odd x then 1 else 0) : serie siguiente where
    siguiente = paso collatz x

collatz :: (Integral a) => a -> a
collatz  = (+1).(*3)

-- Construir el árbol de congruencias

-- como binarios, UP es 1, DN es 0, la potencia mínima está a la
-- izquierda (usar reverse si es necesario)
serieAEntero :: [Op] -> Integer
serieAEntero xs = sum $ zipWith g xs (map (2^) [0..]) where
    g UP k = k
    g DN _ = 0

nivel :: (Integral a) => (a -> a) -> Int -> Array Int Int
nivel f n = array (0,(2^n)-1) ((0,0) : nivel') where
    nivel' = zip (map g [1..((2^n)-1)]) [1..]
    g = fromIntegral . serieAEntero . reverse . take n . serieInfinita collatz

-- requiere entrada de longitud par y no nula (no verificado)
alternos :: [a] -> ([a],[a])
alternos (x:y:xs) = if null xs then ([x],[y]) else (x:m, y:n) where (m, n) = alternos xs

alternosCon :: (a->a->b) -> [a] -> [b]
alternosCon f xs = zipWith f x y where
    (x, y) = alternos xs

torsionesBool :: Int -> [Bool]
torsionesBool n = s where
    tabla = nivel collatz n
    ys = map (tabla !) [0..((2^n)-1)]
    s = alternosCon (>) ys

-- para memoizar
sb16 = torsionesBool 16

recuperar :: Int ->     -- ¿cuantos bits quiero recuperar? max. 16 por sb16
    Integer ->          -- entrada
    [Integer]           -- n bits de la serie collatz de x
recuperar n x = map (`mod` 2) $ resultados 0 paridades where
    congruencias :: [Integer]
    congruencias = take n $ map (\k -> mod x (2^k)) [1..]
    paridades :: [Bool]
    paridades = zipWith (>=) congruencias (map (2^) [0..])
    torsiones :: [Bool]
    torsiones = sb16 -- fuera de la clausura por memoización
    resultados :: Integer -> [Bool] -> [Integer]
    resultados _ [] = []
    resultados y (p:ps) = actual : resultados actual ps where
        t = torsiones !! fromIntegral y
        signo = t /= p
        actual = if signo then y*2 + 1 else y*2