Co-primality dan nomor pi

23

pengantar

Teori bilangan penuh dengan keajaiban, dalam bentuk koneksi yang tidak terduga. Ini salah satunya.

Dua bilangan bulat yang co-prime jika mereka tidak memiliki faktor kesamaan selain 1. Mengingat nomor N , mempertimbangkan semua bilangan bulat dari 1 sampai N . Gambar dua bilangan bulat seperti itu secara acak (semua bilangan bulat memiliki probabilitas yang sama untuk dipilih pada setiap undian; undian bersifat independen dan dengan penggantian). Misalkan p menunjukkan probabilitas bahwa dua bilangan bulat yang dipilih adalah co-prime. Kemudian p cenderung ke 6 / π 2 ≈ 0,6079 ... karena N cenderung tak hingga.

Tantangan

Tujuan dari tantangan ini adalah untuk menghitung p sebagai fungsi dari N .

Sebagai contoh, pertimbangkan N = 4. Ada 16 pasangan yang mungkin diperoleh dari bilangan bulat 1,2,3,4. 11 dari pasangan tersebut adalah co-prime, yaitu (1,1), (1,2), (1,3), (1,4), (2,1), (3,1), (4,1) ), (2,3), (3,2), (3,4), (4,3). Jadi p adalah 11/16 = 0,6875 untuk N = 4.

The tepat nilai p perlu dihitung dengan setidaknya empat desimal. Ini menyiratkan bahwa perhitungan harus bersifat deterministik (berlawanan dengan Monte Carlo). Tapi itu tidak perlu menjadi penghitungan langsung dari semua pasangan seperti di atas; metode apa pun dapat digunakan.

Argumen fungsi atau stdin / stdout dapat digunakan. Jika menampilkan output, trailing zero dapat dihilangkan. Jadi misalnya 0.6300bisa ditampilkan sebagai 0.63. Ini harus ditampilkan sebagai angka desimal, bukan sebagai pecahan (menampilkan string 63/100tidak diperbolehkan).

Kriteria yang menang adalah byte paling sedikit. Tidak ada batasan dalam penggunaan fungsi bawaan.

Uji kasus

Input / output (hanya empat desimal yang wajib, seperti yang ditunjukkan di atas):

1    / 1.000000000000000
2    / 0.750000000000000
4    / 0.687500000000000
10   / 0.630000000000000
100  / 0.608700000000000
1000 / 0.608383000000000
Luis Mendo
sumber
Apakah ada batasan pada kisaran input?
Eric Towers
@EricTowers Program harus bekerja untuk ukuran yang masuk akal hingga batasan memori dan tipe data. Setidaknya 1000
Luis Mendo
Apakah angka rasional sebagai nilai pengembalian (bukan string) diizinkan? Bahasa saya memiliki tipe rasional asli, di mana 63/100merupakan literal yang valid, dapat digunakan dalam perhitungan. (Langs yang saya bicarakan: Factor , Racket )
cat
@cat Tentu, silakan! Namun, perhitungkan ketelitian yang diperlukan
Luis Mendo

Jawaban:

14

Jelly , 12 8 byte

RÆṪSḤ’÷²

Cobalah online!

Kode biner berikut berfungsi dengan versi juru bahasa Jelly ini .

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

Ide

Jelas, jumlah pasangan (j, k) sedemikian rupa sehingga j ≤ k dan j dan k adalah co-prime sama dengan jumlah pasangan (k, j) yang memenuhi kondisi yang sama. Juga, jika j = k , j = 1 = k .

Jadi, untuk menghitung jumlah pasangan co-prime dengan koordinat antara 1 dan n , cukup untuk menghitung jumlah m dari pasangan (j, k) sedemikian sehingga j ≤ k , kemudian menghitung 2m - 1 .

Akhirnya, karena fungsi total Euler φ (k) menghasilkan bilangan bulat antara antara 1 dan k yang co-prime ke k , kita dapat menghitung m sebagai φ (1) + ... + φ (n) .

Kode

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
sumber
2
Oh, Jelly termasuk fungsi totient !? Ide bagus!
Luis Mendo
2
Hitung mundur hingga MATL menyertakan perintah total pada T-1 hari ...
quintopia
@quintopia (saya akhirnya memasukkan fungsi totient) :-D
Luis Mendo
14

Mathematica 43 42 byte

Saya menemukan titik-titik Kisi terlihat dari asal , dari mana gambar di bawah ini diambil, untuk membantu dalam membingkai ulang masalah melalui pertanyaan-pertanyaan berikut mengenai wilayah kuadrat tertentu dari unit kisi:

  • Berapa fraksi poin unit-kisi yang memiliki koordinat co-prime?
  • Apa fraksi poin unit-kisi dapat dilihat dari asal?

kisi


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

Contohnya

Nol trailing dihilangkan.

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}


Pengaturan waktu

Waktunya, dalam detik, mendahului jawabannya.

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

{0.605571, 0.608383}

DavidC
sumber
6

Mathematica, 42 32 byte

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

Fungsi anonim, menggunakan kekuatan kasar sederhana. Kasing terakhir beroperasi sekitar 0,37 detik pada mesin saya. Penjelasan:

                               &   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
sumber
Bisakah Anda memposting contoh lari dan penjelasan, bagi kita yang tidak memiliki Mathematica?
Luis Mendo
2
Ini menyatukan kiriman kami: Count[Array[GCD,{#, #}],1,2]/#^2.& Jadilah tamu saya.
DavidC
4

Dyalog APL, 15 byte

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

Cukup mudah. Ini adalah kereta fungsi monadik. Iota adalah angka dari 1 hingga input, jadi kami mengambil produk luar dengan gcd, lalu menghitung proporsi yang.

lirtosiast
sumber
3

Oktaf, 49 47 byte

Hanya menghitung gcdsemua pasangan dan menghitung.

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

Produk Kroner mengagumkan.

cacat
sumber
kron! Ide bagus!
Luis Mendo
Saya pertama kali menggunakan meshgrid, tetapi kemudian saya perhatikan saya bisa melakukan hal yang sama dengan kroninline! (-> fungsi anonim)
flawr
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`)

Penjelasan

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

Uji

Butuh beberapa saat untuk >100nilai ( ) yang besar dari n.

pengguna81655
sumber
2

Serius, 15 byte

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

Hex Dump:

2c3ba62975317860b1604de4e7442f

Cobalah secara Online

Saya tidak akan repot menjelaskannya karena secara harfiah menggunakan algoritma yang persis sama dengan solusi Dennis's Jelly (meskipun saya mendapatkannya secara mandiri).

kuintopia
sumber
2

J, 19 17 byte

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

Pemakaian:

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

Penjelasan:

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

Solusi Dennis memberikan penjelasan yang bagus bagaimana kita bisa menggunakan fungsi totient.

Cobalah online di sini.

randomra
sumber
2

Mathematica, 35 byte

Menerapkan algoritma @ Dennis.

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

Hitung jumlah totient (fungsi phi Euler) pada rentang dari 1 hingga nilai input. Lipat gandakan dengan floating point dua (dengan presisi empat digit) dan kurangi satu. (Lebih presisi dapat dipertahankan dengan menggunakan sebagai gantinya " 2" dan " 1`4".) Ini adalah jumlah total pasangan koprime, jadi bagi dengan kuadrat input untuk mendapatkan fraksi yang diinginkan. Karena keduanya adalah angka perkiraan, demikian juga hasilnya.

Pengujian (dengan data penghitungan waktu di kolom kiri karena setidaknya salah satu dari kami menganggap itu menarik), dengan fungsi yang ditetapkan fagar garis uji lebih mudah dibaca .:

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

Sunting: Memperlihatkan jangkauan jangkauan input (menukar ketelitian dengan yang satu bukannya yang dua karena sebaliknya hasilnya menjadi agak monoton) dan menantang orang lain untuk melakukan hal yang sama ...

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[]melakukan perhitungan beberapa kali dan mengambil rata-rata kali, berusaha mengabaikan cache dingin dan efek lain yang menyebabkan pencilan waktu. Batas asimptotik adalah

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

jadi untuk argumen> 10 ^ 4, kita bisa mengembalikan "0.6079" dan memenuhi persyaratan presisi dan akurasi yang ditentukan.

Eric Towers
sumber
2

Julia, 95 byte

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

Hanya pendekatan sederhana untuk saat ini; Saya akan mencari solusi yang lebih pendek / lebih elegan segera. Ini adalah fungsi anonim yang menerima integer dan mengembalikan float. Untuk menyebutnya, tetapkan ke variabel.

Tidak Disatukan:

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.
sumber
Sejauh yang saya tahu, Anda tidak perlu collectobjek malas untuk mengambilnya length.
kucing
@ kucing Anda lakukan untuk beberapa tempat yang lengthtidak memiliki metode yang ditentukan, yang merupakan kasus di sini untuk iterator kombinasi yang difilter. Demikian pula endoftidak akan berhasil karena tidak ada metode untuk tipe itu getindex.
Alex A.
@cat rangetidak mengembalikan objek yang sama combinations. Yang terakhir mengembalikan iterator atas kombinasi yang merupakan tipe terpisah tanpa lengthmetode yang ditentukan . Lihat di sini . (Juga :sintaks lebih disukai daripada rangeuntuk membangun rentang;))
Alex A.
2

Sage, 55 byte

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

Berkat Sage menghitung semuanya secara simbolis, masalah epsilon alat berat dan titik apung tidak muncul. Tradeoff adalah, untuk mengikuti aturan format output, diperlukan panggilan tambahan untuk n()(fungsi aproksimal desimal).

Cobalah online

Mego
sumber
Sangat bagus! Anda tampaknya menggunakan Sage cukup sering akhir-akhir ini :-)
Luis Mendo
@LuisMendo Sage hebat dan melakukan semua hal. Ini sangat baik untuk digunakan pada tantangan berbasis matematika karena memiliki perpustakaan builtin besar seperti Mathematica, tetapi sintaks lebih baik (berdasarkan a) tidak menjadi Mathematica, dan b) sedang dibangun di Python).
Mego
2

MATL , 20 17 byte

Ini menggunakan versi bahasa saat ini (5.0.0) .

Pendekatan berdasarkan jawaban @ flawr .

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

Sunting (28 April 2015) : Cobalah online! Setelah jawaban ini diposting, fungsi Y)diubah namanya menjadi X:; tautan termasuk perubahan itu.

Contoh

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

Penjelasan

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

Jawaban lama: 20 byte

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

Penjelasan

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
sumber
Tidak bisakah Anda lebih pendek dengan pendekatan seperti yang saya gunakan dalam oktaf?
flawr
Memang! Terima kasih! 3 byte lebih sedikit. Anda seharusnya melakukannya sendiri di MATL :-)
Luis Mendo
Saya akan mencoba jika itu tidak melewati waktu tidur saya =)
flawr
1

PARI / GP , 25 byte

Membuat fungsi anonim akan menghemat satu byte, tetapi kemudian saya harus menggunakannya selfmembuatnya lebih mahal secara keseluruhan.

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

Factor, 120 113 byte

Saya telah menghabiskan golf kelas ini dan saya tidak bisa membuatnya lebih pendek.

Terjemahan dari: Julia .

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

Contoh berjalan pada 5 kasus uji pertama (nilai 1000 menyebabkannya membekukan editor, dan saya tidak dapat diganggu untuk mengkompilasi yang dapat dieksekusi sekarang):

! 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 }
kucing
sumber
Tambahkan contoh run, mungkin?
Luis Mendo
1
@LuisMendo selesai!
kucing
1

Samau , 12 byte

Penafian: Tidak bersaing karena saya memperbarui bahasa setelah pertanyaan diposting.

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

Hex dump (Samau menggunakan pengkodean CP737):

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

Menggunakan algoritma yang sama dengan jawaban Dennis di Jelly.

alephalpha
sumber
0

Python2 / Pypy, 178 byte

Itu x File:

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

Berlari:

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

Kode ini menghitung pasangan co-prime (n,m) for m<ndua kali ( c+=2). Ini mengabaikan (i,i) for i=1..nyang ok kecuali untuk (1,1), sehingga dikoreksi dengan menginisialisasi counter dengan 1( 1.0untuk mempersiapkan divisi float nanti) untuk mengkompensasi.


sumber