Menu Close

La sucesión ECG

La sucesión ECG estás definida por a(1) = 1, a(2) = 2 y, para n >= 3, a(n) es el menor natural que aún no está en la sucesión tal que a(n) tiene algún divisor común con a(n-1).

Los primeros términos de la sucesión son 1, 2, 4, 6, 3, 9, 12, 8, 10, 5, 15, …

Al dibujar su gráfica, se parece a la de los electrocardiogramas (abreviadamente, ECG). Por ello, la sucesión se conoce como la sucesión ECG.

Definir las funciones

   sucECG :: [Integer]
   graficaSucECG :: Int -> IO ()

tales que

  • sucECG es la lista de los términos de la sucesión ECG. Por ejemplo,
     λ> take 20 sucECG
     [1,2,4,6,3,9,12,8,10,5,15,18,14,7,21,24,16,20,22,11]
     λ> sucECG !! 6000
     6237
  • (graficaSucECG n) dibuja la gráfica de los n primeros términos de la sucesión ECG. Por ejemplo, (graficaSucECG 160) dibuja

Soluciones

import Data.List (delete)
import Graphics.Gnuplot.Simple
 
sucECG :: [Integer]
sucECG = 1 : ecg 2 [2..]
  where ecg x zs = f zs
          where f (y:ys) | gcd x y > 1 = y : ecg y (delete y zs)
                         | otherwise   = f ys
 
graficaSucECG :: Int -> IO ()
graficaSucECG n =
  plotList [ Key Nothing
           , PNG "La_sucesion_ECG.png" 
           ]
           (take n sucECG)

Pensamiento

Algunos desesperados
sólo se curan con soga;
otros, con siete palabras:
la fe se ha puesto de moda.

Antonio Machado

Avanzado

6 soluciones de “La sucesión ECG

  1. frahidzam
    sucECG :: [Integer]
    sucECG = 1 : f 2 [2..]
      where f n xs = (a xs n) : f (a xs n) (delete (a xs n) xs) 
            a xs n= head [x | x <- xs , gcd x n > 1]
     
    graficaSucECG :: Int -> IO ()
    graficaSucECG n = plotList [Key Nothing] (take n sucECG)
  2. alebarmor1

    divisores n = [k | k &lt;- [2..n], mod n k == 0]
    interseccion xs ys = [x | x &lt;- xs, elem x ys]
    divisorComun x y = not ((interseccion (divisores x) (divisores y)) == [])
    genSuc 1 = 1
    genSuc 2 = 2
    genSuc n = head [k | k &lt;- [1..] ,divisorComun k (genSuc (n-1)), and (map (/=k) xs)]
      where xs = [genSuc v | v &lt;- [1..n-1]]
    sucECG :: [Integer]
    sucECG = [genSuc k | k  IO ()
    graficaSucECG n = plotList [Key Nothing] (take n sucECG)
    

  3. luipromor
    sucECG :: [Integer]
    sucECG  = [ f x | x <- [1..]]
                where f 1 = 1
                      f 2 = 2
                      f x = aux  (last ys) ys 3 
                        where  ys = take (x-1) sucECG
                               aux y xs n  | or [ mod n x == 0 && mod y x == 0 | x <- [2..min n y]]                          && notElem n xs = n
                                           | otherwise = aux y xs (n+1) 
    graficaSucECG :: Int -> IO ()
    graficaSucECG n = plotList [] (take n sucECG)
  4. javmarcha1
    sucECG :: [Integer]
    sucECG = [1]++(sucECGaux 2 [2..])
      where sucECGaux x xs = [(d x xs)]++ (sucECGaux (d x xs) (delete (d x xs) xs)) 
            d x xs = head [y | y <- xs , gcd y x > 1]
     
    graficaSucECG :: Int -> IO ()
    graficaSucECG x = plotList [] (take x sucECG)
  5. josejuan

    Usando Set puede reducirse notablemente la constante cuadrática:

    sucECG :: [Int]
    sucECG = 1: f 2 S.empty 3
      where f n s p = n: case filter ((>1).gcd n) (S.toAscList s) of
                           (z:_) -> f z (S.delete z s) p
                           []    -> let w = head $ filter ((>1).gcd n) [p..]
                                    in  f w (S.union s $ S.fromAscList [p..w - 1]) (w + 1)
  6. josejuan

    Usando una criba podemos factorizar rápidamente, pero la búsqueda del siguiente no visitado no veo como hacerla óptima (combinaciones de todos los primos), en todo caso muchísimo más rápida que las anteriores (ej. para n=194561 unos 10seg):

    {-# LANGUAGE BangPatterns #-}
    import System.Environment
    import Data.Vector.Unboxed (Vector, (!))
    import qualified Data.Vector.Unboxed.Mutable as M
    import qualified Data.Vector.Unboxed as V
    import Control.Monad
     
    -- para cada n su mayor divisor primo (biggest prime factor) O(n log log n)
    bpf :: Int -> IO (Vector Int)
    bpf n = do
      v <- M.new (n + 1)
      forM_ [2..n] $ i -> do
        c <- M.read v i
        when (c == 0) $ do
          M.write v i 1
          forM_ [i+i,i+i+i..n] $ j -> M.write v j i
      V.freeze v
     
    -- con lo anterior, podemos factorizar en O(log n)
    primes :: Vector Int -> Int -> [Int]
    primes v = q 0 where q !c !z = let k = v!z
                                   in  if k == c
                                         then q c (z `div` c)
                                         else if k == 1
                                                then if z == c then [] else [z]
                                                else k: q k (z `div` k)
     
    divSup a b = case a `divMod` b of
                  (d, 0) -> d
                  (d, _) -> d + 1
     
    sucECG :: Int -> IO [Int]
    sucECG n = (1:) <$> do
      p <- bpf n -- para factorizar rápido
      b <- M.new (n + 1) -- visitados
      let f !z !a = if z > n then return [] else do
                      M.write b z True
                      -- actualizamos índice menor no visitado
                      let aa a = if a > n then return a else do
                                   u <- M.read b a
                                   if u then aa (a + 1) else return a
                      -- buscamos siguiente
                      z' <- minimum <$> mapM (p -> let mm !i = if i > n then return i else do
                                                                  u <- M.read b i
                                                                  if not u
                                                                    then return i
                                                                    else mm (i + p)
                                                    in  mm (p * (a `divSup` p))) (primes p z)
                      a' <- aa a
                      (z:) <$> f z' a'
      f 2 3
     
    -- para cierta criba hasta n obtendremos el mayor posible
    main = do
      [n] <- map read <$> getArgs
      v <- sucECG n
      let s = length v - 1
      putStrLn $ show s ++ " " ++ show (v!!s)
     
    {-
     
    $ crono ./ecg 300000
    194561 100003
    Mem: 21640 kbytes. Time: 0:10.13
     
    -- usando la solución con Set
    $ crono ./ecg1 2 194561
    100003
    Mem: 21120 kbytes. Time: 1:19.52
     
    -}

Escribe tu solución

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.