PFH: La semana en Exercitium (del 23 al 27 de mayo de 2022)

Esta semana he publicado en Exercitium las soluciones de los siguientes problemas:

1. Densidades de números abundantes, perfectos y deficientes
2. Matriz zigzagueante
3. Numeración con base múltiple
4. El triángulo de Floyd
5. Polinomios cuadráticos generadores de primos

A continuación se muestran las soluciones.

1. Densidades de números abundantes, perfectos y deficientes

La n-ésima densidad de un tipo de número es el cociente entre la cantidad de los números entre 1 y n que son del tipo considerado y n. Por ejemplo, la 7-ésima densidad de los múltiplos de 3 es 2/7 ya que entre los 7 primeros números sólo 2 son múltiplos de 3.

Definir las funciones

   densidades :: Int -> (Double,Double,Double)
   graficas   :: Int -> IO ()

1 2	densidades :: Int -> (Double,Double,Double) graficas :: Int -> IO ()

tales que

(densidades n) es la terna formada por la n-ésima densidad
- de los números abundantes (es decir, para los que la suma de sus divisores propios es mayor que el número),
- de los números perfectos (es decir, para los que la suma de sus divisores propios es mayor que el número) y
- de los números deficientes (es decir, para los que la suma de sus divisores propios es menor que el número).
Por ejemplo,

    densidades 100     ==  (0.22,    2.0e-2, 0.76)
    densidades 1000    ==  (0.246,   3.0e-3, 0.751)
    densidades 10000   ==  (0.2488,  4.0e-4, 0.7508)
    densidades 100000  ==  (0.24795, 4.0e-5, 0.75201)
    densidades 1000000 ==  (0.247545,4.0e-6, 0.752451)

densidades 100 == (0.22, 2.0e-2, 0.76)

densidades 1000 == (0.246, 3.0e-3, 0.751)

densidades 10000 == (0.2488, 4.0e-4, 0.7508)

densidades 100000 == (0.24795, 4.0e-5, 0.75201)

densidades 1000000 == (0.247545,4.0e-6, 0.752451)

(graficas n) dibuja las gráficas de las k-ésimas densidades (para k entre 1 y n) de los números abundantes, de los números perfectos y de los números deficientes. Por ejemplo, (graficas 100) dibuja

y (graficas 400) dibuja

Soluciones

import Data.List (genericLength, group, partition)
import Data.Array (accumArray, assocs)
import Data.Numbers.Primes (primeFactors)
import Graphics.Gnuplot.Simple (plotLists, Attribute (Key))
import Test.QuickCheck (Positive (Positive), quickCheck)

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

densidades1 :: Int -> (Double,Double,Double)
densidades1 n = (f a, f p, f d)
  where a   = nAbundantes n
        p   = nPerfectos n
        d   = n - a - p
        f x = fromIntegral x / fromIntegral n

-- (nAbundantes n) es la cantidad de números abundantes desde 1 hasta
-- n. Por ejemplo,
--    nAbundantes 100 == 22
nAbundantes :: Int -> Int
nAbundantes n = length (filter esAbundante [1..n])

-- (esAbundante n) se verifica si n es un número abundante. Por ejemplo,
--    esAbundante 12 == True
--    esAbundante 22 == False
esAbundante :: Int -> Bool
esAbundante n = sumaDivisores n > n

-- (sumaDivisores n) es la suma de los divisores propios de n. Por
-- ejemplo,
--    sumaDivisores 12 == 16
--    sumaDivisores 22 == 14
sumaDivisores :: Int -> Int
sumaDivisores = sum . divisores

-- (divisores n) es la lista de los divisores propios de n. Por ejemplo,
--    divisores 12== [1,2,3,4,6]
--    divisores 22== [1,2,11]
divisores :: Int -> [Int]
divisores n = [x | x <- [1..n-1],
                   n `mod` x == 0]

-- (nPerfectos n) es la cantidad de números perfectos desde 1 hasta
-- n. Por ejemplo,
--    nPerfectos 100 == 2
nPerfectos :: Int -> Int
nPerfectos n = length (filter esPerfecto [1..n])

-- (esPerfecto n) se verifica si n es un número perfecto. Por ejemplo,
--    esPerfecto 28 == True
--    esPerfecto 38 == False
esPerfecto :: Int -> Bool
esPerfecto n = sumaDivisores n == n

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

densidades2 :: Int -> (Double,Double,Double)
densidades2 n = (f as, f ps, f ds)
  where (as,pds) = partition esAbundante [1..n]
        (ps,ds)  = partition esPerfecto pds
        f xs     = genericLength xs / fromIntegral n

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

densidades3 :: Int -> (Double,Double,Double)
densidades3 n = (f as, f ps, f ds)
  where cs       = map clasificacion [1..n]
        (as,pds) = partition (== Abundante) cs
        (ps,ds)  = partition (== Perfecto) pds
        f xs     = genericLength xs / fromIntegral n

data Clase = Abundante | Perfecto | Deficiente
  deriving (Eq, Show)

-- (clasificacion n) es la clase de número de n. Por ejemplo,
--    clasificacion 12 == Abundante
--    clasificacion 22 == Deficiente
--    clasificacion 28 == Perfecto
clasificacion :: Int -> Clase
clasificacion n
  | sd > n    = Abundante
  | sd < n    = Deficiente
  | otherwise = Perfecto
  where sd = sumaDivisores n

-- 4ª solución
-- ===========

densidades4 :: Int -> (Double,Double,Double)
densidades4 n = (f as, f ps, f ds)
  where cs       = map clasificacion2 [1..n]
        (as,pds) = partition (== Abundante) cs
        (ps,ds)  = partition (== Perfecto) pds
        f xs     = genericLength xs / fromIntegral n

-- 2ª definición de clasificacion
clasificacion2 :: Int -> Clase
clasificacion2 n
  | sd > n    = Abundante
  | sd < n    = Deficiente
  | otherwise = Perfecto
  where sd = sumaDivisores2 n

-- 2ª definición de sumaDivisores
sumaDivisores2 :: Int -> Int
sumaDivisores2 x =
  product [(p^(e+1)-1) `div` (p-1) | (p,e) <- factorizacion x] - x

-- (factorizacion x) es la lista de las bases y exponentes de la
-- descomposición prima de x. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Int -> [(Int,Int)]
factorizacion = map primeroYlongitud . group . primeFactors

-- (primeroYlongitud xs) es el par formado por el primer elemento de xs
-- y la longitud de xs. Por ejemplo,
--    primeroYlongitud [3,2,5,7] == (3,4)
primeroYlongitud :: [a] -> (a,Int)
primeroYlongitud (x:xs) =
  (x, 1 + length xs)

-- 5ª solución
-- ===========

densidades5 :: Int -> (Double,Double,Double)
densidades5 n = (f a, f p, f d)
  where (a,p,d) = distribucion n
        f x = fromIntegral x / fromIntegral n

-- (distribucion n) es la terna (a,p,d) donde a es la cantidad de
-- números abundantes de 1 a n, p la de los perfectos y d la de los
-- deficientes. Por ejemplo,
--    distribucion 100  ==  (22,2,76)
distribucion :: Int -> (Int,Int,Int)
distribucion n = aux (0,0,0) (sumaDivisoresHasta n)
  where aux (a,p,d) [] = (a,p,d)
        aux (a,p,d) ((x,y):xys)
          | x < y     = aux (1+a,p,d) xys
          | x > y     = aux (a,p,1+d) xys
          | otherwise = aux (a,1+p,d) xys

-- (sumaDivisoresHasta n) es la lista de los pares (a,b) tales que a
-- varía entre 1 y n y b es la suma de los divisores propios de a. Por
-- ejemplo,
--    λ> sumaDivisoresHasta 12
--    [(1,0),(2,1),(3,1),(4,3),(5,1),(6,6),(7,1),(8,7),(9,4),(10,8),(11,1),(12,16)]
sumaDivisoresHasta :: Int -> [(Int,Int)]
sumaDivisoresHasta n =
  assocs (accumArray (+) 0 (1,n) (divisoresHasta n))

-- (divisoresHasta n) es la lista de los pares (a,b) tales que a está
-- entre 2 y n y b es un divisor propio e x. Por ejemplo,
--    λ> divisoresHasta 6
--    [(2,1),(3,1),(4,1),(5,1),(6,1),(4,2),(6,2),(6,3)]
--    λ> divisoresHasta 8
--    [(2,1),(3,1),(4,1),(5,1),(6,1),(7,1),(8,1),(4,2),(6,2),(8,2),(6,3),(8,4)]
divisoresHasta :: Int -> [(Int,Int)]
divisoresHasta n = [(a,b) | b <- [1..n `div` 2], a <- [b*2, b*3..n]]

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_densidades :: Positive Int -> Bool
prop_densidades (Positive n) =
  all (== densidades1 n)
      [ densidades2 n
      , densidades3 n
      , densidades4 n
      , densidades5 n
      ]

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

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

-- La comparación es
--    λ> densidades1 2000
--    (0.2465,1.5e-3,0.752)
--    (1.59 secs, 804,883,768 bytes)
--    λ> densidades2 2000
--    (0.2465,1.5e-3,0.752)
--    (1.39 secs, 704,749,552 bytes)
--    λ> densidades3 2000
--    (0.2465,1.5e-3,0.752)
--    (0.81 secs, 403,783,312 bytes)
--    λ> densidades4 2000
--    (0.2465,1.5e-3,0.752)
--    (0.04 secs, 34,364,960 bytes)
--    λ> densidades5 2000
--    (0.2465,1.5e-3,0.752)
--    (0.02 secs, 4,716,720 bytes)
--
--    λ> densidades4 100000
--    (0.24795,4.0e-5,0.75201)
--    (1.61 secs, 4,826,971,624 bytes)
--    λ> densidades5 100000
--    (0.24795,4.0e-5,0.75201)
--    (0.43 secs, 291,989,272 bytes)

-- Gráfica
-- =======

graficas :: Int -> IO ()
graficas n =
  plotLists [Key Nothing]
            [ [x | (x,_,_) <- ts]
            , [y | (_,y,_) <- ts]
            , [z | (_,_,z) <- ts]]
  where ts = [densidades5 k | k <- [1..n]]

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

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

import Data.List (genericLength, group, partition)

import Data.Array (accumArray, assocs)

import Data.Numbers.Primes (primeFactors)

import Graphics.Gnuplot.Simple (plotLists, Attribute (Key))

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

densidades1 :: Int -> (Double,Double,Double)

densidades1 n = (f a, f p, f d)

where a = nAbundantes n

p = nPerfectos n

d = n - a - p

f x = fromIntegral x / fromIntegral n

-- (nAbundantes n) es la cantidad de números abundantes desde 1 hasta

-- n. Por ejemplo,

-- nAbundantes 100 == 22

nAbundantes :: Int -> Int

nAbundantes n = length (filter esAbundante [1..n])

-- (esAbundante n) se verifica si n es un número abundante. Por ejemplo,

-- esAbundante 12 == True

-- esAbundante 22 == False

esAbundante :: Int -> Bool

esAbundante n = sumaDivisores n > n

-- (sumaDivisores n) es la suma de los divisores propios de n. Por

-- ejemplo,

-- sumaDivisores 12 == 16

-- sumaDivisores 22 == 14

sumaDivisores :: Int -> Int

sumaDivisores = sum . divisores

-- (divisores n) es la lista de los divisores propios de n. Por ejemplo,

-- divisores 12== [1,2,3,4,6]

-- divisores 22== [1,2,11]

divisores :: Int -> [Int]

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

n `mod` x == 0]

-- (nPerfectos n) es la cantidad de números perfectos desde 1 hasta

-- n. Por ejemplo,

-- nPerfectos 100 == 2

nPerfectos :: Int -> Int

nPerfectos n = length (filter esPerfecto [1..n])

-- (esPerfecto n) se verifica si n es un número perfecto. Por ejemplo,

-- esPerfecto 28 == True

-- esPerfecto 38 == False

esPerfecto :: Int -> Bool

esPerfecto n = sumaDivisores n == n

-- 2ª solución

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

densidades2 :: Int -> (Double,Double,Double)

densidades2 n = (f as, f ps, f ds)

where (as,pds) = partition esAbundante [1..n]

(ps,ds) = partition esPerfecto pds

f xs = genericLength xs / fromIntegral n

-- 3ª solución

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

densidades3 :: Int -> (Double,Double,Double)

densidades3 n = (f as, f ps, f ds)

where cs = map clasificacion [1..n]

(as,pds) = partition (== Abundante) cs

(ps,ds) = partition (== Perfecto) pds

f xs = genericLength xs / fromIntegral n

data Clase = Abundante | Perfecto | Deficiente

deriving (Eq, Show)

-- (clasificacion n) es la clase de número de n. Por ejemplo,

-- clasificacion 12 == Abundante

-- clasificacion 22 == Deficiente

-- clasificacion 28 == Perfecto

clasificacion :: Int -> Clase

clasificacion n

| sd > n = Abundante

| sd < n = Deficiente

| otherwise = Perfecto

where sd = sumaDivisores n

-- 4ª solución

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

densidades4 :: Int -> (Double,Double,Double)

densidades4 n = (f as, f ps, f ds)

where cs = map clasificacion2 [1..n]

(as,pds) = partition (== Abundante) cs

(ps,ds) = partition (== Perfecto) pds

f xs = genericLength xs / fromIntegral n

-- 2ª definición de clasificacion

clasificacion2 :: Int -> Clase

clasificacion2 n

| sd > n = Abundante

| sd < n = Deficiente

| otherwise = Perfecto

where sd = sumaDivisores2 n

-- 2ª definición de sumaDivisores

sumaDivisores2 :: Int -> Int

sumaDivisores2 x =

product [(p^(e+1)-1) `div` (p-1) | (p,e) <- factorizacion x] - x

-- (factorizacion x) es la lista de las bases y exponentes de la

-- descomposición prima de x. Por ejemplo,

-- factorizacion 600 == [(2,3),(3,1),(5,2)]

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

factorizacion = map primeroYlongitud . group . primeFactors

-- (primeroYlongitud xs) es el par formado por el primer elemento de xs

-- y la longitud de xs. Por ejemplo,

-- primeroYlongitud [3,2,5,7] == (3,4)

primeroYlongitud :: [a] -> (a,Int)

primeroYlongitud (x:xs) =

(x, 1 + length xs)

-- 5ª solución

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

densidades5 :: Int -> (Double,Double,Double)

densidades5 n = (f a, f p, f d)

where (a,p,d) = distribucion n

f x = fromIntegral x / fromIntegral n

-- (distribucion n) es la terna (a,p,d) donde a es la cantidad de

-- números abundantes de 1 a n, p la de los perfectos y d la de los

-- deficientes. Por ejemplo,

-- distribucion 100 == (22,2,76)

distribucion :: Int -> (Int,Int,Int)

distribucion n = aux (0,0,0) (sumaDivisoresHasta n)

where aux (a,p,d) [] = (a,p,d)

aux (a,p,d) ((x,y):xys)

| x < y = aux (1+a,p,d) xys

| x > y = aux (a,p,1+d) xys

| otherwise = aux (a,1+p,d) xys

-- (sumaDivisoresHasta n) es la lista de los pares (a,b) tales que a

-- varía entre 1 y n y b es la suma de los divisores propios de a. Por

-- ejemplo,

-- λ> sumaDivisoresHasta 12

-- [(1,0),(2,1),(3,1),(4,3),(5,1),(6,6),(7,1),(8,7),(9,4),(10,8),(11,1),(12,16)]

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

sumaDivisoresHasta n =

assocs (accumArray (+) 0 (1,n) (divisoresHasta n))

-- (divisoresHasta n) es la lista de los pares (a,b) tales que a está

-- entre 2 y n y b es un divisor propio e x. Por ejemplo,

-- λ> divisoresHasta 6

-- [(2,1),(3,1),(4,1),(5,1),(6,1),(4,2),(6,2),(6,3)]

-- λ> divisoresHasta 8

-- [(2,1),(3,1),(4,1),(5,1),(6,1),(7,1),(8,1),(4,2),(6,2),(8,2),(6,3),(8,4)]

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

divisoresHasta n = [(a,b) | b <- [1..n `div` 2], a <- [b*2, b*3..n]]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_densidades :: Positive Int -> Bool

prop_densidades (Positive n) =

all (== densidades1 n)

[ densidades2 n

, densidades3 n

, densidades4 n

, densidades5 n

]

-- La comprobación es

-- λ> quickCheck prop_densidades

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> densidades1 2000

-- (0.2465,1.5e-3,0.752)

-- (1.59 secs, 804,883,768 bytes)

-- λ> densidades2 2000

-- (0.2465,1.5e-3,0.752)

-- (1.39 secs, 704,749,552 bytes)

-- λ> densidades3 2000

-- (0.2465,1.5e-3,0.752)

-- (0.81 secs, 403,783,312 bytes)

-- λ> densidades4 2000

-- (0.2465,1.5e-3,0.752)

-- (0.04 secs, 34,364,960 bytes)

-- λ> densidades5 2000

-- (0.2465,1.5e-3,0.752)

-- (0.02 secs, 4,716,720 bytes)

-- λ> densidades4 100000

-- (0.24795,4.0e-5,0.75201)

-- (1.61 secs, 4,826,971,624 bytes)

-- λ> densidades5 100000

-- (0.24795,4.0e-5,0.75201)

-- (0.43 secs, 291,989,272 bytes)

-- Gráfica

-- =======

graficas :: Int -> IO ()

graficas n =

plotLists [Key Nothing]

[ [x | (x,_,_) <- ts]

, [y | (_,y,_) <- ts]

, [z | (_,_,z) <- ts]]

where ts = [densidades5 k | k <- [1..n]]

El código se encuentra en GitHub.

2. Matriz zigzagueante

La matriz zizagueante de orden n es la matriz cuadrada con n filas y n columnas y cuyos elementos son los n² primeros números naturales colocados de manera creciente a lo largo de las diagonales secundarias. Por ejemplo, La matriz zigzagueante de orden 5 es

    0  1  5  6 14
    2  4  7 13 15
    3  8 12 16 21
    9 11 17 20 22
   10 18 19 23 24

0 1 5 6 14

2 4 7 13 15

3 8 12 16 21

9 11 17 20 22

10 18 19 23 24

La colocación de los elementos se puede ver gráficamente en esta figura

Definir la función

   zigZag :: Int -> Matrix Int

1	zigZag :: Int -> Matrix Int

tal que (zigZag n) es la matriz zigzagueante de orden n. Por ejemplo,

   λ> zigZag 5
   ┌                ┐
   │  0  1  5  6 14 │
   │  2  4  7 13 15 │
   │  3  8 12 16 21 │
   │  9 11 17 20 22 │
   │ 10 18 19 23 24 │
   └                ┘
   λ> zigZag 4
   ┌             ┐
   │  0  1  5  6 │
   │  2  4  7 12 │
   │  3  8 11 13 │
   │  9 10 14 15 │
   └             ┘
   λ> maximum (zigZag 1500)
   2249999

λ> zigZag 5

┌ ┐

│ 0 1 5 6 14 │

│ 2 4 7 13 15 │

│ 3 8 12 16 21 │

│ 9 11 17 20 22 │

│ 10 18 19 23 24 │

└ ┘

λ> zigZag 4

┌ ┐

│ 0 1 5 6 │

│ 2 4 7 12 │

│ 3 8 11 13 │

│ 9 10 14 15 │

└ ┘

λ> maximum (zigZag 1500)

2249999

Soluciones

import Data.List (sort, sortBy)
import Data.Matrix (Matrix, fromList)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

zigZag1 :: Int -> Matrix Int
zigZag1 n = fromList n n (elementosZigZag n)

-- (elementosZigZag n) es la lista de los elementos de la matriz
-- zizagueante de orden n. Por ejemplo.
--    λ> elementosZigZag 5
--    [0,1,5,6,14,2,4,7,13,15,3,8,12,16,21,9,11,17,20,22,10,18,19,23,24]
elementosZigZag :: Int -> [Int]
elementosZigZag n =
  map snd (sort (zip (ordenZigZag n) [0..]))

-- (ordenZigZag n) es la lista de puntos del cuadrado nxn recorridos en
-- zig-zag por las diagonales secundarias. Por ejemplo,
--    λ> ordenZigZag 4
--    [(1,1), (1,2),(2,1), (3,1),(2,2),(1,3), (1,4),(2,3),(3,2),(4,1),
--     (4,2),(3,3),(2,4), (3,4),(4,3), (4,4)]
ordenZigZag :: Int -> [(Int,Int)]
ordenZigZag n = concat [aux n m | m <- [2..2*n]]
    where aux k m | odd m     = [(x,m-x) | x <- [max 1 (m-k)..min k (m-1)]]
                  | otherwise = [(m-x,x) | x <- [max 1 (m-k)..min k (m-1)]]

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

zigZag2 :: Int -> Matrix Int
zigZag2 n = fromList n n (elementosZigZag2 n)

elementosZigZag2 :: Int -> [Int]
elementosZigZag2 n =
  map snd (sort (zip (ordenZigZag2 n) [0..]))

ordenZigZag2 :: Int -> [(Int,Int)]
ordenZigZag2 n = sortBy comp [(x,y) | x <- [1..n], y <- [1..n]]
    where comp (x1,y1) (x2,y2) | x1+y1 < x2+y2 = LT
                               | x1+y1 > x2+y2 = GT
                               | even (x1+y1)  = compare y1 y2
                               | otherwise     = compare x1 x2

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_zigZag :: Positive Int -> Bool
prop_zigZag (Positive n) =
  zigZag1 n == zigZag2 n

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

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

-- La comparación es
--    λ> length (zigZag1 5000)
--    25000000
--    (2.57 secs, 1,800,683,952 bytes)
--    λ> length (zigZag2 5000)
--    25000000
--    (2.20 secs, 1,800,683,952 bytes)
--
--    λ> maximum (zigZag1 1100)
--    1209999
--    (2.12 secs, 1,840,095,864 bytes)
--    λ> maximum (zigZag2 1100)
--    1209999
--    (21.27 secs, 11,661,088,256 bytes)

import Data.List (sort, sortBy)

import Data.Matrix (Matrix, fromList)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

zigZag1 :: Int -> Matrix Int

zigZag1 n = fromList n n (elementosZigZag n)

-- (elementosZigZag n) es la lista de los elementos de la matriz

-- zizagueante de orden n. Por ejemplo.

-- λ> elementosZigZag 5

-- [0,1,5,6,14,2,4,7,13,15,3,8,12,16,21,9,11,17,20,22,10,18,19,23,24]

elementosZigZag :: Int -> [Int]

elementosZigZag n =

map snd (sort (zip (ordenZigZag n) [0..]))

-- (ordenZigZag n) es la lista de puntos del cuadrado nxn recorridos en

-- zig-zag por las diagonales secundarias. Por ejemplo,

-- λ> ordenZigZag 4

-- [(1,1), (1,2),(2,1), (3,1),(2,2),(1,3), (1,4),(2,3),(3,2),(4,1),

-- (4,2),(3,3),(2,4), (3,4),(4,3), (4,4)]

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

ordenZigZag n = concat [aux n m | m <- [2..2*n]]

where aux k m | odd m = [(x,m-x) | x <- [max 1 (m-k)..min k (m-1)]]

| otherwise = [(m-x,x) | x <- [max 1 (m-k)..min k (m-1)]]

-- 2ª solución

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

zigZag2 :: Int -> Matrix Int

zigZag2 n = fromList n n (elementosZigZag2 n)

elementosZigZag2 :: Int -> [Int]

elementosZigZag2 n =

map snd (sort (zip (ordenZigZag2 n) [0..]))

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

ordenZigZag2 n = sortBy comp [(x,y) | x <- [1..n], y <- [1..n]]

where comp (x1,y1) (x2,y2) | x1+y1 < x2+y2 = LT

| x1+y1 > x2+y2 = GT

| even (x1+y1) = compare y1 y2

| otherwise = compare x1 x2

-- Comprobación de equivalencia

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

-- La propiedad es

prop_zigZag :: Positive Int -> Bool

prop_zigZag (Positive n) =

zigZag1 n == zigZag2 n

-- La comprobación es

-- λ> quickCheck prop_zigZag

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (zigZag1 5000)

-- 25000000

-- (2.57 secs, 1,800,683,952 bytes)

-- λ> length (zigZag2 5000)

-- 25000000

-- (2.20 secs, 1,800,683,952 bytes)

-- λ> maximum (zigZag1 1100)

-- 1209999

-- (2.12 secs, 1,840,095,864 bytes)

-- λ> maximum (zigZag2 1100)

-- 1209999

-- (21.27 secs, 11,661,088,256 bytes)

El código se encuentra en GitHub.

3. Numeración con base múltiple

Sea (b(i) | i ≥ 1) una sucesión infinita de números enteros mayores que 1. Entonces todo entero x mayor que cero se puede escribir de forma única como

   x = x(0) + x(1)b(1) +x(2)b(1)b(2) + ... + x(n)b(1)b(2)...b(n)

1	x = x(0) + x(1)b(1) +x(2)b(1)b(2) + ... + x(n)b(1)b(2)...b(n)

donde cada x(i) satisface la condición 0 ≤ x(i) < b(i+1). Se dice que [x(n),x(n-1),…,x(2),x(1),x(0)] es la representación de x en la base (b(i)). Por ejemplo, la representación de 377 en la base (2, 6, 8, …) es [7,5,0,1] ya que

   377 = 1 + 0*2 + 5*2*4 + 7*2*4*6

1	377 = 1 + 02 + 524 + 7246

y, además, 0 ≤ 1 < 2, 0 ≤ 0 < 4, 0 ≤ 5 < 6 y 0 ≤ 7 < 8.

Definir las funciones

   decimalAmultiple :: [Integer] -> Integer -> [Integer]
   multipleAdecimal :: [Integer] -> [Integer] -> Integer

1 2	decimalAmultiple :: [Integer] -> Integer -> [Integer] multipleAdecimal :: [Integer] -> [Integer] -> Integer

tales que

(decimalAmultiple bs x) es la representación del número x en la base bs. Por ejemplo,

     decimalAmultiple [2,4..] 377                      ==  [7,5,0,1]
     decimalAmultiple [2,5..] 377                      ==  [4,5,3,1]
     decimalAmultiple [2^n | n <- [1..]] 2015          ==  [1,15,3,3,1]
     decimalAmultiple (repeat 10) 2015                 ==  [2,0,1,5]

decimalAmultiple [2,4..] 377 == [7,5,0,1]

decimalAmultiple [2,5..] 377 == [4,5,3,1]

decimalAmultiple [2^n | n <- [1..]] 2015 == [1,15,3,3,1]

decimalAmultiple (repeat 10) 2015 == [2,0,1,5]

(multipleAdecimal bs cs) es el número decimal cuya representación en la base bs es cs. Por ejemplo,

     multipleAdecimal [2,4..] [7,5,0,1]                ==  377
     multipleAdecimal [2,5..] [4,5,3,1]                ==  377
     multipleAdecimal [2^n | n <- [1..]] [1,15,3,3,1]  ==  2015
     multipleAdecimal (repeat 10) [2,0,1,5]            ==  2015

multipleAdecimal [2,4..] [7,5,0,1] == 377

multipleAdecimal [2,5..] [4,5,3,1] == 377

multipleAdecimal [2^n | n <- [1..]] [1,15,3,3,1] == 2015

multipleAdecimal (repeat 10) [2,0,1,5] == 2015

Comprobar con QuickCheck que se verifican las siguientes propiedades

Para cualquier base bs y cualquier entero positivo n,

     multipleAdecimal bs (decimalAmultiple bs x) == x

1	multipleAdecimal bs (decimalAmultiple bs x) == x

Para cualquier base bs y cualquier entero positivo n, el coefiente i-ésimo de la representación múltiple de n en la base bs es un entero no negativo menos que el i-ésimo elemento de bs.

Soluciones

import Test.QuickCheck
import Data.List (unfoldr)

-- 1ª solución de decimalAmultiple
-- ===============================

decimalAmultiple1 :: [Integer] -> Integer -> [Integer]
decimalAmultiple1 bs n = reverse (aux bs n)
  where aux _ 0      = []
        aux (d:ds) m = r : aux ds q
          where (q,r) = quotRem m d

-- 2ª solución de decimalAmultiple
-- ===============================

decimalAmultiple2 :: [Integer] -> Integer -> [Integer]
decimalAmultiple2 bs n = aux bs n []
  where aux _ 0  xs     = xs
        aux (d:ds) m xs = aux ds q (r:xs)
          where (q,r) = quotRem m d

-- 3ª solución de decimalAmultiple
-- ===============================

decimalAmultiple3 :: [Integer] -> Integer -> [Integer]
decimalAmultiple3 xs n = reverse (unfoldr f (xs,n))
  where f (_     ,0) = Nothing
        f ((y:ys),m) = Just (r,(ys,q))
                       where (q,r) = quotRem m y

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_decimalAmultiple :: InfiniteList (Positive Integer) -> Positive Integer -> Bool
prop_decimalAmultiple (InfiniteList xs _) (Positive n) =
  all (== decimalAmultiple1 xs' n)
      [decimalAmultiple2 xs' n,
       decimalAmultiple3 xs' n]
  where xs' = map getPositive xs

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

-- La comparación es
--    λ> length (decimalAmultiple1 [2,7..] (10^(10^5)))
--    21731
--    (0.45 secs, 486,085,256 bytes)
--    λ> length (decimalAmultiple2 [2,7..] (10^(10^5)))
--    21731
--    (0.32 secs, 485,563,664 bytes)
--    λ> length (decimalAmultiple3 [2,7..] (10^(10^5)))
--    21731
--    (0.44 secs, 487,649,768 bytes)

-- 1ª solución de multipleAdecimal
-- ===============================

multipleAdecimal1  :: [Integer] -> [Integer] -> Integer
multipleAdecimal1 xs ns = aux xs (reverse ns)
  where aux (y:ys) (m:ms) = m + y * (aux ys ms)
        aux _ _           = 0

-- 2ª solución de multipleAdecimal
-- ===============================

multipleAdecimal2 :: [Integer] -> [Integer] -> Integer
multipleAdecimal2 bs xs =
  sum (zipWith (*) (reverse xs) (1 : scanl1 (*) bs))

-- Comprobación de equivalencia de multipleAdecimal
-- ================================================

-- La propiedad es
prop_multipleAdecimal :: InfiniteList (Positive Integer) -> [Positive Integer] -> Bool
prop_multipleAdecimal (InfiniteList xs _) ys =
  multipleAdecimal1 xs' ys' == multipleAdecimal2 xs' ys'
  where xs' = map getPositive xs
        ys' = map getPositive ys

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

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

-- La comparación es
--    λ> length (show (multipleAdecimal1 [2,3..] [1..10^4]))
--    35660
--    (0.14 secs, 179,522,152 bytes)
--    λ> length (show (multipleAdecimal2 [2,3..] [1..10^4]))
--    35660
--    (0.22 secs, 243,368,664 bytes)

-- Comprobación de las propiedades
-- ===============================

-- La primera propiedad es
prop_inversas :: InfiniteList (Positive Integer) -> Positive Integer -> Bool
prop_inversas (InfiniteList xs _) (Positive n) =
  multipleAdecimal1 xs' (decimalAmultiple1 xs' n) == n
  where xs' = map getPositive xs

-- Su comprobación es
--    λ> quickCheck prop_inversas
--    +++ OK, passed 100 tests.

-- la 2ª propiedad es
prop_coeficientes :: InfiniteList (Positive Integer) -> Positive Integer -> Bool
prop_coeficientes (InfiniteList xs _) (Positive n) =
  and [0 <= c && c < b | (c,b) <- zip cs xs']
  where xs' = map getPositive xs
        cs = reverse (decimalAmultiple1 xs' n)

-- Su comprobación es
--    λ> quickCheck prop_coeficientes
--    +++ OK, passed 100 tests.

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

import Test.QuickCheck

import Data.List (unfoldr)

-- 1ª solución de decimalAmultiple

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

decimalAmultiple1 :: [Integer] -> Integer -> [Integer]

decimalAmultiple1 bs n = reverse (aux bs n)

where aux _ 0 = []

aux (d:ds) m = r : aux ds q

where (q,r) = quotRem m d

-- 2ª solución de decimalAmultiple

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

decimalAmultiple2 :: [Integer] -> Integer -> [Integer]

decimalAmultiple2 bs n = aux bs n []

where aux _ 0 xs = xs

aux (d:ds) m xs = aux ds q (r:xs)

where (q,r) = quotRem m d

-- 3ª solución de decimalAmultiple

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

decimalAmultiple3 :: [Integer] -> Integer -> [Integer]

decimalAmultiple3 xs n = reverse (unfoldr f (xs,n))

where f (_ ,0) = Nothing

f ((y:ys),m) = Just (r,(ys,q))

where (q,r) = quotRem m y

-- Comprobación de equivalencia

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

-- La propiedad es

prop_decimalAmultiple :: InfiniteList (Positive Integer) -> Positive Integer -> Bool

prop_decimalAmultiple (InfiniteList xs _) (Positive n) =

all (== decimalAmultiple1 xs' n)

[decimalAmultiple2 xs' n,

decimalAmultiple3 xs' n]

where xs' = map getPositive xs

-- Comparación de eficiencia de decimalAmultiple

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

-- La comparación es

-- λ> length (decimalAmultiple1 [2,7..] (10^(10^5)))

-- 21731

-- (0.45 secs, 486,085,256 bytes)

-- λ> length (decimalAmultiple2 [2,7..] (10^(10^5)))

-- 21731

-- (0.32 secs, 485,563,664 bytes)

-- λ> length (decimalAmultiple3 [2,7..] (10^(10^5)))

-- 21731

-- (0.44 secs, 487,649,768 bytes)

-- 1ª solución de multipleAdecimal

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

multipleAdecimal1 :: [Integer] -> [Integer] -> Integer

multipleAdecimal1 xs ns = aux xs (reverse ns)

where aux (y:ys) (m:ms) = m + y * (aux ys ms)

aux _ _ = 0

-- 2ª solución de multipleAdecimal

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

multipleAdecimal2 :: [Integer] -> [Integer] -> Integer

multipleAdecimal2 bs xs =

sum (zipWith (*) (reverse xs) (1 : scanl1 (*) bs))

-- Comprobación de equivalencia de multipleAdecimal

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

-- La propiedad es

prop_multipleAdecimal :: InfiniteList (Positive Integer) -> [Positive Integer] -> Bool

prop_multipleAdecimal (InfiniteList xs _) ys =

multipleAdecimal1 xs' ys' == multipleAdecimal2 xs' ys'

where xs' = map getPositive xs

ys' = map getPositive ys

-- La comprobación es

-- λ> quickCheck prop_multipleAdecimal

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de multipleAdecimal

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

-- La comparación es

-- λ> length (show (multipleAdecimal1 [2,3..] [1..10^4]))

-- 35660

-- (0.14 secs, 179,522,152 bytes)

-- λ> length (show (multipleAdecimal2 [2,3..] [1..10^4]))

-- 35660

-- (0.22 secs, 243,368,664 bytes)

-- Comprobación de las propiedades

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

-- La primera propiedad es

prop_inversas :: InfiniteList (Positive Integer) -> Positive Integer -> Bool

prop_inversas (InfiniteList xs _) (Positive n) =

multipleAdecimal1 xs' (decimalAmultiple1 xs' n) == n

where xs' = map getPositive xs

-- Su comprobación es

-- λ> quickCheck prop_inversas

-- +++ OK, passed 100 tests.

-- la 2ª propiedad es

prop_coeficientes :: InfiniteList (Positive Integer) -> Positive Integer -> Bool

prop_coeficientes (InfiniteList xs _) (Positive n) =

and [0 <= c && c < b | (c,b) <- zip cs xs']

where xs' = map getPositive xs

cs = reverse (decimalAmultiple1 xs' n)

-- Su comprobación es

-- λ> quickCheck prop_coeficientes

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

4. El triángulo de Floyd

El triángulo de Floyd, llamado así en honor a Robert Floyd, es un triángulo rectángulo formado con números naturales. Para crear un triángulo de Floyd, se comienza con un 1 en la esquina superior izquierda, y se continúa escribiendo la secuencia de los números naturales de manera que cada línea contenga un número más que la anterior. Las 5 primeras líneas del triángulo de Floyd son

    1
    2   3
    4   5   6
    7   8   9  10
   11  12  13  14  15

2 3

4 5 6

7 8 9 10

11 12 13 14 15

Definir la función

   trianguloFloyd :: [[Integer]]

1	trianguloFloyd :: [[Integer]]

tal que trianguloFloyd es el triángulo de Floyd. Por ejemplo,

   λ> take 4 trianguloFloyd
   [[1],
    [2,3],
    [4,5,6],
    [7,8,9,10]]
  (trianguloFloyd !! (10^5)) !! 0  ==  5000050001
  (trianguloFloyd !! (10^6)) !! 0  ==  500000500001
  (trianguloFloyd !! (10^7)) !! 0  ==  50000005000001

λ> take 4 trianguloFloyd

[[1],

[2,3],

[4,5,6],

[7,8,9,10]]

(trianguloFloyd !! (10^5)) !! 0 == 5000050001

(trianguloFloyd !! (10^6)) !! 0 == 500000500001

(trianguloFloyd !! (10^7)) !! 0 == 50000005000001

Soluciones

import Data.List (genericLength)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

trianguloFloyd1 :: [[Integer]]
trianguloFloyd1 = floyd 1 [1..]
  where floyd n xs = i : floyd (n+1) r
          where (i,r) = splitAt n xs

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

trianguloFloyd2 :: [[Integer]]
trianguloFloyd2 = iterate siguienteF [1]

-- (siguienteF xs) es la lista de los elementos de la línea xs en el
-- triángulo de Floyd. Por ejemplo,
--    siguienteF [2,3]    ==  [4,5,6]
--    siguienteF [4,5,6]  ==  [7,8,9,10]
siguienteF :: [Integer] -> [Integer]
siguienteF xs = [a..a+n]
    where a = 1 + last xs
          n = genericLength xs

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

trianguloFloyd3 :: [[Integer]]
trianguloFloyd3 =
  [[(n*(n-1) `div` 2) + 1 .. (n*(n+1) `div` 2)] | n <- [1..]]

-- 4ª solución
-- ===========

trianguloFloyd4 :: [[Integer]]
trianguloFloyd4 =
  scanl (\(x:_) y -> [x+y..x+2*y]) [1] [1..]

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_trianguloFloyd :: Positive Int -> Bool
prop_trianguloFloyd (Positive n) =
  all (== (trianguloFloyd1 !! n))
      [trianguloFloyd2 !! n,
       trianguloFloyd3 !! n,
       trianguloFloyd4 !! n]

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

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

-- La comparación es
--    λ> (trianguloFloyd1 !! 5000) !! 5000
--    12507501
--    (1.47 secs, 2,505,005,752 bytes)
--    λ> (trianguloFloyd2 !! 5000) !! 5000
--    12507501
--    (0.79 secs, 2,416,259,176 bytes)
--    λ> (trianguloFloyd3 !! 5000) !! 5000
--    12507501
--    (0.00 secs, 1,809,152 bytes)
--    λ> (trianguloFloyd4 !! 5000) !! 5000
--    12507501
--    (0.01 secs, 3,517,896 bytes)
--
--    λ> (trianguloFloyd3 !! (10^7)) !! 0
--    50000005000001
--    (2.45 secs, 1,656,534,080 bytes)
--    λ> (trianguloFloyd4 !! (10^7)) !! 0
--    50000005000001
--    (10.86 secs, 5,302,760,752 bytes)

import Data.List (genericLength)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

trianguloFloyd1 :: [[Integer]]

trianguloFloyd1 = floyd 1 [1..]

where floyd n xs = i : floyd (n+1) r

where (i,r) = splitAt n xs

-- 2ª solución

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

trianguloFloyd2 :: [[Integer]]

trianguloFloyd2 = iterate siguienteF [1]

-- (siguienteF xs) es la lista de los elementos de la línea xs en el

-- triángulo de Floyd. Por ejemplo,

-- siguienteF [2,3] == [4,5,6]

-- siguienteF [4,5,6] == [7,8,9,10]

siguienteF :: [Integer] -> [Integer]

siguienteF xs = [a..a+n]

where a = 1 + last xs

n = genericLength xs

-- 3ª solución

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

trianguloFloyd3 :: [[Integer]]

trianguloFloyd3 =

[[(n*(n-1) `div` 2) + 1 .. (n*(n+1) `div` 2)] | n <- [1..]]

-- 4ª solución

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

trianguloFloyd4 :: [[Integer]]

trianguloFloyd4 =

scanl (\(x:_) y -> [x+y..x+2*y]) [1] [1..]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_trianguloFloyd :: Positive Int -> Bool

prop_trianguloFloyd (Positive n) =

all (== (trianguloFloyd1 !! n))

[trianguloFloyd2 !! n,

trianguloFloyd3 !! n,

trianguloFloyd4 !! n]

-- La comprobación es

-- λ> quickCheck prop_trianguloFloyd

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> (trianguloFloyd1 !! 5000) !! 5000

-- 12507501

-- (1.47 secs, 2,505,005,752 bytes)

-- λ> (trianguloFloyd2 !! 5000) !! 5000

-- 12507501

-- (0.79 secs, 2,416,259,176 bytes)

-- λ> (trianguloFloyd3 !! 5000) !! 5000

-- 12507501

-- (0.00 secs, 1,809,152 bytes)

-- λ> (trianguloFloyd4 !! 5000) !! 5000

-- 12507501

-- (0.01 secs, 3,517,896 bytes)

-- λ> (trianguloFloyd3 !! (10^7)) !! 0

-- 50000005000001

-- (2.45 secs, 1,656,534,080 bytes)

-- λ> (trianguloFloyd4 !! (10^7)) !! 0

-- 50000005000001

-- (10.86 secs, 5,302,760,752 bytes)

El código se encuentra en GitHub.

5. Polinomios cuadráticos generadores de primos

En 1772, Euler publicó que el polinomio n² + n + 41 genera 40 números primos para todos los valores de n entre 0 y 39. Sin embargo, cuando n = 40, 40²+40+41 = 40(40+1)+41 es divisible por 41.

Usando ordenadores, se descubrió que el polinomio n² – 79n + 1601 genera 80 números primos para todos los valores de n entre 0 y 79.

Definir la función

   generadoresMaximales :: Integer -> (Int,[(Integer,Integer)])

1	generadoresMaximales :: Integer -> (Int,[(Integer,Integer)])

tal que (generadoresMaximales n) es el par (m,xs) donde

xs es la lista de pares (x,y) tales que n²+xn+y es uno de polinomios que genera un número máximo de números primos consecutivos a partir de cero entre todos los polinomios de la forma n²+an+b, con |a| ≤ n y |b| ≤ n y
m es dicho número máximo.

Por ejemplo,

   generadoresMaximales    4  ==  ( 3,[(-2,3),(-1,3),(3,3)])
   generadoresMaximales    6  ==  ( 5,[(-1,5),(5,5)])
   generadoresMaximales   41  ==  (41,[(-1,41)])
   generadoresMaximales   50  ==  (43,[(-5,47)])
   generadoresMaximales  100  ==  (48,[(-15,97)])
   generadoresMaximales  200  ==  (53,[(-25,197)])
   generadoresMaximales 1650  ==  (80,[(-79,1601)])

generadoresMaximales 4 == ( 3,[(-2,3),(-1,3),(3,3)])

generadoresMaximales 6 == ( 5,[(-1,5),(5,5)])

generadoresMaximales 41 == (41,[(-1,41)])

generadoresMaximales 50 == (43,[(-5,47)])

generadoresMaximales 100 == (48,[(-15,97)])

generadoresMaximales 200 == (53,[(-25,197)])

generadoresMaximales 1650 == (80,[(-79,1601)])

Soluciones

import Data.List (sort)
import Data.Numbers.Primes (primes, isPrime)
import I1M.PolOperaciones (valor, consPol, polCero)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

generadoresMaximales1 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales1 n =
  (m,[(a,b) | a <- [-n..n], b <- [-n..n], nPrimos a b == m])
  where m = maximum [nPrimos a b | a <- [-n..n], b <- [-n..n]]

-- (nPrimos a b) es el número de primos consecutivos generados por el
-- polinomio n² + an + b a partir de n=0. Por ejemplo,
--    nPrimos (-1) 41     ==  41
--    nPrimos (-79) 1601  ==  80
nPrimos :: Integer -> Integer -> Int
nPrimos a b =
  length (takeWhile isPrime [n*n+a*n+b | n <- [0..]])

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

-- Notas:
-- 1. Se tiene que b es primo, ya que para n = 0, se tiene que
--    0²+a*0+b = b es primo.
-- 2. Se tiene que 1+a+b es primo, ya que es el valor del polinomio para
--    n = 1.

generadoresMaximales2 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales2 n = (m,map snd zs)
  where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,
                                    a <- [-n..n],
                                    isPrime(1+a+b)]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys

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

generadoresMaximales3 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales3 n = (m,map snd zs)
  where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,
                                    p <- takeWhile (<=1+2*n) primes,
                                    let a = p-b-1]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys

-- 4ª solución (con la librería de polinomios)
-- ===========================================

generadoresMaximales4 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales4 n = (m,map snd zs)
  where xs = [(nPrimos2 a b,(a,b)) | b <- takeWhile (<=n) primes,
                                     p <- takeWhile (<=1+2*n) primes,
                                     let a = p-b-1]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys

-- (nPrimos2 a b) es el número de primos consecutivos generados por el
-- polinomio n² + an + b a partir de n=0. Por ejemplo,
--    nPrimos2 (-1) 41     ==  41
--    nPrimos2 (-79) 1601  ==  80
nPrimos2 :: Integer -> Integer -> Int
nPrimos2 a b =
  length (takeWhile isPrime [valor p n | n <- [0..]])
  where p = consPol 2 1 (consPol 1 a (consPol 0 b polCero))

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_generadoresMaximales :: Positive Integer -> Bool
prop_generadoresMaximales (Positive n) =
  all (equivalentes (generadoresMaximales1 n'))
      [generadoresMaximales2 n',
       generadoresMaximales3 n',
       generadoresMaximales4 n']
  where n' = n+1

equivalentes :: (Int,[(Integer,Integer)]) -> (Int,[(Integer,Integer)]) -> Bool
equivalentes (n,xs) (m,ys) =
  n == m && sort xs == sort ys

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

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

-- La comparación es
--    λ> generadoresMaximales1 300
--    (56,[(-31,281)])
--    (2.10 secs, 2,744,382,760 bytes)
--    λ> generadoresMaximales2 300
--    (56,[(-31,281)])
--    (0.17 secs, 382,103,656 bytes)
--    λ> generadoresMaximales3 300
--    (56,[(-31,281)])
--    (0.19 secs, 346,725,872 bytes)
--    λ> generadoresMaximales4 300
--    (56,[(-31,281)])
--    (0.20 secs, 388,509,808 bytes)

100

101

102

103

104

105

106

107

108

import Data.List (sort)

import Data.Numbers.Primes (primes, isPrime)

import I1M.PolOperaciones (valor, consPol, polCero)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

generadoresMaximales1 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales1 n =

(m,[(a,b) | a <- [-n..n], b <- [-n..n], nPrimos a b == m])

where m = maximum [nPrimos a b | a <- [-n..n], b <- [-n..n]]

-- (nPrimos a b) es el número de primos consecutivos generados por el

-- polinomio n² + an + b a partir de n=0. Por ejemplo,

-- nPrimos (-1) 41 == 41

-- nPrimos (-79) 1601 == 80

nPrimos :: Integer -> Integer -> Int

nPrimos a b =

length (takeWhile isPrime [n*n+a*n+b | n <- [0..]])

-- 2ª solución

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

-- Notas:

-- 1. Se tiene que b es primo, ya que para n = 0, se tiene que

-- 0²+a*0+b = b es primo.

-- 2. Se tiene que 1+a+b es primo, ya que es el valor del polinomio para

-- n = 1.

generadoresMaximales2 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales2 n = (m,map snd zs)

where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,

a <- [-n..n],

isPrime(1+a+b)]

ys = reverse (sort xs)

m = fst (head ys)

zs = takeWhile (\(k,_) -> k == m) ys

-- 3ª solución

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

generadoresMaximales3 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales3 n = (m,map snd zs)

where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,

p <- takeWhile (<=1+2*n) primes,

let a = p-b-1]

ys = reverse (sort xs)

m = fst (head ys)

zs = takeWhile (\(k,_) -> k == m) ys

-- 4ª solución (con la librería de polinomios)

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

generadoresMaximales4 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales4 n = (m,map snd zs)

where xs = [(nPrimos2 a b,(a,b)) | b <- takeWhile (<=n) primes,

p <- takeWhile (<=1+2*n) primes,

let a = p-b-1]

ys = reverse (sort xs)

m = fst (head ys)

zs = takeWhile (\(k,_) -> k == m) ys

-- (nPrimos2 a b) es el número de primos consecutivos generados por el

-- polinomio n² + an + b a partir de n=0. Por ejemplo,

-- nPrimos2 (-1) 41 == 41

-- nPrimos2 (-79) 1601 == 80

nPrimos2 :: Integer -> Integer -> Int

nPrimos2 a b =

length (takeWhile isPrime [valor p n | n <- [0..]])

where p = consPol 2 1 (consPol 1 a (consPol 0 b polCero))

-- Comprobación de equivalencia

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

-- La propiedad es

prop_generadoresMaximales :: Positive Integer -> Bool

prop_generadoresMaximales (Positive n) =

all (equivalentes (generadoresMaximales1 n'))

[generadoresMaximales2 n',

generadoresMaximales3 n',

generadoresMaximales4 n']

where n' = n+1

equivalentes :: (Int,[(Integer,Integer)]) -> (Int,[(Integer,Integer)]) -> Bool

equivalentes (n,xs) (m,ys) =

n == m && sort xs == sort ys

-- La comprobación es

-- λ> quickCheck prop_generadoresMaximales

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> generadoresMaximales1 300

-- (56,[(-31,281)])

-- (2.10 secs, 2,744,382,760 bytes)

-- λ> generadoresMaximales2 300

-- (56,[(-31,281)])

-- (0.17 secs, 382,103,656 bytes)

-- λ> generadoresMaximales3 300

-- (56,[(-31,281)])

-- (0.19 secs, 346,725,872 bytes)

-- λ> generadoresMaximales4 300

-- (56,[(-31,281)])

-- (0.20 secs, 388,509,808 bytes)

El código se encuentra en GitHub.

1. Densidades de números abundantes, perfectos y deficientes

Soluciones

2. Matriz zigzagueante

Soluciones

3. Numeración con base múltiple

Soluciones

4. El triángulo de Floyd

Soluciones

5. Polinomios cuadráticos generadores de primos

Soluciones

Se puede imprimir o compartir con