Co-Primalität und die Zahl pi

23

Einführung

Die Zahlentheorie steckt voller Wunder in Form unerwarteter Zusammenhänge. Hier ist einer von ihnen.

Zwei ganze Zahlen sind Co-Prime , wenn sie keine Faktoren gemeinsam andere als 1. Bei einer Zahl haben N , sollten alle Zahlen von 1 bis N . Ziehe zwei solcher Ganzzahlen nach dem Zufallsprinzip (alle Ganzzahlen haben die gleiche Wahrscheinlichkeit, bei jeder Ziehung ausgewählt zu werden; die Ziehungen sind unabhängig und werden ersetzt). Es sei p die Wahrscheinlichkeit, dass die beiden ausgewählten ganzen Zahlen Co-Primzahlen sind. Dann tendiert p zu 6 / π 2 ≈ 0,6079 ..., während N gegen unendlich tendiert.

Die Herausforderung

Der Zweck dieser Herausforderung besteht darin, p als Funktion von N zu berechnen .

Als Beispiel betrachten wir N = 4. Es gibt 16 mögliche Paare, die aus den ganzen Zahlen 1, 2, 3, 4 erhalten werden. 11 dieser Paare sind Co-Prime, nämlich (1,1), (1,2), (1,3), (1,4), (2,1), (3,1), (4,1 ), (2,3), (3,2), (3,4), (4,3). Somit ist p 11/16 = 0,6875 für N = 4.

Der genaue Wert von p muss mit mindestens vier Dezimalstellen berechnet werden . Dies impliziert, dass die Berechnung deterministisch sein muss (im Gegensatz zu Monte Carlo). Es muss sich aber nicht um eine direkte Aufzählung aller Paare wie oben handeln. Es kann jede Methode angewendet werden.

Funktionsargumente oder stdin / stdout können verwendet werden. Bei der Anzeige der Ausgabe können nachfolgende Nullen weggelassen werden. So kann zum Beispiel 0.6300angezeigt werden als 0.63. Es sollte als Dezimalzahl und nicht als Bruchzahl angezeigt werden (die Anzeige der Zeichenfolge 63/100ist nicht zulässig).

Das Gewinnkriterium sind die wenigsten Bytes. Die Verwendung der integrierten Funktionen unterliegt keinen Einschränkungen.

Testfälle

Eingabe / Ausgabe (nur vier Dezimalstellen sind obligatorisch, wie oben angegeben):

1    / 1.000000000000000
2    / 0.750000000000000
4    / 0.687500000000000
10   / 0.630000000000000
100  / 0.608700000000000
1000 / 0.608383000000000
Luis Mendo
quelle
Gibt es Grenzen für den Eingabebereich?
Eric Towers
@EricTowers Das Programm sollte für jede vernünftige Größe bis hin zu Einschränkungen des Arbeitsspeichers und des Datentyps funktionieren. Mindestens 1000
Luis Mendo
Sind rationale Zahlen als Rückgabewerte (keine Strings) erlaubt? Meine Sprache hat einen nativen rationalen Typ, in dem 63/100es sich um ein gültiges Literal handelt, das sich für die Berechnung eignet. (Langs ich spreche über: Faktor , Schläger )
Katze
@cat Na klar, mach weiter! Beachten Sie jedoch die erforderliche Präzision
Luis Mendo

Antworten:

14

Gelee , 12 8 Bytes

RÆṪSḤ’÷²

Probieren Sie es online!

Der folgende Binärcode funktioniert mit dieser Version des Jelly-Interpreters .

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

Idee

Es ist klar, dass die Anzahl der Paare (j, k), so dass j ≤ k und j und k Co-Primzahlen sind, gleich der Anzahl der Paare (k, j) ist , die die gleichen Bedingungen erfüllen. Auch wenn j = k , ist j = 1 = k .

Um also die Anzahl der Co-Prim-Paare mit Koordinaten zwischen 1 und n zu zählen , genügt es, die Anzahl m der Paare (j, k) so zu berechnen , dass j ≤ k ist , und dann 2m - 1 zu berechnen .

Da schließlich die Eulersche Totientenfunktion φ (k) ergibt die Anzahl ganzer Zahlen zwischen zwischen 1 und K , die co-prim sind k , können wir berechnen , m als φ (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
quelle
2
Oh, Jelly beinhaltet die Totientenfunktion !? Gute Idee!
Luis Mendo
2
Countdown bis MATL ein Totient-Kommando am T-1-Tag enthält ...
Quintopia
@quintopia (ich habe endlich die Totient-Funktion eingefügt) :-D
Luis Mendo
14

Mathematica 43 42 Bytes

Die vom Ursprung aus sichtbaren Gitterpunkte, von denen das folgende Bild stammt, haben sich als hilfreich erwiesen, um das Problem anhand der folgenden Fragen in Bezug auf einen bestimmten quadratischen Bereich des Einheitsgitters neu zu definieren:

  • Welcher Bruchteil der Einheitsgitterpunkte hat Co-Prim-Koordinaten?
  • Welcher Bruchteil von Einheitsgitterpunkten kann vom Ursprung aus gesehen werden?

Gitter


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

Beispiele

Nachgestellte Nullen werden weggelassen.

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}


Zeitliche Koordinierung

Das Timing in Sekunden geht der Antwort voraus.

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

{0,605571, 0,608383}

DavidC
quelle
6

Mathematica, 42 32 Bytes

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

Anonyme Funktion, verwendet einfache rohe Gewalt. Der letzte Fall dauert auf meinem Computer etwa 0,37 Sekunden. Erläuterung:

                               &   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
quelle
Können Sie ein Beispiel und eine Erklärung für diejenigen von uns veröffentlichen, die Mathematica nicht haben?
Luis Mendo
2
Dies vereint unsere Aussagen: Count[Array[GCD,{#, #}],1,2]/#^2.& Seien Sie mein Gast.
DavidC
4

Dyalog APL, 15 Bytes

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

Ziemlich einfach. Es ist ein monadischer Funktionszug. Iota ist die Zahl von 1 bis zur Eingabe, also nehmen wir das äußere Produkt mit gcd und zählen dann den Anteil der Einsen.

Lirtosiast
quelle
3

Oktave, 49 47 Bytes

Einfach die Summe gcdaller Paare berechnen und zählen.

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

Das kronecker Produkt ist fantastisch.

Fehler
quelle
kron! Gute Idee!
Luis Mendo
Ich habe es zuerst benutzt meshgrid, aber dann habe ich gemerkt, dass ich dasselbe mit kronInline machen kann! (-> anonyme Funktion)
Fehler
2

JavaScript (ES6), 88 Byte

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

Erläuterung

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

Prüfung

Dauert eine Weile für große ( >100) Werte von n.

user81655
quelle
2

Im Ernst, 15 Bytes

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

Hex Dump:

2c3ba62975317860b1604de4e7442f

Probieren Sie es online

Ich werde mich nicht darum kümmern, es zu erklären, da es buchstäblich genau den gleichen Algorithmus verwendet wie Dennis 'Jelly-Lösung (obwohl ich es unabhängig abgeleitet habe).

Quintopie
quelle
2

J, 19 Bytes

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

Verwendung:

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

Erläuterung:

*:%~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

Dennis 'Lösung gibt eine nette Erklärung, wie wir die Totientenfunktion verwenden können.

Probieren Sie es hier online aus.

randomra
quelle
2

Mathematica, 35 Bytes

Implementiert den @ Dennis-Algorithmus.

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

Berechnen Sie die Summe des Totienten (Euler-Phi-Funktion) über den Bereich von 1 bis zum Eingabewert. Multiplizieren Sie mit Gleitkomma zwei (mit vier Stellen Genauigkeit) und subtrahieren Sie eine. (Sie können die Genauigkeit erhöhen, indem Sie stattdessen " 2" und " 1`4" verwenden.) Dies ist die Gesamtzahl der Coprime-Paare. Teilen Sie diese durch das Quadrat der Eingabe, um den gewünschten Bruch zu erhalten. Da es sich bei den beiden um eine ungefähre Zahl handelt, gilt dies auch für das Ergebnis.

Testen (mit Timing-Daten in der linken Spalte, da mindestens einer von uns dies für interessant hält), mit der zugewiesenen Funktion, fdamit die Testzeile leichter gelesen werden kann:

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

Bearbeiten: Umfang des Eingabebereichs anzeigen (die Genauigkeit auf eins statt auf zwei ändern, da sonst die Ergebnisse eher eintönig werden) und andere dazu auffordern, dasselbe zu tun ...

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[]Führt die Berechnung mehrmals durch und nimmt einen Mittelwert der Zeiten, wobei versucht wird, Cold-Caches und andere Effekte zu ignorieren, die Timing-Ausreißer verursachen. Die asymptotische Grenze ist

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

Für Argumente> 10 ^ 4 können wir also einfach "0.6079" zurückgeben und die angegebenen Genauigkeits- und Genauigkeitsanforderungen erfüllen.

Eric Towers
quelle
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)))

Nur der einfache Ansatz für den Moment; Ich werde bald auf kürzere / elegantere Lösungen eingehen. Dies ist eine anonyme Funktion, die eine Ganzzahl akzeptiert und einen Gleitkommawert zurückgibt. Um es aufzurufen, weisen Sie es einer Variablen zu.

Ungolfed:

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.
quelle
Soweit ich das beurteilen kann, brauchen Sie kein collectfaules Objekt, um es zu nehmen length.
Katze
@cat Sie tun dies für einige Fälle, in denen lengthkeine Methode definiert ist, wie hier für den gefilterten Kombinationsiterator. Ähnlich endofwürde es nicht funktionieren, da es für diesen Typ keine Methode gibt getindex.
Alex A.
@cat rangegibt nicht dieselbe Art von Objekt zurück wie combinations. Letzterer gibt einen Iterator über Kombinationen zurück, die ein separater Typ ohne definierte lengthMethode sind. Sehen Sie hier . (Auch die :Syntax wird vorgezogen, um rangeBereiche zu konstruieren;)
Alex A.
2

Salbei, 55 Bytes

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

Dank Sage Computing tauchen keine symbolischen Probleme mit dem Maschinen-Epsilon und den Gleitkommazahlen auf. Der Kompromiss ist, um der Ausgabeformatregel zu folgen, ein zusätzlicher Aufruf n()(der Dezimalnäherungsfunktion) erforderlich.

Probieren Sie es online aus

Mego
quelle
Sehr schön! Sie scheinen Sage in letzter Zeit ziemlich oft zu benutzen :-)
Luis Mendo
@ LuisMendo Sage ist großartig und macht alles. Es ist sehr schön, es für mathematische Herausforderungen zu verwenden, da es eine riesige eingebaute Bibliothek wie Mathematica hat, aber die Syntax ist besser (da a) nicht Mathematica ist und b) auf Python aufgebaut ist).
Mego
2

MATL , 20 17 Bytes

Dies verwendet die aktuelle Version (5.0.0) der Sprache.

Ansatz basiert auf der Antwort von @ flawr .

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

Bearbeiten (28. April 2015) : Probieren Sie es online! Nachdem diese Antwort gepostet wurde, wurde die Funktion Y)in umbenannt X:. Der Link enthält diese Änderung.

Beispiel

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

Erläuterung

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

Alte Antwort: 20 Bytes

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

Erläuterung

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
quelle
Könnten Sie mit einem Ansatz wie dem, den ich in der Oktave verwendet habe, nicht noch kürzer sein?
Fehler
Tatsächlich! Vielen Dank! 3 Bytes weniger. Du hättest es selbst in MATL machen sollen :-)
Luis Mendo
Ich hätte es versucht, wenn es nicht weit nach meiner Schlafenszeit gewesen wäre =)
Fehler
1

PARI / GP , 25 Bytes

Das Anonymisieren der Funktion würde ein Byte sparen, aber dann müsste selfes insgesamt teurer werden.

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

Faktor 120 113 Bytes

Ich habe Golf gespielt und kann es nicht kürzer machen.

Übersetzung von: Julia .

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

Das Beispiel wird für die ersten 5 Testfälle ausgeführt (ein Wert von 1000 bewirkt, dass der Editor eingefroren wird, und ich kann mich nicht mehr darum kümmern, eine ausführbare Datei zu kompilieren):

! 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 }
Katze
quelle
Füge vielleicht einen Beispiellauf hinzu?
Luis Mendo
1
@ LuisMendo fertig!
Katze
1

Samau , 12 Bytes

Haftungsausschluss: Kein Wettbewerb, da ich die Sprache nach dem Posten der Frage aktualisiert habe.

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

Hex-Dump (Samau verwendet CP737-Codierung):

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

Verwenden Sie denselben Algorithmus wie Dennis 'Antwort in Jelly.

Alephalpha
quelle
0

Python2 / Pypy, 178 Bytes

Die xDatei:

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

Laufen:

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

Der Code zählt die Co-Prim-Paare (n,m) for m<nzweimal ( c+=2). Dies ignoriert, (i,i) for i=1..nwas bis auf eine Ausnahme in Ordnung ist (1,1), und wird daher durch Initialisieren des Zählers mit 1( 1.0zur Vorbereitung auf die Schwebeteilung später) zur Kompensation korrigiert .


quelle