# Random Golf of the Day #4: The Bertrand Paradox

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

The first post of the series generates a leaderboard.

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

# 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];

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

• I like this answer for using IA32 assembly! Just to say: this is not strictly 386 machine code as 80386 obviously don't have rdrand instruction nor necessary a FP coprocessor – user5572685 Apr 29 '15 at 13:31
• Yeah, IA32 is a better name. Not specific enough but probably more correct than 80386. – anatolyg Apr 29 '15 at 15:24

# Pyth, 2523 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

• Pyth tips: *2_ is the same as y_. Variable Z is initially 0, so you can remove the space in .tO0 4 by writing .tOZ4. – orlp Apr 28 '15 at 22:22
• I'm counting 25 bytes... – Maltysen Apr 28 '15 at 22:31
• Also, you can format the output better with ,+JK-JK – Maltysen Apr 28 '15 at 22:31
• @Maltysen Both fixed. Thanks! – kirbyfan64sos Apr 28 '15 at 22:34
• @orlp I had no clue about y and forgot about Z. Fixed; thanks! – kirbyfan64sos Apr 28 '15 at 22:34

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

# 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
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

• You mispelled my user name... :( – kirbyfan64sos Apr 28 '15 at 22:35
• @kirbyfan64sos derp. Sorry ;) – Maltysen Apr 28 '15 at 22:39

# 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:

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

• That factor M_PI_2 looks suspicious. I think it should be 1 instead. – anatolyg Apr 28 '15 at 19:13
• Yes, completely right, going to fix it now! Thanks a lot! – rorlork Apr 28 '15 at 21:53

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

# 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

• Can't you reduce the size by storing Math somewhere? Or something? (Im not a Java programmer) – Ismael Miguel Apr 28 '15 at 17:18
• @IsmaelMiguel That would cost an extra 2 characters. – TheNumberOne Apr 28 '15 at 17:24
• Sorry :/ It's tempting to try to reduce the number of times that Math shows. What does the meta say about using a code to generate another code to solve the issue? – Ismael Miguel Apr 28 '15 at 17:32
• @IsmaelMiguel That's fair game, although I'll be surprised if you're actually better at metagolfing than at golfing. – Martin Ender Apr 28 '15 at 17:46

# 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)}"}


# 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


# 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

# 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

• You can save about 20 bytes if you use a lambda for the function and return a list comprehension (with an inner generator expression): f=lambda n:[(x,x+2*m.acos(r()))for x in(2*m.pi*r()for _ in range(n))] – Blckknght May 2 '15 at 7:51
• Ah, I had a lambda originally. I guess I didn't think about doubling up on list comprehensions. Thanks! @Blckknght – Zach Gates May 2 '15 at 8:02
• Can be shortened to 109 bytes by fiddling with the imports: tio.run/#python2 – Triggernometry Aug 8 '18 at 21:00

# Perl, 60

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


# R, 605653 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!

• Hi Micky, you can Save 4 bytes by making it a function and not defining x. – JayCe Aug 8 '18 at 17:23
• @JayCe that's a lot better thanks – MickyT Aug 8 '18 at 20:35

# SmileBASIC, 62 bytes

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