Wprowadzenie

Teoria liczb jest pełna cudów w postaci nieoczekiwanych połączeń. Oto jeden z nich.

Dwie liczby całkowite są współ-prime , jeśli nie mają one wspólne czynniki inne niż 1. Biorąc pod uwagę liczbę N , należy rozważyć wszystkie liczby całkowite od 1 do N . Losuj dwie takie liczby całkowite losowo (wszystkie liczby całkowite mają takie samo prawdopodobieństwo, że zostaną wybrane przy każdym losowaniu; losowania są niezależne i zastępują). Niech p oznacza prawdopodobieństwo, że dwie wybrane liczby całkowite są pierwszymi. Następnie p dąży do 6 / π ² ≈ 0,6079 ... tak jak N dąży do nieskończoności.

Wyzwanie

Celem tego wyzwania jest obliczenie s w zależności od N .

Jako przykład rozważmy N = 4. Istnieje 16 możliwych par uzyskanych z liczb całkowitych 1,2,3,4. 11 z tych par jest równoległych, a mianowicie (1,1), (1,2), (1,3), (1,4), (2,1), (3,1), (4,1 ), (2,3), (3,2), (3,4), (4,3). Zatem p wynosi 11/16 = 0,6875 dla N = 4.

Dokładną wartość p musi być obliczana z co najmniej czterech miejsc po przecinku. Oznacza to, że obliczenia muszą być deterministyczne (w przeciwieństwie do Monte Carlo). Ale nie musi to być bezpośrednie wyliczenie wszystkich par jak wyżej; można zastosować dowolną metodę.

Można użyć argumentów funkcyjnych lub stdin / stdout. Jeśli wyświetlasz wynik, zera końcowe można pominąć. Na przykład 0.6300może być wyświetlany jako 0.63. Powinien być wyświetlany jako liczba dziesiętna, a nie jako ułamek (wyświetlanie łańcucha 63/100jest niedozwolone).

Kryterium wygranej to najmniej bajtów. Nie ma żadnych ograniczeń dotyczących korzystania z wbudowanych funkcji.

Przypadki testowe

Wejście / wyjście (tylko cztery miejsca po przecinku są obowiązkowe, jak wskazano powyżej):

1    / 1.000000000000000
2    / 0.750000000000000
4    / 0.687500000000000
10   / 0.630000000000000
100  / 0.608700000000000
1000 / 0.608383000000000

Galaretka , 12 8 bajtów

RÆṪSḤ’÷²

Wypróbuj online!

Poniższy kod binarny działa z tą wersją interpretera Jelly .

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

Pomysł

Oczywiście liczba par (j, k), tak że j ≤ k i j i k są równe, jest równa liczbie par (k, j), które spełniają te same warunki. Ponadto, jeśli j = k , j = 1 = k .

Tak więc, aby policzyć liczbę współpierwotnych par o współrzędnych od 1 do n , wystarczy obliczyć liczbę m par (j, k) taką, że j ≤ k , a następnie obliczyć 2m - 1 .

Wreszcie, ponieważ funkcja sumy Eulera φ (k) daje liczby całkowite od 1 do k, które są równe liczbie pierwotnej do k , możemy obliczyć m jako φ (1) +… + φ (n) .

Kod

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².

Mathematica 43 42 bajty

Zauważyłem, że punkty Kraty widoczne z początku , z którego wykonano poniższe zdjęcie, są pomocne w przeformułowaniu problemu poprzez następujące pytania dotyczące danego kwadratowego obszaru siatki jednostki:

• Jaka część punktów sieci-jednostki ma współrzędne pierwotne?
• Jaką część punktów siatki jednostkowej można zobaczyć na początku?

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

Przykłady

Zera końcowe są pomijane.

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}

wyczucie czasu

Czas w sekundach poprzedza odpowiedź.

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

{0,605571, 0,608383}

Pyth - 13 11 bajtów

.OqR1iM^SQ2

Pakiet testowy .

Mathematica, 42 32 bajty

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

Funkcja anonimowa, wykorzystuje prostą brutalną siłę. Ostatnia sprawa działa na moim komputerze za około 0,37 sekundy. Wyjaśnienie:

                               &   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.

Dyalog APL, 15 bajtów

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

Całkiem proste. Jest to monadyczny ciąg funkcji. Iota to liczby od 1 do wejścia, więc bierzemy zewnętrzny produkt przez gcd, a następnie liczymy proporcję jedności.

Oktawa, 49 47 bajtów

Tylko obliczanie gcdwszystkich par i liczenie.

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

Produkt Kronecker jest niesamowity.

JavaScript (ES6), 88 bajtów

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

Wyjaśnienie

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

Test

Zajmuje trochę czasu dla dużych ( >100) wartości n.

var solution = 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`)

N = <input type="text" id="input" value="100" />
<button onclick="result.textContent=solution(+input.value)">Go</button>
<pre id="result"></pre>

Poważnie, 15 bajtów

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

Hex Dump:

2c3ba62975317860b1604de4e7442f

Wypróbuj online

Nie zamierzam zawracać sobie głowy wyjaśnianiem tego, ponieważ dosłownie używa dokładnie tego samego algorytmu, co rozwiązanie Jelly'ego Jelly'ego (chociaż wyprowadziłem go niezależnie).

J, 19 17 bajtów

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

Stosowanie:

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

Wyjaśnienie:

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

Rozwiązanie Dennisa daje ładne wyjaśnienie, w jaki sposób możemy użyć funkcji totient.

Wypróbuj online tutaj.

Mathematica, 35 bajtów

Implementuje algorytm @ Dennisa.

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

Oblicz sumę totemu (funkcja ph Eulera) w zakresie od 1 do wartości wejściowej. Pomnóż przez liczbę zmiennoprzecinkową dwa (z czterema cyframi dokładności) i odejmij jedną. (Można zachować większą precyzję, używając zamiast tego „ 2” i „ 1`4”.) Jest to całkowita liczba par coprime, więc podziel przez kwadrat danych wejściowych, aby uzyskać żądany ułamek. Ponieważ dwa są liczbą przybliżoną, taki jest też wynik.

Testowanie (z danymi pomiaru czasu w lewej kolumnie, ponieważ przynajmniej jeden z nas uważa to za interesujące), z przypisaną funkcją, aby fłatwiej było odczytać linię testową:

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

Edycja: Pokazywanie zakresu zakresu danych wejściowych (zamiana precyzji na jeden zamiast dwóch, ponieważ w przeciwnym razie wyniki stają się raczej monotonne) i zachęcanie innych do zrobienia tego samego ...

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[]wykonuje obliczenia wiele razy i zajmuje to średnio tyle czasu, próbując zignorować zimne pamięci podręczne i inne efekty powodujące wartości odstające od czasu. Limit asymptotyczny wynosi

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

więc dla argumentów> 10 ^ 4 możemy po prostu zwrócić „0.6079” i spełnić określone wymagania dotyczące precyzji i dokładności.

Julia, 95 bajtów

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

Na razie proste podejście; Wkrótce zajmę się krótszymi / bardziej eleganckimi rozwiązaniami. Jest to anonimowa funkcja, która przyjmuje liczbę całkowitą i zwraca liczbę zmiennoprzecinkową. Aby go wywołać, przypisz go do zmiennej.

Nie golfowany:

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

Szałwia, 55 bajtów

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

Dzięki Sage obliczając wszystko symbolicznie, problemy z epsilon i zmiennoprzecinkowymi maszyn nie pojawiają się. Kompromis polega na tym, że aby zachować zgodność z regułą formatu wyjściowego, konieczne jest dodatkowe wywołanie funkcji n()(funkcja przybliżenia dziesiętnego).

Wypróbuj online

MATL , 20 17 bajtów

Używa bieżącej wersji (5.0.0) języka.

Podejście oparte na odpowiedzi @ flawr .

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

Edytuj (28 kwietnia 2015 r.) : Wypróbuj online! Po opublikowaniu tej odpowiedzi Y)zmieniono nazwę funkcji na X:; link zawiera tę zmianę.

Przykład

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

Wyjaśnienie

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

Stara odpowiedź: 20 bajtów

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

Wyjaśnienie

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

PARI / GP , 25 bajtów

Anonimowość tej funkcji pozwoliłaby zaoszczędzić bajt, ale musiałbym użyć, selfaby ogólnie była droższa.

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

Czynnik, 120 113 bajtów

Spędziłem lekcje gry w golfa i nie mogę tego skrócić.

Tłumaczenie: Julia .

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

Przykład działa na pierwszych 5 testowych przypadkach (wartość 1000 powoduje zawieszenie edytora i nie mogę sobie teraz pozwolić na kompilację pliku wykonywalnego):

! 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 }

Samau , 12 bajtów

Oświadczenie: Nie konkuruje, ponieważ zaktualizowałem język po opublikowaniu pytania.

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

Zrzut szesnastkowy (Samau używa kodowania CP737):

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

Korzystanie z tego samego algorytmu, co odpowiedź Dennisa w Galaretce.

Python2 / Pypy, 178 bajtów

xPliku:

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

Bieganie:

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

Kod zlicza pary równoległe (n,m) for m<ndwa razy ( c+=2). Ignoruje (i,i) for i=1..nto, co jest w porządku, poza (1,1)tym, że jest korygowane przez zainicjowanie licznika za pomocą 1( 1.0aby przygotować się do podziału zmiennoprzecinkowego później) w celu kompensacji.