Co-primalité et le nombre pi

23

introduction

La théorie des nombres regorge de merveilles, sous la forme de connexions inattendues. En voici un.

Deux entiers sont co-prime si elles ne présentent aucun facteur en commun autre que 1. Étant donné un nombre N , considérer tous les entiers de 1 à N . Dessinez deux de ces entiers au hasard (tous les entiers ont la même probabilité d'être sélectionnés à chaque tirage; les tirages sont indépendants et avec remplacement). Soit p la probabilité que les deux entiers sélectionnés soient co-premiers. Alors p tend vers 6 / π 2 ≈ 0,6079 ... comme N tend vers l'infini.

Le défi

Le but de ce défi est de calculer p en fonction de N .

Par exemple, considérons N = 4. Il y a 16 paires possibles obtenues à partir des entiers 1,2,3,4. 11 de ces paires sont co-principales, à savoir (1,1), (1,2), (1,3), (1,4), (2,1), (3,1), (4,1 ), (2,3), (3,2), (3,4), (4,3). Ainsi p est 11/16 = 0,6875 pour N = 4.

La valeur exacte de p doit être calculée avec au moins quatre décimales. Cela implique que le calcul doit être déterministe (par opposition à Monte Carlo). Mais il n'est pas nécessaire que ce soit une énumération directe de toutes les paires comme ci-dessus; n'importe quelle méthode peut être utilisée.

Des arguments de fonction ou stdin / stdout peuvent être utilisés. Si vous affichez la sortie, les zéros de fin peuvent être omis. Ainsi, par exemple, 0.6300peut être affiché comme 0.63. Il doit être affiché sous la forme d'un nombre décimal et non sous forme d'une fraction (l'affichage de la chaîne 63/100n'est pas autorisé).

Le critère de gain est le moins d'octets. Il n'y a aucune restriction sur l'utilisation des fonctions intégrées.

Cas de test

Entrée / sortie (seules quatre décimales sont obligatoires, comme indiqué ci-dessus):

1    / 1.000000000000000
2    / 0.750000000000000
4    / 0.687500000000000
10   / 0.630000000000000
100  / 0.608700000000000
1000 / 0.608383000000000
Luis Mendo
la source
Y a-t-il des limites sur la gamme des entrées?
Eric Towers
@EricTowers Le programme devrait fonctionner pour toute taille raisonnable jusqu'aux limites de mémoire et de type de données. Au moins 1000
Luis Mendo
Les nombres rationnels en tant que valeurs de retour (pas de chaînes) sont-ils autorisés? Ma langue a un type rationnel natif, dans lequel 63/100est un littéral valide, utilisable en calcul. (Langs dont je parle: Factor , Racket )
cat
@cat Bien sûr, allez-y! Tenez compte de la précision requise, cependant
Luis Mendo

Réponses:

14

Gelée , 12 8 octets

RÆṪSḤ’÷²

Essayez-le en ligne!

Le code binaire suivant fonctionne avec cette version de l'interpréteur Jelly .

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

Idée

Il est clair que le nombre de paires (j, k) telles que j ≤ k et j et k sont co-premiers est égal au nombre de paires (k, j) qui remplissent les mêmes conditions. De plus, si j = k , j = 1 = k .

Ainsi, pour compter le nombre de paires co-amorcées avec des coordonnées comprises entre 1 et n , il suffit de calculer le nombre m de paires (j, k) telles que j ≤ k , puis de calculer 2m - 1 .

Enfin, puisque la fonction de totaux d'Euler φ (k) donne le nombre entier compris entre 1 et k qui sont co-premiers à k , nous pouvons calculer m comme φ (1) +… + φ (n) .

Code

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
la source
2
Oh, Jelly inclut la fonction totient!? Bonne idée!
Luis Mendo
2
Compte à rebours jusqu'à ce que MATL inclue une commande de totient à T-1 jour ...
quintopie
@quintopia (j'ai finalement inclus la fonction totient) :-D
Luis Mendo
14

Mathematica 43 42 octets

J'ai trouvé des points de réseau visibles depuis l'origine , à partir desquels la photo ci-dessous est prise, pour aider à recadrer le problème à travers les questions suivantes concernant une région carrée donnée du réseau unitaire:

  • Quelle fraction des points unité-réseau a des coordonnées co-prime?
  • Quelle fraction de points de réseau unitaire peut être vue depuis l'origine?

la grille


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

Exemples

Les zéros de fin sont omis.

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}


Timing

Le timing, en secondes, précède la réponse.

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

{0.605571, 0.608383}

DavidC
la source
6

Mathematica, 42 32 octets

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

Fonction anonyme, utilise une force brute simple. Le dernier cas fonctionne en environ 0,37 seconde sur ma machine. Explication:

                               &   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
la source
Pouvez-vous publier un exemple d'exécution et une explication, pour ceux d'entre nous qui n'ont pas Mathematica?
Luis Mendo
2
Cela unit nos soumissions: Count[Array[GCD,{#, #}],1,2]/#^2.& Soyez mon invité.
DavidC
4

Dyalog APL, 15 octets

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

Assez simple. C'est un train de fonction monadique. Iota est le nombre de 1 à l'entrée, nous prenons donc le produit extérieur par pgcd, puis comptons la proportion de ceux-ci.

lirtosiast
la source
3

Octave, 49 47 octets

Il suffit de calculer la gcdde toutes les paires et de compter.

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

Le produit kronecker est génial.

flawr
la source
kron! Bonne idée!
Luis Mendo
J'ai d'abord utilisé meshgrid, mais j'ai remarqué que je pouvais faire la même chose en kronligne! (-> fonction anonyme)
flawr
2

JavaScript (ES6), 88 octets

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

Explication

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

Tester

Prend un certain temps pour les grandes >100valeurs ( ) de n.

user81655
la source
2

Sérieusement, 15 octets

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

Vidage hexadécimal:

2c3ba62975317860b1604de4e7442f

Essayez-le en ligne

Je ne vais pas m'embêter à l'expliquer car il utilise littéralement exactement le même algorithme que la solution Jelly de Dennis (bien que je l'ait dérivé indépendamment).

quintopie
la source
2

J, 19 17 octets

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

Usage:

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

Explication:

*:%~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 solution de Dennis donne une belle explication sur la façon d'utiliser la fonction de totient.

Essayez-le en ligne ici.

randomra
la source
2

Mathematica, 35 octets

Implémente l'algorithme de @ Dennis.

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

Calculez la somme du total (fonction phi d'Euler) sur la plage de 1 à la valeur d'entrée. Multipliez par virgule flottante deux (avec quatre chiffres de précision) et soustrayez un. (Plus de précision peut être conservée en utilisant à la place " 2" et " 1`4".) Il s'agit du nombre total de paires de coprimes, donc divisez par le carré de l'entrée pour obtenir la fraction souhaitée. Parce que les deux sont un nombre approximatif, le résultat l'est également.

Test (avec les données de chronométrage dans la colonne de gauche car au moins l'un d'entre nous pense que c'est intéressant), avec la fonction assignée à fpour que la ligne de test soit plus facile à lire .:

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} }  *)

Edit: montrant l'étendue de la plage des entrées (permutant la précision à l'une au lieu des deux car sinon les résultats deviennent plutôt monotones) et incitant les autres à faire de même ...

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[]effectue le calcul plusieurs fois et prend une moyenne des temps, essayant d'ignorer les caches froides et autres effets provoquant des valeurs aberrantes de synchronisation. La limite asymptotique est

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

ainsi, pour des arguments> 10 ^ 4, nous pouvons simplement renvoyer "0,6079" et répondre aux exigences de précision et d'exactitude spécifiées.

Eric Towers
la source
2

Julia, 95 octets

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

Juste l'approche simple pour l'instant; J'examinerai bientôt des solutions plus courtes / plus élégantes. Il s'agit d'une fonction anonyme qui accepte un entier et renvoie un flottant. Pour l'appeler, affectez-le à une variable.

Non 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.
la source
Autant que je sache, vous n'avez pas besoin d' collectun objet paresseux pour le prendre length.
cat
@cat Vous le faites pour certains où lengthaucune méthode n'est définie, ce qui est le cas ici pour l'itérateur de combinaisons filtrées. De même endofne fonctionnerait pas car il n'y a pas de méthode pour ce type getindex.
Alex A.
@cat rangene retourne pas le même type d'objet que combinations. Ce dernier retourne un itérateur sur les combinaisons qui est un type séparé sans lengthméthode définie . Voyez ici . (De plus, la :syntaxe est préférée rangeà la construction de plages;))
Alex A.
2

Sauge, 55 octets

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

Grâce à Sage qui calcule tout symboliquement, les problèmes de machine epsilon et de virgule flottante ne surviennent pas. Le compromis est, afin de suivre la règle de format de sortie, un appel supplémentaire à n()(la fonction d'approximation décimale) est nécessaire.

Essayez-le en ligne

Mego
la source
Très agréable! Vous semblez utiliser Sage assez souvent ces derniers temps :-)
Luis Mendo
@LuisMendo Sage est génial et fait tout. C'est très agréable à utiliser pour les défis basés sur les mathématiques car il a une énorme bibliothèque intégrée comme Mathematica, mais la syntaxe est meilleure (en raison de a) ne pas être Mathematica, et b) être construit sur Python).
Mego
2

MATL , 20 17 octets

Cela utilise la version actuelle (5.0.0) de la langue.

Approche basée sur la réponse de @ flawr .

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

Edit (28 avril 2015) : Essayez-le en ligne! Une fois cette réponse publiée, la fonction a Y)été renommée X:; le lien inclut ce changement.

Exemple

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

Explication

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

Ancienne réponse: 20 octets

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

Explication

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
la source
Ne pourriez-vous pas être encore plus court avec une approche comme celle que j'ai utilisée en octave?
flawr
Effectivement! Merci! 3 octets de moins. Vous auriez dû le faire vous-même en MATL :-)
Luis Mendo
J'aurais essayé si ce n'était pas bien passé mon heure de coucher =)
flawr
1

PARI / GP , 25 octets

Rendre la fonction anonyme économiserait un octet, mais je devrais alors l'utiliser pour la selfrendre plus chère dans l'ensemble.

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

Facteur, 120 113 octets

J'ai passé des cours à jouer au golf et je ne peux pas le raccourcir.

Traduction de: Julia .

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

L'exemple s'exécute sur les 5 premiers cas de test (une valeur de 1000 provoque le gel de l'éditeur, et je ne peux pas être dérangé pour compiler un exécutable pour le moment):

! 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 }
chat
la source
Ajouter un exemple de course peut-être?
Luis Mendo
1
@LuisMendo fait!
cat
1

Samau , 12 octets

Avis de non-responsabilité: Pas en compétition car j'ai mis à jour la langue après la publication de la question.

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

Vidage hexadécimal (Samau utilise le codage CP737):

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

En utilisant le même algorithme que la réponse de Dennis dans Jelly.

alephalpha
la source
0

Python2 / Pypy, 178 octets

Le xdossier:

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

Fonctionnement:

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

Le code compte deux (n,m) for m<nfois les paires co-amorces ( c+=2). Ceci ignore (i,i) for i=1..nce qui est correct à l'exception de (1,1), étant ainsi corrigé en initialisant le compteur avec 1( 1.0pour préparer la division flottante plus tard) pour compenser.


la source