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

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 216 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 216 different seeds. If you do, the RNG should be able to return at least 220 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 220. 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.

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

• Wait, are we allowed to use non-integer-range RNG's? Feb 6 '15 at 18:37
• P.S. 2 existing answers seem to use [0, 1) floating point RNG's, is this allowed? Feb 6 '15 at 18:44
• @mniip Yes, floating-point RNG's are allowed, provided they are uniform, seedable and able to return the required number of distinct floats. Feb 6 '15 at 18:46

# Dyalog APL, 33 bytes

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

⎕         ⍝ 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)

• I'm pretty sure this can't beaten by any other language. Feb 8 '15 at 3:53
• This doesn't meet this rule: "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." Feb 8 '15 at 8:03
• @marinus Thanks, fixed. I also changed ⎕rl to be S+1 because ⎕rl←0 has special meaning.
– ngn
Feb 8 '15 at 10:29
• You probably don't actually need the +1, all it says is that you need to support at least 2^16 different values. So working correctly in the range [1..2^16] should be OK. Feb 8 '15 at 18:40
• S=0 would make the computation not repeatable, which violates the rule you quoted above.
– ngn
Feb 8 '15 at 20:26

# 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, • If that's valid syntax in R, you can leave out the f= (the function doesn't necessarily need to be named, if unnamed functions are a thing in your language). Feb 6 '15 at 19:04
• @MartinBüttner: I appreciate the suggestion but to my knowledge R wouldn't know what to do with an unnamed function. Feb 6 '15 at 19:12
• I always get an error message Error: unexpected '}' in "f=fu... besides, are you sure you get the same first numbers if you call f(0,1) and f(0,2)? Feb 6 '15 at 23:33

## 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π⍵)).
• In v14.0 ?0 returns a floating point number between 0 and 1.
– ngn
Feb 7 '15 at 21:22

sub f{srand$_;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$_;map{cos(atan2(1,1)*rand 8)*sqrt-2*log rand}1..pop}print for f(5,3)' -1.59212831801942 0.432167710756345 -0.533673305924252  # 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. • s=0;s++<n; -> ;n-->0; ? Feb 6 '15 at 19:13 • @Geobits Looks like a lambda :D Feb 6 '15 at 20:12 # Commodore 64 Basic, 7670 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). # 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. # 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.

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

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

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

• +1 for using Marsaglia (or rather for not using the straight-forward Box-Muller) ;) Feb 6 '15 at 20:31
• Can this not be a lambda? Something like: (s,n)->{java.util.Random r=new java.util.Random(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*r.nextFloat()-1)*v)v=2*r.nextFloat()-1;} Feb 6 '15 at 23:01
• @Quincunx I could, for one byte. But I don't like hiding my function declarations in uncounted code, whatever the current consensus on meta about it. It's worth the one byte to me ;) Feb 7 '15 at 2:05

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

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

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

• @MartinBüttner I think I've been using too much C these days. Totally forgot about. Feb 7 '15 at 14:00

# 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


# Octave, 9196 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.

• I think this does not meet the conditions, if you call n(0,1) and n(0,2) you get different first numbers, do you? Feb 6 '15 at 23:29
• Crap, you are correct. I've fixed it but it cost me 5 bytes... Feb 9 '15 at 3:23

# 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

• This uses a feature of Pyth (.nZ) that wasn't implemented, when the question was asked. (It was actually implemented today.) Therefore this answer should not be part of the competition (meta.codegolf.stackexchange.com/questions/4867/…). Apr 7 '15 at 20:01
• K, I'll got back to the 32 char solution Apr 7 '15 at 20:03
• Yeah, that would be better. You can still show off you new solution in a seperate section of you answer. But the code your competing should be the one working with the old Pyth. Apr 7 '15 at 20:05
• Btw, I don't think that the 32 solution should be valid either. Since it uses the random seed initialized, which was only added about 5 days ago. Apr 7 '15 at 20:11

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

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

• I'm afraid, the central limit theorem doesn't satisfy "the algorithm you choose must yield a theoretically exact normal distribution". No matter how many uniform variables you add up, as long as the sum is finite, the normal distribution will always be approximate. (While the central limit theorem is a nice idea, I had to rule it out precisely because it's not clear which value should be used for a in your code such that the result is "fair".) Feb 7 '15 at 2:16
• It was worth a shot ;) Feb 7 '15 at 2:18

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

# Perl, 150108 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$_;map{do{$c=($a=-1+2*rand)**2+(-1+2*rand)**2}until$c<1;print$a*sqrt(-2*log($c)/$c)}1..$_}


Removed spare number storage and the definition of $b. 108: sub f{srand$_;map{do{$a=-1+2*rand,$c=$a**2+(-1+2*rand)**2}until$c<1;print$a*sqrt(-2*log($c)/$c)}1..$_}


150:

sub f{srand$_;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..$_}


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

• of course...? it's forEach if you don't want a return value, and I'm pretty sure the _ in is compulsory Jan 22 '19 at 5:14
• what's the /0xffffffff for btw Jan 22 '19 at 5:18

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.

# 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*n]


Try it online!

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

• It is required to have two independent samples for the Box-Muller transformation Jan 22 '19 at 4:56