[Question]
  [
# About the Series
First off, you may treat this like any other code golf challenge, and answer it without worrying about the series at all. However, there is a leaderboard across all challenges. You can find the leaderboard along with some more information about the series [in the first post](https://codegolf.stackexchange.com/q/45302/8478).
Although I have a bunch of ideas lined up for the series, the future challenges are not set in stone yet. If you have any suggestions, please let me know [on the relevant sandbox post](http://meta.codegolf.stackexchange.com/a/4750/8478).
# Hole 4: The Bertrand Paradox
The [Bertrand paradox](http://en.wikipedia.org/wiki/Bertrand_paradox_(probability)#Jaynes.27_solution_using_the_.22maximum_ignorance.22_principle) is an interesting problem, which shows how different methods for picking random chords in a circle can yield different distributions of chords, their midpoints and their lengths.
In this challenge you are supposed to generate random chords of the unit circle, using the "right" method, i.e. one which produces a distribution of chords which is invariant under scaling and translation. In the linked Wikipedia article, "Method 2" is such a method.
Here are the exact rules:
* You should take **one positive integer `N`** which specifies how many chords should be returned. The output should be a list of `N` chords, each specified as two points on the unit circle, given by their polar angle in radians.
* Your code should **be able to return at least 220 different values for each of the two angles**. If your available RNG has a smaller range, you must either first build an RNG with a sufficiently large range on top of the built-in one or you must implement your own [suitable RNG](http://meta.codegolf.stackexchange.com/a/1325/8478). [This page](http://en.wikipedia.org/wiki/Linear_congruential_generator#Parameters_in_common_use) may be helpful for that.
* The distribution of chords must be indistinguishable from the one produced by "Method 2" in the linked Wikipedia article. If you implement a different algorithm to choose chords, please include a proof of correctness. Whatever algorithm you choose to implement, it must **theoretically be able to generate any valid chord in the unit circle** (barring limitations of the underlying PRNG or limited-precision data types).
* Your implementation should use and return either **floating-point numbers (at least 32 bits wide)** or **fixed-point numbers (at least 24 bits wide)** and all arithmetic operations should be accurate within at most 16 [ulp](http://en.wikipedia.org/wiki/Unit_in_the_last_place).
You may write a full program or a function and take input via STDIN (or closest alternative), command-line argument or function argument and produce output via STDOUT (or closest alternative), function return value or function (out) parameter.
Output may be in any convenient list or string format, as long as the individual numbers are clearly distinguishable and their total number is always even.
This is code golf, so the shortest submission (in bytes) wins. And of course, the shortest submission per user will also enter into the overall leaderboard of the series.
## Visualisation
You can use the following snippet to render the generated lines and inspect their distribution. Simply paste a list of pairs of angles into the text area. The snippet should be able to handle almost any list format, as long as the numbers are simple decimal numbers (no scientific notation). I recommend you use at least 1000 lines to get a good idea of the distribution. I've also provided some example data for the different methods presented in the article below.
```
function draw() {
  document.getElementById("output").innerHTML = svg
}
function drawLines() {
  lines = document.getElementById("lines").value;
  var e = prefix;
  //e += '<circle cx="' + offset + '" + cy="' + offset + '" r="' + radius + '" stroke="black" fill="white"/>';
  lines = lines.match(/-?(?:\d*\.\d+|\d+\.\d*|\d+(?!\.))/g);
  for (i = 0; i < lines.length; i+=2)
  {
    t1 = parseFloat(lines[i]);
    t2 = parseFloat(lines[i+1]);
    x1 = Math.cos(t1);
    y1 = Math.sin(t1);
    x2 = Math.cos(t2);
    y2 = Math.sin(t2);
    e += '<path stroke="black" stroke-opacity="0.2" d="M ' + (x1 * radius + offset) + ' ' + (y1 * radius + offset) + ' L ' + (x2 * radius + offset) + ' ' + (y2 * radius + offset) + '"/>';
  }
  e += suffix;
  svg = e;
  draw();
}
var prefix = '<svg height="450" width="450">',
  suffix = "</svg>",
  scale = 0.95,
  offset = 225,
  radius = scale*offset,
  svg = "";
```
```
svg {
  position: absolute;
}
```
```
Paste line specifications here
<br>
<input id="lines" type="textarea" size="60" oninput='drawLines()' />
<br>
<br>
<div id="output"></div>
```
[*Example data generated with method 1.*](http://pastebin.com/raw.php?i=MuW7XU5j)
[*Example data generated with **method 2.***](http://pastebin.com/raw.php?i=B9Amuxrz)
[*Example data generated with method 3.*](http://pastebin.com/raw.php?i=tqPABmt6)
# Leaderboard
**[The first post of the series generates a leaderboard.](https://codegolf.stackexchange.com/q/45302/8478)**
To make sure that your answers show up, please start every answer with a headline, using the following Markdown template:
```
# Language Name, N bytes
```
where `N` is the size of your submission. If you improve your score, you *can* keep old scores in the headline, by striking them through. For instance:
```
# Ruby, <s>104</s> <s>101</s> 96 bytes
```
*(The language is not currently shown, but the snippet does require and parse it, and I may add a by-language leaderboard in the future.)*
  
[Answer]
# IA-32 machine code, 54 bytes
Hexdump of the code:
```
68 00 00 00 4f 0f c7 f0 50 db 04 24 58 d8 34 24
f7 d9 78 f1 d9 c0 dc c8 d9 e8 de e1 d9 fa d9 c9
d9 f3 dc c0 d9 eb de ca d8 c1 dd 1a dd 5a 08 83
c2 10 e2 d1 58 c3
```
It uses a (slightly modified) algorithm that Wikipedia described. In pseudo-code:
```
x = rand_uniform(-1, 1)
y = rand_uniform(-1, 1)
output2 = pi * y
output1 = output2 + 2 * acos(x)
```
I use the range `-1...1` because it's easy to make random numbers in this range: the `rdrand` instruction generates an integer between `-2^31` and `2^31-1`, which can be easily divided by 2^31.
I should have used the range `0...1` for the other random number (x), which is fed into `acos`; however, the negative part is symmetric with the positive part - negative numbers produce chords whose span is greater than pi radians, but for the purpose of illustrating the Bertrand paradox, it doesn't matter.
Since the 80386 (or x87) instruction set doesn't have a dedicated `acos` instruction, I had to express the calculation using only the `atan` instruction:
```
acos(x) = atan(sqrt(1-x^2)/x)
```
Here is the source code that generates the machine code above:
```
__declspec(naked) void __fastcall doit1(int n, std::pair<double, double>* output)
{
    _asm {
        push 0x4f000000;        // [esp] = float representation of 2^32
    myloop:
        rdrand eax;             // generate a random number, -2^31...2^31-1
        push eax;               // convert it
        fild dword ptr [esp];   // to floating-point
        pop eax;                // restore esp
        fdiv dword ptr [esp];   // convert to range -1...1
        neg ecx;
        js myloop;              // do the above 2 times
        // FPU stack contents:  // x           | y
        fld st(0);              // x           | x   | y
        fmul st(0), st;         // x^2         | x   | y
        fld1;                   // 1           | x^2 | x | y
        fsubrp st(1), st;       // 1-x^2       | x   | y
        fsqrt;                  // sqrt(1-x^2) | x   | y
        fxch;                   // x           | sqrt(1-x^2) | y
        fpatan;                 // acos(x)     | y
        fadd st, st(0);         // 2*acos(x)   | y
        fldpi;                  // pi          | 2*acos(x) | y
        fmulp st(2), st;        // 2*acos(x)   | pi*y
        fadd st, st(1);         // output1     | output2
        fstp qword ptr [edx];   // store the numbers
        fstp qword ptr [edx + 8];
        add edx, 16;            // advance the output pointer
        loop myloop;            // loop
        pop eax;                // restore stack pointer
        ret;                    // return
    }
}
```
To generate two random numbers, the code uses a nested loop. To organize the loop, the code takes advantage of the `ecx` register (outer loop counter) being positive. It temporarily negates `ecx`, and then does it again to restore the original value. The `js` instruction repeats the loop when `ecx` is negative (this happens exactly once).
[Answer]
# Pyth, 25 23 22 bytes
A port of rcrmn's C++11 answer. This is my first usage of Pyth, and I've had plenty of fun!
```
VQJ,*y.n0O0.tOZ4,sJ-FJ
```
23-byte version:
```
VQJ*y.n0O0K.tOZ4+JK-JKd
```
Cut a byte by changing the program to use folds+sums and setting J to a tuple, removing K.
Original:
```
VQJ**2.n0O0K.tO0 4+JK-JKd
```
Cut off 2 bytes thanks to @orlp.
Explanation:
```
VQ                         loop as many times as the input number
  J,                       set J to the following tuple expression
    *y.n0O0                2 * .n0 (pi) * O0 (a random number between 0 and 1)
            .tOZ4          .tOZ 4 (acos of OZ (a random number))
                 ,sJ-FJ    print the sum of J and folding J using subtraction in parenthesis, separated by a comma, followed by another newline
```
[Answer]
# Julia, 48 bytes
```
n->(x=2π*rand(n);y=acos(rand(n));hcat(x+y,x-y))
```
This uses the method 2 algorithm, like most of the answers thus far. It creates a lambda function that takes an integer input and returns an n x 2 array. To call it, give it a name, e.g. `f=n->...`.
Ungolfed + explanation:
```
function f(n::Int64)
    # The rand() function returns uniform random numbers using
    # the Mersenne-Twister algorithm
    # Get n random chord angles
    x = 2π*rand(n)
    # Get n random rotations
    y = acos(rand(n))
    # Bind into a 2D array
    hcat(x+y, x-y)
end
```
I really like how the visualizations look, so I'll include one. It's the result of `f(1000)`.
![Circle](https://i.stack.imgur.com/DzxHW.png)
[Answer]
# Pyth, 22 bytes
A port of the C++ answer. I had another 23 byte solution (Now 22!), but it was almost a copy of @kirbyfan64sos's pyth answer with optimizations, so I had to think outside the box a little and creatively (ab)use the fold operator.
```
m,-Fdsdm,y*.nZOZ.tOZ4Q
```
Note that this doesn't work right now because of a bug in the fold operator after the introduction of `reduce2`. I'm putting in a pull request.
```
m             Map    
 ,            Tuple of
  -Fd         Fold subtraction on input
  sd          Fold addition on input
 m      Q     Map over range input
  ,           Tuple           
   y          Double
    *         Product
     .nZ      Pi
     OZ       [0, 1) RNG
  .t  4       Acos
    OZ        [0, 1) RNG
  
```
For refence this was my other solution which works in the same way: [`VQKy*.nZOZJ.tOZ4,+KJ-KJ`](http://pyth.herokuapp.com/?code=VQKy*.nZOZJ.tOZ4%2C%2BKJ-KJ&input=1000)
[Answer]
# IDL, 65 bytes
Obviously this is the same algorithm as @rcrmn, even though I derived it independently before reading their answer.
```
read,n
print,[2,2]#randomu(x,n)*!pi+[-1,1]#acos(randomu(x,n))
end
```
IDL's randomu function uses the Mersenne Twister, which has a period of 219937-1.
EDIT: I ran 1000 chords through the visualizer above, here's a screenshot of the result:
![IDL bertrand](https://i.stack.imgur.com/RG5EQ.png)
[Answer]
# C++11, 214 bytes
```
#include<random>
#include<iostream>
#include<cmath>
int main(){using namespace std;int n;cin>>n;random_device r;uniform_real_distribution<> d;for(;n;--n){float x=2*M_PI*d(r),y=acos(d(r));cout<<x+y<<' '<<x-y<<';';}}
```
So this is a straight out implementation of the *right* algorithm from the wikipedia page. The main problem here in golfing is the oh-so-freaking-long names that the random generator classes have. But, in contrast to good ol' rand, it is at least properly uniform.
Explanation:
```
#include<random>
#include<iostream>
#include<cmath>
int main()
{
    using namespace std;
    int n;
    cin>>n; // Input number
    random_device r; // Get a random number generator
    uniform_real_distribution<> d;   // Get a uniform distribution of 
                                     // floats between 0 and 1
    for(;n;--n)
    {
        float x = 2*M_PI*d(r),       // x: Chosen radius angle
              y = acos(d(r));        // y: Take the distance from the center and 
                                     // apply it an inverse cosine, to get the rotation
        cout<<x+y<<' '<<x-y<<';';    // Print the two numbers: they are the rotation
                                     // of the radius +/- the rotation extracted from
                                     // the distance to the center
    }
}
```
[Answer]
# APL, 46 bytes
```
f←{U←⍵ 2⍴∊{(○2×?0)(¯1○?0)}¨⍳⍵⋄⍉2⍵⍴∊(+/U)(-/U)}
```
My first APL program ever! Surely it can be greatly improved (as my overall understanding of APL is lacking), so any suggestions would be fantastic. This creates a function `f` which takes an integer as input, computes the pairs of chord points using method 2, and prints each pair separated by a newline.
You can [try it online](http://tryapl.org/?a=f%u2190%7BU%u2190%u2375%202%u2374%u220A%7B%28%u25CB2%D7%3F0%29%28%AF1%u25CB%3F0%29%7D%A8%u2373%u2375%u22C4%u23492%u2375%u2374%u220A%28+/U%29%28-/U%29%7D&run)!
Explanation:
```
f←{ ⍝ Create the function f which takes an argument ⍵
    ⍝ Define U to be an ⍵ x 2 array of pairs, where the first
    ⍝ item is 2 times a random uniform float (?0) times pi (○)
    ⍝ and the second is the arccosine (¯1○) of another random
    ⍝ uniform float.
    U ← ⍵ 2 ⍴ ∊{(○2×?0)(¯1○?0)}¨⍳⍵
    ⍝ Create a 2 x ⍵ array filled with U[;1]+U[;2] (+/U) and
    ⍝ U[;1]-U[;2] (-/U). Transpose it into an ⍵ x 2 array of
    ⍝ chord point pairs and return it.
    ⍉ 2 ⍵ ⍴ ∊(+/U)(-/U)
}
```
Note: My previous 19-byte solution was invalid since it returned (x, y) rather than (x+y, x-y). Sadness abounds.
[Answer]
# Java, 114 bytes
```
n->{for(;n-->0;){double a=2*Math.PI*Math.random(),b=Math.acos(Math.random());System.out.println(a+b+" "+(a-b));}};
```
Basic implementation in java. Use as a lambda expression.
[Example usage](http://ideone.com/r0udta)
[Answer]
# Ruby, 72 bytes
My first golf here! I used the same code as everyone, I hope that's okay
```
gets.chomp.to_i.times{puts"#{x=2*Math::PI*rand},#{x+2*Math.acos(rand)}"}
```
[Answer]
# Java, 115 ~~123~~
This is basically the same as most others, but I need a Java score for this hole, so here it goes:
```
void i(int n){for(double x;n-->0;System.out.println(x+2*Math.acos(Math.random())+" "+x))x=2*Math.PI*Math.random();}
```
1000 sample chords can be found at [pastebin](http://pastebin.com/raw.php?i=rpLyzvhj), here are the first five from one run:
```
8.147304676211474 3.772704020731153
8.201346559916786 3.4066194978900106
4.655131524088468 2.887965593766409
4.710707820868578 3.8493686706403984
3.3839198612642423 1.1604092552846672
```
[Answer]
# CJam, ~~24~~ 22 bytes
Similar to other algorithms, here is a version in CJam.
```
{2P*mr_1dmrmC_++]p}ri*
```
An input of 1000 produces a distribution like:
![enter image description here](https://i.snag.gy/QZTkl.jpg)
**How it works**
The algorithm is simply `x = 2 * Pi * rand(); print [x, x + 2 * acos(rand())]`
```
{                 }ri*        e# Run this loop int(input) times
 2P*mr                        e# x := 2 * Pi * rand()
      _                       e# copy x
       1dmr                   e# y := rand()
           mC                 e# z := acos(y)
             _++              e# o := x + z + z
                ]             e# Wrap x and o in an array
                 p            e# Print the array to STDOUT on a new line
```
*Update*: 2 bytes saved thanks to Martin!
[Try it here](http://cjam.aditsu.net/#code=%7B2P*mr_1dmrmC_%2B%2B%5Dp%7Dri*&input=1000)
[Answer]
# Python 3, ~~144~~ 117 bytes
(thanks to Blckknght for the `lambda` pointer)
Using the same method as others:
```
import math as m;from random import random as r;f=lambda n:[(x,x+2*m.acos(r()))for x in(2*m.pi*r()for _ in range(n))]
```
From the Python documentation:
> 
> Python uses the Mersenne Twister as the core generator. It produces 53-bit precision floats and has a period of 219937-1.
> 
> 
> 
## Output
```
>>> f(10)
[(4.8142617617843415, 0.3926824824852387), (3.713855302706769, 1.4014527571152318), (3.0705105305032188, 0.7693910749957577), (1.3583477245841715, 0.9120275474824304), (3.8977143863671646, 1.3309852045392736), (0.9047010644291349, 0.6884780437147916), (3.333698164797664, 1.116653229885653), (3.0027328050516493, 0.6695430795843016), (5.258167740541786, 1.1524381034989306), (4.86435124286598, 1.5676690324824722)]
```
And so on.
## Visualization
![visualization](https://i.stack.imgur.com/S1x3g.png)
[Answer]
# Perl, 60
```
#!perl -l
use Math::Trig;print$x=2*pi*rand,$",$x+2*acos rand for 1..<>
```
[Answer]
# R, ~~60~~ ~~56~~ ~~53~~ 49 bytes
An extra 4 bytes thanks to @JayCe and changing it to a function.
Using the same basic formula as the others. R uses the Mersenne-Twister method by default, but others can be set. Outputs a space separated list.
```
function(n,y=acos(runif(n)))runif(n)*2*pi+c(y,-y)
```
[Try it online!](https://tio.run/##K/qfZvs/rTQvuSQzP08jT6fSNjE5v1ijqDQvM00jT1NTE8bSMtIqyNRO1qjU0a3U/J@mYWig@R8A "R – Try It Online")
[Answer]
# SmileBASIC, 62 bytes
```
INPUT N
FOR I=1TO N
X=PI()*2*RNDF()Y=ACOS(RNDF())?X+Y,X-Y
NEXT
```
]