fromIntegral

Caminos en una retícula

PorJosé A. Alonso 26-mayo-20162-junio-2016

El problema de los caminos en una retícula consiste en, dada una retícula rectangular con m filas y n columnas, determinar todos los caminos para ir desde el vértice inferior izquierdo hasta el vértice superior derecho donde los movimientos permitidos son mover hacia el siguiente vértice a la derecha o arriba.

Por ejemplo, en la siguiente retícula un posible camino es el indicado en rojo.

Para representar los caminos se definen los siguientes tipos de datos:

    data Direccion = D | A deriving (Show, Eq)
    type Camino = [Direccion]

1 2	data Direccion = D \| A deriving (Show, Eq) type Camino = [Direccion]

Por tanto, el camino de la figura anterior se representa por la lista [D,D,A,D,A].

Definir las funciones

   caminos  :: Int -> Int -> [Camino]
   nCaminos :: Int -> Int -> Integer

1 2	caminos :: Int -> Int -> [Camino] nCaminos :: Int -> Int -> Integer

tales que

(caminos m n) es la lista de los caminos en una retícula rectangular con m filas y n columnas. Por ejemplo,

     λ> caminos 2 2
     [[A,A,D,D],[A,D,A,D],[A,D,D,A],[D,A,A,D],[D,A,D,A],[D,D,A,A]]
     λ> caminos 1 3
     [[A,A,A,D],[A,A,D,A],[A,D,A,A],[D,A,A,A]]
     λ> caminos 2 3
     [[A,A,A,D,D],[A,A,D,A,D],[A,A,D,D,A],[A,D,A,A,D],[A,D,A,D,A],[A,D,D,A,A],
      [D,A,A,A,D],[D,A,A,D,A],[D,A,D,A,A],[D,D,A,A,A]]

λ> caminos 2 2

[[A,A,D,D],[A,D,A,D],[A,D,D,A],[D,A,A,D],[D,A,D,A],[D,D,A,A]]

λ> caminos 1 3

[[A,A,A,D],[A,A,D,A],[A,D,A,A],[D,A,A,A]]

λ> caminos 2 3

[[A,A,A,D,D],[A,A,D,A,D],[A,A,D,D,A],[A,D,A,A,D],[A,D,A,D,A],[A,D,D,A,A],

[D,A,A,A,D],[D,A,A,D,A],[D,A,D,A,A],[D,D,A,A,A]]

(nCaminos m n) es el número de los caminos en una retícula rectangular con m filas y n columnas. Por ejemplo,

     nCaminos 2 2  ==  6
     nCaminos 1 3  ==  4
     nCaminos 2 3  ==  10
     length (show (nCaminos 20000 30000))  ==  14612
     length (show (nCaminos 30000 20000))  ==  14612

nCaminos 2 2 == 6

nCaminos 1 3 == 4

nCaminos 2 3 == 10

length (show (nCaminos 20000 30000)) == 14612

length (show (nCaminos 30000 20000)) == 14612

Soluciones

import Data.List (genericLength)

data Direccion = D | A deriving (Show, Eq)

type Camino = [Direccion]

-- Definición de caminos
-- =====================

caminos :: Int -> Int -> [Camino]
caminos 0 n = [replicate n A]
caminos m 0 = [replicate m D]
caminos m n = [A:xs | xs <- caminos m (n-1)] ++
              [D:xs | xs <- caminos (m-1) n]

-- 1ª definición de nCaminos
-- =========================

nCaminos1 :: Int -> Int -> Integer
nCaminos1 m n = genericLength (caminos m n)

-- 2ª definición de nCaminos
-- =========================

nCaminos2 :: Int -> Int -> Integer
nCaminos2 0 n = 1
nCaminos2 m 0 = 1
nCaminos2 m n = nCaminos2 m (n-1) + nCaminos2 (m-1) n

-- 3ª definición de nCaminos
-- =========================

-- Los caminos desde (0,0) a (m,n) son las permutaciones con repetición
-- de m veces la D y n veces la A. Por tanto, su número es
--    (m+n)! / m!*n!

nCaminos3 :: Int -> Int -> Integer
nCaminos3 m n = 
    fact (m+n) `div` (fact m * fact n)

fact :: Int -> Integer
fact n = product [1..fromIntegral n]

-- 4ª solución de nCaminos
-- =======================

nCaminos4 :: Int -> Int -> Integer
nCaminos4 m n = 
    product [a+1..a+b] `div` product [2..b]
    where m' = fromIntegral m
          n' = fromIntegral n
          a  = max m' n'
          b  = min m' n'

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

--    λ> nCaminos1 10 10
--    184756
--    (3.50 secs, 656,909,648 bytes)
--    λ> nCaminos2 10 10
--    184756
--    (0.50 secs, 47,330,272 bytes)
--    λ> nCaminos3 10 10
--    184756
--    (0.01 secs, 0 bytes)
--    λ> nCaminos4 10 10
--    184756
--    (0.01 secs, 0 bytes)
-- 
--    λ> nCaminos2 10 15
--    3268760
--    (8.83 secs, 1,142,623,080 bytes)
--    λ> nCaminos3 10 15
--    3268760
--    (0.01 secs, 0 bytes)
--    λ> nCaminos4 10 15
--    3268760
--    (0.00 secs, 0 bytes)
--
--    λ> let n = 2*10^4 
--    (0.01 secs, 0 bytes)
--    λ> length (show (nCaminos3 n n))
--    12039
--    (8.41 secs, 2,369,767,480 bytes)
--    λ> length (show (nCaminos4 n n))
--    12039
--    (2.98 secs, 833,386,648 bytes)

import Data.List (genericLength)

data Direccion = D | A deriving (Show, Eq)

type Camino = [Direccion]

-- Definición de caminos

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

caminos :: Int -> Int -> [Camino]

caminos 0 n = [replicate n A]

caminos m 0 = [replicate m D]

caminos m n = [A:xs | xs <- caminos m (n-1)] ++

[D:xs | xs <- caminos (m-1) n]

-- 1ª definición de nCaminos

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

nCaminos1 :: Int -> Int -> Integer

nCaminos1 m n = genericLength (caminos m n)

-- 2ª definición de nCaminos

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

nCaminos2 :: Int -> Int -> Integer

nCaminos2 0 n = 1

nCaminos2 m 0 = 1

nCaminos2 m n = nCaminos2 m (n-1) + nCaminos2 (m-1) n

-- 3ª definición de nCaminos

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

-- Los caminos desde (0,0) a (m,n) son las permutaciones con repetición

-- de m veces la D y n veces la A. Por tanto, su número es

-- (m+n)! / m!*n!

nCaminos3 :: Int -> Int -> Integer

nCaminos3 m n =

fact (m+n) `div` (fact m * fact n)

fact :: Int -> Integer

fact n = product [1..fromIntegral n]

-- 4ª solución de nCaminos

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

nCaminos4 :: Int -> Int -> Integer

nCaminos4 m n =

product [a+1..a+b] `div` product [2..b]

where m' = fromIntegral m

n' = fromIntegral n

a = max m' n'

b = min m' n'

-- Comparación de eficiencia

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

-- λ> nCaminos1 10 10

-- 184756

-- (3.50 secs, 656,909,648 bytes)

-- λ> nCaminos2 10 10

-- 184756

-- (0.50 secs, 47,330,272 bytes)

-- λ> nCaminos3 10 10

-- 184756

-- (0.01 secs, 0 bytes)

-- λ> nCaminos4 10 10

-- 184756

-- (0.01 secs, 0 bytes)

-- λ> nCaminos2 10 15

-- 3268760

-- (8.83 secs, 1,142,623,080 bytes)

-- λ> nCaminos3 10 15

-- 3268760

-- (0.01 secs, 0 bytes)

-- λ> nCaminos4 10 15

-- 3268760

-- (0.00 secs, 0 bytes)

-- λ> let n = 2*10^4

-- (0.01 secs, 0 bytes)

-- λ> length (show (nCaminos3 n n))

-- 12039

-- (8.41 secs, 2,369,767,480 bytes)

-- λ> length (show (nCaminos4 n n))

-- 12039

-- (2.98 secs, 833,386,648 bytes)

Densidad de números no monótonos

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

Un número entero positivo se dice que es

creciente si cada uno de sus dígitos es menor o igual que el que está a su derecha; por ejemplo, 134479.
decreciente si cada uno de sus dígitos es menor o igual que el que está a su derecha; por ejemplo, 664210.
no monótono si no es creciente ni decreciente; por ejemplo, 155369.

Para cada entero positivo n, la densidad números no monótonos hasta n es el cociente entre la cantidad de n números no monótonos entre menores o iguales que n y el número n. Por ejemplo, hasta 150 hay 19 números no monótonos (101, 102, 103, 104, 105, 106, 107, 108, 109, 120, 121, 130, 131, 132, 140, 141, 142, 143 y 150); por tanto, la densidad hasta 150 es 19/150 = 0.12666667

Definir las siguientes funciones

   densidad              :: Integer -> Float
   menorConDensidadMayor :: Float -> Integer

1 2	densidad :: Integer -> Float menorConDensidadMayor :: Float -> Integer

tales que

(densidad n) es la densidad de números no monótonos hasta n. Por ejemplo,

     densidad 100  ==  0.0
     densidad 200  ==  0.265
     densidad 400  ==  0.4375
     densidad 800  ==  0.53375

densidad 100 == 0.0

densidad 200 == 0.265

densidad 400 == 0.4375

densidad 800 == 0.53375

(menorConDensidadMayor x) es el menor número n tal que la densidad de números no monótonos hasta n es mayor o igual que x. Por ejemplo,

     menorConDensidadMayor 0.5   ==  538
     densidad 537                ==  0.49906892
     densidad 538                ==  0.5
     menorConDensidadMayor 0.9   ==  21780
     densidad 21779              ==  0.8999954
     densidad 21780              ==  0.9
     menorConDensidadMayor 0.99  ==  1586996

menorConDensidadMayor 0.5 == 538

densidad 537 == 0.49906892

densidad 538 == 0.5

menorConDensidadMayor 0.9 == 21780

densidad 21779 == 0.8999954

densidad 21780 == 0.9

menorConDensidadMayor 0.99 == 1586996

Soluciones

import Data.List (genericLength, sort)

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

densidad :: Integer -> Float
densidad n = 
    genericLength [x | x <- [1..n], esNoMonotono x] / fromIntegral n

esNoMonotono :: Integer -> Bool
esNoMonotono x = not (esCreciente x) && not (esDecreciente x)

esCreciente :: Integer -> Bool
esCreciente x = show x == sort (show x)

esDecreciente :: Integer -> Bool
esDecreciente x = reverse (show x) == sort (show x)

menorConDensidadMayor :: Float -> Integer 
menorConDensidadMayor x  = 
    buscaProporcion x (filter esNoMonotono [1..])
    where buscaProporcion x = 
              snd . head . filter condicion . zip [1..]
              where condicion (a,b) = a >= x * fromIntegral b

import Data.List (genericLength, sort)

-- 1ª definición

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

densidad :: Integer -> Float

densidad n =

genericLength [x | x <- [1..n], esNoMonotono x] / fromIntegral n

esNoMonotono :: Integer -> Bool

esNoMonotono x = not (esCreciente x) && not (esDecreciente x)

esCreciente :: Integer -> Bool

esCreciente x = show x == sort (show x)

esDecreciente :: Integer -> Bool

esDecreciente x = reverse (show x) == sort (show x)

menorConDensidadMayor :: Float -> Integer

menorConDensidadMayor x =

buscaProporcion x (filter esNoMonotono [1..])

where buscaProporcion x =

snd . head . filter condicion . zip [1..]

where condicion (a,b) = a >= x * fromIntegral b

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..]

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)$

Primo suma de dos cuadrados

PorJosé A. Alonso 14-marzo-20161-abril-2016

Definir la sucesión

   primosSumaDe2Cuadrados :: [Integer]

1	primosSumaDe2Cuadrados :: [Integer]

cuyos elementos son los números primos que se pueden escribir como sumas de dos cuadrados. Por ejemplo,

   λ> take 20 primosSumaDe2Cuadrados
   [2,5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181]
   λ> primosSumaDe2Cuadrados !! (2*10^5)
   5803241

λ> take 20 primosSumaDe2Cuadrados

[2,5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181]

λ> primosSumaDe2Cuadrados !! (2*10^5)

5803241

En el ejemplo anterior,

13 está en la sucesión porque es primo y 13 = 2²+3².
11 no está en la sucesión porque no se puede escribir como suma de dos cuadrados (en efecto, 11-1=10, 11-2²=7 y 11-3²=2 no son cuadrados).
20 no está en la sucesión porque, aunque es suma de dos cuadrados (20=4²+2²), no es primo.

Soluciones

import Data.Numbers.Primes (primes)

-- 1ª solución
-- ===========

primosSumaDe2Cuadrados1 :: [Integer]
primosSumaDe2Cuadrados1 =
    [n | n <- primes,
         esSumaDe2Cuadrados n]

esSumaDe2Cuadrados :: Integer -> Bool
esSumaDe2Cuadrados n = not (null (sumas2cuadrados n))

sumas2cuadrados :: Integer -> [(Integer,Integer)]
sumas2cuadrados n = 
    [(x,y) | x <- [1..ceiling (sqrt (fromIntegral n))],
             let z = n - x^2,
             esCuadrado z, 
             let y = raiz z,
             x <= y]
             
-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 26  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x = x == y * y
    where y = raiz x

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz 25  ==  5
--    raiz 24  ==  4
--    raiz 26  ==  5
raiz :: Integer -> Integer
raiz x = floor (sqrt (fromIntegral x))

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

primosSumaDe2Cuadrados2 :: [Integer]
primosSumaDe2Cuadrados2 =
    2 : [n | n <- primes,
             n `mod` 4 == 1]

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

--    λ> primosSumaDe2Cuadrados1 !! (2*10^3)
--    38153
--    (3.03 secs, 568,239,744 bytes)
--    λ> primosSumaDe2Cuadrados2 !! (2*10^3)
--    38153
--    (0.05 secs, 20,017,912 bytes)

import Data.Numbers.Primes (primes)

-- 1ª solución

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

primosSumaDe2Cuadrados1 :: [Integer]

primosSumaDe2Cuadrados1 =

[n | n <- primes,

esSumaDe2Cuadrados n]

esSumaDe2Cuadrados :: Integer -> Bool

esSumaDe2Cuadrados n = not (null (sumas2cuadrados n))

sumas2cuadrados :: Integer -> [(Integer,Integer)]

sumas2cuadrados n =

[(x,y) | x <- [1..ceiling (sqrt (fromIntegral n))],

let z = n - x^2,

esCuadrado z,

let y = raiz z,

x <= y]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por

-- ejemplo,

-- esCuadrado 25 == True

-- esCuadrado 26 == False

esCuadrado :: Integer -> Bool

esCuadrado x = x == y * y

where y = raiz x

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,

-- raiz 25 == 5

-- raiz 24 == 4

-- raiz 26 == 5

raiz :: Integer -> Integer

raiz x = floor (sqrt (fromIntegral x))

-- 2ª solución

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

primosSumaDe2Cuadrados2 :: [Integer]

primosSumaDe2Cuadrados2 =

2 : [n | n <- primes,

n `mod` 4 == 1]

-- Comparación de eficiencia

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

-- λ> primosSumaDe2Cuadrados1 !! (2*10^3)

-- 38153

-- (3.03 secs, 568,239,744 bytes)

-- λ> primosSumaDe2Cuadrados2 !! (2*10^3)

-- 38153

-- (0.05 secs, 20,017,912 bytes)

Referencias

N.J.A. Sloane, Sucesión A002313 de la OEIS.
M. Bates, Primes of the form x² + ny²–
G. Xiao, Two squares.

Sumas de potencias de 3 primos

PorJosé A. Alonso 24-febrero-201625-marzo-2022

Los primeros números de la forma p²+q³+r⁴, con p, q y r primos son

   28 = 2² + 2³ + 2⁴
   33 = 3² + 2³ + 2⁴
   47 = 2² + 3³ + 2⁴
   49 = 5² + 2³ + 2⁴

28 = 2² + 2³ + 2⁴

33 = 3² + 2³ + 2⁴

47 = 2² + 3³ + 2⁴

49 = 5² + 2³ + 2⁴

Definir la sucesión

   sumas3potencias :: [Integer]

1	sumas3potencias :: [Integer]

cuyos elementos son los números que se pueden escribir de la forma p²+q³+r⁴, con p, q y r primos. Por ejemplo,

   λ> take 15 sumas3potencias
   [28,33,47,49,52,68,73,92,93,98,112,114,117,133,138]
   λ> sumas3potencias !! 234567
   8953761

λ> take 15 sumas3potencias

[28,33,47,49,52,68,73,92,93,98,112,114,117,133,138]

λ> sumas3potencias !! 234567

8953761

Soluciones

import Data.Numbers.Primes (primes)

-- 1ª solución
-- ===========

sumas3potencias1 :: [Integer]
sumas3potencias1 = filter esTripleSuma [1..]
 
-- (esTripleSuma n) se verifica si n se puede escribir de la forma
-- p²+q³+r⁴, con p, q y r primos. Por ejemplo, 
--    esTripleSuma 33  ==  True
--    esTripleSuma 45  ==  False
esTripleSuma :: Integer -> Bool
esTripleSuma = not . null . triplesSumas  

-- (triplesSumas n) es la lista de las ternas de primos (p,q,r) tales
-- que p²+q³+r⁴ = n. Por ejemplo,
--    triplesSumas 33   ==  [(3,2,2)]
--    triplesSumas 145  ==  [(11,2,2),(2,5,2)]
--    triplesSumas 45   ==  []
triplesSumas :: Integer -> [(Integer,Integer,Integer)]
triplesSumas n =  
    [(p,q,r) 
    | r <- xs, 
      q <- xs, 
      let p2 = n - q^3 - r^4,
      let p = (floor . sqrt . fromIntegral) p2,
      p*p == p2,
      isPrime p]
    where xs = takeWhile (< (ceiling $ sqrt $ fromIntegral n)) primes

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

sumas3potencias2 :: [Integer]
sumas3potencias2 = sumaOrdenadaInf as (sumaOrdenadaInf bs cs)

sumaOrdenadaInf:: [Integer] -> [Integer] -> [Integer]
sumaOrdenadaInf xs ys = mezclaTodas [map (+x) ys | x <- xs]

mezclaTodas :: Ord a => [[a]] -> [a]
mezclaTodas = foldr1 xmezcla
    where xmezcla (x:xs) ys = x : mezcla xs ys

mezcla :: Ord a => [a] -> [a] -> [a]
mezcla (x:xs) (y:ys) | x < y  = x : mezcla xs (y:ys)
                     | x == y = x : mezcla xs ys
                     | x > y  = y : mezcla (x:xs) ys
 
as = map (^2) primes
bs = map (^3) primes
cs = map (^4) primes

import Data.Numbers.Primes (primes)

-- 1ª solución

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

sumas3potencias1 :: [Integer]

sumas3potencias1 = filter esTripleSuma [1..]

-- (esTripleSuma n) se verifica si n se puede escribir de la forma

-- p²+q³+r⁴, con p, q y r primos. Por ejemplo,

-- esTripleSuma 33 == True

-- esTripleSuma 45 == False

esTripleSuma :: Integer -> Bool

esTripleSuma = not . null . triplesSumas

-- (triplesSumas n) es la lista de las ternas de primos (p,q,r) tales

-- que p²+q³+r⁴ = n. Por ejemplo,

-- triplesSumas 33 == [(3,2,2)]

-- triplesSumas 145 == [(11,2,2),(2,5,2)]

-- triplesSumas 45 == []

triplesSumas :: Integer -> [(Integer,Integer,Integer)]

triplesSumas n =

[(p,q,r)

| r <- xs,

q <- xs,

let p2 = n - q^3 - r^4,

let p = (floor . sqrt . fromIntegral) p2,

p*p == p2,

isPrime p]

where xs = takeWhile (< (ceiling $ sqrt $ fromIntegral n)) primes

-- 2ª solución

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

sumas3potencias2 :: [Integer]

sumas3potencias2 = sumaOrdenadaInf as (sumaOrdenadaInf bs cs)

sumaOrdenadaInf:: [Integer] -> [Integer] -> [Integer]

sumaOrdenadaInf xs ys = mezclaTodas [map (+x) ys | x <- xs]

mezclaTodas :: Ord a => [[a]] -> [a]

mezclaTodas = foldr1 xmezcla

where xmezcla (x:xs) ys = x : mezcla xs ys

mezcla :: Ord a => [a] -> [a] -> [a]

mezcla (x:xs) (y:ys) | x < y = x : mezcla xs (y:ys)

| x == y = x : mezcla xs ys

| x > y = y : mezcla (x:xs) ys

as = map (^2) primes

bs = map (^3) primes

cs = map (^4) primes

Números primos de Hilbert

PorJosé A. Alonso 4-febrero-20161-mayo-2016

Un número de Hilbert es un entero positivo de la forma 4n+1. Los primeros números de Hilbert son 1, 5, 9, 13, 17, 21, 25, 29, 33, 37, 41, 45, 49, 53, 57, 61, 65, 69, …

Un primo de Hilbert es un número de Hilbert n que no es divisible por ningún número de Hilbert menor que n (salvo el 1). Los primeros primos de Hilbert son 5, 9, 13, 17, 21, 29, 33, 37, 41, 49, 53, 57, 61, 69, 73, 77, 89, 93, 97, 101, 109, 113, 121, 129, 133, 137, …

Definir la sucesión

   primosH :: [Integer]

1	primosH :: [Integer]

tal que sus elementos son los primos de Hilbert. Por ejemplo,

   take 15 primosH  == [5,9,13,17,21,29,33,37,41,49,53,57,61,69,73]
   primosH !! 20000 == 203221

1 2	take 15 primosH == [5,9,13,17,21,29,33,37,41,49,53,57,61,69,73] primosH !! 20000 == 203221

Soluciones

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

-- numerosH es la sucesión de los números de Hilbert. Por ejemplo,
--    take 15 numerosH  ==  [1,5,9,13,17,21,25,29,33,37,41,45,49,53,57]
numerosH :: [Integer]
numerosH = [1,5..]

-- (divisoresH n) es la lista de los números de Hilbert que dividen a
-- n. Por ejemplo,
--   divisoresH 117  ==  [1,9,13,117]
--   divisoresH  21  ==  [1,21]
divisoresH :: Integer -> [Integer]
divisoresH n = [x | x <- takeWhile (<=n) numerosH,
                    n `mod` x == 0]

primosH1 :: [Integer]
primosH1 = [n | n <- tail numerosH,
                divisoresH n == [1,n]]

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

primosH2 :: [Integer]
primosH2 = filter esPrimoH (tail numerosH) 
    where esPrimoH n = all noDivideAn [5,9..m]
              where noDivideAn x = n `mod` x /= 0
                    m            = ceiling (sqrt (fromIntegral n))

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

--    λ> primosH1 !! 2000
--    16957
--    (6.93 secs, 750,291,352 bytes)
--    λ> primosH2 !! 2000
--    16957
--    (0.13 secs, 18,066,288 bytes)

-- 1ª definición

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

-- numerosH es la sucesión de los números de Hilbert. Por ejemplo,

-- take 15 numerosH == [1,5,9,13,17,21,25,29,33,37,41,45,49,53,57]

numerosH :: [Integer]

numerosH = [1,5..]

-- (divisoresH n) es la lista de los números de Hilbert que dividen a

-- n. Por ejemplo,

-- divisoresH 117 == [1,9,13,117]

-- divisoresH 21 == [1,21]

divisoresH :: Integer -> [Integer]

divisoresH n = [x | x <- takeWhile (<=n) numerosH,

n `mod` x == 0]

primosH1 :: [Integer]

primosH1 = [n | n <- tail numerosH,

divisoresH n == [1,n]]

-- 2ª definición

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

primosH2 :: [Integer]

primosH2 = filter esPrimoH (tail numerosH)

where esPrimoH n = all noDivideAn [5,9..m]

where noDivideAn x = n `mod` x /= 0

m = ceiling (sqrt (fromIntegral n))

-- Comparación de eficiencia

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

-- λ> primosH1 !! 2000

-- 16957

-- (6.93 secs, 750,291,352 bytes)

-- λ> primosH2 !! 2000

-- 16957

-- (0.13 secs, 18,066,288 bytes)

Suma con redondeos

PorJosé A. Alonso 1-enero-20167-enero-2016

Definir las funciones

   sumaRedondeos       :: Integer -> [Integer]
   limiteSumaRedondeos :: Integer -> Integer

1 2	sumaRedondeos :: Integer -> [Integer] limiteSumaRedondeos :: Integer -> Integer

tales que

(sumaRedondeos n) es la sucesión cuyo k-ésimo término es

     redondeo (n/2) + redondeo (n/4) + ... + redondeo (n/2^k)

1	redondeo (n/2) + redondeo (n/4) + ... + redondeo (n/2^k)

Por ejemplo,

     take 5 (sumaRedondeos 1000)  ==  [500,750,875,937,968]

1	take 5 (sumaRedondeos 1000) == [500,750,875,937,968]

(limiteSumaRedondeos n) es la suma de la serie

     redondeo (n/2) + redondeo (n/4) + redondeo (n/8) + ...

1	redondeo (n/2) + redondeo (n/4) + redondeo (n/8) + ...

Por ejemplo,

     limiteSumaRedondeos 2000                    ==  1999
     limiteSumaRedondeos 2016                    ==  2016
     limiteSumaRedondeos (10^308) `mod` (10^10)  ==  123839487

limiteSumaRedondeos 2000 == 1999

limiteSumaRedondeos 2016 == 2016

limiteSumaRedondeos (10^308) `mod` (10^10) == 123839487

Soluciones

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

sumaRedondeos1 :: Integer -> [Integer]
sumaRedondeos1 n =
    [sum [round (n'/(fromIntegral (2^k))) | k <- [1..m]] | m <- [1..]]
    where n' = fromIntegral n

limiteSumaRedondeos1 :: Integer -> Integer
limiteSumaRedondeos1 = limite . sumaRedondeos1

limite :: [Integer] -> Integer
limite xs = head [x | (x,y) <- zip xs (tail xs), x == y]

-- 2ª definición
sumaRedondeos2 :: Integer -> [Integer]
sumaRedondeos2 n =
    scanl1 (+) [round (n'/(fromIntegral (2^k))) | k <- [1..]]
    where n' = fromIntegral n

limiteSumaRedondeos2 :: Integer -> Integer
limiteSumaRedondeos2 = limite . sumaRedondeos2

-- 3ª definición
sumaRedondeos3 :: Integer -> [Integer]
sumaRedondeos3 n = map fst (iterate f (round (n'/2),4))
    where f (s,d) = (s + round (n'/(fromIntegral d)), 2*d)
          n'      = fromIntegral n

limiteSumaRedondeos3 :: Integer -> Integer
limiteSumaRedondeos3 = limite . sumaRedondeos3

-- 4ª definición
sumaRedondeos4 :: Integer -> [Integer]
sumaRedondeos4 n = xs ++ repeat x
    where n' = fromIntegral n
          m  = round (logBase 2 n')
          xs = scanl1 (+) [round (n'/(fromIntegral (2^k))) | k <- [1..m]]
          x  = last xs

limiteSumaRedondeos4 :: Integer -> Integer
limiteSumaRedondeos4 = limite . sumaRedondeos4

-- Comparación de eficiencia
--    λ> (sumaRedondeos1 4) !! 20000
--    3
--    (0.92 secs, 197,782,232 bytes)
--    λ> (sumaRedondeos1 4) !! 30000
--    3
--    (2.20 secs, 351,084,816 bytes)
--    λ> (sumaRedondeos3 4) !! 20000
--    3
--    (0.30 secs, 53,472,392 bytes)
--    λ> (sumaRedondeos4 4) !! 20000
--    3
--    (0.01 secs, 0 bytes)

-- En lo que sigue, usaremos la 3ª definición
sumaRedondeos :: Integer -> [Integer]
sumaRedondeos = sumaRedondeos3

-- 5ª definición
-- =============

limiteSumaRedondeos5 :: Integer -> Integer
limiteSumaRedondeos5 n = 
    sum [round (n'/(fromIntegral (2^k))) | k <- [1..m]]
    where n' = fromIntegral n
          m  = round (logBase 2 n')

-- 1ª definición

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

sumaRedondeos1 :: Integer -> [Integer]

sumaRedondeos1 n =

[sum [round (n'/(fromIntegral (2^k))) | k <- [1..m]] | m <- [1..]]

where n' = fromIntegral n

limiteSumaRedondeos1 :: Integer -> Integer

limiteSumaRedondeos1 = limite . sumaRedondeos1

limite :: [Integer] -> Integer

limite xs = head [x | (x,y) <- zip xs (tail xs), x == y]

-- 2ª definición

sumaRedondeos2 :: Integer -> [Integer]

sumaRedondeos2 n =

scanl1 (+) [round (n'/(fromIntegral (2^k))) | k <- [1..]]

where n' = fromIntegral n

limiteSumaRedondeos2 :: Integer -> Integer

limiteSumaRedondeos2 = limite . sumaRedondeos2

-- 3ª definición

sumaRedondeos3 :: Integer -> [Integer]

sumaRedondeos3 n = map fst (iterate f (round (n'/2),4))

where f (s,d) = (s + round (n'/(fromIntegral d)), 2*d)

n' = fromIntegral n

limiteSumaRedondeos3 :: Integer -> Integer

limiteSumaRedondeos3 = limite . sumaRedondeos3

-- 4ª definición

sumaRedondeos4 :: Integer -> [Integer]

sumaRedondeos4 n = xs ++ repeat x

where n' = fromIntegral n

m = round (logBase 2 n')

xs = scanl1 (+) [round (n'/(fromIntegral (2^k))) | k <- [1..m]]

x = last xs

limiteSumaRedondeos4 :: Integer -> Integer

limiteSumaRedondeos4 = limite . sumaRedondeos4

-- Comparación de eficiencia

-- λ> (sumaRedondeos1 4) !! 20000

-- 3

-- (0.92 secs, 197,782,232 bytes)

-- λ> (sumaRedondeos1 4) !! 30000

-- 3

-- (2.20 secs, 351,084,816 bytes)

-- λ> (sumaRedondeos3 4) !! 20000

-- 3

-- (0.30 secs, 53,472,392 bytes)

-- λ> (sumaRedondeos4 4) !! 20000

-- 3

-- (0.01 secs, 0 bytes)

-- En lo que sigue, usaremos la 3ª definición

sumaRedondeos :: Integer -> [Integer]

sumaRedondeos = sumaRedondeos3

-- 5ª definición

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

limiteSumaRedondeos5 :: Integer -> Integer

limiteSumaRedondeos5 n =

sum [round (n'/(fromIntegral (2^k))) | k <- [1..m]]

where n' = fromIntegral n

m = round (logBase 2 n')

Año cúbico

PorJosé A. Alonso 14-diciembre-201525-marzo-2022

El año 2016 será un año cúbico porque se puede escribir como la suma de los cubos de 7 números consecutivos; en efecto,

   2016 = 3³+ 4³ +...+ 9³

1	2016 = 3³+ 4³ +...+ 9³

Definir la función

   esCubico :: Integer -> Bool

1	esCubico :: Integer -> Bool

tal que (esCubico x) se verifica si x se puede escribir como la suma de los cubos de 7 números consecutivos. Por ejemplo,

   esCubico 2016                ==  True
   esCubico 2017                ==  False
   esCubico 189005670081900441  ==  True
   esCubico 189005670081900442  ==  False

esCubico 2016 == True

esCubico 2017 == False

esCubico 189005670081900441 == True

esCubico 189005670081900442 == False

Soluciones

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

esCubico1 :: Integer -> Bool
esCubico1 x = pertenece x cubicos

-- cubicos es la lista de los números que se pueden escribir como la
-- suma de los cubos de 7 números consecutivos. Por ejemplo,
--    take 5 cubicos  ==  [784,1295,2016,2989,4256]
cubicos :: [Integer]
cubicos = [sum [x^3 | x <- [y..y+6]] | y <- [1..]]

-- (pertenece x ys) se verifica si x pertenece a la lista ordenada
-- ys. Por ejemplo,
--    pertenece 25 [0,3..]  ==  False
--    pertenece 27 [0,3..]  ==  True
pertenece :: Integer -> [Integer] ->  Bool
pertenece x ys = x == head (dropWhile (<x) ys)

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

esCubico2 :: Integer -> Bool
esCubico2 x = pertenece x (cubicosDesde k)
    where k = floor ((fromIntegral x/7)**(1/3))

cubicosDesde :: Integer -> [Integer]
cubicosDesde k = [sum [x^3 | x <- [y..y+6]] | y <- [k..]]

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

--    λ> esCubico1 189005670081900441
--    True
--    (7.49 secs, 1,917,868,024 bytes)
--    λ> esCubico2 189005670081900441
--    True
--    (0.01 secs, 0 bytes)

-- 1ª definición

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

esCubico1 :: Integer -> Bool

esCubico1 x = pertenece x cubicos

-- cubicos es la lista de los números que se pueden escribir como la

-- suma de los cubos de 7 números consecutivos. Por ejemplo,

-- take 5 cubicos == [784,1295,2016,2989,4256]

cubicos :: [Integer]

cubicos = [sum [x^3 | x <- [y..y+6]] | y <- [1..]]

-- (pertenece x ys) se verifica si x pertenece a la lista ordenada

-- ys. Por ejemplo,

-- pertenece 25 [0,3..] == False

-- pertenece 27 [0,3..] == True

pertenece :: Integer -> [Integer] -> Bool

pertenece x ys = x == head (dropWhile (<x) ys)

-- 2ª definición

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

esCubico2 :: Integer -> Bool

esCubico2 x = pertenece x (cubicosDesde k)

where k = floor ((fromIntegral x/7)**(1/3))

cubicosDesde :: Integer -> [Integer]

cubicosDesde k = [sum [x^3 | x <- [y..y+6]] | y <- [k..]]

-- Comparación de eficiencia

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

-- λ> esCubico1 189005670081900441

-- True

-- (7.49 secs, 1,917,868,024 bytes)

-- λ> esCubico2 189005670081900441

-- True

-- (0.01 secs, 0 bytes)

Raíces enteras de los números primos

PorJosé A. Alonso 2-diciembre-20151-mayo-2016

Definir la sucesión

   raicesEnterasDePrimos :: [Integer]

1	raicesEnterasDePrimos :: [Integer]

cuyos elementos son las partes enteras de las raíces cuadradas de los números primos. Por ejemplo,

   λ> take 30 raicesEnterasDePrimos
   [1,1,2,2,3,3,4,4,4,5,5,6,6,6,6,7,7,7,8,8,8,8,9,9,9,10,10,10,10,10]
   λ> raicesEnterasDePrimos !!  9963
   322
   λ> raicesEnterasDePrimos !!  9964
   323

λ> take 30 raicesEnterasDePrimos

[1,1,2,2,3,3,4,4,4,5,5,6,6,6,6,7,7,7,8,8,8,8,9,9,9,10,10,10,10,10]

λ> raicesEnterasDePrimos !! 9963

322

λ> raicesEnterasDePrimos !! 9964

323

Comprobar con QuickCheck que la diferencia entre dos términos consecutivos de la sucesión es como máximo igual a 1.

Soluciones

import Data.Numbers.Primes (primes)
import Test.QuickCheck

-- 1ª solución
-- ===========

raicesEnterasDePrimos1 :: [Integer]
raicesEnterasDePrimos1 = map raizEntera primes

-- (raizEntera x) es la parte entera de la raíz cuadrada de x. Por
-- ejemplo,
--    raizEntera  8  ==  2
--    raizEntera  9  ==  3
--    raizEntera 10  ==  3
raizEntera :: Integer -> Integer
raizEntera = floor . sqrt . fromIntegral 

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

raicesEnterasDePrimos2 :: [Integer]
raicesEnterasDePrimos2 = map raizEntera2 primes

raizEntera2 :: Integer -> Integer
raizEntera2 n = aux 1
    where aux k | k*k > n   = k-1
                | otherwise = aux (k+1)

-- 3º solución
-- ===========

raicesEnterasDePrimos3 :: [Integer]
raicesEnterasDePrimos3 = aux primes [1..]
    where aux (p:ps) (x:xs) | p > x*x   = aux (p:ps) xs
                            | otherwise = (x-1) : aux ps (x:xs)


-- Comparación de eficiencia
--    ghci> raicesEnterasDePrimos1 !! 400000
--    2408
--    (2.86 secs, 1177922500 bytes)
--    ghci> raicesEnterasDePrimos2 !! 400000
--    2408
--    (3.08 secs, 1177432260 bytes)
--    ghci> raicesEnterasDePrimos3 !! 400000
--    2408
--    (3.88 secs, 1260772112 bytes)

-- En lo sucesivo usaremos la 1ª definición
raicesEnterasDePrimos :: [Integer]
raicesEnterasDePrimos = raicesEnterasDePrimos3

-- La propiedad es
prop_raicesEnterasDePrimos :: Int -> Property
prop_raicesEnterasDePrimos n =
    n >= 0 ==> 
    raicesEnterasDePrimos !! (n+1) - raicesEnterasDePrimos !! n <= 1

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

import Data.Numbers.Primes (primes)

import Test.QuickCheck

-- 1ª solución

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

raicesEnterasDePrimos1 :: [Integer]

raicesEnterasDePrimos1 = map raizEntera primes

-- (raizEntera x) es la parte entera de la raíz cuadrada de x. Por

-- ejemplo,

-- raizEntera 8 == 2

-- raizEntera 9 == 3

-- raizEntera 10 == 3

raizEntera :: Integer -> Integer

raizEntera = floor . sqrt . fromIntegral

-- 2ª solución

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

raicesEnterasDePrimos2 :: [Integer]

raicesEnterasDePrimos2 = map raizEntera2 primes

raizEntera2 :: Integer -> Integer

raizEntera2 n = aux 1

where aux k | k*k > n = k-1

| otherwise = aux (k+1)

-- 3º solución

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

raicesEnterasDePrimos3 :: [Integer]

raicesEnterasDePrimos3 = aux primes [1..]

where aux (p:ps) (x:xs) | p > x*x = aux (p:ps) xs

| otherwise = (x-1) : aux ps (x:xs)

-- Comparación de eficiencia

-- ghci> raicesEnterasDePrimos1 !! 400000

-- 2408

-- (2.86 secs, 1177922500 bytes)

-- ghci> raicesEnterasDePrimos2 !! 400000

-- 2408

-- (3.08 secs, 1177432260 bytes)

-- ghci> raicesEnterasDePrimos3 !! 400000

-- 2408

-- (3.88 secs, 1260772112 bytes)

-- En lo sucesivo usaremos la 1ª definición

raicesEnterasDePrimos :: [Integer]

raicesEnterasDePrimos = raicesEnterasDePrimos3

-- La propiedad es

prop_raicesEnterasDePrimos :: Int -> Property

prop_raicesEnterasDePrimos n =

n >= 0 ==>

raicesEnterasDePrimos !! (n+1) - raicesEnterasDePrimos !! n <= 1

-- La comprobación es

-- λ> quickCheck prop_raicesEnterasDePrimos

-- +++ OK, passed 100 tests.

Productos de N números consecutivos

PorJosé A. Alonso 17-noviembre-201526-marzo-2022

La semana pasada se planteó en Twitter el siguiente problema

Se observa que

      1x2x3x4 = 2x3x4 
      2x3x4x5 = 4x5x6

1 2	1x2x3x4 = 2x3x4 2x3x4x5 = 4x5x6

¿Existen ejemplos de otros productos de cuatro enteros consecutivos iguales a un producto de tres enteros consecutivos?

Definir la función

   esProductoDeNconsecutivos :: Integer -> Integer -> Maybe Integer

1	esProductoDeNconsecutivos :: Integer -> Integer -> Maybe Integer

tal que (esProductoDeNconsecutivos n x) es (Just m) si x es el producto de n enteros consecutivos a partir de m y es Nothing si x no es el producto de n enteros consecutivos. Por ejemplo,

   esProductoDeNconsecutivos 3   6  == Just 1
   esProductoDeNconsecutivos 4   6  == Nothing
   esProductoDeNconsecutivos 4  24  == Just 1
   esProductoDeNconsecutivos 3  24  == Just 2
   esProductoDeNconsecutivos 3 120  == Just 4
   esProductoDeNconsecutivos 4 120  == Just 2

esProductoDeNconsecutivos 3 6 == Just 1

esProductoDeNconsecutivos 4 6 == Nothing

esProductoDeNconsecutivos 4 24 == Just 1

esProductoDeNconsecutivos 3 24 == Just 2

esProductoDeNconsecutivos 3 120 == Just 4

esProductoDeNconsecutivos 4 120 == Just 2

Para ejemplos mayores,

   λ> esProductoDeNconsecutivos 3 (product [10^20..2+10^20])
   Just 100000000000000000000
   λ> esProductoDeNconsecutivos2 4 (product [10^20..2+10^20])
   Nothing
   λ> esProductoDeNconsecutivos2 4 (product [10^20..3+10^20])
   Just 100000000000000000000

λ> esProductoDeNconsecutivos 3 (product [10^20..2+10^20])

Just 100000000000000000000

λ> esProductoDeNconsecutivos2 4 (product [10^20..2+10^20])

Nothing

λ> esProductoDeNconsecutivos2 4 (product [10^20..3+10^20])

Just 100000000000000000000

Usando la función esProductoDeNconsecutivos resolver el problema.

Soluciones

import Data.Maybe

-- 1ª definición
esProductoDeNconsecutivos1 :: Integer -> Integer -> Maybe Integer
esProductoDeNconsecutivos1 n x 
    | null productos = Nothing
    | otherwise      = Just (head productos)
    where productos = [m | m <- [1..x-n], product [m..m+n-1] == x]

-- 2ª definición
esProductoDeNconsecutivos2 :: Integer -> Integer -> Maybe Integer
esProductoDeNconsecutivos2 n x = aux k
    where k = floor (fromIntegral x ** (1/(fromIntegral n))) - (n `div` 2)
          aux m | y == x    = Just m
                | y <  x    = aux (m+1)
                | otherwise = Nothing
                where y = product [m..m+n-1]

-- Comparación de eficiencia
--    λ> esProductoDeNconsecutivos1 3 (product [10^7..2+10^7])
--    Just 10000000
--    (12.37 secs, 5678433692 bytes)
--    λ> esProductoDeNconsecutivos2 3 (product [10^7..2+10^7])
--    Just 10000000
--    (0.00 secs, 1554932 bytes)

-- Solución del problema
-- =====================

soluciones :: [Integer]
soluciones = [x | x <- [121..]
                , isJust (esProductoDeNconsecutivos2 4 x)
                , isJust (esProductoDeNconsecutivos2 3 x)]

-- El cálculo es
--    λ> head soluciones
--    175560
--    λ> esProductoDeNconsecutivos2 4 175560
--    Just 19
--    λ> esProductoDeNconsecutivos2 3 175560
--    Just 55
--    λ> product [19,20,21,22] 
--    175560
--    λ> product [55,56,57]
--    175560
--    λ> product [19,20,21,22] == product [55,56,57]
--    True

-- Se puede definir una función para automatizar el proceso anterior:
soluciones2 :: [(Integer,[Integer],[Integer])]
soluciones2 = [(x,[a..a+3],[b..b+2]) 
               | x <- [121..]
               , let y = esProductoDeNconsecutivos2 4 x
               , isJust y
               , let z = esProductoDeNconsecutivos2 3 x
               , isJust z
               , let a = fromJust y
               , let b = fromJust z
               ]

-- El cálculo es 
--    λ> head soluciones2
--    (175560,[19,20,21,22],[55,56,57])

import Data.Maybe

-- 1ª definición

esProductoDeNconsecutivos1 :: Integer -> Integer -> Maybe Integer

esProductoDeNconsecutivos1 n x

| null productos = Nothing

| otherwise = Just (head productos)

where productos = [m | m <- [1..x-n], product [m..m+n-1] == x]

-- 2ª definición

esProductoDeNconsecutivos2 :: Integer -> Integer -> Maybe Integer

esProductoDeNconsecutivos2 n x = aux k

where k = floor (fromIntegral x ** (1/(fromIntegral n))) - (n `div` 2)

aux m | y == x = Just m

| y < x = aux (m+1)

| otherwise = Nothing

where y = product [m..m+n-1]

-- Comparación de eficiencia

-- λ> esProductoDeNconsecutivos1 3 (product [10^7..2+10^7])

-- Just 10000000

-- (12.37 secs, 5678433692 bytes)

-- λ> esProductoDeNconsecutivos2 3 (product [10^7..2+10^7])

-- Just 10000000

-- (0.00 secs, 1554932 bytes)

-- Solución del problema

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

soluciones :: [Integer]

soluciones = [x | x <- [121..]

, isJust (esProductoDeNconsecutivos2 4 x)

, isJust (esProductoDeNconsecutivos2 3 x)]

-- El cálculo es

-- λ> head soluciones

-- 175560

-- λ> esProductoDeNconsecutivos2 4 175560

-- Just 19

-- λ> esProductoDeNconsecutivos2 3 175560

-- Just 55

-- λ> product [19,20,21,22]

-- 175560

-- λ> product [55,56,57]

-- 175560

-- λ> product [19,20,21,22] == product [55,56,57]

-- True

-- Se puede definir una función para automatizar el proceso anterior:

soluciones2 :: [(Integer,[Integer],[Integer])]

soluciones2 = [(x,[a..a+3],[b..b+2])

| x <- [121..]

, let y = esProductoDeNconsecutivos2 4 x

, isJust y

, let z = esProductoDeNconsecutivos2 3 x

, isJust z

, let a = fromJust y

, let b = fromJust z

]

-- El cálculo es

-- λ> head soluciones2

-- (175560,[19,20,21,22],[55,56,57])

Números libres de cuadrados

PorJosé A. Alonso 3-noviembre-201525-abril-2021

Un número es libre de cuadrados si no es divisible el cuadrado de ningún entero mayor que 1. Por ejemplo, 70 es libre de cuadrado porque sólo es divisible por 1, 2, 5, 7 y 70; en cambio, 40 no es libre de cuadrados porque es divisible por 2².

Definir la función

   libreDeCuadrados :: Integer -> Bool

1	libreDeCuadrados :: Integer -> Bool

tal que (libreDeCuadrados x) se verifica si x es libre de cuadrados. Por ejemplo,

   libreDeCuadrados 70                 ==  True
   libreDeCuadrados 40                 ==  False
   libreDeCuadrados 510510             ==  True
   libreDeCuadrados (((10^10)^10)^10)  ==  False

libreDeCuadrados 70 == True

libreDeCuadrados 40 == False

libreDeCuadrados 510510 == True

libreDeCuadrados (((10^10)^10)^10) == False

Soluciones

-- 1ª definición:
libreDeCuadrados :: Integer -> Bool
libreDeCuadrados x = x == product (divisoresPrimos x)

-- (divisoresPrimos x) es la lista de los divisores primos de x. Por
-- ejemplo,  
--    divisoresPrimos 40  ==  [2,5]
--    divisoresPrimos 70  ==  [2,5,7]
divisoresPrimos :: Integer -> [Integer]
divisoresPrimos x = [n | n <- divisores x, primo n]

-- (divisores n) es la lista de los divisores del número n. Por ejemplo,
--    divisores 30  ==  [1,2,3,5,6,10,15,30]  
divisores :: Integer -> [Integer]
divisores n = [x | x <- [1..n], n `mod` x == 0]

-- (primo n) se verifica si n es primo. Por ejemplo,
--    primo 30  == False
--    primo 31  == True  
primo :: Integer -> Bool
primo n = divisores n == [1, n]

-- 2ª definición
libreDeCuadrados2 :: Integer -> Bool
libreDeCuadrados2 n = 
    null [x | x <- [2..n], rem n (x^2) == 0]

-- 3ª definición
libreDeCuadrados3 :: Integer -> Bool
libreDeCuadrados3 n = 
    null [x | x <- [2..floor (sqrt (fromIntegral n))], 
              rem n (x^2) == 0]

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

--    λ> libreDeCuadrados 510510
--    True
--    (0.76 secs, 89,522,360 bytes)
--    λ> libreDeCuadrados2 510510
--    True
--    (1.78 secs, 371,826,320 bytes)
--    λ> libreDeCuadrados3 510510
--    True
--    (0.01 secs, 0 bytes)

-- 1ª definición:

libreDeCuadrados :: Integer -> Bool

libreDeCuadrados x = x == product (divisoresPrimos x)

-- (divisoresPrimos x) es la lista de los divisores primos de x. Por

-- ejemplo,

-- divisoresPrimos 40 == [2,5]

-- divisoresPrimos 70 == [2,5,7]

divisoresPrimos :: Integer -> [Integer]

divisoresPrimos x = [n | n <- divisores x, primo n]

-- (divisores n) es la lista de los divisores del número n. Por ejemplo,

-- divisores 30 == [1,2,3,5,6,10,15,30]

divisores :: Integer -> [Integer]

divisores n = [x | x <- [1..n], n `mod` x == 0]

-- (primo n) se verifica si n es primo. Por ejemplo,

-- primo 30 == False

-- primo 31 == True

primo :: Integer -> Bool

primo n = divisores n == [1, n]

-- 2ª definición

libreDeCuadrados2 :: Integer -> Bool

libreDeCuadrados2 n =

null [x | x <- [2..n], rem n (x^2) == 0]

-- 3ª definición

libreDeCuadrados3 :: Integer -> Bool

libreDeCuadrados3 n =

null [x | x <- [2..floor (sqrt (fromIntegral n))],

rem n (x^2) == 0]

-- Comparación de eficiencia

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

-- λ> libreDeCuadrados 510510

-- True

-- (0.76 secs, 89,522,360 bytes)

-- λ> libreDeCuadrados2 510510

-- True

-- (1.78 secs, 371,826,320 bytes)

-- λ> libreDeCuadrados3 510510

-- True

-- (0.01 secs, 0 bytes)

Perímetro más frecuente de triángulos rectángulos

PorJosé A. Alonso 7-abril-20151-mayo-2015

El grado perimetral de un número p es la cantidad de tres triángulos rectángulos de lados enteros cuyo perímetro es p. Por ejemplo, el grado perimetral de 120 es 3 ya que sólo hay 3 triángulos rectángulos de lados enteros cuyo perímetro es 120: {20,48,52}, {24,45,51} y {30,40,50}.

Definir la función

   maxGradoPerimetral :: Int -> (Int,[Int])

1	maxGradoPerimetral :: Int -> (Int,[Int])

tal que (maxGradoPerimetral n) es el par (m,ps) tal que m es el máximo grado perimetral de los números menores o iguales que n y ps son los perímetros, menores o iguales que n, cuyo grado perimetral es m. Por ejemplo,

   maxGradoPerimetral   50  ==  (1,[12,24,30,36,40,48])
   maxGradoPerimetral  100  ==  (2,[60,84,90])
   maxGradoPerimetral  200  ==  (3,[120,168,180])
   maxGradoPerimetral  400  ==  (4,[240,360])
   maxGradoPerimetral  500  ==  (5,[420])
   maxGradoPerimetral  750  ==  (6,[720])
   maxGradoPerimetral  839  ==  (6,[720])
   maxGradoPerimetral  840  ==  (8,[840])
   maxGradoPerimetral 1500  ==  (8,[840,1260])
   maxGradoPerimetral 2000  ==  (10,[1680])
   maxGradoPerimetral 3000  ==  (12,[2520])

maxGradoPerimetral 50 == (1,[12,24,30,36,40,48])

maxGradoPerimetral 100 == (2,[60,84,90])

maxGradoPerimetral 200 == (3,[120,168,180])

maxGradoPerimetral 400 == (4,[240,360])

maxGradoPerimetral 500 == (5,[420])

maxGradoPerimetral 750 == (6,[720])

maxGradoPerimetral 839 == (6,[720])

maxGradoPerimetral 840 == (8,[840])

maxGradoPerimetral 1500 == (8,[840,1260])

maxGradoPerimetral 2000 == (10,[1680])

maxGradoPerimetral 3000 == (12,[2520])

Soluciones

import Data.List
import Data.Array (accumArray, assocs)

-- 1ª solución                                                      --
-- ===========

maxGradoPerimetral1 :: Int -> (Int,[Int])
maxGradoPerimetral1 p = (m,[x | (n,x) <- ts, n == m])
    where ts    = [(length (triangulos x),x) | x <- [1..p]] 
          (m,_) = maximum ts 

-- (triangulos p) es el conjunto de triángulos rectángulos de perímetro
-- p. Por ejemplo,
--    triangulos 120  ==  [(20,48,52),(24,45,51),(30,40,50)]
triangulos :: Int -> [(Int,Int,Int)]
triangulos p = 
    [(a,b,c) | a <- [1..q],
               b <- [a..q],
               let c = p-a-b,
               a*a+b*b == c*c]
    where q = p `div` 2

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

maxGradoPerimetral2 :: Int -> (Int,[Int])
maxGradoPerimetral2 p = (m,[x | (n,x) <- ts, n == m])
    where ts    = [(n,x) | (x,n) <- numeroPerimetrosTriangulos p, n > 0]
          (m,_) = maximum ts 

-- (numeroPerimetrosTriangulos p) es la lista formado por los números de
-- 1 a p y la cantidad de triángulos rectángulos enteros cuyo perímetro
-- es dicho número. Por ejemplo,
--    ghci>  [(p,n) | (p,n) <- numeroPerimetrosTriangulos 70, n > 0]
--    [(12,1),(24,1),(30,1),(36,1),(40,1),(48,1),(56,1),(60,2),(70,1)]
numeroPerimetrosTriangulos :: Int -> [(Int,Int)] 
numeroPerimetrosTriangulos p = 
    assocs (accumArray (\x _ -> 1+x) 0 (1,p) (perimetrosTriangulos p))

-- (perimetrosTriangulos p) es la lista formada por los perímetros y los
-- lados de los triángulos rectángulos enteros cuyo perímetro es menor o
-- igual que p. Por ejemplo,
--    ghci> perimetrosTriangulos 70
--    [(12,(3,4,5)),   (30,(5,12,13)),(24,(6,8,10)),  (56,(7,24,25)),
--     (40,(8,15,17)), (36,(9,12,15)),(60,(10,24,26)),(48,(12,16,20)),
--     (60,(15,20,25)),(70,(20,21,29))]
perimetrosTriangulos :: Int -> [(Int,(Int,Int,Int))]
perimetrosTriangulos p =
    [(q,(a,b,c)) | a <- [1..p1],
                   b <- [a..p1],
                   esCuadrado (a*a+b*b), 
                   let c = raizCuadradaE (a*a+b*b), 
                   let q = a+b+c,
                   q <= p]
    where p1 = p `div` 2

-- (esCuadrado n) se verifica si n es un cuadrado. Por ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 27  ==  False
esCuadrado :: Int -> Bool
esCuadrado n = a*a == n
    where a = raizCuadradaE n

-- (raizCuadradaE n) es la raíz cuadrada entera de n. Por ejemplo,
--    raizCuadradaE 25  ==  5
--    raizCuadradaE 27  ==  5
--    raizCuadradaE 35  ==  5
--    raizCuadradaE 36  ==  6
raizCuadradaE :: Int -> Int
raizCuadradaE = floor . sqrt . fromIntegral

-- 3ª solución                                                      --
-- ===========

maxGradoPerimetral3 :: Int -> (Int,[Int])
maxGradoPerimetral3 p = (m,[x | (n,x) <- ts, n == m])
    where ts    = [(n,x) | (x,n) <- numeroPerimetrosTriangulos2 p, n > 0]
          (m,_) = maximum ts 

-- (numeroPerimetrosTriangulos2 p) es la lista formado por los números de
-- 1 a p y la cantidad de triángulos rectángulos enteros cuyo perímetro
-- es dicho número. Por ejemplo,
--    ghci>  [(p,n) | (p,n) <- numeroPerimetrosTriangulos2 70, n > 0]
--    [(12,1),(24,1),(30,1),(36,1),(40,1),(48,1),(56,1),(60,2),(70,1)]
numeroPerimetrosTriangulos2 :: Int -> [(Int,Int)] 
numeroPerimetrosTriangulos2 p = 
    [(head xs, length xs) | xs <- group (sort (perimetrosTriangulos2 p))]

-- (perimetrosTriangulos2 p) es la lista formada por los perímetros de
-- los triángulos rectángulos enteros cuyo perímetro es menor o igual
-- que p. Por ejemplo, 
--    perimetrosTriangulos2 70  ==  [12,30,24,56,40,36,60,48,60,70]
perimetrosTriangulos2 :: Int -> [Int]
perimetrosTriangulos2 p =
    [q | a <- [1..p1],
         b <- [a..p1],
         esCuadrado (a*a+b*b), 
         let c = raizCuadradaE (a*a+b*b), 
         let q = a+b+c,
         q <= p]
    where p1 = p `div` 2

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

-- La comparación es
--    ghci> maxGradoPerimetral1 1000
--    (8,[840])
--    (120.08 secs, 21116625136 bytes)
--    ghci> maxGradoPerimetral2 1000
--    (8,[840])
--    (0.66 secs, 132959056 bytes)
--    ghci> maxGradoPerimetral3 1000
--    (1000,[1])
--    (0.66 secs, 133443816 bytes)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

import Data.List

import Data.Array (accumArray, assocs)

-- 1ª solución --

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

maxGradoPerimetral1 :: Int -> (Int,[Int])

maxGradoPerimetral1 p = (m,[x | (n,x) <- ts, n == m])

where ts = [(length (triangulos x),x) | x <- [1..p]]

(m,_) = maximum ts

-- (triangulos p) es el conjunto de triángulos rectángulos de perímetro

-- p. Por ejemplo,

-- triangulos 120 == [(20,48,52),(24,45,51),(30,40,50)]

triangulos :: Int -> [(Int,Int,Int)]

triangulos p =

[(a,b,c) | a <- [1..q],

b <- [a..q],

let c = p-a-b,

a*a+b*b == c*c]

where q = p `div` 2

-- 2ª solución --

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

maxGradoPerimetral2 :: Int -> (Int,[Int])

maxGradoPerimetral2 p = (m,[x | (n,x) <- ts, n == m])

where ts = [(n,x) | (x,n) <- numeroPerimetrosTriangulos p, n > 0]

(m,_) = maximum ts

-- (numeroPerimetrosTriangulos p) es la lista formado por los números de

-- 1 a p y la cantidad de triángulos rectángulos enteros cuyo perímetro

-- es dicho número. Por ejemplo,

-- ghci> [(p,n) | (p,n) <- numeroPerimetrosTriangulos 70, n > 0]

-- [(12,1),(24,1),(30,1),(36,1),(40,1),(48,1),(56,1),(60,2),(70,1)]

numeroPerimetrosTriangulos :: Int -> [(Int,Int)]

numeroPerimetrosTriangulos p =

assocs (accumArray (\x _ -> 1+x) 0 (1,p) (perimetrosTriangulos p))

-- (perimetrosTriangulos p) es la lista formada por los perímetros y los

-- lados de los triángulos rectángulos enteros cuyo perímetro es menor o

-- igual que p. Por ejemplo,

-- ghci> perimetrosTriangulos 70

-- [(12,(3,4,5)), (30,(5,12,13)),(24,(6,8,10)), (56,(7,24,25)),

-- (40,(8,15,17)), (36,(9,12,15)),(60,(10,24,26)),(48,(12,16,20)),

-- (60,(15,20,25)),(70,(20,21,29))]

perimetrosTriangulos :: Int -> [(Int,(Int,Int,Int))]

perimetrosTriangulos p =

[(q,(a,b,c)) | a <- [1..p1],

b <- [a..p1],

esCuadrado (a*a+b*b),

let c = raizCuadradaE (a*a+b*b),

let q = a+b+c,

q <= p]

where p1 = p `div` 2

-- (esCuadrado n) se verifica si n es un cuadrado. Por ejemplo,

-- esCuadrado 25 == True

-- esCuadrado 27 == False

esCuadrado :: Int -> Bool

esCuadrado n = a*a == n

where a = raizCuadradaE n

-- (raizCuadradaE n) es la raíz cuadrada entera de n. Por ejemplo,

-- raizCuadradaE 25 == 5

-- raizCuadradaE 27 == 5

-- raizCuadradaE 35 == 5

-- raizCuadradaE 36 == 6

raizCuadradaE :: Int -> Int

raizCuadradaE = floor . sqrt . fromIntegral

-- 3ª solución --

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

maxGradoPerimetral3 :: Int -> (Int,[Int])

maxGradoPerimetral3 p = (m,[x | (n,x) <- ts, n == m])

where ts = [(n,x) | (x,n) <- numeroPerimetrosTriangulos2 p, n > 0]

(m,_) = maximum ts

-- (numeroPerimetrosTriangulos2 p) es la lista formado por los números de

-- 1 a p y la cantidad de triángulos rectángulos enteros cuyo perímetro

-- es dicho número. Por ejemplo,

-- ghci> [(p,n) | (p,n) <- numeroPerimetrosTriangulos2 70, n > 0]

-- [(12,1),(24,1),(30,1),(36,1),(40,1),(48,1),(56,1),(60,2),(70,1)]

numeroPerimetrosTriangulos2 :: Int -> [(Int,Int)]

numeroPerimetrosTriangulos2 p =

[(head xs, length xs) | xs <- group (sort (perimetrosTriangulos2 p))]

-- (perimetrosTriangulos2 p) es la lista formada por los perímetros de

-- los triángulos rectángulos enteros cuyo perímetro es menor o igual

-- que p. Por ejemplo,

-- perimetrosTriangulos2 70 == [12,30,24,56,40,36,60,48,60,70]

perimetrosTriangulos2 :: Int -> [Int]

perimetrosTriangulos2 p =

[q | a <- [1..p1],

b <- [a..p1],

esCuadrado (a*a+b*b),

let c = raizCuadradaE (a*a+b*b),

let q = a+b+c,

q <= p]

where p1 = p `div` 2

-- Comparación de eficiencia

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

-- La comparación es

-- ghci> maxGradoPerimetral1 1000

-- (8,[840])

-- (120.08 secs, 21116625136 bytes)

-- ghci> maxGradoPerimetral2 1000

-- (8,[840])

-- (0.66 secs, 132959056 bytes)

-- ghci> maxGradoPerimetral3 1000

-- (1000,[1])

-- (0.66 secs, 133443816 bytes)

Números de suma prima hereditarios por la derecha

PorJosé A. Alonso 27-marzo-20151-mayo-2015

Decimos que un número es de suma prima si la suma de todos sus dígitos es un número primo. Por ejemplo el número 562 es de suma prima pues la suma de sus dígitos es el número primo 13; sin embargo, el número 514 no es de suma prima pues la suma de sus dígitos es 10, que no es primo.

Decimos que un número es de suma prima hereditario por la derecha si es de suma prima y los números que se obtienen eliminando sus últimas cifras también son de suma prima. Por ejemplo 7426 es de suma prima hereditario por la derecha pues 7426, 742, 74 y 7 son todos números de suma prima.

Definir la constante

   listaSumaPrimaHD :: [Integer]

1	listaSumaPrimaHD :: [Integer]

cuyo valor es la lista infinita de los números de suma prima hereditarios por la derecha. Por ejemplo,

   take 10 listaSumaPrimaHD     ==  [2,3,5,7,20,21,23,25,29,30]
   listaSumaPrimaHD !! 2000000  ==  3800024668046

1 2	take 10 listaSumaPrimaHD == [2,3,5,7,20,21,23,25,29,30] listaSumaPrimaHD !! 2000000 == 3800024668046

Soluciones

import Data.Char (digitToInt)
import Data.Array
import Data.Numbers.Primes (isPrime)

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

listaSumaPrimaHD1 :: [Integer]
listaSumaPrimaHD1 = filter sumaPrimaHD [1..]

-- (sumaPrimaHD n) se verifica si n es de suma prima hereditario por la
-- derecha. Por ejemplo,
--    sumaPrimaHD 7426  ==  True
--    sumaPrimaHD 7427  ==  False
sumaPrimaHD n
    | n < 10    = isPrime n
    | otherwise = sumaPrima n && sumaPrimaHD (n `div` 10)

-- (sumaPrima n) se verifica si n es un número de suma prima. Por
-- ejemplo, 
--    sumaPrima 562  ==  True
--    sumaPrima 514  ==  False
sumaPrima :: Integer -> Bool
sumaPrima = isPrime . sum . digitos

-- (digitos n) es la lista de los dígitos de n. Por ejemplo,
--    digitos 325  ==  [3,2,5]
digitos :: Integer -> [Integer]
digitos = map (fromIntegral . digitToInt) . show


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

listaSumaPrimaHD2 :: [Integer]
listaSumaPrimaHD2 = map fst paresSumaPrimaHDDigitos

paresSumaPrimaHDDigitos :: [(Integer, Integer)]
paresSumaPrimaHDDigitos =
    paresSumaPrimaHDDigitosAux 1 [(2,2),(3,3),(5,5),(7,7)]

paresSumaPrimaHDDigitosAux :: Integer -> [(Integer,Integer)] -> 
                              [(Integer,Integer)]
paresSumaPrimaHDDigitosAux n ac =
    ac ++ paresSumaPrimaHDDigitosAux (n+1) 
                                     (concatMap extiendeSumaPrimaHD ac)

extiendeSumaPrimaHD :: (Integer,Integer) -> [(Integer,Integer)]
extiendeSumaPrimaHD (n,s) = [(n*10+k,s+k) | k <- [0..9], isPrime (s+k)]

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

listaSumaPrimaHD3 :: [Integer]
listaSumaPrimaHD3 = 
    map fst (concat (iterate (concatMap extiendeSumaPrimaHD3) 
                             [(2,2),(3,3),(5,5),(7,7)]))

extiendeSumaPrimaHD3 :: (Integer,Integer) -> [(Integer,Integer)]
extiendeSumaPrimaHD3 (n,s) = [(n*10+k,s+k) | k <- extensiones ! s]

extensiones :: Array Integer [Integer]
extensiones = array (1,1000) 
              [(n,[k | k <- [0..9], isPrime (n+k)]) | n <- [1..1000]]

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

--    ghci> listaSumaPrimaHD1 !! 600
--    34004
--    (2.47 secs, 1565301720 bytes)
--    ghci> listaSumaPrimaHD2 !! 600
--    34004
--    (0.02 secs, 7209000 bytes)
--    ghci> listaSumaPrimaHD3 !! 600
--    34004
--    (0.01 secs, 1579920 bytes)
-- 
--    ghci> listaSumaPrimaHD2 !! 2000000
--    3800024668046
--    (45.41 secs, 29056613824 bytes)
--    ghci> listaSumaPrimaHD3 !! 2000000
--    3800024668046
--    (4.29 secs, 973265400 bytes)

import Data.Char (digitToInt)

import Data.Array

import Data.Numbers.Primes (isPrime)

-- 1ª definición

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

listaSumaPrimaHD1 :: [Integer]

listaSumaPrimaHD1 = filter sumaPrimaHD [1..]

-- (sumaPrimaHD n) se verifica si n es de suma prima hereditario por la

-- derecha. Por ejemplo,

-- sumaPrimaHD 7426 == True

-- sumaPrimaHD 7427 == False

sumaPrimaHD n

| n < 10 = isPrime n

| otherwise = sumaPrima n && sumaPrimaHD (n `div` 10)

-- (sumaPrima n) se verifica si n es un número de suma prima. Por

-- ejemplo,

-- sumaPrima 562 == True

-- sumaPrima 514 == False

sumaPrima :: Integer -> Bool

sumaPrima = isPrime . sum . digitos

-- (digitos n) es la lista de los dígitos de n. Por ejemplo,

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

digitos :: Integer -> [Integer]

digitos = map (fromIntegral . digitToInt) . show

-- 2ª definición

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

listaSumaPrimaHD2 :: [Integer]

listaSumaPrimaHD2 = map fst paresSumaPrimaHDDigitos

paresSumaPrimaHDDigitos :: [(Integer, Integer)]

paresSumaPrimaHDDigitos =

paresSumaPrimaHDDigitosAux 1 [(2,2),(3,3),(5,5),(7,7)]

paresSumaPrimaHDDigitosAux :: Integer -> [(Integer,Integer)] ->

[(Integer,Integer)]

paresSumaPrimaHDDigitosAux n ac =

ac ++ paresSumaPrimaHDDigitosAux (n+1)

(concatMap extiendeSumaPrimaHD ac)

extiendeSumaPrimaHD :: (Integer,Integer) -> [(Integer,Integer)]

extiendeSumaPrimaHD (n,s) = [(n*10+k,s+k) | k <- [0..9], isPrime (s+k)]

-- 3ª definición

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

listaSumaPrimaHD3 :: [Integer]

listaSumaPrimaHD3 =

map fst (concat (iterate (concatMap extiendeSumaPrimaHD3)

[(2,2),(3,3),(5,5),(7,7)]))

extiendeSumaPrimaHD3 :: (Integer,Integer) -> [(Integer,Integer)]

extiendeSumaPrimaHD3 (n,s) = [(n*10+k,s+k) | k <- extensiones ! s]

extensiones :: Array Integer [Integer]

extensiones = array (1,1000)

[(n,[k | k <- [0..9], isPrime (n+k)]) | n <- [1..1000]]

-- Comparación de eficiencia

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

-- ghci> listaSumaPrimaHD1 !! 600

-- 34004

-- (2.47 secs, 1565301720 bytes)

-- ghci> listaSumaPrimaHD2 !! 600

-- 34004

-- (0.02 secs, 7209000 bytes)

-- ghci> listaSumaPrimaHD3 !! 600

-- 34004

-- (0.01 secs, 1579920 bytes)

-- ghci> listaSumaPrimaHD2 !! 2000000

-- 3800024668046

-- (45.41 secs, 29056613824 bytes)

-- ghci> listaSumaPrimaHD3 !! 2000000

-- 3800024668046

-- (4.29 secs, 973265400 bytes)

Listas con los ceros emparejados

PorJosé A. Alonso 4-marzo-201511-marzo-2015

Sea S un conjunto de números. Las listas de ceros emparejados de S son las listas formadas con los elementos de S y en las cuales los ceros aparecen en sublistas de longitud par. Por ejemplo, si S = {0,1,2} entonces [1], [2], [2,1], [2,0,0,2,0,0,1] y [0,0,0,0,1,2] son listas de ceros emparejados de S; pero [0,0,0,2,1,0,0] y [0,0,1,0,1] no lo son.

Definir las funciones

   cerosEmparejados  :: Integer -> Integer -> [[Integer]]
   nCerosEmparejados :: Integer -> Integer -> Integer

1 2	cerosEmparejados :: Integer -> Integer -> [[Integer]] nCerosEmparejados :: Integer -> Integer -> Integer

tales que
+ (cerosEmparejados m n) es la lista de las listas de longitud n de ceros emparejados con los números 0, 1, 2,…, m. Por ejemplo,

     ghci> cerosEmparejados 2 0
     [[]]
     ghci> cerosEmparejados 2 1
     [[1],[2]]
     ghci> cerosEmparejados 3 1
     [[1],[2],[3]]
     ghci> cerosEmparejados 2 2
     [[1,1],[1,2],[2,1],[2,2],[0,0]]
     ghci> cerosEmparejados 2 3
     [[1,1,1],[1,1,2],[1,2,1],[1,2,2],[1,0,0],[2,1,1],[2,1,2],
      [2,2,1],[2,2,2],[2,0,0],[0,0,1],[0,0,2]]
     ghci> cerosEmparejados 2 4
     [[1,1,1,1],[1,1,1,2],[1,1,2,1],[1,1,2,2],[1,1,0,0],[1,2,1,1],
      [1,2,1,2],[1,2,2,1],[1,2,2,2],[1,2,0,0],[1,0,0,1],[1,0,0,2],
      [2,1,1,1],[2,1,1,2],[2,1,2,1],[2,1,2,2],[2,1,0,0],[2,2,1,1],
      [2,2,1,2],[2,2,2,1],[2,2,2,2],[2,2,0,0],[2,0,0,1],[2,0,0,2],
      [0,0,1,1],[0,0,1,2],[0,0,2,1],[0,0,2,2],[0,0,0,0]]

ghci> cerosEmparejados 2 0

[[]]

ghci> cerosEmparejados 2 1

[[1],[2]]

ghci> cerosEmparejados 3 1

[[1],[2],[3]]

ghci> cerosEmparejados 2 2

[[1,1],[1,2],[2,1],[2,2],[0,0]]

ghci> cerosEmparejados 2 3

[[1,1,1],[1,1,2],[1,2,1],[1,2,2],[1,0,0],[2,1,1],[2,1,2],

[2,2,1],[2,2,2],[2,0,0],[0,0,1],[0,0,2]]

ghci> cerosEmparejados 2 4

[[1,1,1,1],[1,1,1,2],[1,1,2,1],[1,1,2,2],[1,1,0,0],[1,2,1,1],

[1,2,1,2],[1,2,2,1],[1,2,2,2],[1,2,0,0],[1,0,0,1],[1,0,0,2],

[2,1,1,1],[2,1,1,2],[2,1,2,1],[2,1,2,2],[2,1,0,0],[2,2,1,1],

[2,2,1,2],[2,2,2,1],[2,2,2,2],[2,2,0,0],[2,0,0,1],[2,0,0,2],

[0,0,1,1],[0,0,1,2],[0,0,2,1],[0,0,2,2],[0,0,0,0]]

(nCerosEmparejados m n) es el número de listas de longitud n de ceros emparejados con los números 0, 1, 2,…, m. Por ejemplo,

     nCerosEmparejados 2 2   ==  5
     nCerosEmparejados 2 3   ==  12
     nCerosEmparejados 2 4   ==  29
     nCerosEmparejados 9 27  ==  79707842493701635611689499

nCerosEmparejados 2 2 == 5

nCerosEmparejados 2 3 == 12

nCerosEmparejados 2 4 == 29

nCerosEmparejados 9 27 == 79707842493701635611689499

Soluciones

import Data.List (genericIndex)
import Data.List (genericLength) 

cerosEmparejados :: Integer -> Integer -> [[Integer]]
cerosEmparejados m 0 = [[]]
cerosEmparejados m 1 = [[k] | k <- [1..m]]
cerosEmparejados m n = 
    [x:ys | x <- [1..m], ys <- cerosEmparejados m (n-1)] ++
    [0:0:ys | ys <- cerosEmparejados m (n-2)]

-- 1ª definición de nCerosEmparejados: 
nCerosEmparejados1 :: Integer -> Integer -> Integer
nCerosEmparejados1 m n = fromIntegral (length (cerosEmparejados m n))

-- 2ª definición de nCerosEmparejados: 
nCerosEmparejados2 :: Integer -> Integer -> Integer
nCerosEmparejados2 _ 0 = 1
nCerosEmparejados2 m 1 = m
nCerosEmparejados2 m n = 
    m * nCerosEmparejados2 m (n-1) + nCerosEmparejados2 m (n-2)

-- 3ª definición de nCerosEmparejados: 
nCerosEmparejados3 :: Integer -> Integer -> Integer
nCerosEmparejados3 m n = aux `genericIndex` n
    where aux = 1 : m : zipWith (\x y -> x+m*y) aux (tail aux)

-- Comparación de eficiencia
--    ghci> nCerosEmparejados1 9 7
--    5144589
--    (22.30 secs, 6556279384 bytes)
--    ghci> nCerosEmparejados2 9 7
--    5144589
--    (0.01 secs, 515464 bytes)
--    ghci> nCerosEmparejados3 9 7
--    5144589
--    (0.01 secs, 518104 bytes)
--
--    ghci> nCerosEmparejados2 9 33
--    45556060883025783396845717812863
--    (21.59 secs, 2676070480 bytes)
--    ghci> nCerosEmparejados3 9 33
--    45556060883025783396845717812863
--    (0.00 secs, 1017240 bytes)

import Data.List (genericIndex)

import Data.List (genericLength)

cerosEmparejados :: Integer -> Integer -> [[Integer]]

cerosEmparejados m 0 = [[]]

cerosEmparejados m 1 = [[k] | k <- [1..m]]

cerosEmparejados m n =

[x:ys | x <- [1..m], ys <- cerosEmparejados m (n-1)] ++

[0:0:ys | ys <- cerosEmparejados m (n-2)]

-- 1ª definición de nCerosEmparejados:

nCerosEmparejados1 :: Integer -> Integer -> Integer

nCerosEmparejados1 m n = fromIntegral (length (cerosEmparejados m n))

-- 2ª definición de nCerosEmparejados:

nCerosEmparejados2 :: Integer -> Integer -> Integer

nCerosEmparejados2 _ 0 = 1

nCerosEmparejados2 m 1 = m

nCerosEmparejados2 m n =

m * nCerosEmparejados2 m (n-1) + nCerosEmparejados2 m (n-2)

-- 3ª definición de nCerosEmparejados:

nCerosEmparejados3 :: Integer -> Integer -> Integer

nCerosEmparejados3 m n = aux `genericIndex` n

where aux = 1 : m : zipWith (\x y -> x+m*y) aux (tail aux)

-- Comparación de eficiencia

-- ghci> nCerosEmparejados1 9 7

-- 5144589

-- (22.30 secs, 6556279384 bytes)

-- ghci> nCerosEmparejados2 9 7

-- 5144589

-- (0.01 secs, 515464 bytes)

-- ghci> nCerosEmparejados3 9 7

-- 5144589

-- (0.01 secs, 518104 bytes)

-- ghci> nCerosEmparejados2 9 33

-- 45556060883025783396845717812863

-- (21.59 secs, 2676070480 bytes)

-- ghci> nCerosEmparejados3 9 33

-- 45556060883025783396845717812863

-- (0.00 secs, 1017240 bytes)

2015 como raíz cuadrada de la suma de tres cubos

PorJosé A. Alonso 19-diciembre-20141-mayo-2016

Todos los años, en las proximidades del final de año suelen aparecer cuestiones con propiedades del número del nuevo año. Una sobre el 2015 es la publicada el martes en la entrada 2015 como raíz de la suma de tres cubos del blog Números y algo más en la que se pide calcular tres números tales que 2015 sea igual a la raíz cuadrada de la suma de dichos tres números.

A partir de dicha entrada, se propone el siguiente problema: Definir la sucesión

   raicesCuadradasDeSumasDe3Cubos :: [Integer]

1	raicesCuadradasDeSumasDe3Cubos :: [Integer]

cuyos elementos son los números que se pueden escribir como raíces cuadradas de sumas de tres cubos. Por ejemplo,

   take 9 raicesCuadradasDeSumasDe3Cubos = [6,9,15,27,48,72,53,59,78]

1	take 9 raicesCuadradasDeSumasDe3Cubos = [6,9,15,27,48,72,53,59,78]

El 6 está en la sucesión porque 1³+2³+3³ = 36 y la raíz cuadrada de36 es 6 y el 9 está porque 3³+3³+3³ = 81 y la raíz cuadrada de 81 es 9. Algunos números tienen varias descomposiones como raíz cuadrada de suma de tres cubos; por ejemplo, el 71 se puede escribir como la raíz cuadrada de la suma de los cubos de 6, 9 y 16 y también como la de 4, 4, y 17.

A partir de la sucesión se plantean las siguientes cuestiones:

¿Qué lugar ocupa el 2015 en la sucesión?
¿Cuál será el próximo año que se podrá escribir como la raíz cuadrada de suma de tres cubos?
¿Cuáles son las descomposiciones de 2015 como raíz cuadrada de suma de tres cubos?
¿Cuáles son los años hasta el 2015 que se pueden escribir como raíz cuadrada de suma de tres cubos de más formas distintas?

Soluciones

import Data.List (sort)

raicesCuadradasDeSumasDe3Cubos :: [Integer]
raicesCuadradasDeSumasDe3Cubos =
    [n | n <- [1..], not (null (descomposiciones n))]

-- (descomposiciones n) es la lista de ternas de números tales que n es la
-- raíz cuadrada de la suma de los cubos de los tres números de la
-- terna. Por ejemplo,
--    descomposiciones  6  ==  [(1,2,3)]
--    descomposiciones  9  ==  [(3,3,3)]
--    descomposiciones 71  ==  [(6,9,16),(4,4,17)]
descomposiciones :: Integer -> [(Integer,Integer,Integer)]
descomposiciones n = 
    [(a,b,c) | c <- [1..floor ((fromIntegral (n*n))**(1/3))],
               b <- [1..c],
               let d = n^2 - c^3 - b^3,
               d > 0,
               let a = round ((fromIntegral d)**(1/3)),
               a <= b,
               a^3 == d]    

-- El cálculo de la posición de 2015 es
--    ghci> length (takeWhile (<=2015) raicesCuadradasDeSumasDe3Cubos)
--    343

-- El cálculo del próximo año expresable como la raízcuadrada de la suma
-- de tres cubos es
--    ghci> head (dropWhile (<=2015) raicesCuadradasDeSumasDe3Cubos)
--    2022

-- (masDescomponibles xs) es la lista de elementos de xs que se pueden
-- escribir de más formas como raíz cuadrada de suma de tres cubos.
masDescomponibles :: [Integer] -> [Integer]
masDescomponibles xs =
    [y | (x,y) <- takeWhile (\p -> fst p == u) zs]
    where zs = reverse (sort [(length (descomposiciones x),x) | x <- xs])
          u  = fst (head zs)

-- El cálculo de los años hasta el 2015 con mayor número de
-- descomposiciones es
--    ghci> masDescomponibles (takeWhile (<=2015) raicesCuadradasDeSumasDe3Cubos)
--    [1728]

import Data.List (sort)

raicesCuadradasDeSumasDe3Cubos :: [Integer]

raicesCuadradasDeSumasDe3Cubos =

[n | n <- [1..], not (null (descomposiciones n))]

-- (descomposiciones n) es la lista de ternas de números tales que n es la

-- raíz cuadrada de la suma de los cubos de los tres números de la

-- terna. Por ejemplo,

-- descomposiciones 6 == [(1,2,3)]

-- descomposiciones 9 == [(3,3,3)]

-- descomposiciones 71 == [(6,9,16),(4,4,17)]

descomposiciones :: Integer -> [(Integer,Integer,Integer)]

descomposiciones n =

[(a,b,c) | c <- [1..floor ((fromIntegral (n*n))**(1/3))],

b <- [1..c],

let d = n^2 - c^3 - b^3,

d > 0,

let a = round ((fromIntegral d)**(1/3)),

a <= b,

a^3 == d]

-- El cálculo de la posición de 2015 es

-- ghci> length (takeWhile (<=2015) raicesCuadradasDeSumasDe3Cubos)

-- 343

-- El cálculo del próximo año expresable como la raízcuadrada de la suma

-- de tres cubos es

-- ghci> head (dropWhile (<=2015) raicesCuadradasDeSumasDe3Cubos)

-- 2022

-- (masDescomponibles xs) es la lista de elementos de xs que se pueden

-- escribir de más formas como raíz cuadrada de suma de tres cubos.

masDescomponibles :: [Integer] -> [Integer]

masDescomponibles xs =

[y | (x,y) <- takeWhile (\p -> fst p == u) zs]

where zs = reverse (sort [(length (descomposiciones x),x) | x <- xs])

u = fst (head zs)

-- El cálculo de los años hasta el 2015 con mayor número de

-- descomposiciones es

-- ghci> masDescomponibles (takeWhile (<=2015) raicesCuadradasDeSumasDe3Cubos)

-- [1728]

Caminos en una retícula

Soluciones

Densidad de números no monótonos

Soluciones

Operaciones con series de potencias

Soluciones

Raíz entera

Soluciones

Soluciones en Maxima

Primo suma de dos cuadrados

Soluciones

Referencias

Sumas de potencias de 3 primos

Soluciones

Números primos de Hilbert

Soluciones

Suma con redondeos

Soluciones

Año cúbico

Soluciones

Raíces enteras de los números primos

Soluciones

Productos de N números consecutivos

Soluciones

Números libres de cuadrados

Soluciones

Números de suma prima hereditarios por la derecha

Soluciones

Listas con los ceros emparejados

Soluciones

2015 como raíz cuadrada de la suma de tres cubos

Soluciones

Entradas recientes

Buscar

Correo electrónico