Random Golf of the Day #4: The Bertrand Paradox

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.

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.

Hole 4: The Bertrand Paradox

The Bertrand paradox 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 2²⁰ 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. This page 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.

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.

Example data generated with method 2.

Example data generated with method 3.

Leaderboard

The first post of the series generates a leaderboard.

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 1

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 2

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 3

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

Answer 4

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

Answer 5

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 2¹⁹⁹³⁷-1.

EDIT: I ran 1000 chords through the visualizer above, here's a screenshot of the result:

Answer 6

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 7

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!

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 8

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

Answer 9

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 10

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, 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 11

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:

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

Answer 12

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 2¹⁹⁹³⁷-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

Answer 13

Perl, 60

#!perl -l
use Math::Trig;print$x=2*pi*rand,$",$x+2*acos rand for 1..<>

Answer 14

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!

Answer 15

SmileBASIC, 62 bytes

INPUT N
FOR I=1TO N
X=PI()*2*RNDF()Y=ACOS(RNDF())?X+Y,X-Y
NEXT