Random Golf of the Day #2: Numbers from a Normal Distribution

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 2: Numbers from a Normal Distribution

I can't believe this hasn't been done yet! You're to generate random numbers, drawing from a normal distribution. Some rules (the majority of them are probably automatically covered by most submission, but some of them are in place to ensure consistency of results between vastly different languages):

• You should take two non-negative integers as input: a seed S and the amount N of numbers to return. The output should be a list of N floating point numbers, drawn from a normal distribution with mean 0 and variance 1. Whenever your submission is given the same seed S it should produce the same number. In particular, if it is called once with (S, N1) and once with (S, N2), the first min(N1, N2) entries of the two outputs should be identical. In addition, at least 2¹⁶ different values of S should produce different sequences.

• You may use any built-in random number generator that is documented to draw numbers from an (approximately) uniform distribution, provided you can pass S on to it and it supports at least 2¹⁶ different seeds. If you do, the RNG should be able to return at least 2²⁰ different values for any given number you request from it.

• If your available uniform RNG has a smaller range, is not seedable, or supports too few seeds, you must either first build a uniform RNG with a sufficiently large range on top of the built-in one or you must implement your own suitable RNG using the seed. This page may be helpful for that.
• If you don't implement an established algorithm for generating normal distributions, please include a proof of correctness. In either case, the algorithm you choose must yield a theoretically exact normal distribution (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 make use of the full width of the chosen type.
• You must not use any built-in functions directly related to normal distribution or Gaussian integrals, like the Error function or its inverse.

You may write a full program or a function and take input via STDIN, command-line argument, function argument or prompt and produce output via return value or by printing to STDOUT (or closest alternative).

S and N will be non-negative integers, each less than 2²⁰. Output may be in any convenient, unambiguous list or string format.

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.

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

Dyalog APL, 33 bytes

{(.5*⍨¯2×⍟?0)×1○○2×?0}¨⍳⎕⊣⎕rl←1+⎕

Box-Muller:

⎕         ⍝ evaluated input
⎕rl←1+⎕   ⍝ set APL's random seed to 1+⎕ (S)
          ⍝   add 1 because ⎕rl←0 has special meaning: sets the seed randomly
{ }¨⍳N    ⍝ do the thing in braces N times
?0        ⍝ random number 0≤x<1
1○○2×A    ⍝ sin(2πA)
.5*⍨¯2×⍟B ⍝ sqrt(-2lnB)

Answer 2

R, 68 bytes

function(S,N){set.seed(S);sqrt(-2*log(runif(N)))*cos(2*pi*runif(N))}

This uses the runif() function, which generates random deviates from a uniform distribution. The seed for the random number generation is specified using set.seed(), which by default uses the Mersenne-Twister algorithm with a period of 2^19937-1.

The result is an R vector of length N containing the computed standard normal deviates.

This uses the Box-Muller method: For two independent uniform random variables U and V,

Answer 3

Dyalog APL, 42 34

{⎕RL←⍺⋄{(.5*⍨¯2×⍟⍺)×1○⍵×○2}/?⍵2⍴0}

This is a function that takes S as its left argument and N as its right argument.

     5{⎕RL←⍺⋄{(.5*⍨¯2×⍟⍺)×1○⍵×○2}/?⍵2⍴0}10
3.019132549 ¯0.2903143175 ¯0.7353414637 1.421417015 2.327544764 ¯0.00005019747711 ¯0.9582127248 ¯0.2764568462
      ¯0.1602736853 ¯0.9912352616
     5{⎕RL←⍺⋄{(.5*⍨¯2×⍟⍺)×1○⍵×○2}/?⍵2⍴0}20
3.019132549 ¯0.2903143175 ¯0.7353414637 1.421417015 2.327544764 ¯0.00005019747711 ¯0.9582127248 ¯0.2764568462
      ¯0.1602736853 ¯0.9912352616 0.642585109 ¯0.2450019151 ¯0.415034463 0.03481768503 ¯0.4621212815 ¯0.760925979
      0.2592913013 1.884867889 ¯0.9621252731 0.3062560446

It's an implementation of the Box-Muller transform, using Dyalog APL's built-in random operator ?, which by default is a Mersenne twister that returns 64-bit values, which should be enough.

Explanation:

• ⎕RL←⍺: set the random seed to ⍺.
• ?⍵2⍴0: generate ⍵ pairs of random numbers between 0 and 1.
• {...}/: apply the following function to each pair:
  • (.5*⍨¯2×⍟⍺)×1○⍵×○2: calculate the Z0 value (sqrt(-2 ln ⍺)×cos(2π⍵)).

Answer 4

Perl, 67

sub f{srand$_[0];map{cos(atan2(1,1)*rand 8)*sqrt-2*log rand}1..pop}

Box-Muller as in other entries. f takes parameters in order S, N.

Use:

$ perl -le 'sub f{srand$_[0];map{cos(atan2(1,1)*rand 8)*sqrt-2*log rand}1..pop}print for f(5,3)'
-1.59212831801942
0.432167710756345
-0.533673305924252

Answer 5

Java, 164 161 bytes

class B extends java.util.Random{B(int s,int n){super(s);for(;n-->0;System.out.println(Math.sqrt(-2*Math.log(nextDouble()))*Math.cos(2*Math.PI*nextDouble())));}}

This takes in input via function and output via stdout. It uses the Box-Muller method.

Answer 6

Commodore 64 Basic, 76 70 63 bytes

1INPUTS,N:S=R/(-S):F┌I=1TON:?S●(-2*LOG(R/(1)))*S╮(2*π*R/(1)):N─

Since the PETSCII character set contains some symbols not present in Unicode, I've made substitutions: / = SHIFT+N, ┌ = SHIFT+O, ● = SHIFT+Q, ╮ = SHIFT+I, ─ = SHIFT+E

This implements the standard Box-Muller transform to generate the numbers; I picked the sin(x) half of the transform because Commodore 64 Basic has a two-character shortcut for sin(), but not for cos().

Although the manual states otherwise, the value of the argument to RND does matter: if a negative number is passed, the random number generator is not merely re-seeded, it is re-seeded with that number. This makes seeding much simpler: instead of needing to POKE five memory locations, I merely need to make a do-nothing call to RND, which reduces the code from two lines/121 bytes to 1 line/76 bytes.

Edit: Golfed six bytes off by realizing I could combine the two INPUT statements, and that the space after TO was optional.

Edit: Golfed another seven off: Commodore Basic does, in fact, have Pi as a built-in constant, and it's even typeable on a modern keyboard (SHIFT+PgDn, in case you were wondering).

Answer 7

80386 machine code, 72 bytes

Hexdump of the code:

60 8b 7c 24 24 33 ed 8d 75 fb 8d 46 79 f7 e2 f7
f6 8b da b3 7f 0f cb d1 eb 89 1f d9 07 d9 e8 de
e9 33 ee 75 e5 d9 ed d9 c9 d9 f1 dc c0 d9 e0 d9
fa d9 c9 d9 eb de c9 dc c0 d9 ff de c9 d9 1f 83
c7 04 e2 c6 61 c2 04 00

Here is the source code (can be compiled by Visual Studio):

__declspec(naked) void __fastcall doit(int count, unsigned seed, float* output)
{
    _asm {
                                // ecx = count
                                // edx = seed
        // save registers
        pushad;
        mov edi, [esp + 0x24];  // edi = pointer to output
        xor ebp, ebp;           // ebp = 0
        lea esi, [ebp - 5];     // esi = 4294967291 (a prime number)

    myloop:
        // Calculate the next random number
        lea eax, [esi + 121];   // eax = 116
        mul edx;
        div esi;
        mov ebx, edx;

        // Convert it to a float in the range 1...2
        mov bl, 0x7f;
        bswap ebx;
        shr ebx, 1;

        // Convert to range 0...1 and push onto the FPU stack
        mov [edi], ebx;
        fld dword ptr [edi];
        fld1;
        fsubp st(1), st;

        // Make 2 such random numbers
        xor ebp, esi;
        jnz myloop;

        // Calculate sqrt(-2*ln(x))
        fldln2;
        fxch;
        fyl2x;
        fadd st, st(0);
        fchs;
        fsqrt;

        // Calculate cos(2*pi*y)
        fxch st(1);
        fldpi;
        fmul;
        fadd st, st(0);
        fcos;

        // Calculate and write output
        fmulp st(1), st;
        fstp dword ptr [edi];
        add edi, 4;

        // Repeat
        loop myloop

        // Return
        popad;
        ret 4;
    }
}

---

Here I use a a Lehmer random number generator. It uses the following algorithm:

x(k+1) = 116 * x(k) mod 4294967291

Here 4294967291 is a big (2^32-5) prime number, and 116 is a small (less than 128; see below) number that is its primitive root. I chose a primitive root that has a more-or-less random distribution of zeros and ones in binary representation (01110100). This RNG has the maximum possible period of 4294967290, if the seed is nonzero.

---

The relatively small numbers I used here (116 and 4294967291, which can be represented also as -5) let me take advantage of the lea instruction encoding:

8d 46 79     lea eax, [esi+121]

It is assembled to 3 bytes if the numbers can fit into 1 byte.

---

The multiplication and division use edx and eax as their working registers, which is why I made seed the second parameter to the function (fastcall calling convention uses edx to pass the second parameter). In addition, the first parameter is passed in ecx, which is a good place to hold a counter: a loop can be organized in 1 instruction!

e2 c6        loop myloop

---

To convert an integer into a floating-point number, I exploited the representation of single-precision floating-point numbers: if I set the high 9 bits (exponent) to the bit pattern 001111111, and leave the 23 low bits random, I'll get a random number in the range 1...2. I took the idea from here. To set the high 9 bits, I used some bit-fiddling on ebx:

mov ebx, edx;    xxxxxxxx|yyyyyyyy|zzzzzzzz|aaaaaaaa
mov bl, 0x7f;    xxxxxxxx|yyyyyyyy|zzzzzzzz|01111111
bswap ebx;       01111111|zzzzzzzz|yyyyyyyy|xxxxxxxx
shr ebx, 1;      00111111|1zzzzzzz|zyyyyyyy|yxxxxxxx

---

To generate two random numbers, I used a nested loop of 2 iterations. I organized it with xor:

xor ebp, esi;    first time, the result is -5
jnz myloop;      second time, the result is 0 - exit loop

---

The floating-point code implements the Box-Muller transform.

Answer 8

Haskell, 118  144 

import System.Random;h z=let(u,r)=random z in(cos$2*pi*fst(random r)::Float)*sqrt(-2*log u):h r;g=(.(h.mkStdGen)).take

Example usage:

*Main> g 3 0x6AE4A92CAFA8A742
[0.50378895,-0.20593005,-0.16684927]
*Main> g 6 0x6AE4A92CAFA8A742
[0.50378895,-0.20593005,-0.16684927,1.1229043,-0.10026576,0.4279402]
*Main> g 6 0xE09B1088DF461F7D
[2.5723906,-0.5177805,-1.3535261,0.7400385,3.5619608e-3,-8.246434e-2]

The return type of random is constrained to Float, which makes random generate a uniform float in [0, 1). From then on it's a simlpe box-muller formula with some pointless magic for list generation.

Answer 9

Golflua, 63 70

Golflua info and instructions.

\g(n,s)`_ENV,b=M,{}rs(s)~@i=1,n b[i]=q(l(r()^-2))*c(r()*pi)$~b$

Returns a table containing the values. In the example I'm using ~T.u( ), which is the same as return table.unpack( ) in lua.

> ~T.u(g(3,1234567))
0.89302672974232 0.36330401643578 -0.64762161593981
> ~T.u(g(5,1234567))
0.89302672974232 0.36330401643578 -0.64762161593981 -0.70654636393063 -0.65662878785425
> ~T.u(g(5,7654321))
0.3867923683064 -0.31758512485963 -0.58059120409317 1.2079459300077 1.1500121921242

A lot of chars were saved by setting the function's environment to M (aka math).

Answer 10

SAS, 108

I already posted an answer in R that's shorter than this, but there are very few SAS answers on PPCG, so why not add another?

%macro f(s,n);data;do i=1 to &n;x=sqrt(-2*log(ranuni(&s)))*cos(8*atan2(1,1)*ranuni(&s));put x;end;run;%mend;

With some white space:

%macro f(s, n);
    data;
        do i = 1 to &n;
            x = sqrt(-2 * log(ranuni(&s))) * cos(8 * atan2(1, 1) * ranuni(&s));
            put x;
        end;
    run;
%mend;

This defines a macro which can be called like %f(5, 3). The macro executes a data step which loops through the integers 1 to N, and at each iteration it computes a random normal deviate using Box-Muller and prints it to the log using the put statement.

SAS has no built-in for pi, so the best we can do is approximate it with arctangent.

The ranuni() function (which is deprecated but requires a couple fewer characters than the newer function) returns a random number from a uniform distribution. The SAS documentation doesn't give a whole lot of detail on the RNG implementation other than it has a period of 2^31-2.

In SAS macros, macro variables are referenced with a preceding & and resolve to their values at run time.

As you've probably witnessed, SAS is rarely an actual contender in a code-golf contest.

Answer 11

Java, 193 bytes

While this doesn't beat the current Java leader, I decided to post anyway to show a different method of calculation. It's a golfed version of OpenJDK's nextGaussian().

class N extends java.util.Random{N(int s,int n){super(s);for(float a,v;n-->0;System.out.println(v*Math.sqrt(-2*Math.log(a)/a)))for(a=0;a>=1|a==0;a=v*v+(v=2*nextFloat()-1)*v)v=2*nextFloat()-1;}}

With line breaks:

class N extends java.util.Random{
    N(int s,int n){
        super(s);
        for(float a,v;
            n-->0;
            System.out.println(v*Math.sqrt(-2*Math.log(a)/a)))
                for(a=0;
                    a>=1|a==0;
                    a=v*v+(v=2*nextFloat()-1)*v)v=2*nextFloat()-1;
    }
}

Answer 12

T-SQL, 155 bytes

CREATE PROC R(@S BIGINT,@N INT)AS
DECLARE @ INT=0,@K INT=8388607WHILE @<@N
BEGIN
SELECT SQRT(-2*LOG(RAND(@S*@%@K)))*COS(2*PI()*RAND(@S*9*@%@K))SET @+=1
END

Use with EXEC R S,N because there is no STD_IN in T-SQL where S and N are the seed and N respectively. S will produce "random" (RAND(seed) is a really bad random number implementation) sequences when S>2^16 (possibly before that, but I won't guarantee it). Uses Box-Muller like most solutions so far. 8388607 is 2^23-1, which should hopefully generate 2^20 different values.

Answer 13

Powershell, 164 bytes

Param($s,$n)$q=2147483647
$a=GET-RANDOM -SETSEED $s
FOR(;$n---gt0;){$a=GET-RANDOM
$b=GET-RANDOM
[math]::SQRT(-2*[math]::LOG($a/$q))*[math]::COS(2*[math]::PI*$b/$q)}

Same as most answers with Box-Muller. Not very experienced with Powershell, so any help golfing would be appreciated.

Answer 14

Ruby, 72 bytes

->n,s{include Math;srand s;n.times{p sqrt(-2*log(rand))*sin(2*rand*PI)}}

Input (as lambda function):

f.(6, 12353405)

Output:

-1.1565142460805273
0.9352802655317097
1.3566720571574993
-0.9683973210257978
0.9851210877202192
0.14709635752306677

PS: I would like to know if this can be golfed further. I'm just a beginner.

Answer 15

Matlab, 77

The first input should be n, the second one s.

a=input('');
rand('seed',a(2));
for i=1:a;
    (-2*log(rand))^.5*cos(2*pi*rand)
end

Answer 16

Octave, 91 96 88 bytes

function r=n(s,n)rand("seed",s);x=rand(2,n);r=cos(2*pi*x(1,:)).*sqrt(-2*log(x(2,:)));end

Or, with whitespace:

function r=n(s,n)
  rand("seed",s);
  x=rand(2,n);
  r=cos(2*pi*x(1,:)).*sqrt(-2*log(x(2,:)));
end

Set the seed up front, and use the Box-Mueller method.

NB: Octave allows for generation of arrays of random numbers, and can use standard operations on these arrays which produce array outputs. The .* operator is element-by-element multiplication of two arrays to produce the result.

Answer 17

Pyth, 32 bytes

No Python being used now in super-quotes because of the new functions Pyth has now. Yet another Box-Mueller.

 .xvzVQ*@_y.lOZ2.71 2.ty*3.14OZ1

That space in the beginning is important.

.xvz             Seed RNG with evaluated input
VQ               For N in second input
*                Multiplication
 @       2       Square root
   _y            Times negative 2
    .l )         Natural log
     OZ          Of random float (RNG with zero give [0, 1) float)
 .t       1      Cosine from list of trig functions
  y              Double
   *             Multiplication
    .nZ          Pi from constants list
    OZ           Random Float

The seeding doesn't seem to work in the online interpreter, but it works fine in the local version. The online interpreter seems to be fixed, so here's a permalink: permalink

Answer 18

STATA, 85 bytes

di _r(s)_r(n)
set se $s
set ob $n
g a=sqrt(-2*ln(runiform()))*cos(2*runiform()*_pi)
l

Takes input via standard in (first number is S, then N). Sets the seed to S. Sets the number of observations to N. Makes a variable and sets its value to the Box Muller transform value (thanks to @Alex for showing it). Then lists all of the observations in a table with column header a and observation numbers next to them. If these are not okay, let me know, and I can remove headers and/or observation numbers.

Answer 19

R, 89 Bytes

I know R has been done before but I wanted to show a different approach than the Box-Muller that everyone else used. My solution uses the Central Limit Theorem .

f=function(S,N){set.seed(S);a=1000;for(i in(1:N)){print(sqrt(12/a)*(sum(runif(a))-a/2))}}

Answer 20

TI-Basic, 74 bytes

Prompt S,N:S→rand:For(X,1,N:0→A:0→V:0→W:While A≥1 or A=0:2rand-1→V:2rand-1→W:V²+W²→A:End:Disp VW√(Aֿ¹-2log(A:End

1      1111111   11   1111111111111111111     1111   111111   1111111   11111111111111  11    111111111   111

The ¹ is actually the inverse operator.

Answer 21

Perl, 150 108 107 bytes

This uses the Marsaglia Polar Method. Called with f(S, N).

Moved the assignment of $a into the calculation of $c.

107:

sub f{srand$_[0];map{do{$c=($a=-1+2*rand)**2+(-1+2*rand)**2}until$c<1;print$a*sqrt(-2*log($c)/$c)}1..$_[1]}

Removed spare number storage and the definition of $b.

108:

sub f{srand$_[0];map{do{$a=-1+2*rand,$c=$a**2+(-1+2*rand)**2}until$c<1;print$a*sqrt(-2*log($c)/$c)}1..$_[1]}

150:

sub f{srand$_[0];map{$h?$h=!print$s:do{do{$a=-1+2*rand,$b=-1+2*rand,$c=$a*$a+$b*$b}until$c<1;$d=sqrt(-2*log($c)/$c);$s=$b*$d;$h=print$a*$d;}}1..$_[1]}

Answer 22

Swift, 144 142

Nothing clever, just seeing how Swift works.

import Foundation;func r(s:UInt32,n:Int){srand(s);for i in 0..<n{println(sqrt(-2*log(Double(rand())/0xffffffff))*sin(2*Double(rand())*M_PI))}}

I was hoping I could use (0...n).map{} but it the compiler doesn't seem to recognise map{} unless you use a parameter.

Answer 23

Haskell, 97 bytes

import System.Random
h(a:b:c)=sqrt(-2*log a::Float)*cos(2*pi*b):h c
f a=take a.h.randoms.mkStdGen

Try it online!

Just your basic Box-Muller transformation, on an infinite list of random numbers.

Answer 24

Python 3, 118 bytes

from random import*
from math import*
def f(s,n):seed(s);return[sqrt(-2*log(random()))*cos(tau*random())for _ in[0]*n]

Try it online!

Answer 25

SmileBASIC, 81 bytes

Well, now that I've answered the first question, I have to do all the rest...

Generating random numbers is cheap, but seeding the RNG uses the longest builtin function in the language, RANDOMIZE.

DEF N S,N
RANDOMIZE.,S
FOR I=1TO N?SQR(-2*LOG(RNDF()))*COS(PI()*2*RNDF())NEXT
END

Maybe there's some way to optimize the formula. I don't see how it's required to use two RNG calls.