Co-primalidad y el número pi

23

Introducción

La teoría de números está llena de maravillas, en forma de conexiones inesperadas. Aquí hay uno de ellos.

Dos enteros son co-prime si no tienen factores en común distinto de 1. Dado un número N , tenga en cuenta todos los números enteros de 1 a N . Dibuje dos de estos enteros al azar (todos los enteros tienen la misma probabilidad de ser seleccionados en cada sorteo; los sorteos son independientes y con reemplazo). Supongamos que p denota la probabilidad de que los dos enteros seleccionados sean coprimos. Entonces p tiende a 6 / π 2 ≈ 0.6079 ... como N tiende a infinito.

El reto

El propósito de este reto es calcular p como una función de N .

Como ejemplo, considere N = 4. Hay 16 pares posibles obtenidos de los enteros 1,2,3,4. 11 de esos pares son primos, a saber (1,1), (1,2), (1,3), (1,4), (2,1), (3,1), (4,1 ), (2,3), (3,2), (3,4), (4,3). Por lo tanto, p es 16/11 = 0.6875 para N = 4.

El valor exacto de p debe calcularse con al menos cuatro decimales. Esto implica que el cálculo tiene que ser determinista (a diferencia de Monte Carlo). Pero no es necesario que sea una enumeración directa de todos los pares como se indicó anteriormente; Se puede utilizar cualquier método.

Se pueden usar argumentos de función o stdin / stdout. Si se muestra la salida, se pueden omitir los ceros finales. Entonces, por ejemplo, 0.6300se puede mostrar como 0.63. Debe mostrarse como un número decimal, no como una fracción ( 63/100no se permite mostrar la cadena ).

El criterio ganador es la menor cantidad de bytes. No hay restricciones en el uso de funciones integradas.

Casos de prueba

Entrada / salida (solo cuatro decimales son obligatorios, como se indicó anteriormente):

1    / 1.000000000000000
2    / 0.750000000000000
4    / 0.687500000000000
10   / 0.630000000000000
100  / 0.608700000000000
1000 / 0.608383000000000
Luis Mendo
fuente
¿Hay algún límite en el rango de entradas?
Eric Towers
@EricTowers El programa debería funcionar para cualquier tamaño razonable hasta limitaciones de memoria y tipo de datos. Al menos 1000
Luis Mendo
¿Están permitidos los números racionales como valores de retorno (no cadenas)? Mi lenguaje tiene un tipo racional nativo, en el cual 63/100es un literal válido, utilizable en computación. (Langs estoy hablando de: Factor , Racket )
gato
@cat ¡Claro, adelante! Sin embargo, tenga en cuenta la precisión requerida
Luis Mendo

Respuestas:

14

Jalea , 12 8 bytes

RÆṪSḤ’÷²

Pruébalo en línea!

El siguiente código binario funciona con esta versión del intérprete Jelly .

0000000: 52 91 b0 53 aa b7 9a 8a  R..S....

Idea

Claramente, el número de pares (j, k) tal que j ≤ k y j y k son primos iguales es igual al número de pares (k, j) que satisfacen las mismas condiciones. Además, si j = k , j = 1 = k .

Por lo tanto, para contar el número de pares de co-prime con coordenadas entre 1 y n , es suficiente para calcular la cantidad m de pares (j, k) que tal j ≤ k , entonces compute 2m - 1 .

Por último, ya que la función totient de Euler φ (k) se obtienen los números enteros entre números entre 1 y k que son co-prime a k , podemos calcular m como φ (1) + ... + φ (n) .

Código

RÆṪSḤ’÷²    Input: n

R           Yield [1, ..., n].
 ÆṪ         Apply Euler's totient function to each k in [1, ..., n].
   S        Compute the sum of all results.
    Ḥ       Multiply the result by 2.
     ’      Subtract 1.
      ÷²    Divide the result by n².
Dennis
fuente
2
Oh, ¿Jelly incluye la función totient? ¡Buena idea!
Luis Mendo
2
La cuenta regresiva hasta que MATL incluye un comando totient en el día T-1 ...
quintopia
@quintopia (finalmente incluí la función totient) :-D
Luis Mendo
14

Mathematica 43 42 bytes

Encontré que los puntos de celosía visibles desde el origen , de los cuales se toma la siguiente imagen, son útiles para replantear el problema a través de las siguientes preguntas con respecto a una región cuadrada dada de la red de unidades:

  • ¿Qué fracción de los puntos de celosía unitaria tienen coordenadas co-primas?
  • ¿Qué fracción de puntos de celosía unitaria puede verse desde el origen?

cuadrícula


N@Mean[Mean/@Boole@Array[CoprimeQ,{#,#}]]&

Ejemplos

Los ceros finales se omiten.

N@Mean[Mean/@Boole@Array[CoprimeQ,{#,#}]]&/@Range@10

{1., 0.75, 0.777778, 0.6875, 0.76, 0.638889, 0.714286, 0.671875, 0.679012, 0.63}


Sincronización

El tiempo, en segundos, precede a la respuesta.

N@Mean[Mean/@Boole@Array[CoprimeQ,{#,#}]]&[1000]// AbsoluteTiming

{0.605571, 0.608383}

DavidC
fuente
6

Mathematica, 42 32 bytes

Count[GCD~Array~{#,#},1,2]/#^2.&

Función anónima, utiliza fuerza bruta simple. El último caso se ejecuta en aproximadamente .37 segundos en mi máquina. Explicación:

                               &   A function taking N and returning
Count[               , , ]           the number of
                      1               ones
                                     in the
                        2             second
                                     level of
         ~Array~                      a matrix of
      GCD                              the greatest common denominators of
                {#,#}                 all 2-tuples in [1..N]
                          /         divided by
                           #          N
                            ^2.      squared.
LegionMammal978
fuente
¿Puede publicar un ejemplo de ejecución y explicación para aquellos de nosotros que no tenemos Mathematica?
Luis Mendo
2
Esto une nuestras presentaciones: Count[Array[GCD,{#, #}],1,2]/#^2.& Sé mi invitado.
DavidC
4

Dyalog APL, 15 bytes

(+/÷⍴)∘∊1=⍳∘.∨⍳

Muy claro. Es un tren de funciones monádicas. Iota son los números del 1 a la entrada, por lo que tomamos el producto externo por mcd, luego contamos la proporción de unos.

lirtosiast
fuente
3

Octava, 49 47 bytes

Simplemente calculando el gcdde todos los pares y contando.

@(n)mean(mean(gcd(c=kron(ones(n,1),1:n),c')<2))

El producto kronecker es impresionante.

falla
fuente
kron! ¡Buena idea!
Luis Mendo
Lo utilicé por primera vez meshgrid, ¡pero luego noté que podía hacer lo mismo con kroninline! (-> función anónima)
flawr
2

JavaScript (ES6), 88 bytes

n=>eval(`p=0;for(i=n;i;i--)for(j=n;j;j--,p+=a)for(a=1,k=j;k>1;k--)a=i%k||j%k?a:0;p/n/n`)

Explicación

n=>
  eval(`                     // use eval to enable for loop without {} or return
    p=0;                     // p = number of pairs
    for(i=n;i;i--)           // i = 1 to n
      for(j=n;j;j--,p+=a)    // j = i to n, a will equal 1 if i and j are coprime, else 0
        for(a=1,k=j;k>1;k--) // for each number from 0 to j
          a=i%k||j%k?a:0;    // if i%k and j%k are both 0, this pair is not coprime
    p/n/n                    // return result (equivalent to p/(n*n))
  `)

Prueba

Toma un tiempo para >100valores grandes ( ) de n.

usuario81655
fuente
2

En serio, 15 bytes

,;ª)u1x`▒`MΣτD/

Hex Dump:

2c3ba62975317860b1604de4e7442f

Pruébalo en línea

No me voy a molestar en explicarlo, ya que literalmente usa exactamente el mismo algoritmo que la solución Jelly de Dennis (aunque lo deduje de forma independiente).

quintapia
fuente
2

J, 19 17 bytes

*:%~1-~5+/@p:|@i:

Uso:

   (*:%~1-~5+/@p:|@i:) 4
0.6875

Explicación:

*:%~1-~5+/@p:|@i:
               i: [-n..n]
             |@   absolute value of each element ([n..1,0,1,..n])
       5+/@p:     sum of the totient function for each element
    1-~           decreased by one, giving the number of co-prime pairs
*:%~              divided by N^2

La solución de Dennis da una buena explicación de cómo podemos usar la función totient.

Pruébelo en línea aquí.

randomra
fuente
2

Mathematica, 35 bytes

Implementa el algoritmo de @ Dennis.

(2`4Plus@@EulerPhi@Range[#]-1)/#^2&

Calcule la suma del totient (función phi de Euler) en el rango de 1 al valor de entrada. Multiplique por coma flotante dos (con cuatro dígitos de precisión) y reste uno. (Se puede retener más precisión usando en su lugar " 2" y " 1`4".) Este es el número total de pares coprimos, así que divídalo por el cuadrado de la entrada para obtener la fracción deseada. Debido a que los dos son un número aproximado, también lo es el resultado.

Pruebas (con datos de tiempo en la columna de la izquierda ya que al menos uno de nosotros piensa que es interesante), con la función asignada para fque la línea de prueba se lea más fácilmente .:

f=(2`4Plus@@EulerPhi@Range[#]-1)/#^2&
RepeatedTiming[f[#]] & /@ {1, 2, 4, 10, 100, 1000}
(* {{5.71*10^-6, 1.000}, 
    {5.98*10^-6, 0.750}, 
    {0.000010  , 0.6875}, 
    {0.0000235 , 0.6300}, 
    {0.00028   , 0.6087}, 
    {0.0033    , 0.6084} }  *)

Editar: Mostrar la extensión del rango de entradas (intercambiando la precisión a una en lugar de las dos porque de lo contrario los resultados se vuelven bastante monótonos) y desafiando a otros a hacer lo mismo ...

f = (2 Plus @@ EulerPhi@Range[#] - 1`4)/#^2 &
{#}~Join~RepeatedTiming[f[#]] & /@ {1, 2, 4, 10, 100, 1000, 10^4, 10^5, 10^6, 10^7}
(*  Results are {input, wall time, output}
    {{       1,  5.3*10^-6, 1.000}, 
     {       2,  6.0*10^-6, 0.7500}, 
     {       4,  0.0000102, 0.68750}, 
     {      10,  0.000023 , 0.63000}, 
     {     100,  0.00028  , 0.6087000}, 
     {    1000,  0.0035   , 0.608383000}, 
     {   10000,  0.040    , 0.60794971000}, 
     {  100000,  0.438    , 0.6079301507000}, 
     { 1000000,  4.811    , 0.607927104783000}, 
     {10000000, 64.0      , 0.60792712854483000}}  *)

RepeatedTiming[]realiza el cálculo varias veces y toma un promedio de las veces, intentando ignorar los cachés fríos y otros efectos que causan valores atípicos de tiempo. El límite asintótico es

N[6/Pi^2,30]
(*  0.607927101854026628663276779258  *)

entonces, para argumentos> 10 ^ 4, podemos simplemente devolver "0.6079" y cumplir con los requisitos de precisión y exactitud especificados.

Eric Towers
fuente
2

Julia, 95 bytes

n->(c=combinations([1:n;1:n],2);((L=length)(collect(filter(x->gcd(x...)<2,c)))÷2+1)/L(∪(c)))

Solo el enfoque directo por ahora; Buscaré soluciones más cortas / más elegantes pronto. Esta es una función anónima que acepta un entero y devuelve un flotante. Para llamarlo, asígnelo a una variable.

Sin golf:

function f(n::Integer)
    # Get all pairs of the integers from 1 to n
    c = combinations([1:n; 1:n], 2)

    # Get the coprime pairs
    p = filter(x -> gcd(x...) == 1, c)

    # Compute the probability
    return (length(collect(p)) ÷ 2 + 1) / length(unique(c))
end
Alex A.
fuente
Por lo que puedo decir, no necesitas collectun objeto perezoso para tomarlo length.
gato
@cat Lo haces para algunos lugares donde lengthno tiene un método definido, que es el caso aquí para el iterador de combinaciones filtradas. Del mismo modo endof, no funcionaría porque no hay un método para ese tipo getindex.
Alex A.
@cat rangeno devuelve el mismo tipo de objeto que combinations. Este último devuelve un iterador sobre combinaciones que es un tipo separado sin un lengthmétodo definido . Ver aquí . (También :se prefiere la sintaxis sobre la rangeconstrucción de rangos;))
Alex A.
2

Sabio, 55 bytes

lambda a:n((sum(map(euler_phi,range(1,a+1)))*2-1)/a**2)

Gracias a Sage que computa todo simbólicamente, los problemas de épsilon de máquina y punto flotante no surgen. La compensación es, para seguir la regla del formato de salida, n()se necesita una llamada adicional a (la función de aproximación decimal).

Pruébalo en línea

Mego
fuente
¡Muy agradable! Parece que estás usando Sage con bastante frecuencia últimamente :-)
Luis Mendo
@LuisMendo Sage es genial y hace todas las cosas. Es muy agradable de usar en desafíos matemáticos porque tiene una gran biblioteca integrada como Mathematica, pero la sintaxis es mejor (en virtud de que a) no es Mathematica yb) está construida en Python).
Mego
2

MATL , 20 17 bytes

Esto usa la versión actual (5.0.0) del idioma.

Enfoque basado en la respuesta de @ flawr .

i:tg!X*t!Zd1=Y)Ym

Editar (28 de abril de 2015) : ¡ Pruébelo en línea! Después de que se publicó esta respuesta, Y)se cambió el nombre de la función a X:; El enlace incluye ese cambio.

Ejemplo

>> matl i:tg!X*t!Zd1=Y)Ym
> 100
0.6087

Explicación

i:         % vector 1, 2, ... up to input number
tg!        % copy, convert into ones, transpose
X*         % Kronecker product. Produces a matrix
t!         % copy, transpose
Zd         % gcd of all pairs
1=         % is it equal to 1?
Y)         % linearize into column array
Ym         % compute mean

Respuesta anterior: 20 bytes

Oi:t"t@Zd1=sb+w]n2^/

Explicación

O             % produce literal 0. Initiallizes count of co-prime pairs.
i             % input number, say N
:             % create vector 1, 2, ... N
t             % duplicate of vector
"             % for loop
    t         % duplicate of vector
    @         % loop variable: each element from vector
    Zd        % gcd of vector and loop variable. Produces a vector
    1=s       % number of "1" in result. Each "1" indicates co-primality
    b+w       % accumulate into count of co-prime pairs
]             % end
n2^/          % divide by N^2
Luis Mendo
fuente
¿No podrías ser aún más corto con un enfoque como el que usé en octava?
defecto
¡En efecto! ¡Gracias! 3 bytes menos. Deberías haberlo hecho tú mismo en MATL :-)
Luis Mendo
Lo hubiera intentado si no hubiera pasado mi hora de dormir =)
error
1

PARI / GP , 25 bytes

Hacer que la función sea anónima ahorraría un byte, pero luego tendría que usar para selfhacerla más costosa en general.

f(n)=n^2-sum(j=2,n,f(n\j))
Charles
fuente
1

Factor, 120 113 bytes

He pasado clase jugando al golf y no puedo acortarlo.

Traducción de: Julia .

[ [a,b] dup append 2 <combinations> [ members ] keep [ first2 coprime? ] filter [ length ] bi@ 2 /i 1 + swap /f ]

El ejemplo se ejecuta en los primeros 5 casos de prueba (un valor de 1000 hace que congele el editor, y no puedo molestarme en compilar un ejecutable en este momento):

! with floating point division
IN: scratchpad auto-use {
      1    
      2    
      4    
      10   
      100  
    }
    [ 
      [1,b] dup append 2 <combinations> [ members ] keep 
      [ first2 coprime? ] filter [ length ] bi@ 2 /i 1 + swap /f 
    ]
    map

--- Data stack:
{ 1.0 0.75 0.6875 0.63 0.6087 }
! with rational division
IN: scratchpad auto-use {
      1    
      2    
      4    
      10   
      100  
    }
    [ 
      [1,b] dup append 2 <combinations> [ members ] keep 
      [ first2 coprime? ] filter [ length ] bi@ 2 /i 1 + swap / 
    ]
    map

--- Data stack:
{ 1.0 0.75 0.6875 0.63 0.6087 }
{ 1 3/4 11/16 63/100 6087/10000 }
gato
fuente
Añadir un ejemplo ejecutar tal vez?
Luis Mendo
1
@LuisMendo hecho!
gato
1

Samau , 12 bytes

Descargo de responsabilidad: No estoy compitiendo porque actualicé el idioma después de que se publicó la pregunta.

▌;\φΣ2*($2^/

Volcado hexadecimal (Samau utiliza la codificación CP737):

dd 3b 5c ad 91 32 2a 28 24 32 5e 2f

Usando el mismo algoritmo que la respuesta de Dennis en Jelly.

alephalpha
fuente
0

Python2 / Pypy, 178 bytes

El xarchivo:

N={1:set([1])}
n=0
c=1.0
e=input()
while n<e:
 n+=1
 for d in N[n]:
  m=n+d
  try:N[m].add(d)
  except:N[m]=set([d,m])
 for m in range(1,n):
  if N[m]&N[n]==N[1]:c+=2
print c/n/n

Corriendo:

$ pypy x <<<1
1.0
$ pypy x <<<10
0.63
$ pypy x <<<100
0.6087
$ pypy x <<<1000
0.608383

El código cuenta los pares co-prime dos (n,m) for m<nveces ( c+=2). Esto ignora (i,i) for i=1..ncuál está bien, excepto por (1,1)lo que se corrige inicializando el contador con 1( 1.0para prepararse para la división de flotación más adelante) para compensar.


fuente