abril 2016 – Exercitium

Operaciones con series de potencias

PorJosé A. Alonso 29-abril-201625-marzo-2022

Una serie de potencias es una serie de la forma

   a₀ + a₁x + a₂x² + a₃x³ + ...

1	a₀ + a₁x + a₂x² + a₃x³ + ...

Las series de potencias se pueden representar mediante listas
infinitas. Por ejemplo, la serie de la función exponencial es

   e^x = 1 + x + x²/2! + x³/3! + ...

1	e^x = 1 + x + x²/2! + x³/3! + ...

y se puede representar por [1, 1, 1/2, 1/6, 1/24, 1/120, …]

Las operaciones con series se pueden ver como una generalización de las de los polinomios.

En lo que sigue, usaremos el tipo (Serie a) para representar las series de potencias con coeficientes en a y su definición es

   type Serie a = [a]

1	type Serie a = [a]

Definir las siguientes funciones

   opuesta      :: Num a => Serie a -> Serie a
   suma         :: Num a => Serie a -> Serie a -> Serie a
   resta        :: Num a => Serie a -> Serie a -> Serie a
   producto     :: Num a => Serie a -> Serie a -> Serie a
   cociente     :: Fractional a => Serie a -> Serie a -> Serie a
   derivada     :: (Num a, Enum a) => Serie a -> Serie a
   integral     :: (Fractional a, Enum a) => Serie a -> Serie a
   expx         :: Serie Rational

opuesta :: Num a => Serie a -> Serie a

suma :: Num a => Serie a -> Serie a -> Serie a

resta :: Num a => Serie a -> Serie a -> Serie a

producto :: Num a => Serie a -> Serie a -> Serie a

cociente :: Fractional a => Serie a -> Serie a -> Serie a

derivada :: (Num a, Enum a) => Serie a -> Serie a

integral :: (Fractional a, Enum a) => Serie a -> Serie a

expx :: Serie Rational

tales que

(opuesta xs) es la opuesta de la serie xs. Por ejemplo,

     λ> take 7 (opuesta [-6,-4..])
     [6,4,2,0,-2,-4,-6]

1 2	λ> take 7 (opuesta [-6,-4..]) [6,4,2,0,-2,-4,-6]

(suma xs ys) es la suma de las series xs e ys. Por ejemplo,

     λ> take 7 (suma [1,3..] [2,4..])
     [3,7,11,15,19,23,27]

1 2	λ> take 7 (suma [1,3..] [2,4..]) [3,7,11,15,19,23,27]

(resta xs ys) es la resta de las series xs es ys. Por ejemplo,

     λ> take 7 (resta [3,5..] [2,4..])
     [1,1,1,1,1,1,1]
     λ> take 7 (resta ([3,7,11,15,19,23,27] ++ repeat 0) [1,3..])
     [2,4,6,8,10,12,14]

λ> take 7 (resta [3,5..] [2,4..])

[1,1,1,1,1,1,1]

λ> take 7 (resta ([3,7,11,15,19,23,27] ++ repeat 0) [1,3..])

[2,4,6,8,10,12,14]

(producto xs ys) es el producto de las series xs e ys. Por ejemplo,

     λ> take 7 (producto [3,5..] [2,4..])
     [6,22,52,100,170,266,392]

1 2	λ> take 7 (producto [3,5..] [2,4..]) [6,22,52,100,170,266,392]

(cociente xs ys) es el cociente de las series xs e ys. Por ejemplo,

     λ> take 7 (cociente ([6,22,52,100,170,266,392] ++ repeat 0) [3,5..])
     [2.0,4.0,6.0,8.0,10.0,12.0,14.0]

1 2	λ> take 7 (cociente ([6,22,52,100,170,266,392] ++ repeat 0) [3,5..]) [2.0,4.0,6.0,8.0,10.0,12.0,14.0]

(derivada xs) es la derivada de la serie xs. Por ejemplo,

     λ> take 7 (derivada [2,4..])
     [4,12,24,40,60,84,112]

1 2	λ> take 7 (derivada [2,4..]) [4,12,24,40,60,84,112]

(integral xs) es la integral de la serie xs. Por ejemplo,

     λ> take 7 (integral ([4,12,24,40,60,84,112] ++ repeat 0))
     [0.0,4.0,6.0,8.0,10.0,12.0,14.0]

1 2	λ> take 7 (integral ([4,12,24,40,60,84,112] ++ repeat 0)) [0.0,4.0,6.0,8.0,10.0,12.0,14.0]

expx es la serie de la función exponencial. Por ejemplo,

     λ> take 8 expx
     [1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]
     λ> take 8 (derivada expx)
     [1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]
     λ> take 8 (integral expx)
     [0 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

λ> take 8 expx

[1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

λ> take 8 (derivada expx)

[1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

λ> take 8 (integral expx)

[0 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

Soluciones

type Serie a = [a] 

opuesta :: Num a => Serie a -> Serie a
opuesta = map negate

suma :: Num a => Serie a -> Serie a -> Serie a
suma = zipWith (+)

resta :: Num a => Serie a -> Serie a -> Serie a
resta xs ys = suma xs (opuesta ys)

producto :: Num a => Serie a -> Serie a -> Serie a
producto (x:xs) zs@(y:ys) = 
    x*y : suma (producto xs zs) (map (x*) ys)

cociente :: Fractional a => Serie a -> Serie a -> Serie a
cociente (x:xs) (y:ys) = zs 
    where zs = x/y : map (/y) (resta xs (producto zs ys))  

derivada :: (Num a, Enum a) => Serie a -> Serie a
derivada (_:xs) = zipWith (*) xs [1..]

integral :: (Fractional a, Enum a) => Serie a -> Serie a
integral xs = 0 : zipWith (/) xs [1..]

expx :: Serie Rational
expx = map (1/) (map fromIntegral factoriales)

-- factoriales es la lista de los factoriales. Por ejemplo, 
--    take 7 factoriales  ==  [1,1,2,6,24,120,720]
factoriales :: [Integer]
factoriales = 1 : scanl1 (*) [1..]

type Serie a = [a]

opuesta :: Num a => Serie a -> Serie a

opuesta = map negate

suma :: Num a => Serie a -> Serie a -> Serie a

suma = zipWith (+)

resta :: Num a => Serie a -> Serie a -> Serie a

resta xs ys = suma xs (opuesta ys)

producto :: Num a => Serie a -> Serie a -> Serie a

producto (x:xs) zs@(y:ys) =

x*y : suma (producto xs zs) (map (x*) ys)

cociente :: Fractional a => Serie a -> Serie a -> Serie a

cociente (x:xs) (y:ys) = zs

where zs = x/y : map (/y) (resta xs (producto zs ys))

derivada :: (Num a, Enum a) => Serie a -> Serie a

derivada (_:xs) = zipWith (*) xs [1..]

integral :: (Fractional a, Enum a) => Serie a -> Serie a

integral xs = 0 : zipWith (/) xs [1..]

expx :: Serie Rational

expx = map (1/) (map fromIntegral factoriales)

-- factoriales es la lista de los factoriales. Por ejemplo,

-- take 7 factoriales == [1,1,2,6,24,120,720]

factoriales :: [Integer]

factoriales = 1 : scanl1 (*) [1..]

Mínimo número de cambios para igualar una lista

PorJosé A. Alonso 27-abril-20164-mayo-2016

Definir la función

   nMinimoCambios :: Ord a => [a] -> Int

1	nMinimoCambios :: Ord a => [a] -> Int

tal que (nMinimoCambios xs) es el menor número de elementos de xs que hay que cambiar para que todos sean iguales. Por ejemplo,

   nMinimoCambios [3,5,3,7,9,6]      ==  4
   nMinimoCambios [3,5,3,7,3,3]      ==  2
   nMinimoCambios "Salamanca"        ==  5
   nMinimoCambios (4 : [1..500000])  ==  499999

nMinimoCambios [3,5,3,7,9,6] == 4

nMinimoCambios [3,5,3,7,3,3] == 2

nMinimoCambios "Salamanca" == 5

nMinimoCambios (4 : [1..500000]) == 499999

En el primer ejemplo, los elementos que hay que cambiar son 5, 7, 9 y 6.

Soluciones

import Data.List (group, nub, sort)

-- 1ª definición
-- =============

nMinimoCambios1 :: Ord a => [a] -> Int
nMinimoCambios1 xs = 
    length xs - fst (last (sort (frecuencias xs)))

-- (frecuencias xs) es la lista de los pares de los elementos de xs y el
-- número de veces que ocurren en xs. Por ejemplo,
--    frecuencias [3,5,3,7,9,6]  ==  [(2,3),(1,5),(1,7),(1,9),(1,6)]
--    frecuencias [3,5,3,7,5,5]  ==  [(2,3),(3,5),(1,7)]
frecuencias :: Ord a => [a] -> [(Int,a)]
frecuencias xs = [(cuenta x xs,x) | x <- nub xs]

-- (cuenta x ys) es el número de veces que ocurre x en ys. Por ejemplo,
--    cuenta 3 [3,5,3,7,9,6]  ==  2
cuenta :: Ord a => a -> [a] -> Int
cuenta x ys = length [y | y <- ys, y == x]

-- 2ª definición
-- =============

nMinimoCambios2 :: Ord a => [a] -> Int
nMinimoCambios2 xs = 
    length xs - fst (last (sort (frecuencias2 xs)))

-- (frecuencias2 xs) es la lista de los pares de los elementos de xs y el
-- número de veces que ocurren en xs. Por ejemplo,
--    frecuencias2 [3,5,3,7,9,6]  ==  [(2,3),(1,5),(1,7),(1,9),(1,6)]
--    frecuencias2 [3,5,3,7,5,5]  ==  [(2,3),(3,5),(1,7)]
frecuencias2 :: Ord a  => [a] -> [(Int,a)]
frecuencias2 xs = 
    [(1 + length ys, y) | (y:ys) <- group (sort xs)] 

-- 3ª definición
-- =============

nMinimoCambios3 :: Ord a => [a] -> Int
nMinimoCambios3 xs = sum ys - maximum ys
    where ys = [length ys | ys <- group (sort xs)]

import Data.List (group, nub, sort)

-- 1ª definición

-- =============

nMinimoCambios1 :: Ord a => [a] -> Int

nMinimoCambios1 xs =

length xs - fst (last (sort (frecuencias xs)))

-- (frecuencias xs) es la lista de los pares de los elementos de xs y el

-- número de veces que ocurren en xs. Por ejemplo,

-- frecuencias [3,5,3,7,9,6] == [(2,3),(1,5),(1,7),(1,9),(1,6)]

-- frecuencias [3,5,3,7,5,5] == [(2,3),(3,5),(1,7)]

frecuencias :: Ord a => [a] -> [(Int,a)]

frecuencias xs = [(cuenta x xs,x) | x <- nub xs]

-- (cuenta x ys) es el número de veces que ocurre x en ys. Por ejemplo,

-- cuenta 3 [3,5,3,7,9,6] == 2

cuenta :: Ord a => a -> [a] -> Int

cuenta x ys = length [y | y <- ys, y == x]

-- 2ª definición

-- =============

nMinimoCambios2 :: Ord a => [a] -> Int

nMinimoCambios2 xs =

length xs - fst (last (sort (frecuencias2 xs)))

-- (frecuencias2 xs) es la lista de los pares de los elementos de xs y el

-- número de veces que ocurren en xs. Por ejemplo,

-- frecuencias2 [3,5,3,7,9,6] == [(2,3),(1,5),(1,7),(1,9),(1,6)]

-- frecuencias2 [3,5,3,7,5,5] == [(2,3),(3,5),(1,7)]

frecuencias2 :: Ord a => [a] -> [(Int,a)]

frecuencias2 xs =

[(1 + length ys, y) | (y:ys) <- group (sort xs)]

-- 3ª definición

-- =============

nMinimoCambios3 :: Ord a => [a] -> Int

nMinimoCambios3 xs = sum ys - maximum ys

where ys = [length ys | ys <- group (sort xs)]

Productos simultáneos de dos y tres números consecutivos

PorJosé A. Alonso 25-abril-201625-abril-2021

Definir la función

   productos :: Integer -> Integer -> [[Integer]]

1	productos :: Integer -> Integer -> [[Integer]]

tal que (productos n x) es las listas de n elementos consecutivos cuyo producto es x. Por ejemplo,

   productos 2 6     ==  [[2,3]]
   productos 3 6     ==  [[1,2,3]]
   productos 4 1680  ==  [[5,6,7,8]]
   productos 2 5     ==  []

productos 2 6 == [[2,3]]

productos 3 6 == [[1,2,3]]

productos 4 1680 == [[5,6,7,8]]

productos 2 5 == []

Comprobar con QuickCheck que si n > 0 y x > 0, entonces

   productos n (product [x..x+n-1]) == [[x..x+n-1]]

1	productos n (product [x..x+n-1]) == [[x..x+n-1]]

Usando productos, definir la función

   productosDe2y3consecutivos :: [Integer]

1	productosDe2y3consecutivos :: [Integer]

cuyos elementos son los números naturales (no nulos) que pueden expresarse simultáneamente como producto de dos y tres números consecutivos. Por ejemplo,

   head productosDe2y3consecutivos  ==  6

1	head productosDe2y3consecutivos == 6

Nota. Según demostró Mordell en 1962, productosDe2y3consecutivos sólo tiene dos elementos.

Soluciones

import Test.QuickCheck

-- 1ª definición
productos1 :: Integer -> Integer -> [[Integer]]
productos1 n x = [[y..y+n-1] | y <- [1..x],
                               product [y..y+n-1] == x]

-- 2ª definición
-- =============

-- Se puede reducir el intervalo de búsqueda teniendo en cuenta las
-- siguientes desigualdades
--    y*(y+1)* ... (y+n-1) = x
--    y^n <= x <= (y+n-1)^n
--    y <= x^(1/n), x^(1/n)-n+1 <= y
--    x^(1/n)-n+1 <= y <= x^(1/n)

productos2 :: Integer -> Integer -> [[Integer]]
productos2 n x = [[z..z+n-1] | z <- [y-n+1..y],
                               product [z..z+n-1] == x]
    where y = floor ((fromIntegral x)**(1/(fromIntegral n)))

productos :: Integer -> Integer -> [[Integer]]
productos = productos2

prop_productos n x =
    n > 0 && x > 0 ==> productos n (product [x..x+n-1]) == [[x..x+n-1]]

-- La comprobación es
--    ghci> quickCheck prop_productos
--    +++ OK, passed 100 tests.
--    (0.10 secs, 26409644 bytes)

productosDe2y3consecutivos :: [Integer]
productosDe2y3consecutivos = [x| x <- [1..], 
                                 let ys = productos 2 x,
                                 not (null ys), 
                                 let zs = productos 3 x,
                                 not (null zs)] 

-- El cálculo es
--    ghci> take 2 productosDe2y3consecutivos
--    [6,210]
--    ghci> productos 2 210
--    [[14,15]]
--    ghci> productos 3 210
--    [[5,6,7]]

import Test.QuickCheck

-- 1ª definición

productos1 :: Integer -> Integer -> [[Integer]]

productos1 n x = [[y..y+n-1] | y <- [1..x],

product [y..y+n-1] == x]

-- 2ª definición

-- =============

-- Se puede reducir el intervalo de búsqueda teniendo en cuenta las

-- siguientes desigualdades

-- y*(y+1)* ... (y+n-1) = x

-- y^n <= x <= (y+n-1)^n

-- y <= x^(1/n), x^(1/n)-n+1 <= y

-- x^(1/n)-n+1 <= y <= x^(1/n)

productos2 :: Integer -> Integer -> [[Integer]]

productos2 n x = [[z..z+n-1] | z <- [y-n+1..y],

product [z..z+n-1] == x]

where y = floor ((fromIntegral x)**(1/(fromIntegral n)))

productos :: Integer -> Integer -> [[Integer]]

productos = productos2

prop_productos n x =

n > 0 && x > 0 ==> productos n (product [x..x+n-1]) == [[x..x+n-1]]

-- La comprobación es

-- ghci> quickCheck prop_productos

-- +++ OK, passed 100 tests.

-- (0.10 secs, 26409644 bytes)

productosDe2y3consecutivos :: [Integer]

productosDe2y3consecutivos = [x| x <- [1..],

let ys = productos 2 x,

not (null ys),

let zs = productos 3 x,

not (null zs)]

-- El cálculo es

-- ghci> take 2 productosDe2y3consecutivos

-- [6,210]

-- ghci> productos 2 210

-- [[14,15]]

-- ghci> productos 3 210

-- [[5,6,7]]

Conflictos de horarios

PorJosé A. Alonso 22-abril-201629-abril-2016

Los horarios de los cursos se pueden representar mediante matrices donde las filas indican los curso, las columnas las horas de clase y el valor correspondiente al curso i y la hora j es verdadero indica que i tiene clase a la hora j.

En Haskell, podemos usar la matrices de la librería Data.Matrix y definir el tipo de los horarios por

   type Horario = Matrix Bool

1	type Horario = Matrix Bool

Un ejemplo de horario es

   ejHorarios1 :: Horario
   ejHorarios1 = fromLists [[True,  True,  False, False],
                            [False, True,  True,  False],
                            [False, False, True,  True]]

ejHorarios1 :: Horario

ejHorarios1 = fromLists [[True, True, False, False],

[False, True, True, False],

[False, False, True, True]]

en el que el 1º curso tiene clase a la 1ª y 2ª hora, el 2º a la 2ª y a la 3ª y el 3º a la 3ª y a la 4ª.

Definir la función

   cursosConflictivos :: Horario -> [Int] -> Bool

1	cursosConflictivos :: Horario -> [Int] -> Bool

tal que (cursosConflictivos h is) se verifica para si los cursos de la lista is hay alguna hora en la que más de uno tiene clase a dicha hora. Por ejemplo,

   cursosConflictivos ejHorarios1 [1,2]  ==  True
   cursosConflictivos ejHorarios1 [1,3]  ==  False

1 2	cursosConflictivos ejHorarios1 [1,2] == True cursosConflictivos ejHorarios1 [1,3] == False

Soluciones

import Data.Matrix

type Horario = Matrix Bool

ejHorarios1 :: Horario
ejHorarios1 = fromLists [[True,  True,  False, False],
                         [False, True,  True,  False],
                         [False, False, True,  True]]

--    horaConflictiva ejHorarios1 [1,2]   2  ==  True
--    horaConflictiva ejHorarios1 [1,3]   2  ==  False
--    horaConflictiva ejHorarios1 [1,2,3] 1  ==  False
horaConflictiva :: Horario -> [Int] -> Int -> Bool
horaConflictiva h is j =
  length [i | i <- is, h ! (i,j)] > 1

cursosConflictivos :: Horario -> [Int] -> Bool
cursosConflictivos h is =
  or [horaConflictiva h is j | j <- [1..n]]
  where n = ncols h

import Data.Matrix

type Horario = Matrix Bool

ejHorarios1 :: Horario

ejHorarios1 = fromLists [[True, True, False, False],

[False, True, True, False],

[False, False, True, True]]

-- horaConflictiva ejHorarios1 [1,2] 2 == True

-- horaConflictiva ejHorarios1 [1,3] 2 == False

-- horaConflictiva ejHorarios1 [1,2,3] 1 == False

horaConflictiva :: Horario -> [Int] -> Int -> Bool

horaConflictiva h is j =

length [i | i <- is, h ! (i,j)] > 1

cursosConflictivos :: Horario -> [Int] -> Bool

cursosConflictivos h is =

or [horaConflictiva h is j | j <- [1..n]]

where n = ncols h

Mínima diferencia entre elementos de una lista

PorJosé A. Alonso 18-abril-201625-abril-2016

Definir la función

   minimaDiferencia :: (Num a, Ord a) => [a] -> a

1	minimaDiferencia :: (Num a, Ord a) => [a] -> a

tal que (minimaDiferencia xs) es el menor valor absoluto de las diferencias entre todos los pares de elementos de xs (que se supone que tiene al menos 2 elementos). Por ejemplo,

   minimaDiferencia [1,5,3,19,18,25]  ==  1
   minimaDiferencia [30,5,20,9]       ==  4
   minimaDiferencia [30,5,20,9,5]     ==  0
   minimaDiferencia [1..10^6]         ==  1

minimaDiferencia [1,5,3,19,18,25] == 1

minimaDiferencia [30,5,20,9] == 4

minimaDiferencia [30,5,20,9,5] == 0

minimaDiferencia [1..10^6] == 1

En el primer ejemplo la menor diferencia es 1 y se da entre los elementos 19 y 18; en el 2ª es 4 entre los elementos 5 y 9 y en la 3ª es 0 porque el elemento 5 está repetido.

Soluciones

import Data.List (delete, sort)

-- 1ª definición    
minimaDiferencia1 :: (Num a, Ord a) => [a] -> a
minimaDiferencia1 xs =
    minimum [abs (x - y) | x <- xs, y <- delete x xs]

-- 2ª definición    
minimaDiferencia2 :: (Num a, Ord a) => [a] -> a
minimaDiferencia2 [x,y]  = abs (x - y)
minimaDiferencia2 (x:xs) = min (minimum [abs (x - y) | y <- xs])
                               (minimaDiferencia2 xs)      

-- 3ª definición    
minimaDiferencia3 :: (Num a, Ord a) => [a] -> a
minimaDiferencia3 xs =
    minimum [y - x | (x,y) <- zip ys (tail ys)]
    where ys = sort xs

-- Comparación de eficiencia
--    λ> minimaDiferencia1 [1..2*10^3]
--    1
--    (3.47 secs, 1,419,790,688 bytes)
--    λ> minimaDiferencia2 [1..2*10^3]
--    1
--    (0.78 secs, 684,030,712 bytes)
--    λ> minimaDiferencia3 [1..2*10^3]
--    1
--    (0.02 secs, 0 bytes)
--    
--    λ> minimaDiferencia2 [1..5*10^3]
--    1
--    (5.14 secs, 4,233,952,776 bytes)
--    λ> minimaDiferencia3 [1..5*10^3]
--    1
--    (0.02 secs, 0 bytes)

import Data.List (delete, sort)

-- 1ª definición

minimaDiferencia1 :: (Num a, Ord a) => [a] -> a

minimaDiferencia1 xs =

minimum [abs (x - y) | x <- xs, y <- delete x xs]

-- 2ª definición

minimaDiferencia2 :: (Num a, Ord a) => [a] -> a

minimaDiferencia2 [x,y] = abs (x - y)

minimaDiferencia2 (x:xs) = min (minimum [abs (x - y) | y <- xs])

(minimaDiferencia2 xs)

-- 3ª definición

minimaDiferencia3 :: (Num a, Ord a) => [a] -> a

minimaDiferencia3 xs =

minimum [y - x | (x,y) <- zip ys (tail ys)]

where ys = sort xs

-- Comparación de eficiencia

-- λ> minimaDiferencia1 [1..2*10^3]

-- 1

-- (3.47 secs, 1,419,790,688 bytes)

-- λ> minimaDiferencia2 [1..2*10^3]

-- 1

-- (0.78 secs, 684,030,712 bytes)

-- λ> minimaDiferencia3 [1..2*10^3]

-- 1

-- (0.02 secs, 0 bytes)

-- λ> minimaDiferencia2 [1..5*10^3]

-- 1

-- (5.14 secs, 4,233,952,776 bytes)

-- λ> minimaDiferencia3 [1..5*10^3]

-- 1

-- (0.02 secs, 0 bytes)

Números como sumas de primos consecutivos

PorJosé A. Alonso 15-abril-201625-abril-2021

En el artículo Integers as a sum of consecutive primes in 2,3,4,.. ways se presentan números que se pueden escribir como sumas de primos consecutivos de varias formas. Por ejemplo, el 41 se puede escribir de dos formas distintas

   41 =  2 +  3 +  5 + 7 + 11 + 13
   41 = 11 + 13 + 17

1 2	41 = 2 + 3 + 5 + 7 + 11 + 13 41 = 11 + 13 + 17

el 240 se puede escribir de tres formas

   240 =  17 +  19 + 23 + 29 + 31 + 37 + 41 + 43
   240 =  53 +  59 + 61 + 67
   240 = 113 + 127

240 = 17 + 19 + 23 + 29 + 31 + 37 + 41 + 43

240 = 53 + 59 + 61 + 67

240 = 113 + 127

y el 311 se puede escribir de 4 formas

   311 =  11 +  13 +  17 + 19 + 23 + 29 + 31 + 37 + 41 + 43 + 47
   311 =  31 +  37 +  41 + 43 + 47 + 53 + 59
   311 =  53 +  59 +  61 + 67 + 71
   311 = 101 + 103 + 107

311 = 11 + 13 + 17 + 19 + 23 + 29 + 31 + 37 + 41 + 43 + 47

311 = 31 + 37 + 41 + 43 + 47 + 53 + 59

311 = 53 + 59 + 61 + 67 + 71

311 = 101 + 103 + 107

Definir la función

   sumas :: Integer -> [[Integer]]

1	sumas :: Integer -> [[Integer]]

tal que (sumas x) es la lista de las formas de escribir x como suma
de dos o más números primos consecutivos. Por ejemplo,

   ghci> sumas 41
   [[2,3,5,7,11,13],[11,13,17]]
   ghci> sumas 240
   [[17,19,23,29,31,37,41,43],[53,59,61,67],[113,127]]
   ghci> sumas 311
   [[11,13,17,19,23,29,31,37,41,43,47],[31,37,41,43,47,53,59],[53,59,61,67,71],[101,103,107]]
   ghci> maximum [length (sumas n) | n <- [1..600]]
   4

ghci> sumas 41

[[2,3,5,7,11,13],[11,13,17]]

ghci> sumas 240

[[17,19,23,29,31,37,41,43],[53,59,61,67],[113,127]]

ghci> sumas 311

[[11,13,17,19,23,29,31,37,41,43,47],[31,37,41,43,47,53,59],[53,59,61,67,71],[101,103,107]]

ghci> maximum [length (sumas n) | n <- [1..600]]

Soluciones

import Data.Numbers.Primes (primes)
import Data.List (span)

sumas :: Integer -> [[Integer]]
sumas x = [ys | n <- takeWhile (< x) primes, 
                let ys = sumaDesde x n,
                not (null ys)]

-- (sumaDesde x n) es la lista de al menos dos números primos
-- consecutivos a partir del número primo n cuya suma es x, si existen y
-- la lista vacía en caso contrario. Por ejemplo,
--    sumaDesde 15 3  ==  [3,5,7]
--    sumaDesde  7 3  ==  []
sumaDesde :: Integer -> Integer -> [Integer]
sumaDesde x n | x == y    = take (1 + length us) ys
              | otherwise = []
    where ys       = dropWhile (<n) primes
          (us,y:_) = span (<x) (scanl1 (+) ys)

import Data.Numbers.Primes (primes)

import Data.List (span)

sumas :: Integer -> [[Integer]]

sumas x = [ys | n <- takeWhile (< x) primes,

let ys = sumaDesde x n,

not (null ys)]

-- (sumaDesde x n) es la lista de al menos dos números primos

-- consecutivos a partir del número primo n cuya suma es x, si existen y

-- la lista vacía en caso contrario. Por ejemplo,

-- sumaDesde 15 3 == [3,5,7]

-- sumaDesde 7 3 == []

sumaDesde :: Integer -> Integer -> [Integer]

sumaDesde x n | x == y = take (1 + length us) ys

| otherwise = []

where ys = dropWhile (<n) primes

(us,y:_) = span (<x) (scanl1 (+) ys)

Caminos reducidos

PorJosé A. Alonso 12-abril-201625-abril-2021

Un camino es una sucesión de pasos en una de las cuatros direcciones Norte, Sur, Este, Oeste. Ir en una dirección y a continuación en la opuesta es un esfuerzo que se puede reducir, Por ejemplo, el camino [Norte,Sur,Este,Sur] se puede reducir a [Este,Sur].

Un camino se dice que es reducido si no tiene dos pasos consecutivos en direcciones opuesta. Por ejemplo, [Este,Sur] es reducido y [Norte,Sur,Este,Sur] no lo es.

En Haskell, las direcciones y los caminos se pueden definir por

   data Direccion = N | S | E | O deriving (Show, Eq)
   type Camino = [Direccion]

1 2	data Direccion = N \| S \| E \| O deriving (Show, Eq) type Camino = [Direccion]

Definir la función

   reducido :: Camino -> Camino

1	reducido :: Camino -> Camino

tal que (reducido ds) es el camino reducido equivalente al camino ds. Por ejemplo,

   reducido []                              ==  []
   reducido [N]                             ==  [N]
   reducido [N,O]                           ==  [N,O]
   reducido [N,O,E]                         ==  [N]
   reducido [N,O,E,S]                       ==  [] 
   reducido [N,O,S,E]                       ==  [N,O,S,E]
   reducido [S,S,S,N,N,N]                   ==  []
   reducido [N,S,S,E,O,N]                   ==  []
   reducido [N,S,S,E,O,N,O]                 ==  [O]
   reducido (take (10^7) (cycle [N,E,O,S])) ==  []

reducido [] == []

reducido [N] == [N]

reducido [N,O] == [N,O]

reducido [N,O,E] == [N]

reducido [N,O,E,S] == []

reducido [N,O,S,E] == [N,O,S,E]

reducido [S,S,S,N,N,N] == []

reducido [N,S,S,E,O,N] == []

reducido [N,S,S,E,O,N,O] == [O]

reducido (take (10^7) (cycle [N,E,O,S])) == []

Nótese que en el penúltimo ejemplo las reducciones son

       [N,S,S,E,O,N,O]  
   --> [S,E,O,N,O]  
   --> [S,N,O]  
   --> [O]

[N,S,S,E,O,N,O]

--> [S,E,O,N,O]

--> [S,N,O]

--> [O]

Soluciones

data Direccion = N | S | E | O deriving (Show, Eq)

type Camino = [Direccion]

-- 1ª solución (por recursión):
reducido1 :: Camino -> Camino
reducido1 [] = []
reducido1 (d:ds) | null ds'                = [d]
                 | d == opuesta (head ds') = tail ds'
                 | otherwise               = d:ds'
    where ds' = reducido1 ds

opuesta :: Direccion -> Direccion
opuesta N = S
opuesta S = N
opuesta E = O
opuesta O = E

-- 2ª solución (por plegado)
reducido2 :: Camino -> Camino
reducido2 = foldr aux []
    where aux N (S:xs) = xs
          aux S (N:xs) = xs
          aux E (O:xs) = xs
          aux O (E:xs) = xs
          aux x xs     = x:xs

-- 3ª solución 
reducido3 :: Camino -> Camino
reducido3 []       = []
reducido3 (N:S:ds) = reducido3 ds
reducido3 (S:N:ds) = reducido3 ds
reducido3 (E:O:ds) = reducido3 ds
reducido3 (O:E:ds) = reducido3 ds
reducido3 (d:ds) | null ds'                = [d]
                 | d == opuesta (head ds') = tail ds'
                 | otherwise               = d:ds'
    where ds' = reducido3 ds

-- 4ª solución
reducido4 :: Camino -> Camino
reducido4 ds = reverse (aux ([],ds)) where 
    aux (N:xs, S:ys) = aux (xs,ys)
    aux (S:xs, N:ys) = aux (xs,ys)
    aux (E:xs, O:ys) = aux (xs,ys)
    aux (O:xs, E:ys) = aux (xs,ys)
    aux (  xs, y:ys) = aux (y:xs,ys)
    aux (  xs,   []) = xs

-- Comparación de eficiencia
--    ghci> reducido1 (take (10^6) (cycle [N,E,O,S]))
--    []
--    (3.87 secs, 460160736 bytes)
--    ghci> reducido2 (take (10^6) (cycle [N,E,O,S]))
--    []
--    (1.16 secs, 216582880 bytes)
--    ghci> reducido3 (take (10^6) (cycle [N,E,O,S]))
--    []
--    (0.58 secs, 98561872 bytes)
--    ghci> reducido4 (take (10^6) (cycle [N,E,O,S]))
--    []
--    (0.64 secs, 176154640 bytes)
--    
--    ghci> reducido3 (take (10^7) (cycle [N,E,O,S]))
--    []
--    (5.43 secs, 962694784 bytes)
--    ghci> reducido4 (take (10^7) (cycle [N,E,O,S]))
--    []
--    (9.29 secs, 1722601528 bytes)
-- 
--    ghci> length $ reducido3 (take 2000000 $ cycle [N,O,N,S,E,N,S,O,S,S])
--    400002
--    (4.52 secs, 547004960 bytes)
--    ghci> length $ reducido4 (take 2000000 $ cycle [N,O,N,S,E,N,S,O,S,S])
--    400002
--    (2.17 secs, 379049224 bytes)

data Direccion = N | S | E | O deriving (Show, Eq)

type Camino = [Direccion]

-- 1ª solución (por recursión):

reducido1 :: Camino -> Camino

reducido1 [] = []

reducido1 (d:ds) | null ds' = [d]

| d == opuesta (head ds') = tail ds'

| otherwise = d:ds'

where ds' = reducido1 ds

opuesta :: Direccion -> Direccion

opuesta N = S

opuesta S = N

opuesta E = O

opuesta O = E

-- 2ª solución (por plegado)

reducido2 :: Camino -> Camino

reducido2 = foldr aux []

where aux N (S:xs) = xs

aux S (N:xs) = xs

aux E (O:xs) = xs

aux O (E:xs) = xs

aux x xs = x:xs

-- 3ª solución

reducido3 :: Camino -> Camino

reducido3 [] = []

reducido3 (N:S:ds) = reducido3 ds

reducido3 (S:N:ds) = reducido3 ds

reducido3 (E:O:ds) = reducido3 ds

reducido3 (O:E:ds) = reducido3 ds

reducido3 (d:ds) | null ds' = [d]

| d == opuesta (head ds') = tail ds'

| otherwise = d:ds'

where ds' = reducido3 ds

-- 4ª solución

reducido4 :: Camino -> Camino

reducido4 ds = reverse (aux ([],ds)) where

aux (N:xs, S:ys) = aux (xs,ys)

aux (S:xs, N:ys) = aux (xs,ys)

aux (E:xs, O:ys) = aux (xs,ys)

aux (O:xs, E:ys) = aux (xs,ys)

aux ( xs, y:ys) = aux (y:xs,ys)

aux ( xs, []) = xs

-- Comparación de eficiencia

-- ghci> reducido1 (take (10^6) (cycle [N,E,O,S]))

-- []

-- (3.87 secs, 460160736 bytes)

-- ghci> reducido2 (take (10^6) (cycle [N,E,O,S]))

-- []

-- (1.16 secs, 216582880 bytes)

-- ghci> reducido3 (take (10^6) (cycle [N,E,O,S]))

-- []

-- (0.58 secs, 98561872 bytes)

-- ghci> reducido4 (take (10^6) (cycle [N,E,O,S]))

-- []

-- (0.64 secs, 176154640 bytes)

-- ghci> reducido3 (take (10^7) (cycle [N,E,O,S]))

-- []

-- (5.43 secs, 962694784 bytes)

-- ghci> reducido4 (take (10^7) (cycle [N,E,O,S]))

-- []

-- (9.29 secs, 1722601528 bytes)

-- ghci> length $ reducido3 (take 2000000 $ cycle [N,O,N,S,E,N,S,O,S,S])

-- 400002

-- (4.52 secs, 547004960 bytes)

-- ghci> length $ reducido4 (take 2000000 $ cycle [N,O,N,S,E,N,S,O,S,S])

-- 400002

-- (2.17 secs, 379049224 bytes)

Primos permutables

PorJosé A. Alonso 11-abril-201618-abril-2016

Un primo permutable es un número primo tal que todos los números obtenidos permutando sus cifras son primos. Por ejemplo, 337 es un primo permutable ya que 337, 373 y 733 son primos.

Definir las funciones

   esPrimoPermutable :: Integer -> Bool
   primosPermutables :: [Integer]

1 2	esPrimoPermutable :: Integer -> Bool primosPermutables :: [Integer]

tales que

(esPrimoPermutable x) se verifica si x es un primo permutable. Por ejemplo,

     esPrimoPermutable 97  ==  True
     esPrimoPermutable 337 ==  True
     esPrimoPermutable 23  ==  False

esPrimoPermutable 97 == True

esPrimoPermutable 337 == True

esPrimoPermutable 23 == False

primosPermutables es la lista de los primos permutables. Por ejemplo,

     λ> take 20 primosPermutables
     [2,3,5,7,11,13,17,31,37,71,73,79,97,113,131,199,311,337,373,733]

1 2	λ> take 20 primosPermutables [2,3,5,7,11,13,17,31,37,71,73,79,97,113,131,199,311,337,373,733]

Soluciones

import Data.List (nub, permutations)
import Data.Numbers.Primes (primes, isPrime)

-- 1ª solución
-- ===========
    
esPrimoPermutable :: Integer -> Bool
esPrimoPermutable x =
    all isPrime (permutacionesN x)

permutacionesN :: Integer -> [Integer]
permutacionesN x =
    [read ys | ys <- nub (permutations xs)]
    where xs = show x
          n  = length xs
        
primosPermutables :: [Integer]
primosPermutables =
    [x | x <- primes, esPrimoPermutable x]

-- 2ª solución
-- ===========

esPrimoPermutable2 :: Integer -> Bool
esPrimoPermutable2 = all isPrime . map read . nub . permutations . show
 
primosPermutables2 :: [Integer]
primosPermutables2 = filter esPrimoPermutable2 primes

import Data.List (nub, permutations)

import Data.Numbers.Primes (primes, isPrime)

-- 1ª solución

-- ===========

esPrimoPermutable :: Integer -> Bool

esPrimoPermutable x =

all isPrime (permutacionesN x)

permutacionesN :: Integer -> [Integer]

permutacionesN x =

[read ys | ys <- nub (permutations xs)]

where xs = show x

n = length xs

primosPermutables :: [Integer]

primosPermutables =

[x | x <- primes, esPrimoPermutable x]

-- 2ª solución

-- ===========

esPrimoPermutable2 :: Integer -> Bool

esPrimoPermutable2 = all isPrime . map read . nub . permutations . show

primosPermutables2 :: [Integer]

primosPermutables2 = filter esPrimoPermutable2 primes

Referencias

Permutable prime en Wikipedia.
Sucesión A003459 de la OEIS.

Raíz entera

PorJosé A. Alonso 8-abril-20161-mayo-2016

Definir la función

   raizEnt :: Integer -> Integer -> Integer

1	raizEnt :: Integer -> Integer -> Integer

tal que (raizEnt x n) es la raíz entera n-ésima de x; es decir, el mayor número entero y tal que y^n <= x. Por ejemplo,

   raizEnt  8 3      ==  2
   raizEnt  9 3      ==  2
   raizEnt 26 3      ==  2
   raizEnt 27 3      ==  3
   raizEnt (10^50) 2 ==  10000000000000000000000000

raizEnt 8 3 == 2

raizEnt 9 3 == 2

raizEnt 26 3 == 2

raizEnt 27 3 == 3

raizEnt (10^50) 2 == 10000000000000000000000000

Comprobar con QuickCheck que para todo número natural n,

    raizEnt (10^(2*n)) 2 == 10^n

1	raizEnt (10^(2*n)) 2 == 10^n

Soluciones

import Test.QuickCheck

-- 1ª definición
raizEnt1 :: Integer -> Integer -> Integer
raizEnt1 x n =
    last (takeWhile (\y -> y^n <= x) [0..])

-- 2ª definición         
raizEnt2 :: Integer -> Integer -> Integer
raizEnt2 x n =
    floor ((fromIntegral x)**(1 / fromIntegral n))

-- Nota. La definición anterior falla para números grandes. Por ejemplo,
--    λ> raizEnt2 (10^50) 2 == 10^25
--    False
          
-- 3ª definición          
raizEnt3 :: Integer -> Integer -> Integer
raizEnt3 x n = aux (1,x)
    where aux (a,b) | d == x    = c
                    | c == a    = c
                    | d < x     = aux (c,b)
                    | otherwise = aux (a,c) 
              where c = (a+b) `div` 2
                    d = c^n

-- Comparación de eficiencia
--    λ> raizEnt1 (10^14) 2
--    10000000
--    (6.15 secs, 6,539,367,976 bytes)
--    λ> raizEnt2 (10^14) 2
--    10000000
--    (0.00 secs, 0 bytes)
--    λ> raizEnt3 (10^14) 2
--    10000000
--    (0.00 secs, 25,871,944 bytes)
--    
--    λ> raizEnt2 (10^50) 2
--    9999999999999998758486016
--    (0.00 secs, 0 bytes)
--    λ> raizEnt3 (10^50) 2
--    10000000000000000000000000
--    (0.00 secs, 0 bytes)
                        
-- La propiedad es                        
prop_raizEnt :: Integer -> Bool
prop_raizEnt n =
    raizEnt3 (10^(2*m)) 2 == 10^m
    where m = abs n

-- La comprobación es              
--    λ> quickCheck prop_raizEnt
--    +++ OK, passed 100 tests.

import Test.QuickCheck

-- 1ª definición

raizEnt1 :: Integer -> Integer -> Integer

raizEnt1 x n =

last (takeWhile (\y -> y^n <= x) [0..])

-- 2ª definición

raizEnt2 :: Integer -> Integer -> Integer

raizEnt2 x n =

floor ((fromIntegral x)**(1 / fromIntegral n))

-- Nota. La definición anterior falla para números grandes. Por ejemplo,

-- λ> raizEnt2 (10^50) 2 == 10^25

-- False

-- 3ª definición

raizEnt3 :: Integer -> Integer -> Integer

raizEnt3 x n = aux (1,x)

where aux (a,b) | d == x = c

| c == a = c

| d < x = aux (c,b)

| otherwise = aux (a,c)

where c = (a+b) `div` 2

d = c^n

-- Comparación de eficiencia

-- λ> raizEnt1 (10^14) 2

-- 10000000

-- (6.15 secs, 6,539,367,976 bytes)

-- λ> raizEnt2 (10^14) 2

-- 10000000

-- (0.00 secs, 0 bytes)

-- λ> raizEnt3 (10^14) 2

-- 10000000

-- (0.00 secs, 25,871,944 bytes)

-- λ> raizEnt2 (10^50) 2

-- 9999999999999998758486016

-- (0.00 secs, 0 bytes)

-- λ> raizEnt3 (10^50) 2

-- 10000000000000000000000000

-- (0.00 secs, 0 bytes)

-- La propiedad es

prop_raizEnt :: Integer -> Bool

prop_raizEnt n =

raizEnt3 (10^(2*m)) 2 == 10^m

where m = abs n

-- La comprobación es

-- λ> quickCheck prop_raizEnt

-- +++ OK, passed 100 tests.

Soluciones en Maxima

raizEnt (x,n) := inrt (x,n)$

1	raizEnt (x,n) := inrt (x,n)$

Números de Lucas

PorJosé A. Alonso 7-abril-201614-abril-2016

Los números de Lucas son los elementos de la sucesión L(n) definida por

   L(0) = 2
   L(1) = 1
   L(n) = L(n-1) + L(n-2), si n > 1.

L(0) = 2

L(1) = 1

L(n) = L(n-1) + L(n-2), si n > 1.

Los primeros números de Lucas son

   2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, ...

1	2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, ...

Definir las funciones

   nLucas :: Integer -> Integer 
   lucas  :: [Integer]

1 2	nLucas :: Integer -> Integer lucas :: [Integer]

tales que

(nLucas n) es el n-ésimo número de Lucas. Por ejemplo,

   nLucas 5                       ==  11
   nLucas 32                      ==  4870847
   length (show (nLucas (10^5)))  ==  20899

nLucas 5 == 11

nLucas 32 == 4870847

length (show (nLucas (10^5))) == 20899

lucas es la lista de los números de Lucas. Por ejemplo,

   take 11 lucas ==  [2,1,3,4,7,11,18,29,47,76,123]

1	take 11 lucas == [2,1,3,4,7,11,18,29,47,76,123]

Soluciones

import Data.List (genericIndex)

-- 1ª definición
-- =============

nLucas1 :: Integer -> Integer
nLucas1 0 = 2
nLucas1 1 = 1
nLucas1 n = nLucas1 (n-1) + nLucas1 (n-2)

lucas1 :: [Integer]
lucas1 = [nLucas1 n | n <- [0..]]

-- 2ª definición
-- =============

lucas2 :: [Integer]
lucas2 = 2 : 1 : zipWith (+) lucas2 (tail lucas2)

nLucas2 :: Integer -> Integer
nLucas2 n = lucas2 `genericIndex` n

-- 3ª definición
-- =============

lucas3  :: [Integer]
lucas3 = 2 : scanl (+) 1 lucas3
 
nLucas3 :: Integer -> Integer
nLucas3 = genericIndex lucas3
            
-- Comparación de eficiencia
-- =========================

--    λ> nLucas1 32
--    4870847
--    (3.22 secs, 1,467,677,208 bytes)
--    λ> nLucas2 32
--    4870847
--    (0.00 secs, 0 bytes)
--    λ> nLucas3 32
--    4870847
--    (0.00 secs, 0 bytes)
--    
--    λ> length (show (nLucas2 230000))
--    48068
--    (1.57 secs, 2,415,008,392 bytes)
--    λ> length (show (nLucas3 230000))
--    48068
--    (2.08 secs, 2,411,107,352 bytes)

import Data.List (genericIndex)

-- 1ª definición

-- =============

nLucas1 :: Integer -> Integer

nLucas1 0 = 2

nLucas1 1 = 1

nLucas1 n = nLucas1 (n-1) + nLucas1 (n-2)

lucas1 :: [Integer]

lucas1 = [nLucas1 n | n <- [0..]]

-- 2ª definición

-- =============

lucas2 :: [Integer]

lucas2 = 2 : 1 : zipWith (+) lucas2 (tail lucas2)

nLucas2 :: Integer -> Integer

nLucas2 n = lucas2 `genericIndex` n

-- 3ª definición

-- =============

lucas3 :: [Integer]

lucas3 = 2 : scanl (+) 1 lucas3

nLucas3 :: Integer -> Integer

nLucas3 = genericIndex lucas3

-- Comparación de eficiencia

-- =========================

-- λ> nLucas1 32

-- 4870847

-- (3.22 secs, 1,467,677,208 bytes)

-- λ> nLucas2 32

-- 4870847

-- (0.00 secs, 0 bytes)

-- λ> nLucas3 32

-- 4870847

-- (0.00 secs, 0 bytes)

-- λ> length (show (nLucas2 230000))

-- 48068

-- (1.57 secs, 2,415,008,392 bytes)

-- λ> length (show (nLucas3 230000))

-- 48068

-- (2.08 secs, 2,411,107,352 bytes)

Soluciones en Maxima

nLucas (n) := lucas (n)$

1	nLucas (n) := lucas (n)$

La evaluación de los ejemplos es

(%i5) nLucas (5);
(%o5) 11
(%i6) nLucas (32);
(%o6) 4870847

(%i5) nLucas (5);

(%o5) 11

(%i6) nLucas (32);

(%o6) 4870847

Inverso multiplicativo modular

PorJosé A. Alonso 6-abril-201626-marzo-2022

El inverso multiplicativo modular de un entero n módulo p es el número m, entre 1 y p-1, tal que

   mn = 1 (mod p)

1	mn = 1 (mod p)

Por ejemplo, el inverso multiplicativo de 2 módulo 5 es 3, ya que 1 <= 3 <= 4 y 2×3 = 1 (mod 5).

El inverso multipicativo de n módulo p existe si y sólo si n y p son coprimos; es decir, si mcd(n,p) = 1.

Definir la función

   invMod :: Integer -> Integer -> Maybe Integer

1	invMod :: Integer -> Integer -> Maybe Integer

tal que (invMod n p) es justo el inverso multiplicativo de n módulo p, si existe y Nothing en caso contrario. Por ejemplo,

   λ> invMod 2 5
   Just 3
   λ> invMod 2 6
   Nothing
   λ> [(x,invMod x 5) | x <- [0..4]]
   [(0,Nothing),(1,Just 1),(2,Just 3),(3,Just 2),(4,Just 4)]
   λ> [(x,invMod x 6) | x <- [0..5]]
   [(0,Nothing),(1,Just 1),(2,Nothing),(3,Nothing),(4,Nothing),(5,Just 5)]
   λ> let n = 10^7 in invMod (10^n) (1+10^n) == Just (10^n)
   True

λ> invMod 2 5

Just 3

λ> invMod 2 6

Nothing

λ> [(x,invMod x 5) | x <- [0..4]]

[(0,Nothing),(1,Just 1),(2,Just 3),(3,Just 2),(4,Just 4)]

λ> [(x,invMod x 6) | x <- [0..5]]

[(0,Nothing),(1,Just 1),(2,Nothing),(3,Nothing),(4,Nothing),(5,Just 5)]

λ> let n = 10^7 in invMod (10^n) (1+10^n) == Just (10^n)

True

Soluciones

-- 1ª definición
-- =============

invMod1 :: Integer -> Integer -> Maybe Integer
invMod1 n p | gcd n p == 1 = Just (head [m | m <- [1..p-1], m*n `mod` p == 1])
            | otherwise    = Nothing

-- 2ª definición
-- =============

invMod2 :: Integer -> Integer -> Maybe Integer
invMod2 n p | m /= 1    = Nothing
            | x < 0     = Just (x + p)
            | otherwise = Just x
    where (x,_,m) = mcdExt n p

-- [Algoritmo extendido de Euclides]
-- (mcd a b) es la terna (x,y,g) tal que g es el máximo común divisor
-- de a y b y se cumple que ax + by = g. Por ejemplo,
--    mcdExt  2  5  ==  (-2,1,1)
--    mcdExt  2  6  ==  ( 1,0,2)
--    mcdExt 12 15  ==  (-1,1,3)
mcdExt :: Integer -> Integer -> (Integer,Integer,Integer)
mcdExt a 0 = (1, 0, a)
mcdExt a b = (t, s - q * t, g)
  where (q, r)    = a `quotRem` b
        (s, t, g) = mcdExt b r

-- 3ª definición
-- =============

invMod3 :: Integer -> Integer -> Maybe Integer
invMod3 n p | gcd n p == 1 = Just (aux n p)
            | otherwise    = Nothing
  where aux 1 p = 1
        aux n p = (x * p + 1) `div` n
          where x = n - aux (p `mod` n) n                    

-- Comparación de eficiencia
-- =========================

--    λ> invMod1 (10^7) (1+10^7)
--    Just 10000000
--    (5.05 secs, 2,872,350,896 bytes)
--    λ> invMod2 (10^7) (1+10^7)
--    Just 10000000
--    (0.00 secs, 0 bytes)
--    λ> invMod3 (10^7) (1+10^7)
--    Just 10000000
--    (0.00 secs, 0 bytes)
--    
--    λ> let n = 10^7 in invMod2 (10^n) (1+10^n) == Just (10^n)
--    True
--    (0.88 secs, 55,728,120 bytes)
--    λ> let n = 10^7 in invMod3 (10^n) (1+10^n) == Just (10^n)
--    True
--    (1.99 secs, 80,644,352 bytes)

-- 1ª definición

-- =============

invMod1 :: Integer -> Integer -> Maybe Integer

invMod1 n p | gcd n p == 1 = Just (head [m | m <- [1..p-1], m*n `mod` p == 1])

| otherwise = Nothing

-- 2ª definición

-- =============

invMod2 :: Integer -> Integer -> Maybe Integer

invMod2 n p | m /= 1 = Nothing

| x < 0 = Just (x + p)

| otherwise = Just x

where (x,_,m) = mcdExt n p

-- [Algoritmo extendido de Euclides]

-- (mcd a b) es la terna (x,y,g) tal que g es el máximo común divisor

-- de a y b y se cumple que ax + by = g. Por ejemplo,

-- mcdExt 2 5 == (-2,1,1)

-- mcdExt 2 6 == ( 1,0,2)

-- mcdExt 12 15 == (-1,1,3)

mcdExt :: Integer -> Integer -> (Integer,Integer,Integer)

mcdExt a 0 = (1, 0, a)

mcdExt a b = (t, s - q * t, g)

where (q, r) = a `quotRem` b

(s, t, g) = mcdExt b r

-- 3ª definición

-- =============

invMod3 :: Integer -> Integer -> Maybe Integer

invMod3 n p | gcd n p == 1 = Just (aux n p)

| otherwise = Nothing

where aux 1 p = 1

aux n p = (x * p + 1) `div` n

where x = n - aux (p `mod` n) n

-- Comparación de eficiencia

-- =========================

-- λ> invMod1 (10^7) (1+10^7)

-- Just 10000000

-- (5.05 secs, 2,872,350,896 bytes)

-- λ> invMod2 (10^7) (1+10^7)

-- Just 10000000

-- (0.00 secs, 0 bytes)

-- λ> invMod3 (10^7) (1+10^7)

-- Just 10000000

-- (0.00 secs, 0 bytes)

-- λ> let n = 10^7 in invMod2 (10^n) (1+10^n) == Just (10^n)

-- True

-- (0.88 secs, 55,728,120 bytes)

-- λ> let n = 10^7 in invMod3 (10^n) (1+10^n) == Just (10^n)

-- True

-- (1.99 secs, 80,644,352 bytes)

Solución en Maxima

invMod (n,p) := block ([r],
  r : inv_mod (n,p),
  if integerp (r)
  then Just (r)
  else Nothing)$

invMod (n,p) := block ([r],

r : inv_mod (n,p),

if integerp (r)

then Just (r)

else Nothing)$

La evaluación de los ejemplos es

(%i2) invMod (2,5);
(%o2) Just(3)
(%i3) invMod (2,6);
(%o3) Nothing
(%i4) makelist([x,invMod(x,5)],x,0,4);
(%o4) [[0,Nothing],[1,Just(1)],[2,Just(3)],[3,Just(2)],[4,Just(4)]]
(%i5) makelist([x,invMod(x,6)],x,0,5);
(%o5) [[0,Nothing],[1,Just(1)],[2,Nothing],[3,Nothing],[4,Nothing],[5,Just(5)]]
(%i6) (n : 10^7, is (invMod (10^n,1+10^n) = Just (10^n)));
(%o6) true

(%i2) invMod (2,5);

(%o2) Just(3)

(%i3) invMod (2,6);

(%o3) Nothing

(%i4) makelist([x,invMod(x,5)],x,0,4);

(%o4) [[0,Nothing],[1,Just(1)],[2,Just(3)],[3,Just(2)],[4,Just(4)]]

(%i5) makelist([x,invMod(x,6)],x,0,5);

(%o5) [[0,Nothing],[1,Just(1)],[2,Nothing],[3,Nothing],[4,Nothing],[5,Just(5)]]

(%i6) (n : 10^7, is (invMod (10^n,1+10^n) = Just (10^n)));

(%o6) true

Referencia

Calculating multiplicative inverses in modular arithmetic

Clases de equivalencia

PorJosé A. Alonso 5-abril-201625-abril-2021

Definir la función

   clasesEquivalencia :: [a] -> (a -> a -> Bool) -> [[a]]

1	clasesEquivalencia :: [a] -> (a -> a -> Bool) -> [[a]]

tal que (clasesEquivalencia xs r) es la lista de las clases de equivalencia de xs respecto de la relación de equivalencia r. Por ejemplo,

   λ> clasesEquivalencia [1..20] (\x y -> x `mod` 3 == y `mod` 3)
   [[1,4,7,10,13,16,19],[2,5,8,11,14,17,20],[3,6,9,12,15,18]]
   λ> clasesEquivalencia [1..20] (\x y -> (x - y) `mod` 5 == 0)
   [[1,6,11,16],[2,7,12,17],[3,8,13,18],[4,9,14,19],[5,10,15,20]]
   λ> clasesEquivalencia [-4..4] (\x y -> abs x == abs y)
   [[-4,4],[-3,3],[-2,2],[-1,1],[0]]

λ> clasesEquivalencia [1..20] (\x y -> x `mod` 3 == y `mod` 3)

[[1,4,7,10,13,16,19],[2,5,8,11,14,17,20],[3,6,9,12,15,18]]

λ> clasesEquivalencia [1..20] (\x y -> (x - y) `mod` 5 == 0)

[[1,6,11,16],[2,7,12,17],[3,8,13,18],[4,9,14,19],[5,10,15,20]]

λ> clasesEquivalencia [-4..4] (\x y -> abs x == abs y)

[[-4,4],[-3,3],[-2,2],[-1,1],[0]]

Soluciones

import Data.List (partition)

clasesEquivalencia :: [a] -> (a -> a -> Bool) -> [[a]]
clasesEquivalencia [] r = []
clasesEquivalencia ys@(x:xs) r =
    us : clasesEquivalencia vs r
    where (us,vs) = partition (r x) ys

import Data.List (partition)

clasesEquivalencia :: [a] -> (a -> a -> Bool) -> [[a]]

clasesEquivalencia [] r = []

clasesEquivalencia ys@(x:xs) r =

us : clasesEquivalencia vs r

where (us,vs) = partition (r x) ys

Solución en Maxima

clasesEquivalencia (xs,r) :=
  equiv_classes (setify (xs),r)$

1 2	clasesEquivalencia (xs,r) := equiv_classes (setify (xs),r)$

La evaluación de los ejemplos es

(%i6) xs : makelist (k,k,1,10)$
(%i7) clasesEquivalencia (xs, lambda ([x,y], is (remainder (x-y,3) = 0)));
(%o7) {{1, 4, 7, 10}, {2, 5, 8}, {3, 6, 9}}
(%i8) clasesEquivalencia (xs, lambda ([x,y], is (remainder (x-y,5) = 0)));
(%o8) {{1, 6}, {2, 7}, {3, 8}, {4, 9}, {5, 10}}
(%i9) clasesEquivalencia (makelist(k,k,-4,4), lambda ([x,y], is (abs (x) = abs (y))));
(%o9) {{- 4, 4}, {- 3, 3}, {- 2, 2}, {- 1, 1}, {0}}

(%i6) xs : makelist (k,k,1,10)$

(%i7) clasesEquivalencia (xs, lambda ([x,y], is (remainder (x-y,3) = 0)));

(%o7) {{1, 4, 7, 10}, {2, 5, 8}, {3, 6, 9}}

(%i8) clasesEquivalencia (xs, lambda ([x,y], is (remainder (x-y,5) = 0)));

(%o8) {{1, 6}, {2, 7}, {3, 8}, {4, 9}, {5, 10}}

(%i9) clasesEquivalencia (makelist(k,k,-4,4), lambda ([x,y], is (abs (x) = abs (y))));

(%o9) {{- 4, 4}, {- 3, 3}, {- 2, 2}, {- 1, 1}, {0}}

Factorial generalizado

PorJosé A. Alonso 4-abril-201611-abril-2016

El factorial generalizado de x respecto de y y z es el producto x(x-z)(x-2z) … (x-(y-1)z). Por ejemplo, el factorial generalizado de 7 respecto de 3 y 2 es 7x5x3 = 105 y el de 7 respecto de 2 y 3 es 7×4 = 28

Definir la función

   factGen :: Integer -> Integer -> Integer -> Integer

1	factGen :: Integer -> Integer -> Integer -> Integer

tal que (factGen x y z) es el factorial generalizado de x respecto de y y z. Por ejemplo,

   factGen 7 3 2  ==  105
   factGen 7 2 3  ==  28

1 2	factGen 7 3 2 == 105 factGen 7 2 3 == 28

Nota: Se supone que x, y y z son positivos y z < x.

Comprobar con QuickCheck que (factGen x x 1) es el factorial de x.

Soluciones

import Test.QuickCheck

-- 1ª definición (por comprensión):
factGen :: Integer -> Integer -> Integer -> Integer
factGen x y z = product [x,x-z..x-(y-1)*z]

-- 2ª definición (por recursión):
factGen2 :: Integer -> Integer -> Integer -> Integer
factGen2 x 1 _ = x
factGen2 x y z = x * factGen2 (x-z) (y-1) z

-- Propiedad de equivalencia
prop_equiv :: Positive Integer
           -> Positive Integer -> Positive Integer -> Property
prop_equiv (Positive x) (Positive y) (Positive z) =
    z < x ==> factGen x y z == factGen2 x y z

-- La comprobación es
--    λ> quickCheck prop_equiv
--    +++ OK, passed 100 tests.
      
-- La propiedad es       
prop_factGen :: Positive Integer -> Bool
prop_factGen (Positive x) =
    factGen x x 1 == product [1..x]

-- La comprobación es      
--    λ> quickCheck prop_factGen
--    +++ OK, passed 100 tests.

import Test.QuickCheck

-- 1ª definición (por comprensión):

factGen :: Integer -> Integer -> Integer -> Integer

factGen x y z = product [x,x-z..x-(y-1)*z]

-- 2ª definición (por recursión):

factGen2 :: Integer -> Integer -> Integer -> Integer

factGen2 x 1 _ = x

factGen2 x y z = x * factGen2 (x-z) (y-1) z

-- Propiedad de equivalencia

prop_equiv :: Positive Integer

-> Positive Integer -> Positive Integer -> Property

prop_equiv (Positive x) (Positive y) (Positive z) =

z < x ==> factGen x y z == factGen2 x y z

-- La comprobación es

-- λ> quickCheck prop_equiv

-- +++ OK, passed 100 tests.

-- La propiedad es

prop_factGen :: Positive Integer -> Bool

prop_factGen (Positive x) =

factGen x x 1 == product [1..x]

-- La comprobación es

-- λ> quickCheck prop_factGen

-- +++ OK, passed 100 tests.

Solución en Maxima

factGen (x,y,z) := genfact (x,y,z)$

1	factGen (x,y,z) := genfact (x,y,z)$

Representación decimal de números racionales

PorJosé A. Alonso 1-abril-20168-abril-2016

Los números decimales se representan por ternas, donde el primer elemento es la parte entera, el segundo es el anteperíodo y el tercero es el período. Por ejemplo,

    6/2  = 3                  se representa por (3,[],[])
    1/2  = 0.5                se representa por (0,[5],[])
    1/3  = 0.333333...        se representa por (0,[],[3])  
   23/14 = 1.6428571428571... se representa por (1,[6],[4,2,8,5,7,1])

6/2 = 3 se representa por (3,[],[])

1/2 = 0.5 se representa por (0,[5],[])

1/3 = 0.333333... se representa por (0,[],[3])

23/14 = 1.6428571428571... se representa por (1,[6],[4,2,8,5,7,1])

Su tipo es

   type Decimal = (Integer,[Integer],[Integer])

1	type Decimal = (Integer,[Integer],[Integer])

Los números racionales se representan por un par de enteros, donde el primer elemento es el numerador y el segundo el denominador. Por ejemplo, el número 2/3 se representa por (2,3). Su tipo es

   type Racional = (Integer,Integer)

1	type Racional = (Integer,Integer)

Definir las funciones

   decimal  :: Racional -> Decimal
   racional :: Decimal -> Racional

1 2	decimal :: Racional -> Decimal racional :: Decimal -> Racional

tales que

(decimal r) es la representación decimal del número racional r. Por ejemplo,

        decimal (1,4)    ==  (0,[2,5],[])
        decimal (1,3)    ==  (0,[],[3])
        decimal (23,14)  ==  (1,[6],[4,2,8,5,7,1])

decimal (1,4) == (0,[2,5],[])

decimal (1,3) == (0,[],[3])

decimal (23,14) == (1,[6],[4,2,8,5,7,1])

(racional d) es el número racional cuya representación decimal es d. Por ejemplo,

        racional (0,[2,5],[])           ==  (1,4)
        racional (0,[],[3])             ==  (1,3)
        racional (1,[6],[4,2,8,5,7,1])  ==  (23,14)

racional (0,[2,5],[]) == (1,4)

racional (0,[],[3]) == (1,3)

racional (1,[6],[4,2,8,5,7,1]) == (23,14)

Con la función decimal se puede calcular los períodos de los números racionales. Por ejemplo,

   ghci> let (_,_,p) = decimal (23,14) in concatMap show p
   "428571"
   ghci> let (_,_,p) = decimal (1,47) in concatMap show p
   "0212765957446808510638297872340425531914893617"
   ghci> let (_,_,p) = decimal (1,541) in length (concatMap show p)
   540

ghci> let (_,_,p) = decimal (23,14) in concatMap show p

"428571"

ghci> let (_,_,p) = decimal (1,47) in concatMap show p

"0212765957446808510638297872340425531914893617"

ghci> let (_,_,p) = decimal (1,541) in length (concatMap show p)

540

Comprobar con QuickCheck si las funciones decimal y racional son inversas.

Soluciones

-- ---------------------------------------------------------------------
-- § Librerías auxiliares                                             --
-- ---------------------------------------------------------------------

import Data.List
import Test.QuickCheck

-- ---------------------------------------------------------------------
-- § Tipos                                                            --
-- ---------------------------------------------------------------------

-- Un número racional se representa por un par de enteros, donde el
-- primer elemento es el numerador y el segundo el denominador.
type Racional = (Integer,Integer)

-- Un número decimal es una terna donde el primer elemento es la parte
-- entera, el segundo es el anteperíodo y el tercero es el período.
type Decimal = (Integer,[Integer],[Integer])

-- ---------------------------------------------------------------------
-- § De racional a decimal                                            --
-- ---------------------------------------------------------------------

-- (mayorExponente a x) es el mayor n tal que a^n divide a x. Por ejemplo,
--    mayorExponente 2 40  ==  3
--    mayorExponente 5 40  ==  1
--    mayorExponente 5 41  ==  0
mayorExponente :: Integer -> Integer -> Integer
mayorExponente a x = head [n | n <- [1..], mod x (a^n) /= 0] - 1

-- (longitudAnteperiodo n) es la longitud del anteperíodo de 1/n (y de
-- cualquier fracción irreducible de denominador n). Por ejemplo, 
--    longitudAnteperiodo 2   ==  1
--    longitudAnteperiodo 3   ==  0
--    longitudAnteperiodo 14  ==  1
longitudAnteperiodo :: Integer -> Integer
longitudAnteperiodo n = max (mayorExponente 2 n) (mayorExponente 5 n)

-- (longitudPeriodo n) es la longitud del período de 1/n (y de
-- cualquier fracción irreducible de denominador n). Por ejemplo, 
--    longitudPeriodo 2  ==  0
--    longitudPeriodo 3  ==  1
--    longitudPeriodo 7  ==  6
longitudPeriodo :: Integer -> Integer
longitudPeriodo n
    | m == 1    = 0 
    | otherwise = head [k | k <- [1..], (10^k) `mod` m == 1]
    where m = n `div` (2^(mayorExponente 2 n) * 5^(mayorExponente 5 n))

-- (expansionDec (x,y)) es la expansión decimal de x/y. Por ejemplo, 
--    take 10 (expansionDec (1,4))    ==  [0,2,5]
--    take 10 (expansionDec (1,7))    ==  [0,1,4,2,8,5,7,1,4,2]
--    take 12 (expansionDec (90,7))   ==  [12,8,5,7,1,4,2,8,5,7,1,4]
--    take 12 (expansionDec (23,14))  ==  [1,6,4,2,8,5,7,1,4,2,8,5]
expansionDec :: Racional -> [Integer]
expansionDec (x,y) 
    | r == 0    = [q] 
    | otherwise = q : expansionDec (r*10,y)
    where (q,r) = quotRem x y

-- (parteEntera (a,b)) es la parte entera de a/b. Por ejemplo,
--    parteEntera (125,3)  ==  41
parteEntera :: Racional -> Integer
parteEntera (a,b) = a `div` b

-- (antePeriodo (a,b)) es el anteperíodo de a/b; es decir, la lista de
-- cifras que se encuentran entre la parte entera y el primer período de
-- a/b. Por ejemplo,
--    antePeriodo (23,14)  ==  [6]
--    antePeriodo (1,5)    ==  [2]
antePeriodo :: Racional -> [Integer]
antePeriodo (a,b) = genericTake s (tail xs)
    where xs = expansionDec (a,b)
          s  = longitudAnteperiodo b

-- (periodo (a,b)) es el período de a/b; es decir, la lista de cifras que
-- se encuentran entre la parte entera y el primer período de a/b. Por
-- ejemplo, 
--    periodo (1,3)                   ==  [3]
--    periodo (1,5)                   ==  []
--    periodo (1,7)                   ==  [1,4,2,8,5,7]
--    periodo (23,14)                 ==  [4,2,8,5,7,1]
--    concatMap show $ periodo (1,29) == "0344827586206896551724137931"
periodo :: Racional -> [Integer]
periodo (a,b) = genericTake t (genericDrop (1+s) xs)
    where xs = expansionDec (a,b)
          s  = longitudAnteperiodo b
          t  = longitudPeriodo b

-- (reducido (a,b)) es la forma reducida del número racional a/b. Por
-- ejemplo, 
--    reducido (40,80)  ==  (1,2)
reducido :: Racional -> Racional
reducido (a,b) = (a `div` m, b `div` m) 
    where m = gcd a b

-- (decimal (x,y)) es la forma decimal de x/y; es decir, la terna
-- formada por la parte entera, la parte decimal pura y la parte decimal
-- periódica. Por ejemplo, 
--    decimal (1,4)    ==  (0,[2,5],[])
--    decimal (1,3)    ==  (0,[],[3])
--    decimal (23,14)  ==  (1,[6],[4,2,8,5,7,1])
decimal :: Racional -> Decimal
decimal (a,b) = (parteEntera r, antePeriodo r, periodo r)
    where r = reducido (a,b)

-- ---------------------------------------------------------------------
-- § De decimal a racional                                            --
-- ---------------------------------------------------------------------

-- (digitosAnumero xs) es el número correspondiente a la lista de
-- dígitos xs. Por ejemplo,
--    digitosAnumero [3,2,5]  ==  325
digitosAnumero :: [Integer] -> Integer
digitosAnumero = foldl (\a y -> 10*a+y) 0

-- 2ª definición de digitosAnumero (por comprensión)
digitosAnumero2 :: [Integer] -> Integer
digitosAnumero2 xs = sum [x*(10^i) | (x,i) <- zip xs [n-1,n-2..0]]
    where n = length xs

-- (racional (x,ys,zs)) es el número racional cuya representación
-- decimal es (x,ys,zs). Por ejemplo,
--    racional (0,[2,5],[])           ==  (1,4)
--    racional (0,[],[3])             ==  (1,3)
--    racional (1,[6],[4,2,8,5,7,1])  ==  (23,14)
racional :: Decimal -> Racional
racional (x,ys,[]) = reducido (a,b) where
    a = digitosAnumero (x:ys)
    b = 10^(length ys)
racional (x,ys,zs) = reducido (a,b) where 
    a = digitosAnumero (x:ys++zs) - digitosAnumero (x:ys)
    b = digitosAnumero (replicate (length zs) 9 ++ replicate (length ys) 0)

-- ---------------------------------------------------------------------
-- § Propiedades                                                      --
-- ---------------------------------------------------------------------

-- La 1ª propiedad es
prop_RacionalDecimal :: Racional -> Property
prop_RacionalDecimal (a,b) =
    a >= 0 && b > 0  ==>
    racional (decimal (a,b)) == reducido (a,b)

-- La comprobación es 
--    ghci> quickCheck prop_RacionalDecimal
--    +++ OK, passed 100 tests.

-- En lugar de reducido se puede usar un generador de números racionales
numeroRacional :: Gen Racional
numeroRacional = do a <- arbitrary
                    b <- arbitrary
                    return (reducido (abs a, 1+ abs b))

-- La propiedad es
prop_RacionalDecimal2 :: Property
prop_RacionalDecimal2 =
    forAll numeroRacional
           (\(a,b) -> racional (decimal (a,b)) == (a,b))

-- La comprobación es
--   ghci> quickCheck prop_RacionalDecimal2
--   +++ OK, passed 100 tests.

-- Para la 2ª propiedad se define un generador de números decimales
numeroDecimal :: Gen Decimal
numeroDecimal = do a <- arbitrary
                   b <- arbitrary
                   return (decimal (abs a, 1+ abs b))

-- La 2ª propiedad es
prop_DecimalRacional :: Property
prop_DecimalRacional =
    forAll numeroDecimal 
           (\(x,ys,zs) -> decimal (racional (x,ys,zs)) == (x,ys,zs))

-- La comprobación es
--   ghci> quickCheck prop_DecimalRacional
--   +++ OK, passed 100 tests.

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

-- ---------------------------------------------------------------------

-- § Librerías auxiliares --

-- ---------------------------------------------------------------------

import Data.List

import Test.QuickCheck

-- ---------------------------------------------------------------------

-- § Tipos --

-- ---------------------------------------------------------------------

-- Un número racional se representa por un par de enteros, donde el

-- primer elemento es el numerador y el segundo el denominador.

type Racional = (Integer,Integer)

-- Un número decimal es una terna donde el primer elemento es la parte

-- entera, el segundo es el anteperíodo y el tercero es el período.

type Decimal = (Integer,[Integer],[Integer])

-- ---------------------------------------------------------------------

-- § De racional a decimal --

-- ---------------------------------------------------------------------

-- (mayorExponente a x) es el mayor n tal que a^n divide a x. Por ejemplo,

-- mayorExponente 2 40 == 3

-- mayorExponente 5 40 == 1

-- mayorExponente 5 41 == 0

mayorExponente :: Integer -> Integer -> Integer

mayorExponente a x = head [n | n <- [1..], mod x (a^n) /= 0] - 1

-- (longitudAnteperiodo n) es la longitud del anteperíodo de 1/n (y de

-- cualquier fracción irreducible de denominador n). Por ejemplo,

-- longitudAnteperiodo 2 == 1

-- longitudAnteperiodo 3 == 0

-- longitudAnteperiodo 14 == 1

longitudAnteperiodo :: Integer -> Integer

longitudAnteperiodo n = max (mayorExponente 2 n) (mayorExponente 5 n)

-- (longitudPeriodo n) es la longitud del período de 1/n (y de

-- cualquier fracción irreducible de denominador n). Por ejemplo,

-- longitudPeriodo 2 == 0

-- longitudPeriodo 3 == 1

-- longitudPeriodo 7 == 6

longitudPeriodo :: Integer -> Integer

longitudPeriodo n

| m == 1 = 0

| otherwise = head [k | k <- [1..], (10^k) `mod` m == 1]

where m = n `div` (2^(mayorExponente 2 n) * 5^(mayorExponente 5 n))

-- (expansionDec (x,y)) es la expansión decimal de x/y. Por ejemplo,

-- take 10 (expansionDec (1,4)) == [0,2,5]

-- take 10 (expansionDec (1,7)) == [0,1,4,2,8,5,7,1,4,2]

-- take 12 (expansionDec (90,7)) == [12,8,5,7,1,4,2,8,5,7,1,4]

-- take 12 (expansionDec (23,14)) == [1,6,4,2,8,5,7,1,4,2,8,5]

expansionDec :: Racional -> [Integer]

expansionDec (x,y)

| r == 0 = [q]

| otherwise = q : expansionDec (r*10,y)

where (q,r) = quotRem x y

-- (parteEntera (a,b)) es la parte entera de a/b. Por ejemplo,

-- parteEntera (125,3) == 41

parteEntera :: Racional -> Integer

parteEntera (a,b) = a `div` b

-- (antePeriodo (a,b)) es el anteperíodo de a/b; es decir, la lista de

-- cifras que se encuentran entre la parte entera y el primer período de

-- a/b. Por ejemplo,

-- antePeriodo (23,14) == [6]

-- antePeriodo (1,5) == [2]

antePeriodo :: Racional -> [Integer]

antePeriodo (a,b) = genericTake s (tail xs)

where xs = expansionDec (a,b)

s = longitudAnteperiodo b

-- (periodo (a,b)) es el período de a/b; es decir, la lista de cifras que

-- se encuentran entre la parte entera y el primer período de a/b. Por

-- ejemplo,

-- periodo (1,3) == [3]

-- periodo (1,5) == []

-- periodo (1,7) == [1,4,2,8,5,7]

-- periodo (23,14) == [4,2,8,5,7,1]

-- concatMap show $ periodo (1,29) == "0344827586206896551724137931"

periodo :: Racional -> [Integer]

periodo (a,b) = genericTake t (genericDrop (1+s) xs)

where xs = expansionDec (a,b)

s = longitudAnteperiodo b

t = longitudPeriodo b

-- (reducido (a,b)) es la forma reducida del número racional a/b. Por

-- ejemplo,

-- reducido (40,80) == (1,2)

reducido :: Racional -> Racional

reducido (a,b) = (a `div` m, b `div` m)

where m = gcd a b

-- (decimal (x,y)) es la forma decimal de x/y; es decir, la terna

-- formada por la parte entera, la parte decimal pura y la parte decimal

-- periódica. Por ejemplo,

-- decimal (1,4) == (0,[2,5],[])

-- decimal (1,3) == (0,[],[3])

-- decimal (23,14) == (1,[6],[4,2,8,5,7,1])

decimal :: Racional -> Decimal

decimal (a,b) = (parteEntera r, antePeriodo r, periodo r)

where r = reducido (a,b)

-- ---------------------------------------------------------------------

-- § De decimal a racional --

-- ---------------------------------------------------------------------

-- (digitosAnumero xs) es el número correspondiente a la lista de

-- dígitos xs. Por ejemplo,

-- digitosAnumero [3,2,5] == 325

digitosAnumero :: [Integer] -> Integer

digitosAnumero = foldl (\a y -> 10*a+y) 0

-- 2ª definición de digitosAnumero (por comprensión)

digitosAnumero2 :: [Integer] -> Integer

digitosAnumero2 xs = sum [x*(10^i) | (x,i) <- zip xs [n-1,n-2..0]]

where n = length xs

-- (racional (x,ys,zs)) es el número racional cuya representación

-- decimal es (x,ys,zs). Por ejemplo,

-- racional (0,[2,5],[]) == (1,4)

-- racional (0,[],[3]) == (1,3)

-- racional (1,[6],[4,2,8,5,7,1]) == (23,14)

racional :: Decimal -> Racional

racional (x,ys,[]) = reducido (a,b) where

a = digitosAnumero (x:ys)

b = 10^(length ys)

racional (x,ys,zs) = reducido (a,b) where

a = digitosAnumero (x:ys++zs) - digitosAnumero (x:ys)

b = digitosAnumero (replicate (length zs) 9 ++ replicate (length ys) 0)

-- ---------------------------------------------------------------------

-- § Propiedades --

-- ---------------------------------------------------------------------

-- La 1ª propiedad es

prop_RacionalDecimal :: Racional -> Property

prop_RacionalDecimal (a,b) =

a >= 0 && b > 0 ==>

racional (decimal (a,b)) == reducido (a,b)

-- La comprobación es

-- ghci> quickCheck prop_RacionalDecimal

-- +++ OK, passed 100 tests.

-- En lugar de reducido se puede usar un generador de números racionales

numeroRacional :: Gen Racional

numeroRacional = do a <- arbitrary

b <- arbitrary

return (reducido (abs a, 1+ abs b))

-- La propiedad es

prop_RacionalDecimal2 :: Property

prop_RacionalDecimal2 =

forAll numeroRacional

(\(a,b) -> racional (decimal (a,b)) == (a,b))

-- La comprobación es

-- ghci> quickCheck prop_RacionalDecimal2

-- +++ OK, passed 100 tests.

-- Para la 2ª propiedad se define un generador de números decimales

numeroDecimal :: Gen Decimal

numeroDecimal = do a <- arbitrary

b <- arbitrary

return (decimal (abs a, 1+ abs b))

-- La 2ª propiedad es

prop_DecimalRacional :: Property

prop_DecimalRacional =

forAll numeroDecimal

(\(x,ys,zs) -> decimal (racional (x,ys,zs)) == (x,ys,zs))

-- La comprobación es

-- ghci> quickCheck prop_DecimalRacional

-- +++ OK, passed 100 tests.

Meses: abril 2016

Operaciones con series de potencias

Soluciones

Productos simultáneos de dos y tres números consecutivos

Soluciones

Conflictos de horarios

Soluciones

Mínima diferencia entre elementos de una lista

Soluciones

Números como sumas de primos consecutivos

Soluciones

Caminos reducidos

Soluciones

Primos permutables

Soluciones

Referencias

Raíz entera

Soluciones

Soluciones en Maxima

Números de Lucas

Soluciones

Soluciones en Maxima

Inverso multiplicativo modular

Soluciones

Solución en Maxima

Referencia

Clases de equivalencia

Soluciones

Solución en Maxima

Factorial generalizado

Soluciones

Solución en Maxima

Representación decimal de números racionales

Soluciones

Entradas recientes

Buscar

Correo electrónico