group

Las conjeturas de Catalan y de Pillai

PorJosé A. Alonso 15-enero-202022-enero-2020

La conjetura de Catalan, enunciada en 1844 por Eugène Charles Catalan y demostrada 2002 por Preda Mihăilescu1, afirma que

Las únicas dos potencias de números enteros consecutivos son 8 y 9 (que son respectivamente 2³ y 3²).

En otras palabras, la única solución entera de la ecuación

   x^a - y^b = 1

1	x^a - y^b = 1

para x, a, y, b > 1 es x = 3, a = 2, y = 2, b = 3.

La conjetura de Pillai, propuesta por S.S. Pillai en 1942, generaliza este resultado y es un problema abierto. Afirma que cada entero se puede escribir sólo un número finito de veces como una diferencia de dos potencias perfectas. En otras palabras, para todo entero positivo n, el conjunto de soluciones de

   x^a - y^b = n

1	x^a - y^b = n

para x, a, y, b > 1 es finito.

Por ejemplo, para n = 4, hay 3 soluciones

   (2,3, 2,2) ya que 2³ -  2² =   8 -   4 = 4
   (6,2, 2,5) ya que 6² -  2⁵ =  36 -  32 = 4
   (5,3,11,2) ya que 5³ - 11² = 125 - 121 = 4

(2,3, 2,2) ya que 2³ - 2² = 8 - 4 = 4

(6,2, 2,5) ya que 6² - 2⁵ = 36 - 32 = 4

(5,3,11,2) ya que 5³ - 11² = 125 - 121 = 4

Las soluciones se pueden representar por la menor potencia (en el caso anterior, por 4, 32 y 121) ya que dado n (en el caso anterior es 4), la potencia mayor es la menor más n.

Definir las funciones

   potenciasPerfectas :: [Integer]
   solucionesPillati :: Integer -> [Integer]
   solucionesPillatiAcotadas :: Integer -> Integer -> [Integer]

potenciasPerfectas :: [Integer]

solucionesPillati :: Integer -> [Integer]

solucionesPillatiAcotadas :: Integer -> Integer -> [Integer]

tales que

potenciasPerfectas es la lista de las potencias perfectas (es decir, de los números de la forma x^a con x y a mayores que 1). Por ejemplo,

     take 10 potenciasPerfectas  ==  [4,8,9,16,25,27,32,36,49,64]
     potenciasPerfectas !! 200   ==  28224

1 2	take 10 potenciasPerfectas == [4,8,9,16,25,27,32,36,49,64] potenciasPerfectas !! 200 == 28224

(solucionesPillati n) es la lista de las menores potencias de las soluciones de la ecuación de Pillati x^a – y^b = n; es decir, es la lista de los u tales que u y u+n son potencias perfectas. Por ejemplo,

     take 3 (solucionesPillati 4)  ==  [4,32,121]
     take 2 (solucionesPillati 5)  ==  [4,27]
     take 4 (solucionesPillati 7)  ==  [9,25,121,32761]

take 3 (solucionesPillati 4) == [4,32,121]

take 2 (solucionesPillati 5) == [4,27]

take 4 (solucionesPillati 7) == [9,25,121,32761]

(solucionesPillatiAcotadas c n) es la lista de elementos de (solucionesPillati n) menores que n. Por ejemplo,

     solucionesPillatiAcotadas (10^3) 1  ==  [8]
     solucionesPillatiAcotadas (10^3) 2  ==  [25]
     solucionesPillatiAcotadas (10^3) 3  ==  [125]
     solucionesPillatiAcotadas (10^3) 4  ==  [4,32,121]
     solucionesPillatiAcotadas (10^3) 5  ==  [4,27]
     solucionesPillatiAcotadas (10^3) 6  ==  []
     solucionesPillatiAcotadas (10^3) 7  ==  [9,25,121]
     solucionesPillatiAcotadas (10^5) 7  ==  [9,25,121,32761]

solucionesPillatiAcotadas (10^3) 1 == [8]

solucionesPillatiAcotadas (10^3) 2 == [25]

solucionesPillatiAcotadas (10^3) 3 == [125]

solucionesPillatiAcotadas (10^3) 4 == [4,32,121]

solucionesPillatiAcotadas (10^3) 5 == [4,27]

solucionesPillatiAcotadas (10^3) 6 == []

solucionesPillatiAcotadas (10^3) 7 == [9,25,121]

solucionesPillatiAcotadas (10^5) 7 == [9,25,121,32761]

Soluciones

import Data.List (group)
import Data.Numbers.Primes (primeFactors)

-- Definiciones de potenciasPerfectas
-- ==================================

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

potenciasPerfectas1 :: [Integer]
potenciasPerfectas1 = filter esPotenciaPerfecta [4..]

-- (esPotenciaPerfecta x) se verifica si x es una potencia perfecta. Por
-- ejemplo, 
--    esPotenciaPerfecta 36  ==  True
--    esPotenciaPerfecta 72  ==  False
esPotenciaPerfecta :: Integer -> Bool
esPotenciaPerfecta = not . null. potenciasPerfectasDe 

-- (potenciasPerfectasDe x) es la lista de pares (a,b) tales que 
-- x = a^b. Por ejemplo,
--    potenciasPerfectasDe 64  ==  [(2,6),(4,3),(8,2)]
--    potenciasPerfectasDe 72  ==  []
potenciasPerfectasDe :: Integer -> [(Integer,Integer)]
potenciasPerfectasDe n = 
    [(m,k) | m <- takeWhile (\x -> x*x <= n) [2..]
           , k <- takeWhile (\x -> m^x <= n) [2..]
           , m^k == n]

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

potenciasPerfectas2 :: [Integer]
potenciasPerfectas2 = [x | x <- [4..], esPotenciaPerfecta2 x]

-- (esPotenciaPerfecta2 x) se verifica si x es una potencia perfecta. Por
-- ejemplo, 
--    esPotenciaPerfecta2 36  ==  True
--    esPotenciaPerfecta2 72  ==  False
esPotenciaPerfecta2 :: Integer -> Bool
esPotenciaPerfecta2 x = mcd (exponentes x) > 1

-- (exponentes x) es la lista de los exponentes de l factorización prima
-- de x. Por ejemplos,
--    exponentes 36  ==  [2,2]
--    exponentes 72  ==  [3,2]
exponentes :: Integer -> [Int]
exponentes x = [length ys | ys <- group (primeFactors x)] 

-- (mcd xs) es el máximo común divisor de la lista xs. Por ejemplo,
--    mcd [4,6,10]  ==  2
--    mcd [4,5,10]  ==  1
mcd :: [Int] -> Int
mcd = foldl1 gcd

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

potenciasPerfectas3 :: [Integer]
potenciasPerfectas3 = mezclaTodas potencias

-- potencias es la lista las listas de potencias de todos los números
-- mayores que 1 con exponentes mayores que 1. Por ejemplo,
--    λ> map (take 3) (take 4 potencias)
--    [[4,8,16],[9,27,81],[16,64,256],[25,125,625]]
potencias :: [[Integer]]
potencias = [[n^k | k <- [2..]] | n <- [2..]]

-- (mezclaTodas xss) es la mezcla ordenada sin repeticiones de las
-- listas ordenadas xss. Por ejemplo,
--    take 7 (mezclaTodas potencias)  ==  [4,8,9,16,25,27,32]
mezclaTodas :: Ord a => [[a]] -> [a]
mezclaTodas = foldr1 xmezcla
  where xmezcla (x:xs) ys = x : mezcla xs ys

-- (mezcla xs ys) es la mezcla ordenada sin repeticiones de las
-- listas ordenadas xs e ys. Por ejemplo,
--    take 7 (mezcla [2,5..] [4,6..])  ==  [2,4,5,6,8,10,11]
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

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

--    λ> potenciasPerfectas1 !! 200
--    28224
--    (7.24 secs, 9,245,991,160 bytes)
--    λ> potenciasPerfectas2 !! 200
--    28224
--    (0.30 secs, 814,597,152 bytes)
--    λ> potenciasPerfectas3 !! 200
--    28224
--    (0.01 secs, 7,061,120 bytes)

-- En lo que sigue se usa la 3ª definición
potenciasPerfectas :: [Integer]
potenciasPerfectas = potenciasPerfectas3

-- Definición de solucionesPillati
-- ===============================

solucionesPillati :: Integer -> [Integer]
solucionesPillati n =
  [x | x <- potenciasPerfectas
     , esPotenciaPerfecta2 (x+n)]

-- Definición de solucionesPillatiAcotadas
-- =======================================

solucionesPillatiAcotadas :: Integer -> Integer -> [Integer]
solucionesPillatiAcotadas c n =
  [x | x <- takeWhile (< (c-n)) potenciasPerfectas
     , esPotenciaPerfecta2 (x+n)]

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

import Data.List (group)

import Data.Numbers.Primes (primeFactors)

-- Definiciones de potenciasPerfectas

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

-- 1ª definición

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

potenciasPerfectas1 :: [Integer]

potenciasPerfectas1 = filter esPotenciaPerfecta [4..]

-- (esPotenciaPerfecta x) se verifica si x es una potencia perfecta. Por

-- ejemplo,

-- esPotenciaPerfecta 36 == True

-- esPotenciaPerfecta 72 == False

esPotenciaPerfecta :: Integer -> Bool

esPotenciaPerfecta = not . null. potenciasPerfectasDe

-- (potenciasPerfectasDe x) es la lista de pares (a,b) tales que

-- x = a^b. Por ejemplo,

-- potenciasPerfectasDe 64 == [(2,6),(4,3),(8,2)]

-- potenciasPerfectasDe 72 == []

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

potenciasPerfectasDe n =

[(m,k) | m <- takeWhile (\x -> x*x <= n) [2..]

, k <- takeWhile (\x -> m^x <= n) [2..]

, m^k == n]

-- 2ª definición

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

potenciasPerfectas2 :: [Integer]

potenciasPerfectas2 = [x | x <- [4..], esPotenciaPerfecta2 x]

-- (esPotenciaPerfecta2 x) se verifica si x es una potencia perfecta. Por

-- ejemplo,

-- esPotenciaPerfecta2 36 == True

-- esPotenciaPerfecta2 72 == False

esPotenciaPerfecta2 :: Integer -> Bool

esPotenciaPerfecta2 x = mcd (exponentes x) > 1

-- (exponentes x) es la lista de los exponentes de l factorización prima

-- de x. Por ejemplos,

-- exponentes 36 == [2,2]

-- exponentes 72 == [3,2]

exponentes :: Integer -> [Int]

exponentes x = [length ys | ys <- group (primeFactors x)]

-- (mcd xs) es el máximo común divisor de la lista xs. Por ejemplo,

-- mcd [4,6,10] == 2

-- mcd [4,5,10] == 1

mcd :: [Int] -> Int

mcd = foldl1 gcd

-- 3ª definición

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

potenciasPerfectas3 :: [Integer]

potenciasPerfectas3 = mezclaTodas potencias

-- potencias es la lista las listas de potencias de todos los números

-- mayores que 1 con exponentes mayores que 1. Por ejemplo,

-- λ> map (take 3) (take 4 potencias)

-- [[4,8,16],[9,27,81],[16,64,256],[25,125,625]]

potencias :: [[Integer]]

potencias = [[n^k | k <- [2..]] | n <- [2..]]

-- (mezclaTodas xss) es la mezcla ordenada sin repeticiones de las

-- listas ordenadas xss. Por ejemplo,

-- take 7 (mezclaTodas potencias) == [4,8,9,16,25,27,32]

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

mezclaTodas = foldr1 xmezcla

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

-- (mezcla xs ys) es la mezcla ordenada sin repeticiones de las

-- listas ordenadas xs e ys. Por ejemplo,

-- take 7 (mezcla [2,5..] [4,6..]) == [2,4,5,6,8,10,11]

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

-- Comparación de eficiencia

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

-- λ> potenciasPerfectas1 !! 200

-- 28224

-- (7.24 secs, 9,245,991,160 bytes)

-- λ> potenciasPerfectas2 !! 200

-- 28224

-- (0.30 secs, 814,597,152 bytes)

-- λ> potenciasPerfectas3 !! 200

-- 28224

-- (0.01 secs, 7,061,120 bytes)

-- En lo que sigue se usa la 3ª definición

potenciasPerfectas :: [Integer]

potenciasPerfectas = potenciasPerfectas3

-- Definición de solucionesPillati

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

solucionesPillati :: Integer -> [Integer]

solucionesPillati n =

[x | x <- potenciasPerfectas

, esPotenciaPerfecta2 (x+n)]

-- Definición de solucionesPillatiAcotadas

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

solucionesPillatiAcotadas :: Integer -> Integer -> [Integer]

solucionesPillatiAcotadas c n =

[x | x <- takeWhile (< (c-n)) potenciasPerfectas

, esPotenciaPerfecta2 (x+n)]

Referencia

Catalan’s conjecture en Wikipedia.

Pensamiento

Y te enviaré mi canción:
«Se canta lo que se pierde»,
con un papagayo verde
que la diga en tu balcón.

Antonio Machado

Teorema de Liouville sobre listas CuCu

PorJosé A. Alonso 3-enero-202010-enero-2020

Una lista CuCu es una lista de números enteros positivos tales que la suma de sus Cubos es igual al Cuadrado de su suma. Por ejemplo, [1, 2, 3, 2, 4, 6] es una lista CuCu ya que

   1³ + 2³ + 3³ + 2³ + 4³ + 6³ = (1 + 2 + 3 + 2 + 4 + 6)²

1	1³ + 2³ + 3³ + 2³ + 4³ + 6³ = (1 + 2 + 3 + 2 + 4 + 6)²

La lista de Liouville correspondiente al número entero positivo n es la lista formada por el número de divisores de cada divisor de n. Por ejemplo, para el número 20 se tiene que sus divisores son

   1, 2, 4, 5, 10, 20

1	1, 2, 4, 5, 10, 20

puesto que el número de sus divisores es

El 1 tiene 1 divisor (el 1 solamente).
El 2 tiene 2 divisores (el 1 y el 2).
El 4 tiene 3 divisores (el 1, el 2 y el 4).
El 5 tiene 2 divisores (el 1 y el 5).
El 10 tiene 4 divisores (el 1, el 2, el 5 y el 10).
El 20 tiene 6 divisores (el 1, el 2, el 4, el 5, el 10 y el 20).

la lista de Liouville de 20 es [1, 2, 3, 2, 4, 6] que, como se comentó anteriormente, es una lista CuCu.

El teorema de Lioville afirma que todas las lista de Lioville son CuCu.

Definir las funciones

   esCuCu :: [Integer] -> Bool
   liouville :: Integer -> [Integer]

1 2	esCuCu :: [Integer] -> Bool liouville :: Integer -> [Integer]

tales que

(esCuCu xs) se verifica si la lista xs es CuCu; es decir, la suma de los cubos de sus elementos es igual al cuadrado de su suma. Por ejemplo,

     esCuCu [1,2,3]        ==  True
     esCuCu [1,2,3,2]      ==  False
     esCuCu [1,2,3,2,4,6]  ==  True

esCuCu [1,2,3] == True

esCuCu [1,2,3,2] == False

esCuCu [1,2,3,2,4,6] == True

(liouville n) es la lista de Lioville correspondiente al número n. Por ejemplo,

     liouville 20  ==  [1,2,3,2,4,6]
     liouville 60  ==  [1,2,2,3,2,4,4,6,4,6,8,12]
     length (liouville (product [1..25]))  ==  340032

liouville 20 == [1,2,3,2,4,6]

liouville 60 == [1,2,2,3,2,4,4,6,4,6,8,12]

length (liouville (product [1..25])) == 340032

Comprobar con QuickCheck

que para todo entero positivo n, (liouville (2^n)) es la lista [1,2,3,…,n+1] y
el teorema de Lioville; es decir, para todo entero positivo n, (liouville n) es una lista CuCu.

Nota: Este ejercicio está basado en Cómo generar conjuntos CuCu de Gaussianos.

Soluciones

import Data.List (genericLength, group, inits, sort)
import Data.Numbers.Primes (primeFactors)
import Test.QuickCheck

esCuCu :: [Integer] -> Bool
esCuCu xs = sum (map (^3) xs) == (sum xs)^2

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

liouville :: Integer -> [Integer]
liouville n = map numeroDivisores (divisores n)

-- (divisores x) es el conjunto de divisores de los x. 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]

-- (numeroDivisores x) es el número de divisores de x. Por ejemplo, 
--    numeroDivisores 12  ==  6
--    numeroDivisores 25  ==  3
numeroDivisores :: Integer -> Integer
numeroDivisores n = genericLength (divisores n) 

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

liouville2 :: Integer -> [Integer]
liouville2 n = map numeroDivisores2 (divisores2 n)

-- Se usan las funciones
-- + divisores de "Conjunto de divisores" http://bit.ly/2OtbFIj
-- + numeroDivisores de "Número de divisores" http://bit.ly/2DgVh74

-- (divisores2 x) es el conjunto de divisores de los x. Por ejemplo, 
--   divisores2 30  ==  [1,2,3,5,6,10,15,30]
divisores2 :: Integer -> [Integer]
divisores2 = sort
           . map (product . concat)
           . sequence
           . map inits
           . group
           . primeFactors

-- (numeroDivisores2 x) es el número de divisores de x. Por ejemplo, 
--    numeroDivisores2 12  ==  6
--    numeroDivisores2 25  ==  3
numeroDivisores2 :: Integer -> Integer
numeroDivisores2 =
  product . map ((+1) . genericLength) . group . primeFactors

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

-- La comparación es
--    λ> length (liouville (product [1..11]))
--    540
--    (13.66 secs, 7,983,550,640 bytes)
--    λ> length (liouville2 (product [1..11]))
--    540
--    (0.01 secs, 1,255,328 bytes)

-- Propiedad
-- =========

-- La propiedad es
prop_Liouville :: Integer -> Property
prop_Liouville n =
  n > 0 ==> liouville2 (2^n) == [1..n+1]

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

-- Teorema de Liouville
-- ====================

-- La propiedad es
teorema_Liouville :: Integer -> Property
teorema_Liouville n =
  n > 0 ==> esCuCu (liouville n)

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

import Data.List (genericLength, group, inits, sort)

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

esCuCu :: [Integer] -> Bool

esCuCu xs = sum (map (^3) xs) == (sum xs)^2

-- 1ª definición de liouville

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

liouville :: Integer -> [Integer]

liouville n = map numeroDivisores (divisores n)

-- (divisores x) es el conjunto de divisores de los x. 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]

-- (numeroDivisores x) es el número de divisores de x. Por ejemplo,

-- numeroDivisores 12 == 6

-- numeroDivisores 25 == 3

numeroDivisores :: Integer -> Integer

numeroDivisores n = genericLength (divisores n)

-- 2ª definición de liouville

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

liouville2 :: Integer -> [Integer]

liouville2 n = map numeroDivisores2 (divisores2 n)

-- Se usan las funciones

-- + divisores de "Conjunto de divisores" http://bit.ly/2OtbFIj

-- + numeroDivisores de "Número de divisores" http://bit.ly/2DgVh74

-- (divisores2 x) es el conjunto de divisores de los x. Por ejemplo,

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

divisores2 :: Integer -> [Integer]

divisores2 = sort

. map (product . concat)

. sequence

. map inits

. group

. primeFactors

-- (numeroDivisores2 x) es el número de divisores de x. Por ejemplo,

-- numeroDivisores2 12 == 6

-- numeroDivisores2 25 == 3

numeroDivisores2 :: Integer -> Integer

numeroDivisores2 =

product . map ((+1) . genericLength) . group . primeFactors

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (liouville (product [1..11]))

-- 540

-- (13.66 secs, 7,983,550,640 bytes)

-- λ> length (liouville2 (product [1..11]))

-- 540

-- (0.01 secs, 1,255,328 bytes)

-- Propiedad

-- =========

-- La propiedad es

prop_Liouville :: Integer -> Property

prop_Liouville n =

n > 0 ==> liouville2 (2^n) == [1..n+1]

-- La comprobación es

-- λ> quickCheck prop_Liouville

-- +++ OK, passed 100 tests.

-- Teorema de Liouville

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

-- La propiedad es

teorema_Liouville :: Integer -> Property

teorema_Liouville n =

n > 0 ==> esCuCu (liouville n)

-- La comprobación es

-- λ> quickCheck teorema_Liouville

-- +++ OK, passed 100 tests.

Pensamiento

¡Oh, tarde viva y quieta
que opuso al panta rhei su nada corre.

Antonio Machado

Derivada aritmética

PorJosé A. Alonso 31-diciembre-20197-enero-2020

La derivada aritmética es una función definida sobre los números naturales por analogía con la regla del producto para el cálculo de las derivadas usada en análisis.

Para un número natural n su derivada D(n) se define por

   D(0)  = 0
   D(1)  = 0
   D(p)  = 1, si p es primo
   D(ab) = D(a)b + aD(b) (regla de Leibniz para derivar productos)

D(0) = 0

D(1) = 0

D(p) = 1, si p es primo

D(ab) = D(a)b + aD(b) (regla de Leibniz para derivar productos)

Por ejemplo,

   D(6)  = D(2*3) = D(2)*3 + 2*D(3) = 1*3 + 2*1 =  5
   D(12) = D(2*6) = D(2)*6 + 2*D(6) = 1*6 + 2*5 = 16

1 2	D(6) = D(23) = D(2)3 + 2D(3) = 13 + 21 = 5 D(12) = D(26) = D(2)6 + 2D(6) = 16 + 25 = 16

Definir la función

   derivada :: Integer -> Integer

1	derivada :: Integer -> Integer

tal que (derivada n) es la derivada aritmética de n. Por ejemplo,

   derivada  6  ==  5
   derivada 12  ==  16
   maximum [derivada n | n <- [1..60000]]  ==  380928

derivada 6 == 5

derivada 12 == 16

maximum [derivada n | n <- [1..60000]] == 380928

Comprobar con QuickCheck que si x es un número entero positivo y su descomposición en factores primos es

   x = p(1)^e(1) + p(2)^e(2) +...+ p(n)^e(n)

1	x = p(1)^e(1) + p(2)^e(2) +...+ p(n)^e(n)

entonces la derivada de x es

   x * [e(1)/p(1) + e(2)/p(2) +...+ e(n)/p(n)]

1	x * [e(1)/p(1) + e(2)/p(2) +...+ e(n)/p(n)]

Nota: No usar en la definición la propiedad que hay que comprobar.

Soluciones

import Data.List (genericLength, group)
import Data.Numbers.Primes (isPrime, primeFactors)
import Test.QuickCheck

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

derivada :: Integer -> Integer
derivada 0 = 0
derivada 1 = 0
derivada n | esPrimo n = 1
           | otherwise = (derivada a) * b + a * (derivada b)
  where a = menorFactor n
        b = n `div` a

-- (esPrimo n) se verifica si n es primo. Por ejemplo,
--    esPrimo 5  ==  True
--    esPrimo 6  ==  False
esPrimo :: Integer -> Bool
esPrimo 0 = False
esPrimo 1 = False
esPrimo n = n == menorFactor n

-- (menorFactor n) es el menor divisor primo de n (con n >= 2). Por
-- ejemplo, 
--    menorFactor 6   ==  2
--    menorFactor 7   ==  7
--    menorFactor 15  ==  3
menorFactor :: Integer -> Integer
menorFactor n
  | even n = 2
  | otherwise = head [x | x <- [3,5..]
                        , n `mod` x == 0]

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

derivada2 :: Integer -> Integer
derivada2 0 = 0
derivada2 1 = 0
derivada2 n | isPrime n = 1
            | otherwise = (derivada2 a) * b + a * (derivada2 b)
  where (a:_) = primeFactors n
        b     = n `div` a

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

--    λ> maximum [derivada n | n <- [1..10000]]
--    53248
--    (1.59 secs, 1,091,452,552 bytes)
--    λ> maximum [derivada2 n | n <- [1..10000]]
--    53248
--    (0.17 secs, 457,819,120 bytes)

-- Propiedad
-- =========

-- La propiedad es
prop_derivada :: Integer -> Property
prop_derivada x =
  x > 0 ==>
  derivada x == sum [(x * e) `div` p | (p,e) <- factorizacion 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 :: Integer -> [(Integer,Integer)]
factorizacion n =
  [(head xs,genericLength xs) | xs <- group (primeFactors n)]

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

import Data.List (genericLength, group)

import Data.Numbers.Primes (isPrime, primeFactors)

import Test.QuickCheck

-- 1ª solución

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

derivada :: Integer -> Integer

derivada 0 = 0

derivada 1 = 0

derivada n | esPrimo n = 1

| otherwise = (derivada a) * b + a * (derivada b)

where a = menorFactor n

b = n `div` a

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

-- esPrimo 5 == True

-- esPrimo 6 == False

esPrimo :: Integer -> Bool

esPrimo 0 = False

esPrimo 1 = False

esPrimo n = n == menorFactor n

-- (menorFactor n) es el menor divisor primo de n (con n >= 2). Por

-- ejemplo,

-- menorFactor 6 == 2

-- menorFactor 7 == 7

-- menorFactor 15 == 3

menorFactor :: Integer -> Integer

menorFactor n

| even n = 2

| otherwise = head [x | x <- [3,5..]

, n `mod` x == 0]

-- 2ª solución

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

derivada2 :: Integer -> Integer

derivada2 0 = 0

derivada2 1 = 0

derivada2 n | isPrime n = 1

| otherwise = (derivada2 a) * b + a * (derivada2 b)

where (a:_) = primeFactors n

b = n `div` a

-- Comparación de eficiencia

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

-- λ> maximum [derivada n | n <- [1..10000]]

-- 53248

-- (1.59 secs, 1,091,452,552 bytes)

-- λ> maximum [derivada2 n | n <- [1..10000]]

-- 53248

-- (0.17 secs, 457,819,120 bytes)

-- Propiedad

-- =========

-- La propiedad es

prop_derivada :: Integer -> Property

prop_derivada x =

x > 0 ==>

derivada x == sum [(x * e) `div` p | (p,e) <- factorizacion 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 :: Integer -> [(Integer,Integer)]

factorizacion n =

[(head xs,genericLength xs) | xs <- group (primeFactors n)]

-- Su comprobación es

-- λ> quickCheck prop_derivada

-- +++ OK, passed 100 tests.

Referencias

Arithmetic derivative en Wikipedia.
Sucesión A003415 de OEIS.

Pensamiento

En ese jardín, Guiomar,
el mutuo jardín que inventan
dos corazones al par,
se funden y complementan
nuestras horas.

Antonio Machado

Pares definidos por su MCD y su MCM

PorJosé A. Alonso 13-diciembre-201920-diciembre-2019

Definir las siguientes funciones

   pares  :: Integer -> Integer -> [(Integer,Integer)]
   nPares :: Integer -> Integer -> Integer

1 2	pares :: Integer -> Integer -> [(Integer,Integer)] nPares :: Integer -> Integer -> Integer

tales que

(pares a b) es la lista de los pares de números enteros positivos tales que su máximo común divisor es a y su mínimo común múltiplo es b. Por ejemplo,

     pares 3 3  == [(3,3)]
     pares 4 12 == [(4,12),(12,4)]
     pares 2 12 == [(2,12),(4,6),(6,4),(12,2)]
     pares 2 60 == [(2,60),(4,30),(6,20),(10,12),(12,10),(20,6),(30,4),(60,2)]
     pares 2 7  == []
     pares 12 3  ==  []
     length (pares 3 (product [3,5..91]))  ==  8388608

pares 3 3 == [(3,3)]

pares 4 12 == [(4,12),(12,4)]

pares 2 12 == [(2,12),(4,6),(6,4),(12,2)]

pares 2 60 == [(2,60),(4,30),(6,20),(10,12),(12,10),(20,6),(30,4),(60,2)]

pares 2 7 == []

pares 12 3 == []

length (pares 3 (product [3,5..91])) == 8388608

(nPares a b) es el número de pares de enteros positivos tales que su máximo común divisor es a y su mínimo común múltiplo es b. Por ejemplo,

     nPares 3 3   ==  1
     nPares 4 12  ==  2
     nPares 2 12  ==  4
     nPares 2 60  ==  8
     nPares 2 7   ==  0
     nPares 12 3  ==  0
     nPares 3 (product [3..3*10^4]) `mod` (10^12)  ==  477999992832
     length (show (nPares 3 (product [3..3*10^4])))  ==  977

nPares 3 3 == 1

nPares 4 12 == 2

nPares 2 12 == 4

nPares 2 60 == 8

nPares 2 7 == 0

nPares 12 3 == 0

nPares 3 (product [3..3*10^4]) `mod` (10^12) == 477999992832

length (show (nPares 3 (product [3..3*10^4]))) == 977

Soluciones

import Data.Numbers.Primes (primeFactors)
import Data.List (genericLength, group, nub, sort, subsequences)
import Test.QuickCheck

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

pares1 :: Integer -> Integer -> [(Integer,Integer)]
pares1 a b = [(x,y) | x <- [1..b]
                    , y <- [1..b]
                    , gcd x y == a
                    , lcm x y == b]

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

pares2 :: Integer -> Integer -> [(Integer,Integer)]
pares2 a b = [(x,y) | x <- [a,a+a..b]
                    , y <- [a,a+a..b]
                    , gcd x y == a
                    , lcm x y == b]

-- Comparación de eficiencia
--    λ> length (pares1 3 (product [3,5..11]))
--    16
--    (95.12 secs, 86,534,165,528 bytes)
--    λ> length (pares2 3 (product [3,5..11]))
--    16
--    (15.80 secs, 14,808,762,128 bytes)

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

pares3 :: Integer -> Integer -> [(Integer,Integer)]
pares3 a b = [(x,y) | x <- [a,a+a..b]
                    , c `rem` x == 0
                    , let y = c `div` x
                    , gcd x y == a
                    ]
  where c = a * b

-- Comparacioń de eficiencia
--    λ> length (pares2 3 (product [3,5..11]))
--    16
--    (15.80 secs, 14,808,762,128 bytes)
--    λ> length (pares3 3 (product [3,5..11]))
--    16
--    (0.01 secs, 878,104 bytes)

-- 4ª definición de pares
-- ======================

-- Para la cuarta definición de pares se observa la relación con los
-- factores primos
--    λ> [(primeFactors x, primeFactors y) | (x,y) <- pares1 2 12]
--    [([2],[2,2,3]),([2,2],[2,3]),([2,3],[2,2]),([2,2,3],[2])]
--    λ> [primeFactors x | (x,y) <- pares1 2 12]
--    [[2],[2,2],[2,3],[2,2,3]]
--    λ> [primeFactors x | (x,y) <- pares1 2 60]
--    [[2],[2,2],[2,3],[2,5],[2,2,3],[2,2,5],[2,3,5],[2,2,3,5]]
--    λ> [primeFactors x | (x,y) <- pares1 6 60]
--    [[2,3],[2,2,3],[2,3,5],[2,2,3,5]]
--    λ> [primeFactors x | (x,y) <- pares1 2 24]
--    [[2],[2,3],[2,2,2],[2,2,2,3]]
-- Se observa que cada pares se obtiene de uno de los subconjuntos de los
-- divisores primos de b/a. Por ejemplo,
--    λ> (a,b) = (2,24)
--    λ> b `div` a
--    12
--    λ> primeFactors it
--    [2,2,3]
--    λ> group it
--    [[2,2],[3]]
--    λ> subsequences it
--    [[],[[2,2]],[[3]],[[2,2],[3]]]
--    λ> map concat it
--    [[],[2,2],[3],[2,2,3]]
--    λ> map product it
--    [1,4,3,12]
--    λ> [(a * x, b `div` x) | x <- it]
--    [(2,24),(8,6),(6,8),(24,2)]
-- A partir de la observación se construye la siguiente definición

pares4 :: Integer -> Integer -> [(Integer,Integer)]
pares4 a b
  | b `mod` a /= 0 = []
  | otherwise =
    [(a * x, b `div` x)
    | x <- map (product . concat)
               ((subsequences . group . primeFactors) (b `div` a))]

-- Nota. La función pares4 calcula el mismo conjunto que las anteriores,
-- pero no necesariamente en el mismo orden. Por ejemplo,
--    λ> pares3 2 60 
--    [(2,60),(4,30),(6,20),(10,12),(12,10),(20,6),(30,4),(60,2)]
--    λ> pares4 2 60 
--    [(2,60),(4,30),(6,20),(12,10),(10,12),(20,6),(30,4),(60,2)]
--    λ> pares3 2 60 == sort (pares4 2 60)
--    True

-- Comparacioń de eficiencia
--    λ> length (pares3 3 (product [3,5..17]))
--    64
--    (4.44 secs, 2,389,486,440 bytes)
--    λ> length (pares4 3 (product [3,5..17]))
--    64
--    (0.00 secs, 177,704 bytes)
  
-- Propiedades de equivalencia de pares
-- ====================================

prop_pares :: Integer -> Integer -> Property
prop_pares a b =
  a > 0 && b > 0 ==>
  all (== pares1 a b)
      [sort (f a b) | f <- [ pares2
                           , pares3
                           , pares4
                           ]]
  
prop_pares2 :: Integer -> Integer -> Property
prop_pares2 a b =
  a > 0 && b > 0 ==>
  all (== pares1 a (a * b))
      [sort (f a (a * b)) | f <- [ pares2
                                 , pares3
                                 , pares4
                                 ]]
  
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares
--    +++ OK, passed 100 tests.
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares
--    +++ OK, passed 100 tests.
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares2
--    +++ OK, passed 100 tests.
             
-- 1ª definición de nPares
-- =======================

nPares1 :: Integer -> Integer -> Integer
nPares1 a b = genericLength (pares4 a b)

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

nPares2 :: Integer -> Integer -> Integer
nPares2 a b = 2^(length (nub (primeFactors (b `div` a))))

-- Comparación de eficiencia
--    λ> nPares1 3 (product [3,5..91])
--    8388608
--    (4.68 secs, 4,178,295,920 bytes)
--    λ> nPares2 3 (product [3,5..91])
--    8388608
--    (0.00 secs, 234,688 bytes)

-- Propiedad de equivalencia de nPares
-- ===================================

prop_nPares :: Integer -> Integer -> Property
prop_nPares a b =
  a > 0 && b > 0 ==>
  nPares1 a (a * b) == nPares2 a (a * b)
  
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_nPares
--    +++ OK, passed 100 tests.

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

import Data.Numbers.Primes (primeFactors)

import Data.List (genericLength, group, nub, sort, subsequences)

import Test.QuickCheck

-- 1ª definición de pares

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

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

pares1 a b = [(x,y) | x <- [1..b]

, y <- [1..b]

, gcd x y == a

, lcm x y == b]

-- 2ª definición de pares

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

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

pares2 a b = [(x,y) | x <- [a,a+a..b]

, y <- [a,a+a..b]

, gcd x y == a

, lcm x y == b]

-- Comparación de eficiencia

-- λ> length (pares1 3 (product [3,5..11]))

-- 16

-- (95.12 secs, 86,534,165,528 bytes)

-- λ> length (pares2 3 (product [3,5..11]))

-- 16

-- (15.80 secs, 14,808,762,128 bytes)

-- 3ª definición de pares

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

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

pares3 a b = [(x,y) | x <- [a,a+a..b]

, c `rem` x == 0

, let y = c `div` x

, gcd x y == a

]

where c = a * b

-- Comparacioń de eficiencia

-- λ> length (pares2 3 (product [3,5..11]))

-- 16

-- (15.80 secs, 14,808,762,128 bytes)

-- λ> length (pares3 3 (product [3,5..11]))

-- 16

-- (0.01 secs, 878,104 bytes)

-- 4ª definición de pares

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

-- Para la cuarta definición de pares se observa la relación con los

-- factores primos

-- λ> [(primeFactors x, primeFactors y) | (x,y) <- pares1 2 12]

-- [([2],[2,2,3]),([2,2],[2,3]),([2,3],[2,2]),([2,2,3],[2])]

-- λ> [primeFactors x | (x,y) <- pares1 2 12]

-- [[2],[2,2],[2,3],[2,2,3]]

-- λ> [primeFactors x | (x,y) <- pares1 2 60]

-- [[2],[2,2],[2,3],[2,5],[2,2,3],[2,2,5],[2,3,5],[2,2,3,5]]

-- λ> [primeFactors x | (x,y) <- pares1 6 60]

-- [[2,3],[2,2,3],[2,3,5],[2,2,3,5]]

-- λ> [primeFactors x | (x,y) <- pares1 2 24]

-- [[2],[2,3],[2,2,2],[2,2,2,3]]

-- Se observa que cada pares se obtiene de uno de los subconjuntos de los

-- divisores primos de b/a. Por ejemplo,

-- λ> (a,b) = (2,24)

-- λ> b `div` a

-- 12

-- λ> primeFactors it

-- [2,2,3]

-- λ> group it

-- [[2,2],[3]]

-- λ> subsequences it

-- [[],[[2,2]],[[3]],[[2,2],[3]]]

-- λ> map concat it

-- [[],[2,2],[3],[2,2,3]]

-- λ> map product it

-- [1,4,3,12]

-- λ> [(a * x, b `div` x) | x <- it]

-- [(2,24),(8,6),(6,8),(24,2)]

-- A partir de la observación se construye la siguiente definición

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

pares4 a b

| b `mod` a /= 0 = []

| otherwise =

[(a * x, b `div` x)

| x <- map (product . concat)

((subsequences . group . primeFactors) (b `div` a))]

-- Nota. La función pares4 calcula el mismo conjunto que las anteriores,

-- pero no necesariamente en el mismo orden. Por ejemplo,

-- λ> pares3 2 60

-- [(2,60),(4,30),(6,20),(10,12),(12,10),(20,6),(30,4),(60,2)]

-- λ> pares4 2 60

-- [(2,60),(4,30),(6,20),(12,10),(10,12),(20,6),(30,4),(60,2)]

-- λ> pares3 2 60 == sort (pares4 2 60)

-- True

-- Comparacioń de eficiencia

-- λ> length (pares3 3 (product [3,5..17]))

-- 64

-- (4.44 secs, 2,389,486,440 bytes)

-- λ> length (pares4 3 (product [3,5..17]))

-- 64

-- (0.00 secs, 177,704 bytes)

-- Propiedades de equivalencia de pares

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

prop_pares :: Integer -> Integer -> Property

prop_pares a b =

a > 0 && b > 0 ==>

all (== pares1 a b)

[sort (f a b) | f <- [ pares2

, pares3

, pares4

]]

prop_pares2 :: Integer -> Integer -> Property

prop_pares2 a b =

a > 0 && b > 0 ==>

all (== pares1 a (a * b))

[sort (f a (a * b)) | f <- [ pares2

, pares3

, pares4

]]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares

-- +++ OK, passed 100 tests.

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares

-- +++ OK, passed 100 tests.

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares2

-- +++ OK, passed 100 tests.

-- 1ª definición de nPares

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

nPares1 :: Integer -> Integer -> Integer

nPares1 a b = genericLength (pares4 a b)

-- 2ª definición de nPares

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

nPares2 :: Integer -> Integer -> Integer

nPares2 a b = 2^(length (nub (primeFactors (b `div` a))))

-- Comparación de eficiencia

-- λ> nPares1 3 (product [3,5..91])

-- 8388608

-- (4.68 secs, 4,178,295,920 bytes)

-- λ> nPares2 3 (product [3,5..91])

-- 8388608

-- (0.00 secs, 234,688 bytes)

-- Propiedad de equivalencia de nPares

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

prop_nPares :: Integer -> Integer -> Property

prop_nPares a b =

a > 0 && b > 0 ==>

nPares1 a (a * b) == nPares2 a (a * b)

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_nPares

-- +++ OK, passed 100 tests.

Pensamiento

Largo es el camino de la enseñanza por medio de teorías; breve y eficaz por medio de ejemplos. ~ Séneca

Menor número triangular con más de n divisores

PorJosé A. Alonso 27-noviembre-20194-diciembre-2019

La sucesión de los números triangulares se obtiene sumando los números naturales.

   *     *      *        *         *   
        * *    * *      * *       * *  
              * * *    * * *     * * * 
                      * * * *   * * * *
                               * * * * * 
   1     3      6        10        15

* * * * *

* * * * * * * *

* * * * * * * * *

* * * * * * * *

* * * * *

1 3 6 10 15

Así, el 7º número triangular es

   1 + 2 + 3 + 4 + 5 + 6 + 7 = 28.

1	1 + 2 + 3 + 4 + 5 + 6 + 7 = 28.

Los primeros 10 números triangulares son

   1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

1	1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

Los divisores de los primeros 7 números triangulares son:

    1: 1
    3: 1,3
    6: 1,2,3,6
   10: 1,2,5,10
   15: 1,3,5,15
   21: 1,3,7,21
   28: 1,2,4,7,14,28

1: 1

3: 1,3

6: 1,2,3,6

10: 1,2,5,10

15: 1,3,5,15

21: 1,3,7,21

28: 1,2,4,7,14,28

Como se puede observar, 28 es el menor número triangular con más de 5 divisores.

Definir la función

   menorTriangularConAlMenosNDivisores :: Int -> Integer

1	menorTriangularConAlMenosNDivisores :: Int -> Integer

tal que (menorTriangularConAlMenosNDivisores n) es el menor número triangular que tiene al menos n divisores. Por ejemplo,

   menorTriangularConAlMenosNDivisores 5    ==  28
   menorTriangularConAlMenosNDivisores 50   ==  25200
   menorTriangularConAlMenosNDivisores 500  ==  76576500

menorTriangularConAlMenosNDivisores 5 == 28

menorTriangularConAlMenosNDivisores 50 == 25200

menorTriangularConAlMenosNDivisores 500 == 76576500

Nota: Este ejercicio está basado en el problema 12 del Proyecto Euler

Soluciones

import Data.List (group)
import Data.Numbers.Primes (primeFactors)

menorTriangularConAlMenosNDivisores :: Int -> Integer
menorTriangularConAlMenosNDivisores n = 
  head [x | x <- triangulares, nDivisores x >= n]

-- Nota: Se usarán las funciones
-- + triangulares definida en [Números triangulares](http://bit.ly/2rtr6a3) y
-- + nDivisores definida en [Número de divisores](http://bit.ly/2DgVh74)

-- triangulares es la sucesión de los números triangulares. Por ejemplo,
--    take 10 triangulares  ==  [1,3,6,10,15,21,28,36,45,55]
triangulares :: [Integer]
triangulares = scanl1 (+) [1..]

-- (nDivisores x) es el número de divisores de x. Por ejemplo,
--    nDivisores 28  ==  6
nDivisores :: Integer -> Int
nDivisores = product . map ((+1) . length) . group . primeFactors

import Data.List (group)

import Data.Numbers.Primes (primeFactors)

menorTriangularConAlMenosNDivisores :: Int -> Integer

menorTriangularConAlMenosNDivisores n =

head [x | x <- triangulares, nDivisores x >= n]

-- Nota: Se usarán las funciones

-- + triangulares definida en [Números triangulares](http://bit.ly/2rtr6a3) y

-- + nDivisores definida en [Número de divisores](http://bit.ly/2DgVh74)

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

-- take 10 triangulares == [1,3,6,10,15,21,28,36,45,55]

triangulares :: [Integer]

triangulares = scanl1 (+) [1..]

-- (nDivisores x) es el número de divisores de x. Por ejemplo,

-- nDivisores 28 == 6

nDivisores :: Integer -> Int

nDivisores = product . map ((+1) . length) . group . primeFactors

Pensamiento

«La Matemática es una ciencia experimental y la computación es el experimento.» ~ Rivin

Entradas recientes

Buscar

Correo electrónico