# Random numbers

(Redirected from Random Numbers)
Random numbers
You are encouraged to solve this task according to the task description, using any language you may know.

Generate a collection filled with   1000   normally distributed random (or pseudo-random) numbers with a mean of   1.0   and a   standard deviation   of   0.5

Many libraries only generate uniformly distributed random numbers. If so, you may use one of these algorithms.

```with Ada.Numerics;                       use Ada.Numerics;

procedure Normal_Random is
function Normal_Distribution
(  Seed  : Generator;
Mu    : Float := 1.0;
Sigma : Float := 0.5
)  return Float is
begin
return
Mu + (Sigma * Sqrt (-2.0 * Log (Random (Seed), 10.0)) * Cos (2.0 * Pi * Random (Seed)));
end Normal_Distribution;

Seed         : Generator;
Distribution : array (1..1_000) of Float;
begin
Reset (Seed);
for I in Distribution'Range loop
Distribution (I) := Normal_Distribution (Seed);
end loop;
end Normal_Random;
```

## ALGOL 68

Translation of: C
Works with: ALGOL 68 version Revision 1 - no extensions to language used
Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny
```PROC random normal = REAL:  # normal distribution, centered on 0, std dev 1 #
(
sqrt(-2*log(random)) * cos(2*pi*random)
);

test:(
[1000]REAL rands;
FOR i TO UPB rands DO
rands[i] := 1 + random normal/2
OD;
INT limit=10;
printf((\$"("n(limit-1)(-d.6d",")-d.5d" ... )"\$, rands[:limit]))
)```
Output:
```( 0.693461, 0.948424, 0.482261, 1.045939, 0.890818, 1.467935, 0.604153, 0.804811, 0.690227, 0.83462 ... )
```

## Arturo

```rnd: function []-> (random 0 10000)//10000

rands: map 1..1000 'x [
1 + (sqrt neg 2 * ln rnd) * (cos 2 * pi * rnd)
]

print rands
```
Output:
`0.6219537961087694 1.279396486161406 1.619019280815647 2.119538294228789 -0.1598383851981044 2.67797211803156 0.9304703037226587 1.629254364659528 -0.4171704717398712 0.9082342931486092 -0.5929704390625219 2.117000897984871 -0.1981633787460266 0.01132471973856153 2.102359263212924 0.2408823232884222 2.046195035792376 0.6374831627030295 0.000839808324124558 1.117061838266626 0.7413355299469649 0.4485598815755762 2.999434800016997 2.560580541932842 1.703197984879731 2.889159248353575 -1.800157205708138 1.756810020187321 0.7136708180852145 0.5929151678321705 0.332519993787973 2.660212054362758 0.5835660585480075 0.8527946892567934 1.640573993747053 0.09471843345263908 1.051402997891346 1.116149156137905 -0.7400139019343499 1.782572831979232 2.531779039786426 0.5240268064639871 0.07099232630526586 -0.854892656700071 1.54381929430469 -0.4416899008614745 0.4274356035015117 0.7350027625573482 2.153583935076981 1.461215281535983 -1.041723064151266 2.338060763553139 -0.1492967916030414 0.3799517724040202 0.4577924541353815 0.673317567666373 -2.27731583876462 1.28480355806061 -0.6925811023772748 -0.2642224122781984 0.6590513830891744 2.55537133425143 1.67933335469247 0.8659013395355968 -1.211026941441126 0.9524579534222226 -0.1931750631835656 -0.5119479869237693 0.1814749003063878 3.03139579963414 ...`

## AutoHotkey

contributed by Laszlo on the ahk forum

```Loop 40
R .= RandN(1,0.5) "`n"  ; mean = 1.0, standard deviation = 0.5
MsgBox %R%

RandN(m,s) { ; Normally distributed random numbers of mean = m, std.dev = s by Box-Muller method
Static i, Y
If (i := !i) { ; every other call
Random U, 0, 1.0
Random V, 0, 6.2831853071795862
U := sqrt(-2*ln(U))*s
Y := m + U*sin(V)
Return m + U*cos(V)
}
Return Y
}
```

## Avail

```Method "U(_,_)" is
[
lower : number,
upper : number
|
divisor ::= ((1<<32)) ÷ (upper - lower)→double;
map a pRNG through [i : integer | (i ÷ divisor) + lower]
];

Method "a Marsaglia polar sampler" is
[
generator for
[
yield : [double]→⊤
|
source ::= U(-1, 1);
Repeat [
x ::= take 1 from source[1];
y ::= take 1 from source[1];
s ::= x^2 + y^2;
If 0 < s < 1 then
[
factor ::= ((-2 × ln s) ÷ s) ^ 0.5;
yield(x × factor);
yield(y × factor);
];
]
]
];

// the default distribution has mean 0 and std dev 1.0, so we scale the values
sampler ::= map a Marsaglia polar sampler through [d : double | d ÷ 2.0 + 1.0];
values ::= take 1000 from sampler;```

## AWK

One-liner:

```\$ awk 'func r(){return sqrt(-2*log(rand()))*cos(6.2831853*rand())}BEGIN{for(i=0;i<1000;i++)s=s" "1+0.5*r();print s}'
```

```function r() {
return sqrt( -2*log( rand() ) ) * cos(6.2831853*rand() )
}

BEGIN {
n=1000
for(i=0;i<n;i++) {
x = 1 + 0.5*r()
s = s" "x
}
print s
}
```
Output:

first few values only

`0.783753 1.16682 1.17989 1.14975 1.34784 0.29296 0.979227 1.04402 0.567835 1.58812 0.465559 1.27186 0.324533 0.725827 -0.0626549 0.632273 1.0145 1.3387 0.861667 1.04147 1.2576 1.02497 0.58453 0.9619 1.26902 0.851048 -0.126259 0.863256 `

...

## BASIC

### ANSI BASIC

Translation of: FreeBASIC
Works with: Decimal BASIC
```100 REM Random numbers
110 RANDOMIZE
120 DEF RandomNormal = COS(2 * PI * RND) * SQR(-2 * LOG(RND))
130 DIM R(0 TO 999)
140 LET Sum = 0
150 FOR I = 0 TO 999
160    LET R(I) = 1 + RandomNormal / 2
170    LET Sum = Sum + R(I)
180 NEXT I
190 LET Mean = Sum / 1000
200 LET Sum = 0
210 FOR I = 0 TO 999
220    LET Sum = Sum + (R(I) - Mean) ^ 2
230 NEXT I
240 LET SD = SQR(Sum / 1000)
250 PRINT "Mean is              "; Mean
260 PRINT "Standard Deviation is"; SD
270 PRINT
280 END
```
Output:

Two runs.

```Mean is               1.00216454061435
Standard Deviation is .504515904812839
```
```Mean is               .995781408878628
Standard Deviation is .499307289407576
```

### Applesoft BASIC

The Commodore BASIC code works in Applesoft BASIC.

### BASIC256

Translation of: FreeBASIC
```# Generates normally distributed random numbers with mean 0 and standard deviation 1
function randomNormal()
return cos(2.0 * pi * rand) * sqr(-2.0 * log(rand))
end function

dim r(1000)
sum = 0.0
# Generate 1000 normally distributed random numbers
# with mean 1 and standard deviation 0.5
# and calculate their sum
for i = 0 to 999
r[i] = 1.0 + randomNormal() / 2.0
sum += r[i]
next i

mean = sum / 1000.0
sum = 0.0
# Now calculate their standard deviation
for i = 0 to 999
sum += (r[i] - mean) ^ 2.0
next i
sd = sqr(sum/1000.0)

print "Mean is               "; mean
print "Standard Deviation is "; sd
end```
Output:
```Mean is               1.002092
Standard Deviation is 0.4838570687```

### BBC BASIC

```      DIM array(999)
FOR number% = 0 TO 999
array(number%) = 1.0 + 0.5 * SQR(-2*LN(RND(1))) * COS(2*PI*RND(1))
NEXT

mean = SUM(array()) / (DIM(array(),1) + 1)
array() -= mean
stdev = MOD(array()) / SQR(DIM(array(),1) + 1)

PRINT "Mean = " ; mean
PRINT "Standard deviation = " ; stdev
```
Output:
```Mean = 1.01848064
Standard deviation = 0.503551814```

### Chipmunk Basic

Translation of: FreeBASIC
```10 ' Random numbers
20 randomize timer
30 dim r(999)
40 sum = 0
50 for i = 0 to 999
60   r(i) = 1+randomnormal()/2
70   sum = sum+r(i)
80 next
90 mean = sum/1000
100 sum = 0
110 for i = 0 to 999
120   sum = sum+(r(i)-mean)^2
130 next
140 sd = sqr(sum/1000)
150 print "Mean is               ";mean
160 print "Standard Deviation is ";sd
170 print
180 end

500 sub randomnormal()
510   randomnormal = cos(2*pi*rnd(1))*sqr(-2*log(rnd(1)))
520 end sub
```
Output:

Two runs.

```Mean is               1.007087
Standard Deviation is 0.496848
```
```Mean is               0.9781
Standard Deviation is 0.508147
```

### Commodore BASIC

```10 DIM AR(999): DIM DE(999)
20 FOR N = 0 TO 999
30 AR(N)= 0 + SQR(-1.3*LOG(RND(1))) * COS(1.2*PI*RND(1))
40 NEXT N
50 :
60 REM SUM
70 LET SU = 0
80 FOR N = 0 TO 999
90 LET SU = SU + AR(N)
100 NEXT N
110 :
120 REM MEAN
130 LET ME= 0
140 LET ME = SU/1000
150 :
160 REM DEVIATION
170 FOR N = 0 TO 999
180 T = AR(N)-ME: REM SUBTRACT MEAN FROM NUMBER
190 T = T * T: REM SQUARE THE RESULT
200 DE(N) = T : REM STORE IN DEVIATION ARRAY
210 NEXT N
220 LET DS=0: REM SUM OF DEVIATION ARRAY
230 FOR N = 0 TO 999
240 LET DS = DS + DE(N)
250 NEXT N
260 LET DM=0: REM MEAN OF DEVIATION ARRAY
270 LET DM = DS / 1000
280 LET DE = 0:
290 LET DE = SQR(DM)
300 :
310 PRINT "MEAN = "ME
320 PRINT "STANDARD DEVIATION ="DE
330 END
```

### FreeBASIC

```' FB 1.05.0 Win64

Const pi As Double = 3.141592653589793
Randomize

' Generates normally distributed random numbers with mean 0 and standard deviation 1
Function randomNormal() As Double
Return Cos(2.0 * pi * Rnd) * Sqr(-2.0 * Log(Rnd))
End Function

Dim r(0 To 999) As Double
Dim sum As Double = 0.0

' Generate 1000 normally distributed random numbers
' with mean 1 and standard deviation 0.5
' and calculate their sum
For i As Integer = 0 To 999
r(i) = 1.0 + randomNormal/2.0
sum += r(i)
Next

Dim mean As Double = sum / 1000.0

Dim sd As Double
sum = 0.0
' Now calculate their standard deviation
For i As Integer = 0 To 999
sum += (r(i) - mean) ^ 2.0
Next
sd  = Sqr(sum/1000.0)

Print "Mean is              "; mean
Print "Standard Deviation is"; sd
Print
Print "Press any key to quit"
Sleep```

Sample result:

Output:
```Mean is               1.000763573902885
Standard Deviation is 0.500653063426955
```

### FutureBasic

Note: To generate the random number, rather than using FB's native "rnd" function, this code wraps C code into the RandomZeroToOne function.

```window 1

local fn RandomZeroToOne as double
double result
cln result = (double)( (rand() % 100000 ) * 0.00001 );
end fn = result

local fn RandomGaussian as double
double r = fn RandomZeroToOne
end fn = 1 + .5 * ( sqr( -2 * log(r) ) * cos( 2 * pi * r ) )

long i
double mean, std, a(1000)

for i = 1 to 1000
a(i) = fn RandomGaussian
mean += a(i)
next
mean = mean / 1000

for i = 1 to 1000
std += ( a(i) - mean )^2
next
std = std / 1000

print "           Average:  "; mean
print "Standard Deviation:  "; std

HandleEvents```
Output:
```           Average: 1.053724951604593
Standard Deviation: 0.2897370762627166
```

### GW-BASIC

The Commodore BASIC code works in GW-BASIC.

### Liberty BASIC

```dim a(1000)
mean =1
sd =0.5
for i = 1 to 1000   '   throw 1000 normal variates
a( i)  =mean +sd *( sqr( -2 * log( rnd( 0))) * cos( 2 * pi * rnd( 0)))
next i```

### Minimal BASIC

Translation of: FreeBASIC
```10 REM Random numbers
20 LET P = 4*ATN(1)
30 RANDOMIZE
40 DEF FNN = COS(2*P*RND)*SQR(-2*LOG(RND))
50 DIM R(999)
60 LET S = 0
70 FOR I = 0 TO 999
80 LET R(I) = 1+FNN/2
90 LET S = S+R(I)
100 NEXT I
110 LET M = S/1000
120 LET S = 0
130 FOR I = 0 TO 999
140 LET S = S+(R(I)-M)^2
150 NEXT I
160 LET D = SQR(S/1000)
170 PRINT "Mean is              "; M
180 PRINT "Standard Deviation is"; D
190 PRINT
200 END
```

### PureBasic

```Procedure.f RandomNormal()
; This procedure can return any real number.
Protected.f x1, x2

; random numbers from the open interval ]0, 1[
x1 = (Random(999998)+1) / 1000000       ; must be > 0 because of Log(x1)
x2 = (Random(999998)+1) / 1000000

ProcedureReturn Sqr(-2*Log(x1)) * Cos(2*#PI*x2)
EndProcedure

Define i, n=1000

Dim a.q(n-1)
For i = 0 To n-1
a(i) = 1 + 0.5 * RandomNormal()
Next```

### QuickBASIC

Works with: QuickBasic version 4.5
```RANDOMIZE TIMER 'seeds random number generator with the system time
pi = 3.141592653589793#
DIM a(1 TO 1000) AS DOUBLE
CLS
FOR i = 1 TO 1000
a(i) = 1 + SQR(-2 * LOG(RND)) * COS(2 * pi * RND)
NEXT i
```

### Run BASIC

```dim a(1000)
pi = 22/7
for i = 1 to 1000
a( i)  = 1 + .5 * (sqr(-2 * log(rnd(0))) * cos(2 * pi * rnd(0)))
next i```

### TI-83 BASIC

Built-in function: randNorm()

``` randNorm(1,.5)
```

Or by a program:

Calculator symbol translations:

"STO" arrow: →

Square root sign: √

```ClrList L1
For(A,1,1000)
√(-2*ln(rand))*cos(2*π*A)→L1(A)
End
```

### ZX Spectrum Basic

Here we have converted the QBasic code to suit the ZX Spectrum:

```10 RANDOMIZE 0 : REM seeds random number generator based on uptime
20 DIM a(1000)
30 CLS
40 FOR i = 1 TO 1000
50 LET a(i) = 1 + SQR(-2 * LN(RND)) * COS(2 * PI * RND)
60 NEXT i```

## C

```#include <stdlib.h>
#include <math.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

double drand()   /* uniform distribution, (0..1] */
{
return (rand()+1.0)/(RAND_MAX+1.0);
}
double random_normal()  /* normal distribution, centered on 0, std dev 1 */
{
return sqrt(-2*log(drand())) * cos(2*M_PI*drand());
}
int main()
{
int i;
double rands[1000];
for (i=0; i<1000; i++)
rands[i] = 1.0 + 0.5*random_normal();
return 0;
}
```

## C#

Translation of: JavaScript
```private static double randomNormal()
{
return Math.Cos(2 * Math.PI * tRand.NextDouble()) * Math.Sqrt(-2 * Math.Log(tRand.NextDouble()));
}
```

Then the methods in Random numbers#Metafont are used to calculate the average and the Standard Deviation:

```static Random tRand = new Random();

static void Main(string[] args)
{
double[] a = new double[1000];

double tAvg = 0;
for (int x = 0; x < a.Length; x++)
{
a[x] = randomNormal() / 2 + 1;
tAvg += a[x];
}

tAvg /= a.Length;
Console.WriteLine("Average: " + tAvg.ToString());

double s = 0;
for (int x = 0; x < a.Length; x++)
{
s += Math.Pow((a[x] - tAvg), 2);
}
s = Math.Sqrt(s / 1000);

Console.WriteLine("Standard Deviation: " + s.ToString());

}
```

An example result:

```Average: 1,00510073053613
Standard Deviation: 0,502540443430955
```

## C++

Works with: C++11

The new C++ standard looks very similar to the Boost library example below.

```#include <random>
#include <functional>
#include <vector>
#include <algorithm>
using namespace std;

int main()
{
random_device seed;
mt19937 engine(seed());
normal_distribution<double> dist(1.0, 0.5);
auto rnd = bind(dist, engine);

vector<double> v(1000);
generate(v.begin(), v.end(), rnd);
return 0;
}
```
Works with: C++03
```#include <cstdlib>   // for rand
#include <cmath>     // for atan, sqrt, log, cos
#include <algorithm> // for generate_n

double const pi = 4*std::atan(1.0);

// simple functor for normal distribution
class normal_distribution
{
public:
normal_distribution(double m, double s): mu(m), sigma(s) {}
double operator() const // returns a single normally distributed number
{
double r1 = (std::rand() + 1.0)/(RAND_MAX + 1.0); // gives equal distribution in (0, 1]
double r2 = (std::rand() + 1.0)/(RAND_MAX + 1.0);
return mu + sigma * std::sqrt(-2*std::log(r1))*std::cos(2*pi*r2);
}
private:
const double mu, sigma;
};

int main()
{
double array[1000];
std::generate_n(array, 1000, normal_distribution(1.0, 0.5));
return 0;
}
```
Library: Boost

This example used Mersenne Twister generator. It can be changed by changing the typedef.

```#include <vector>
#include "boost/random.hpp"
#include "boost/generator_iterator.hpp"
#include <boost/random/normal_distribution.hpp>
#include <algorithm>

typedef boost::mt19937 RNGType; ///< mersenne twister generator

int main() {
RNGType rng;
boost::normal_distribution<> rdist(1.0,0.5); /**< normal distribution
with mean of 1.0 and standard deviation of 0.5 */

boost::variate_generator< RNGType, boost::normal_distribution<> >
get_rand(rng, rdist);

std::vector<double> v(1000);
generate(v.begin(),v.end(),get_rand);
return 0;
}
```

## Clojure

```(import '(java.util Random))
(def normals
(let [r (Random.)]
(take 1000 (repeatedly #(-> r .nextGaussian (* 0.5) (+ 1.0))))))
```

## COBOL

```       IDENTIFICATION DIVISION.
PROGRAM-ID. RANDOM.
AUTHOR.  Bill Gunshannon
INSTALLATION.  Home.
DATE-WRITTEN.  14 January 2022.
************************************************************
** Program Abstract:
**   Able to get the Mean to be really close to 1.0 but
**     couldn't get the Standard Deviation any closer than
**     .3 to .4.
************************************************************

DATA DIVISION.

WORKING-STORAGE SECTION.

01  Sample-Size          PIC 9(5)         VALUE 1000.
01  Total                PIC 9(10)V9(5)  VALUE 0.0.
01  Arith-Mean           PIC 999V999  VALUE 0.0.
01  Std-Dev              PIC 999V999  VALUE 0.0.
01  Seed                 PIC 999V999.
01  TI                   PIC 9(8).

01  Idx                  PIC 99999     VALUE 0.
01  Intermediate         PIC 9(10)V9(5)  VALUE 0.0.
01  Rnd-Work.
05  Rnd-Tbl
OCCURS 1 TO 99999 TIMES DEPENDING ON Sample-Size.
10  Rnd              PIC 9V9999999  VALUE 0.0.

PROCEDURE DIVISION.

Main-Program.
ACCEPT TI FROM TIME.
MOVE FUNCTION RANDOM(TI) TO Seed.
PERFORM WITH TEST AFTER VARYING Idx
FROM 1 BY 1
UNTIL Idx = Sample-Size
COMPUTE Intermediate =
(FUNCTION RANDOM() * 2.01)
MOVE Intermediate TO Rnd(Idx)
END-PERFORM.
PERFORM WITH TEST AFTER VARYING Idx
FROM 1 BY 1
UNTIL Idx = Sample-Size
COMPUTE Total = Total + Rnd(Idx)
END-PERFORM.

COMPUTE Arith-Mean = Total / Sample-Size.
DISPLAY "Mean: " Arith-Mean.

PERFORM WITH TEST AFTER VARYING Idx
FROM 1 BY 1
UNTIL Idx = Sample-Size
COMPUTE Intermediate =
Intermediate + (Rnd(Idx) - Arith-Mean) ** 2
END-PERFORM.
COMPUTE Std-Dev = Intermediate / Sample-Size.

DISPLAY "Std-Dev: " Std-Dev.

STOP RUN.

END PROGRAM RANDOM.
```

## Common Lisp

```(loop for i from 1 to 1000
collect (1+ (* (sqrt (* -2 (log (random 1.0)))) (cos (* 2 pi (random 1.0))) 0.5)))
```

## Crystal

Translation of: Phix
```n, mean, sd, tau = 1000, 1, 0.5, (2 * Math::PI)
array = Array.new(n) { mean + sd * Math.sqrt(-2 * Math.log(rand)) * Math.cos(tau * rand) }

mean = array.sum / array.size
standev = Math.sqrt( array.sum{ |x| (x - mean) ** 2 } / array.size )
puts "mean = #{mean}, standard deviation = #{standev}"
```
Output:
`mean = 1.0093442539237896, standard deviation = 0.504694489463623`

## D

```import std.stdio, std.random, std.math;

struct NormalRandom {
double mean, stdDev;

// Necessary because it also defines an opCall.
this(in double mean_, in double stdDev_) pure nothrow {
this.mean = mean_;
this.stdDev = stdDev_;
}

double opCall() const /*nothrow*/ {
immutable r1 = uniform01, r2 = uniform01; // Not nothrow.
return mean + stdDev * sqrt(-2 * r1.log) * cos(2 * PI * r2);
}
}

void main() {
double[1000] array;
auto nRnd = NormalRandom(1.0, 0.5);
foreach (ref x; array)
//x = nRnd;
x = nRnd();
}
```

### Alternative Version

(Untested)

Library: tango
```import tango.math.random.Random;

void main() {
double[1000] list;
auto r = new Random();
foreach (ref l; list) {
r.normalSource!(double)()(l);
l = 1.0 + 0.5 * l;
}
}
```

## Delphi

Delphi has RandG function which generates random numbers with normal distribution using Marsaglia-Bray algorithm:

```program Randoms;

{\$APPTYPE CONSOLE}

uses
Math;

var
Values: array[0..999] of Double;
I: Integer;

begin
//  Randomize;   Commented to obtain reproducible results
for I:= Low(Values) to High(Values) do
Values[I]:= RandG(1.0, 0.5);  // Mean = 1.0, StdDev = 0.5
Writeln('Mean          = ', Mean(Values):6:4);
Writeln('Std Deviation = ', StdDev(Values):6:4);
end.
```
Output:
```Mean          = 1.0098
Std deviation = 0.5016```

## DWScript

```var values : array [0..999] of Float;
var i : Integer;

for i := values.Low to values.High do
values[i] := RandG(1, 0.5);
```

## E

`accum [] for _ in 1..1000 { _.with(entropy.nextGaussian()) }`

## EasyLang

```numfmt 5 0
e = 2.7182818284590452354
for i = 1 to 1000
a[] &= 1 + 0.5 * sqrt (-2 * log10 randomf / log10 e) * cos (360 * randomf)
.
for v in a[]
avg += v / len a[]
.
print "Average: " & avg
for v in a[]
s += pow (v - avg) 2
.
s = sqrt (s / len a[])
print "Std deviation: " & s
```

## Eiffel

```class
APPLICATION

inherit
ARGUMENTS

create
make

feature {NONE} -- Initialization

l_time: TIME
l_seed: INTEGER
math:DOUBLE_MATH
rnd:RANDOM
Size:INTEGER
once
Result:= 1000
end

make
-- Run application.
local
ergebnis:ARRAY[DOUBLE]
tavg: DOUBLE
x: INTEGER
tmp: DOUBLE
text : STRING

do
-- initialize random generator
create l_time.make_now
l_seed := l_time.hour
l_seed := l_seed * 60 + l_time.minute
l_seed := l_seed * 60 + l_time.second
l_seed := l_seed * 1000 + l_time.milli_second
create rnd.set_seed (l_seed)

-- initialize random number container and math
create ergebnis.make_filled (0.0, 1, size)
tavg := 0;
create math

from
x := 1
until
x > ergebnis.count
loop
tmp := randomNormal / 2 + 1
tavg := tavg + tmp
ergebnis.enter (tmp , x)
x := x + 1
end

tavg := tavg / ergebnis.count
text := "Average: "
text.append_double (tavg)
text.append ("%N")
print(text)

tmp := 0
from
x:= 1
until
x > ergebnis.count
loop
tmp := tmp + (ergebnis.item (x) - tavg)^2
x := x + 1
end

tmp := math.sqrt (tmp / ergebnis.count)
text := "Standard Deviation: "
text.append_double (tmp)
text.append ("%N")
print(text)

end

randomNormal:DOUBLE

local

first: DOUBLE
second: DOUBLE

do
rnd.forth
first := rnd.double_item
rnd.forth
second := rnd.double_item

Result := math.cosine (2 * math.pi * first) * math.sqrt (-2 * math.log (second))

end
end
```

Example Result

```Average: 1.0079398405028137
Standard Deviation: 0.49042787564453988
```

## Elena

Translation of: C#

ELENA 6.x :

```import extensions;
import extensions'math;

randomNormal()
{
^ cos(2 * Pi_value * randomGenerator.nextReal())
* sqrt(-2 * ln(randomGenerator.nextReal()))
}

public program()
{
real[] a := new real[](1000);

real tAvg := 0;
for (int x := 0; x < a.Length; x += 1)
{
a[x] := (randomNormal()) / 2 + 1;
tAvg += a[x]
};

tAvg /= a.Length;
console.printLine("Average: ", tAvg);

real s := 0;
for (int x := 0; x < a.Length; x += 1)
{
s += power(a[x] - tAvg, 2)
};

s := sqrt(s / 1000);

console.printLine("Standard Deviation: ", s);

}```
Output:
```Average: 0.9842420481571
Standard Deviation: 0.5109070975558
```

## Elixir

```defmodule Random do
def normal(mean, sd) do
{a, b} = {:rand.uniform, :rand.uniform}
mean + sd * (:math.sqrt(-2 * :math.log(a)) * :math.cos(2 * :math.pi * b))
end
end

std_dev = fn (list) ->
mean = Enum.sum(list) / length(list)
sd = Enum.reduce(list, 0, fn x,acc -> acc + (x-mean)*(x-mean) end) / length(list)
|> :math.sqrt
IO.puts "Mean: #{mean},\tStdDev: #{sd}"
end

xs = for _ <- 1..1000, do: Random.normal(1.0, 0.5)
std_dev.(xs)
```
Output:
```Mean: 1.009079383094275,        StdDev: 0.4991894476975088
```

used Erlang function `:rand.normal`

```xs = for _ <- 1..1000, do: 1.0 + :rand.normal * 0.5
std_dev.(xs)
```
Output:
```Mean: 0.9955701150615597,       StdDev: 0.5036412260426065
```

## Erlang

Works with: Erlang
```mean(Values) ->
mean(tl(Values), hd(Values), 1).

mean([], Acc, Length) ->
Acc / Length;
mean(Values, Acc, Length) ->
mean(tl(Values), hd(Values)+Acc, Length+1).

variance(Values) ->
Mean = mean(Values),
variance(Values, Mean, 0) / length(Values).

variance([], _, Acc) ->
Acc;
variance(Values, Mean, Acc) ->
Diff = hd(Values) - Mean,
DiffSqr = Diff * Diff,
variance(tl(Values), Mean, Acc + DiffSqr).

stddev(Values) ->
math:sqrt(variance(Values)).

normal(Mean, StdDev) ->
U = random:uniform(),
V = random:uniform(),
Mean + StdDev * ( math:sqrt(-2 * math:log(U)) * math:cos(2 * math:pi() * V) ).  % Erlang's math:log is the natural logarithm.

main(_) ->
X = [ normal(1.0, 0.5) || _ <- lists:seq(1, 1000) ],
io:format("mean = ~w\n", [mean(X)]),
io:format("stddev = ~w\n", [stddev(X)]).
```
Output:
```mean = 1.0118289913718608
stddev = 0.5021636849524854
```

## ERRE

```PROGRAM DISTRIBUTION

!
! for rosettacode.org
!

! formulas taken from TI-59 Master Library manual

CONST NUM_ITEM=1000

!VAR SUMX#,SUMX2#,R1#,R2#,Z#,I%

DIM A#[1000]

BEGIN
! seeds random number generator with system time
RANDOMIZE(TIMER)

PRINT(CHR\$(12);)  !CLS
SUMX#=0  SUMX2#=0

FOR I%=1 TO NUM_ITEM DO
R1#=RND(1)  R2#=RND(1)
Z#=SQR(-2*LOG(R1#))*COS(2*π*R2#)
A#[I%]=Z#/2+1   ! I want a normal distribution with
!      mean=1 and std.dev=0.5
SUMX#+=A#[I%]  SUMX2#+=A#[I%]*A#[I%]
END FOR

Z#=SUMX#/NUM_ITEM

PRINT("Average is";Z#)
PRINT("Standard dev. is";SQR(SUMX2#/NUM_ITEM-Z#*Z#))

END PROGRAM
```

## Euler Math Toolbox

```>v=normal(1,1000)*0.5+1;
>mean(v), dev(v)
1.00291801071
0.498226876528```

## Euphoria

Translation of: PureBasic
```include misc.e

function RandomNormal()
atom x1, x2
x1 = rand(999999) / 1000000
x2 = rand(999999) / 1000000
return sqrt(-2*log(x1)) * cos(2*PI*x2)
end function

constant n = 1000
sequence s
s = repeat(0,n)
for i = 1 to n do
s[i] = 1 + 0.5 * RandomNormal()
end for```

## F#

```let n = MathNet.Numerics.Distributions.Normal(1.0,0.5)
List.init 1000 (fun _->n.Sample())
```
Output:
```  [0.734433576; 1.54225304; 0.4407528678; 1.177675412; 0.4318617021;
0.6026656337; 0.769764924; 1.104693934; 0.6297500925; 0.9594598077;
1.684736389; 1.160376323; 0.883354356; 0.9513968363; 0.9727698268;
0.5315570949; 0.9599239266; 1.564976755; 0.7232002879; 1.084139442;
1.220914517; 0.3553085946; 1.112549824; 1.989443553; 0.5752307543;
1.156682549; 0.7886670467; 0.02050745923; 1.532060208; 1.18789591;
1.408946777; 1.038812004; 1.724679503; 1.671565045; 1.266831442;
1.363611654; 1.705819067; 0.5772366328; 0.4503488498; 1.496891481;
0.9831877282; 0.3845460366; 0.8253240671; 1.298969969; 0.4265904553;
0.9303696876; 0.445003361; 0.753175816; 0.6143534043; 1.059982235;
0.7143206784; 0.2233328038; 1.005178481; 0.7697392436; 0.5904948577;
0.5127953044; 0.6467346747; 0.7929387604; -0.1501790761; 0.8750780903;
0.941704369; 1.37941579; 0.4739006145; 1.998886344; 1.219428519;
0.06270791476; 1.097739804; 0.7584232803; 1.042177217; 1.166561247;
1.502357164; 1.171525776; 0.1528807432; 0.2289389756; 1.36208422;
0.3714421124; 1.299571092; 1.171553369; 1.317807265; 1.616662281;
1.724223246; 1.059580642; 1.270520918; -0.1827677907; 1.938593232;
1.420362143; 1.888357595; 0.7851629936; 0.7080554899; 0.7747215818;
1.403719877; 0.5765950249; 1.275206565; 0.6292054813; 1.525562798;
0.6224640457; 0.8524078517; 0.7646595627; 0.6799834691; 0.773111053; ...]
```

==F#==

```let gaussianRand count =
let o = new System.Random()
let pi = System.Math.PI
let gaussrnd =
(fun _ -> 1. + 0.5 * sqrt(-2. * log(o.NextDouble())) * cos(2. * pi * o.NextDouble()))
[ for i in {0 .. (int count)} -> gaussrnd() ]
```

## Factor

```USING: random ;
1000 [ 1.0 0.5 normal-random-float ] replicate
```

## Falcon

```a = []
for i in [0:1000] : a+= norm_rand_num()

function norm_rand_num()
pi = 2*acos(0)
return 1 + (cos(2 * pi * random()) * pow(-2 * log(random()) ,1/2)) /2
end```

## Fantom

Two solutions. The first uses Fantom's random-number generator, which produces a uniform distribution. So, convert to a normal distribution using a formula:

```class Main
{
static const Float PI := 0.0f.acos * 2  // we need to precompute PI

static Float randomNormal ()
{
return (Float.random * PI * 2).cos * (Float.random.log * -2).sqrt
}

public static Void main ()
{
mean := 1.0f
sd := 0.5f
Float[] values := [,] // this is the collection to fill with random numbers
1000.times { values.add (randomNormal * sd + mean) }
}
}```

The second calls out to Java's Gaussian random-number generator:

```using [java] java.util::Random

class Main
{
Random generator := Random()

Float randomNormal ()
{
return generator.nextGaussian
}

public static Void main ()
{
rnd := Main()  // create an instance of Main class, which holds the generator
mean := 1.0f
sd := 0.5f
Float[] values := [,] // this is the collection to fill with random numbers
1000.times { values.add (rnd.randomNormal * sd + mean) }
}
}```

## Forth

Works with: gforth version 0.6.2
```require random.fs
here to seed

-1. 1 rshift 2constant MAX-D	\ or s" MAX-D" ENVIRONMENT? drop

: frnd ( -- f )			\ uniform distribution 0..1
rnd rnd dabs d>f MAX-D d>f f/ ;

: frnd-normal ( -- f )		\ centered on 0, std dev 1
frnd pi f* 2e f* fcos
frnd fln -2e f* fsqrt f* ;

: ,normals ( n -- )		\ store many, centered on 1, std dev 0.5
0 do frnd-normal 0.5e f* 1e f+ f, loop ;

create rnd-array 1000 ,normals
```

For newer versions of gforth (tested on 0.7.3), it seems you need to use HERE SEED ! instead of HERE TO SEED, because SEED has been made a variable instead of a value.

```rnd rnd dabs d>f
```

is necessary, but surprising and definitely not well documented / perhaps not compliant.

## Fortran

Works with: Fortran version 90 and later
```PROGRAM Random

INTEGER, PARAMETER :: n = 1000
INTEGER :: i
REAL :: array(n), pi, temp, mean = 1.0, sd = 0.5

pi = 4.0*ATAN(1.0)
CALL RANDOM_NUMBER(array) ! Uniform distribution

! Now convert to normal distribution
DO i = 1, n-1, 2
temp = sd * SQRT(-2.0*LOG(array(i))) * COS(2*pi*array(i+1)) + mean
array(i+1) = sd * SQRT(-2.0*LOG(array(i))) * SIN(2*pi*array(i+1)) + mean
array(i) = temp
END DO

! Check mean and standard deviation
mean = SUM(array)/n
sd = SQRT(SUM((array - mean)**2)/n)

WRITE(*, "(A,F8.6)") "Mean = ", mean
WRITE(*, "(A,F8.6)") "Standard Deviation = ", sd

END PROGRAM Random
```
Output:
``` Mean = 0.995112
Standard Deviation = 0.503373
```

## Free Pascal

Free Pascal provides the randg function in the RTL math unit that produces Gaussian-distributed random numbers with the Box-Müller algorithm.

```function randg(mean,stddev: float): float;
```

## Frink

`a = new array[[1000], {|x| randomGaussian[1, 0.5]}]`

## Go

This solution uses math/rand package in the standard library. See also though the subrepository rand package at https://godoc.org/golang.org/x/exp/rand, which also has a NormFloat64 and has a rand source with a number of advantages over the one in standard library.

```package main

import (
"fmt"
"math"
"math/rand"
"strings"
"time"
)

const mean = 1.0
const stdv = .5
const n = 1000

func main() {
var list [n]float64
rand.Seed(time.Now().UnixNano())
for i := range list {
list[i] = mean + stdv*rand.NormFloat64()
}
// show computed mean and stdv of list
var s, sq float64
for _, v := range list {
s += v
}
cm := s / n
for _, v := range list {
d := v - cm
sq += d * d
}
fmt.Printf("mean %.3f, stdv %.3f\n", cm, math.Sqrt(sq/(n-1)))
// show histogram by hdiv divisions per stdv over +/-hrange stdv
const hdiv = 3
const hrange = 2
var h [1 + 2*hrange*hdiv]int
for _, v := range list {
bin := hrange*hdiv + int(math.Floor((v-mean)/stdv*hdiv+.5))
if bin >= 0 && bin < len(h) {
h[bin]++
}
}
const hscale = 10
for _, c := range h {
fmt.Println(strings.Repeat("*", (c+hscale/2)/hscale))
}
}
```
Output:
```mean 0.995, stdv 0.503
**
****
******
********
************
************
*************
************
**********
********
*****
***
**
```

## Groovy

```rnd = new Random()
result = (1..1000).inject([]) { r, i -> r << rnd.nextGaussian() }
```

```import System.Random

pairs :: [a] -> [(a,a)]
pairs (x:y:zs) = (x,y):pairs zs
pairs _        = []

gauss mu sigma (r1,r2) =
mu + sigma * sqrt (-2 * log r1) * cos (2 * pi * r2)

gaussians :: (RandomGen g, Random a, Floating a) => Int -> g -> [a]
gaussians n g = take n \$ map (gauss 1.0 0.5) \$ pairs \$ randoms g

result :: IO [Double]
result = getStdGen >>= \g -> return \$ gaussians 1000 g
```

Or using Data.Random from random-fu package:

```replicateM 1000 \$ normal 1 0.5
```

To print them:

```import  Data.Random

thousandRandomNumbers :: RVar [Double]
thousandRandomNumbers =  replicateM 1000 \$ normal 1 0.5

main = do
x <- sample thousandRandomNumbers
print x
```

## HicEst

```REAL :: n=1000, m=1, s=0.5, array(n)

pi = 4 * ATAN(1)
array = s * (-2*LOG(RAN(1)))^0.5  * COS(2*pi*RAN(1)) + m```

## Icon and Unicon

The seed &random may be assigned in either language; either to randomly seed or to pick a fixed starting point. ?i is the random number generator, returning an integer from 0 to i - 1 for non-zero integer i. As a special case, ?0 yields a random floating point number from 0.0 <= r < 1.0

Note that Unicon randomly seeds it's generator.

```procedure main()
local L
L := list(1000)
every L[1 to 1000] := 1.0 + 0.5 * sqrt(-2.0 * log(?0)) * cos(2.0 * &pi * ?0)

every write(!L)
end
```

## IDL

```result = 1.0 + 0.5*randomn(seed,1000)
```

## J

Solution:

```urand=: ?@\$ 0:
zrand=: (2 o. 2p1 * urand) * [: %: _2 * [: ^. urand

1 + 0.5 * zrand 100
```

Alternative Solution:
Using the normal script from the stats/distribs addon.

```   require 'stats/distribs/normal'
1 0.5 rnorm 1000
1.44868803 1.21548637 0.812460657 1.54295452 1.2470606 ...
```

## Java

```double[] list = new double[1000];
double mean = 1.0, std = 0.5;
Random rng = new Random();
for(int i = 0;i<list.length;i++) {
list[i] = mean + std * rng.nextGaussian();
}
```

## JavaScript

```function randomNormal() {
return Math.cos(2 * Math.PI * Math.random()) * Math.sqrt(-2 * Math.log(Math.random()))
}

var a = []
for (var i=0; i < 1000; i++){
a[i] = randomNormal() / 2 + 1
}
```

## jq

Works with: jq version 1.4

Since jq is a purely functional language, it is convenient to define the pseudo-random number generator functions as filters whose inputs and outputs are arrays containing a "seed".

The following uses the same pseudo-random number generator as the Microsoft C Runtime (see Linear congruential generator).

'A Pseudo-Random Number Generator'

```# 15-bit integers generated using the same formula as rand() from the Microsoft C Runtime.
# The random numbers are in [0 -- 32767] inclusive.
# Input: an array of length at least 2 interpreted as [count, state, ...]
# Output: [count+1, newstate, r] where r is the next pseudo-random number.
def next_rand_Microsoft:
.[0] as \$count | .[1] as \$state
| ( (214013 * \$state) + 2531011) % 2147483648 # mod 2^31
| [\$count+1 , ., (. / 65536 | floor) ] ;```

'Box-Muller Method'

```# Generate a single number following the normal distribution with mean 0, variance 1,
# using the Box-Muller method: X = sqrt(-2 ln U) * cos(2 pi V) where U and V are uniform on [0,1].
# Input: [n, state]
# Output [n+1, nextstate, r]
def next_rand_normal:
def u: next_rand_Microsoft | .[2] /= 32767;
u as \$u1
| (\$u1 | u) as \$u2
| ((( (8*(1|atan)) * \$u1[2]) | cos)
* ((-2 * ((\$u2[2]) | log)) | sqrt)) as \$r
| [ (.[0]+1), \$u2[1], \$r] ;

# Generate "count" arrays, each containing a random normal variate with the given mean and standard deviation.
# Input: [count, state]
# Output: [updatedcount, updatedstate, rnv]
# where "state" is a seed and "updatedstate" can be used as a seed.
def random_normal_variate(mean; sd; count):
next_rand_normal
| recurse( if .[0] < count then next_rand_normal else empty end)
| .[2] = (.[2] * sd) + mean;```

Example The task can be completed using: [0,1] | random_normal_variate(1; 0.5; 1000) | .[2]

We show just the sample average and standard deviation:

```def summary:
length as \$l | add as \$sum | (\$sum/\$l) as \$a
| reduce .[] as \$x (0; . + ( (\$x - \$a) | .*. ))
| [ \$a, (./\$l | sqrt)] ;

[ [0,1] | random_normal_variate(1; 0.5; 1000) | .[2] ] | summary```
Output:
```\$ jq -n -c -f Random_numbers.jq
[0.9932830741018853,0.4977760644490579]
```

## Julia

Julia's standard library provides a `randn` function to generate normally distributed random numbers (with mean 0 and standard deviation 0.5, which can be easily rescaled to any desired values):

```randn(1000) * 0.5 + 1
```

## Kotlin

```// version 1.0.6

import java.util.Random

fun main(args: Array<String>) {
val r = Random()
val da = DoubleArray(1000)
for (i in 0 until 1000)  da[i] = 1.0 + 0.5 * r.nextGaussian()
// now check actual mean and SD
val mean = da.average()
val sd = Math.sqrt(da.map { (it - mean) * (it - mean) }.average())
println("Mean is \$mean")
println("S.D. is \$sd")
}
```

Sample output:

Output:
```Mean is 1.0071784073168768
S.D. is 0.48567118114896807
```

## LabVIEW

Works with: LabVIEW version 8.6

## Lingo

```-- Returns a random float value in range 0..1
on randf ()
n = random(the maxinteger)-1
return n / float(the maxinteger-1)
end```
```normal = []
repeat with i = 1 to 1000
normal.add(1 + sqrt(-2 * log(randf())) * cos(2 * PI * randf()) / 2)
end repeat```

## Lobster

Uses built-in `rnd_gaussian`

```let mean = 1.0
let stdv = 0.5
let count = 1000

// stats computes a running mean and variance
// See Knuth TAOCP vol 2, 3rd edition, page 232

def stats(xs: [float]) -> float, float: // variance, mean
var M = xs[0]
var S = 0.0
var n = 1.0
for(xs.length - 1) i:
let x = xs[i + 1]
n = n + 1.0
let mm = (x - M)
M += mm / n
S += mm * (x - M)
return (if n > 0.0: S / n else: 0.0), M

def test_random_normal() -> [float]:
rnd_seed(floor(seconds_elapsed() * 1000000))
let r = vector_reserve(typeof return, count)
for (count):
r.push(rnd_gaussian() * stdv + mean)
let cvar, cmean = stats(r)
let cstdv = sqrt(cvar)
print concat_string(["Mean: ", string(cmean), ", Std.Deviation: ", string(cstdv)], "")

test_random_normal()```

## Logo

Works with: UCB Logo

The earliest Logos only have a RANDOM function for picking a random non-negative integer. Many modern Logos have floating point random generators built-in.

```to random.float   ; 0..1
localmake "max.int lshift -1 -1
output quotient random :max.int :max.int
end

to random.gaussian
output product cos random 360  sqrt -2 / ln random.float
end

make "randoms cascade 1000 [fput random.gaussian / 2 + 1 ?] []```

## Lua

```local list = {}
for i = 1, 1000 do
list[i] = 1 + math.sqrt(-2 * math.log(math.random())) * math.cos(2 * math.pi * math.random()) / 2
end
```

## M2000 Interpreter

M2000 use a Wichmann - Hill Pseudo Random Number Generator.

```Module CheckIt {
Function StdDev (A()) {
\\ A()  has a copy of values
N=Len(A())
if N<1 then Error "Empty Array"
M=Each(A())
k=0
\\ make sum, dev same type as A(k)
sum=A(k)-A(k)
dev=sum
\\ find mean
While M {
sum+=Array(M)
}
Mean=sum/N
\\ make a pointet to A()
P=A()
\\ subtruct from each item
P-=Mean

M=Each(P)
While M {
dev+=Array(M)*Array(M)
}
\\ as pointer to arrray
=(if(dev>0->Sqrt(dev/N), 0), Mean)
}
Function randomNormal {
\\ by default all numbers are double
\\ cos() get degrees
=1+Cos(360 * rnd) * Sqrt(-2 * Ln(rnd)) /2
}
\\ fill array calling  randomNormal() for each item
Dim A(1000)<<randomNormal()
\\ we can pass a pointer to array and place it to stack of values
DisplayMeanAndStdDeviation(A())  ' mean ~ 1 std deviation ~0.5
\\ check M2000 rnd only
Dim B(1000)<<rnd
DisplayMeanAndStdDeviation(B())  ' mean ~ 0.5 std deviation ~0.28

DisplayMeanAndStdDeviation((0,0,14,14))  ' mean = 7 std deviation = 7
DisplayMeanAndStdDeviation((0,6,8,14))  ' mean = 7 std deviation = 5
DisplayMeanAndStdDeviation((6,6,8,8))  ' mean = 7 std deviation = 1

Sub DisplayMeanAndStdDeviation(A)
\\ push to stack all items of an array (need an array pointer)
Push  ! StdDev(A)
\\ read from strack two numbers
Print "Mean is               "; Number
Print "Standard Deviation is "; Number
End Sub
}
Checkit```

## Maple

```with(Statistics):
Sample(Normal(1, 0.5), 1000);```

or

`1+0.5*ArrayTools[RandomArray](1000,1,distribution=normal);`

## Mathematica/Wolfram Language

Built-in function RandomReal with built-in distribution NormalDistribution as an argument:

```RandomReal[NormalDistribution[1, 1/2], 1000]
```

## MATLAB

Native support :

```    mu = 1; sd = 0.5;
x = randn(1000,1) * sd + mu;
```

The statistics toolbox provides this function

```   x = normrnd(mu, sd, [1000,1]);
```

This script uses the Box-Mueller Transform to transform a number from the uniform distribution to a normal distribution of mean = mu0 and standard deviation = chi2.

```function randNum = randNorm(mu0,chi2, sz)

u = ( 2 * rand(sz) ) - 1;
v = ( 2 * rand(sz) ) - 1;

end

randNum = (v .* scaleFactor .* chi2) + mu0;

end
```

Output:

```>> randNorm(1,.5, [1000,1])

ans =

0.693984121077029
```

## Maxima

```load(distrib)\$

random_normal(1.0, 0.5, 1000);
```

## MAXScript

```arr = #()
for i in 1 to 1000 do
(
a = random 0.0 1.0
b = random 0.0 1.0
c = 1.0 + 0.5 * sqrt (-2*log a) * cos (360*b) -- Maxscript cos takes degrees
append arr c
)```

## Metafont

Metafont has `normaldeviate` which produces pseudorandom normal distributed numbers with mean 0 and variance one. So the following complete the task:

```numeric col[];

m := 0;               % m holds the mean, for testing purposes
for i = 1 upto 1000:
col[i] := 1 + .5normaldeviate;
m := m + col[i];
endfor

% testing
m := m / 1000;       % finalize the computation of the mean

s := 0;              % in s we compute the standard deviation
for i = 1 upto 1000:
s := s + (col[i] - m)**2;
endfor
s := sqrt(s / 1000);

show m, s;    % and let's show that really they get what we wanted
end```

A run gave

```>> 0.99947
>> 0.50533```

Assigning a value to the special variable randomseed will allow to have always the same sequence of pseudorandom numbers

## MiniScript

```randNormal = function(mean=0, stddev=1)
return mean + sqrt(-2 * log(rnd,2.7182818284)) * cos(2*pi*rnd) * stddev
end function

x = []
for i in range(1,1000)
x.push randNormal(1, 0.5)
end for
```

## Mirah

```import java.util.Random

list = double[999]
mean = 1.0
std = 0.5
rng = Random.new
0.upto(998) do | i |
list[i] = mean + std * rng.nextGaussian
end```

## МК-61/52

```П7	<->	П8	1/x	П6	ИП6	П9	СЧ	П6	1/x
ln	ИП8	*	2	*	КвКор	ИП9	2	*	пи
*	sin	*	ИП7	+	С/П	БП	05
```

Input: РY - variance, РX - expectation.

Or:

```3	10^x	П0	ПП	13	2	/	1	+	С/П	L0	03	С/П
СЧ	lg	2	/-/	*	КвКор	2	пи	^	СЧ	*	*	cos	*	В/О
```

to generate 1000 numbers with a mean of 1.0 and a standard deviation of 0.5.

## Modula-3

Translation of: C
```MODULE Rand EXPORTS Main;

IMPORT Random;
FROM Math IMPORT log, cos, sqrt, Pi;

VAR rands: ARRAY [1..1000] OF LONGREAL;

(* Normal distribution. *)
PROCEDURE RandNorm(): LONGREAL =
BEGIN
WITH rand = NEW(Random.Default).init() DO
RETURN
sqrt(-2.0D0 * log(rand.longreal())) * cos(2.0D0 * Pi * rand.longreal());
END;
END RandNorm;

BEGIN
FOR i := FIRST(rands) TO LAST(rands) DO
rands[i] := 1.0D0 + 0.5D0 * RandNorm();
END;
END Rand.```

## Nanoquery

Translation of: Java
```list = {0} * 1000
mean = 1.0; std = 0.5
rng  = new(Nanoquery.Util.Random)

for i in range(0, len(list) - 1)
list[i] = mean + std * rng.getGaussian()
end```

## NetRexx

```/* NetRexx */
options replace format comments java crossref symbols nobinary

import java.math.BigDecimal
import java.math.MathContext

-- prologue
numeric digits 20

-- get input, set defaults
parse arg dp mu sigma ec .
if mu    = '' | mu    = '.' then mean             =  1.0; else mean             = mu
if sigma = '' | sigma = '.' then stdDeviation     =  0.5; else stdDeviation     = sigma
if dp    = '' | dp    = '.' then displayPrecision =    1; else displayPrecision = dp
if ec    = '' | ec    = '.' then elements         = 1000; else elements         = ec

-- set up
RNG = Random()
numberList = java.util.List
numberList = ArrayList()

-- generate list of random numbers
loop for elements
rn = mean + stdDeviation * RNG.nextGaussian()
end

-- report
say "Mean:              " mean
say "Standard Deviation:" stdDeviation
say "Precision:         " displayPrecision
say
drawBellCurve(numberList, displayPrecision)

return

-- -----------------------------------------------------------------------------
method drawBellCurve(numberList = java.util.List, precision) static
Collections.sort(numberList)
val = BigDecimal
lastN = ''
nextN = ''
loop val over numberList
nextN = Rexx(val.toPlainString()).format(5, precision)
select
when lastN = '' then nop
when lastN \= nextN then say lastN
otherwise nop
end
say '*\-'
lastN = nextN
end val
say lastN

return```
Output:
```Mean:               1.0
Standard Deviation: 0.5
Precision:          1

*    2.7
**    2.5
*    2.4
***    2.3
*****    2.2
*******    2.1
*************    2.0
*************    1.9
*****************************    1.8
*************************    1.7
*************************************    1.6
******************************************************    1.5
********************************************    1.4
********************************************************************    1.3
*****************************************************************    1.2
**************************************************************************    1.1
*********************************************************************************************    1.0
*************************************************************    0.9
**********************************************************************    0.8
**************************************************************    0.7
***********************************************************************    0.6
**************************************************************    0.5
******************************************    0.4
*******************************    0.3
***************************    0.2
***************    0.1
*********    0.0
******   -0.1
***   -0.2
***   -0.3
*   -0.4
*   -0.6
**   -0.7
```

## NewLISP

```(normal 1 .5 1000)
```

## Nim

```import random, stats, strformat

var rs: RunningStat

randomize()

for _ in 1..5:
for _ in 1..1000: rs.push gauss(1.0, 0.5)
echo &"mean: {rs.mean:.5f}    stdDev: {rs.standardDeviation:.5f}"
```
Output:
```mean: 1.01294    stdDev: 0.49692
mean: 1.00262    stdDev: 0.50028
mean: 0.99878    stdDev: 0.49662
mean: 0.99830    stdDev: 0.49820
mean: 1.00658    stdDev: 0.49703```

## Objeck

```bundle Default {
class RandomNumbers {
function : Main(args : String[]) ~ Nil {
rands := Float->New[1000];
for(i := 0; i < rands->Size(); i += 1;) {
rands[i] := 1.0 + 0.5 * RandomNormal();
};

each(i : rands) {
rands[i]->PrintLine();
};
}

function : native : RandomNormal() ~ Float {
return (2 * Float->Pi() * Float->Random())->Cos() * (-2 * (Float->Random()->Log()))->SquareRoot();
}
}
}```

## OCaml

```let pi = 4. *. atan 1.;;
let random_gaussian () =
1. +. sqrt (-2. *. log (Random.float 1.)) *. cos (2. *. pi *. Random.float 1.);;
let a = Array.init 1000 (fun _ -> random_gaussian ());;
```

## Octave

```p = normrnd(1.0, 0.5, 1000, 1);
disp(mean(p));
disp(sqrt(sum((p - mean(p)).^2)/numel(p)));
```
Output:
```1.0209
0.51048```

## ooRexx

Translation of: REXX

### version 1

```/*REXX pgm gens 1,000 normally distributed #s: mean=1, standard dev.=0.5*/
pi=RxCalcPi()                     /* get value of pi                */
Parse Arg n seed .                /* allow specification of N & seed*/
If n==''|n==',' Then
n=1000                          /* N  is the size of the array.   */
If seed\=='' Then
Call random,,seed               /* use seed for repeatable RANDOM#*/
mean=1                            /* desired new mean (arith. avg.) */
sd=1/2                            /* desired new standard deviation.*/
Do g=1 For n                      /* generate N uniform random nums.*/
n.g=random(0,1e5)/1e5           /* REXX gens uniform rand integers*/
End

Say '              old mean=' mean()
Say 'old standard deviation=' stddev()
Say
Do j=1 To n-1 By 2
m=j+1
/*use Box-Muller method           */
_=sd*RxCalcPower(-2*RxCalcLog(n.j),.5)*RxCalcCos(2*pi*n.m,,'R')+mean
n.m=sd*RxCalcpower(-2*RxCalcLog(n.j),.5)*RxCalcSin(2*pi*n.m,,'R')+,
mean                              /* rand # must be 0???1.          */
n.j=_
End                             /* j                              */
Say '              new mean=' mean()
Say 'new standard deviation=' stddev()
Exit
mean:
_=0
Do k=1 For n
_=_+n.k
End
Return _/n
stddev:
_avg=mean()
_=0
Do k=1 For n
_=_+(n.k-_avg)**2
End
Return RxCalcPower(_/n,.5)

:: requires rxmath library```
Output:
```              old mean= 0.49830002
old standard deviation= 0.283199568

new mean= 1.00377404
new standard deviation= 0.501444536  ```

### version 2

Using the nice function names in the algorithm.

```/*REXX pgm gens 1,000 normally distributed #s: mean=1, standard dev.=0.5*/
pi=RxCalcPi()                     /* get value of pi                */
Parse Arg n seed .                /* allow specification of N & seed*/
If n==''|n==',' Then
n=1000                          /* N  is the size of the array.   */
If seed\=='' Then
Call random,,seed               /* use seed for repeatable RANDOM#*/
mean=1                            /* desired new mean (arith. avg.) */
sd=1/2                            /* desired new standard deviation.*/
Do g=1 For n                      /* generate N uniform random nums.*/
n.g=random(0,1e5)/1e5           /* REXX gens uniform rand integers*/
End

Say '              old mean=' mean()
Say 'old standard deviation=' stddev()
Say
Do j=1 To n-1 By 2
m=j+1
/*use Box-Muller method           */
_=sd*sqrt(-2*ln(n.j))*cos(2*pi*n.m)+mean
n.m=sd*sqrt(-2*ln(n.j))*sin(2*pi*n.m)+mean
n.j=_
End
Say '              new mean=' mean()
Say 'new standard deviation=' stddev()
Exit
mean:
_=0
Do k=1 For n
_=_+n.k
End
Return _/n
stddev:
_avg=mean()
_=0
Do k=1 For n
_=_+(n.k-_avg)**2
End
Return sqrt(_/n)

sqrt: Return RxCalcSqrt(arg(1))
ln:   Return RxCalcLog(arg(1))
cos:  Return RxCalcCos(arg(1),,'R')
sin:  Return RxCalcSin(arg(1),,'R')

:: requires rxmath library```

## PARI/GP

```rnormal()={
my(pr=32*ceil(default(realprecision)*log(10)/log(4294967296)),u1=random(2^pr)*1.>>pr,u2=random(2^pr)*1.>>pr);
sqrt(-2*log(u1))*cos(2*Pi*u2)  \\ in previous version "u1" instead of "u2" was used --> has given crap distribution
\\ Could easily be extended with a second normal at very little cost.
};
vector(1000,unused,rnormal()/2+1)```

## Pascal

The following function calculates Gaussian-distributed random numbers with the Box-Müller algorithm:

```function rnorm (mean, sd: real): real;
{Calculates Gaussian random numbers according to the Box-Müller approach}
var
u1, u2: real;
begin
u1 := random;
u2 := random;
rnorm := mean * abs(1 + sqrt(-2 * (ln(u1))) * cos(2 * pi * u2) * sd);
/* error !?! Shouldn't it be "mean +" instead of "mean *" ? */
end;
```

Delphi and Free Pascal support implement a randg function that delivers Gaussian-distributed random numbers.

## Perl

```my \$PI = 2 * atan2 1, 0;

my @nums = map {
1 + 0.5 * sqrt(-2 * log rand) * cos(2 * \$PI * rand)
} 1..1000;
```

## Phix

Translation of: Euphoria
```function RandomNormal()
return sqrt(-2*log(rnd())) * cos(2*PI*rnd())
end function

sequence s = repeat(0,1000)
for i=1 to length(s) do
s[i] = 1 + 0.5 * RandomNormal()
end for
```

## Phixmonti

```include ..\Utilitys.pmt

def RandomNormal
drop rand log -2 * sqrt 2 pi * rand * cos * 0.5 * 1 +
enddef

1000 var n
0 n repeat

getid RandomNormal map

dup
sum n / var mean
"Mean: " print mean print nl

0 swap n for
get mean - 2 power rot + swap
endfor
swap n / sqrt "Standard deviation: " print print```

## PHP

```function random() {
return mt_rand() / mt_getrandmax();
}

\$pi 	= pi();          // Set PI

\$a = array();
for (\$i = 0; \$i < 1000; \$i++) {
\$a[\$i] = 1.0 + ((sqrt(-2 * log(random())) * cos(2 * \$pi * random())) * 0.5);

}
```

## Picat

```main =>
_ = random2(), % random seed
G = [gaussian_dist(1,0.5) : _ in 1..1000],
println(first_10=G[1..10]),
println([mean=avg(G),stdev=stdev(G)]),
nl.

% Gaussian (Normal) distribution, Box-Muller algorithm
gaussian01() = Y =>
U = frand(0,1),
V = frand(0,1),
Y = sqrt(-2*log(U))*sin(2*math.pi*V).

gaussian_dist(Mean,Stdev) = Mean + (gaussian01() * Stdev).

% Variance of Xs
variance(Xs) = Variance =>
Mu = avg(Xs),
N  = Xs.len,
Variance = sum([ (X-Mu)**2 : X in Xs ]) / N.

% Standard deviation
stdev(Xs) = sqrt(variance(Xs)).```
Output:
```first_10 = [1.639965415776091,0.705425965005482,0.981532402477848,0.309148743347499,1.252800181962738,0.098829881195179,0.74888084504147,0.181494956495445,1.304931340021904,0.595939453660087]
[mean = 0.99223677282248,stdev = 0.510336641737154]```

## PicoLisp

Translation of: C
```(load "@lib/math.l")

(de randomNormal ()  # Normal distribution, centered on 0, std dev 1
(*/
(sqrt (* -2.0 (log (rand 0 1.0))))
(cos (*/ 2.0 pi (rand 0 1.0) `(* 1.0 1.0)))
1.0 ) )

(seed (time))                                      # Randomize

(let Result
(make                                           # Build list
(do 1000                                     # of 1000 elements
(link (+ 1.0 (/ (randomNormal) 2))) ) )
(for N (head 7 Result)                          # Print first 7 results
(prin (format N *Scl) " ") ) )```
Output:
`1.500334 1.212931 1.095283 0.433122 0.459116 1.302446 0.402477`

## PL/I

```/* CONVERTED FROM WIKI FORTRAN */
Normal_Random: procedure options (main);
declare (array(1000), pi, temp,
mean initial (1.0), sd initial (0.5)) float (18);
declare (i, n) fixed binary;

n = hbound(array, 1);
pi = 4.0*ATAN(1.0);
array = random(); /* Uniform distribution */
/* Now convert to normal distribution */
DO i = 1 to n-1 by 2;
temp = sd * SQRT(-2.0*LOG(array(i))) * COS(2*pi*array(i+1)) + mean;
array(i+1) = sd * SQRT(-2.0*LOG(array(i))) * SIN(2*pi*array(i+1)) + mean;
array(i) = temp;
END;
/* Check mean and standard deviation */
mean = SUM(array)/n;
sd = SQRT(SUM((array - mean)**2)/n);
put skip edit ( "Mean = ", mean ) (a, F(18,16) );
put skip edit ( "Standard Deviation = ", sd) (a, F(18,16));
END Normal_Random;```
Output:
```Mean = 1.0125630677913652  Standard Deviation = 0.5067289784535284
3 runs with different seeds to random():
Mean = 1.0008390411168471  Standard Deviation = 0.5095810511317908
Mean = 0.9754351286894838  Standard Deviation = 0.4804376530558166
Mean = 1.0177411222687990  Standard Deviation = 0.5165899662493400
```

## PL/SQL

```DECLARE
--The desired collection
type t_coll is table of number index by binary_integer;
l_coll t_coll;

c_max pls_integer := 1000;
BEGIN
FOR l_counter IN 1 .. c_max
LOOP
-- dbms_random.normal delivers normal distributed random numbers with a mean of 0 and a variance of 1
-- We just adjust the values and get the desired result:
l_coll(l_counter) := DBMS_RANDOM.normal * 0.5 + 1;
DBMS_OUTPUT.put_line (l_coll(l_counter));
END LOOP;
END;```

## Pop11

```;;; Choose radians as arguments to trigonometic functions

;;; procedure generating standard normal distribution
define random_normal() -> result;
lvars r1 = random0(1.0), r2 = random0(1.0);
cos(2*pi*r1)*sqrt(-2*log(r2)) -> result
enddefine;

lvars array, i;

;;; Put numbers on the stack
for i from 1 to 1000 do 1.0+0.5*random_normal() endfor;
;;; collect them into array
consvector(1000) -> array;```

## PowerShell

```function Get-RandomNormal
{
[CmdletBinding()]
Param ( [double]\$Mean, [double]\$StandardDeviation )

\$RandomNormal = \$Mean + \$StandardDeviation * [math]::Sqrt( -2 * [math]::Log( ( Get-Random -Minimum 0.0 -Maximum 1.0 ) ) ) * [math]::Cos( 2 * [math]::PI * ( Get-Random -Minimum 0.0 -Maximum 1.0 ) )

return \$RandomNormal
}

#  Standard deviation function for testing
function Get-StandardDeviation
{
[CmdletBinding()]
param ( [double[]]\$Numbers )

\$Measure = \$Numbers | Measure-Object -Average
\$PopulationDeviation = 0
ForEach (\$Number in \$Numbers) { \$PopulationDeviation += [math]::Pow( ( \$Number - \$Measure.Average ), 2 ) }
\$StandardDeviation = [math]::Sqrt( \$PopulationDeviation / ( \$Measure.Count - 1 ) )
return \$StandardDeviation
}

#  Test
\$RandomNormalNumbers = 1..1000 | ForEach { Get-RandomNormal -Mean 1 -StandardDeviation 0.5 }

\$Measure = \$RandomNormalNumbers | Measure-Object -Average

\$Stats = [PSCustomObject]@{
Count             = \$Measure.Count
Average           = \$Measure.Average
StandardDeviation = Get-StandardDeviation -Numbers \$RandomNormalNumbers
}

\$Stats | Format-List
```
Output:
```Count             : 1000
Average           : 1.01206560135809
StandardDeviation : 0.489099623426272
```

## Python

Using random.gauss
```>>> import random
>>> values = [random.gauss(1, .5) for i in range(1000)]
>>>
```
Quick check of distribution
```>>> def quick_check(numbers):
count = len(numbers)
mean = sum(numbers) / count
sdeviation = (sum((i - mean)**2 for i in numbers) / count)**0.5
return mean, sdeviation

>>> quick_check(values)
(1.0140373306786599, 0.49943411329234066)
>>>
```

Note that the random module in the Python standard library supports a number of statistical distribution methods.

Alternatively using random.normalvariate
```>>> values = [ random.normalvariate(1, 0.5) for i in range(1000)]
>>> quick_check(values)
(0.990099111944864, 0.5029847005836282)
>>>
```

## R

```# For reproducibility, set the seed:
set.seed(12345L)

result <- rnorm(1000, mean = 1, sd = 0.5)```

## Racket

```#lang racket
(for/list ([i 1000])
(add1 (* (sqrt (* -2 (log (random)))) (cos (* 2 pi (random))) 0.5)))
```

Alternative:

```#lang racket
(require math/distributions)
(sample (normal-dist 1.0 0.5) 1000)
```

## Raku

(formerly Perl 6)

Works with: Rakudo version #22 "Thousand Oaks"
```sub randnorm (\$mean, \$stddev) {
\$mean + \$stddev * sqrt(-2 * log rand) * cos(2 * pi * rand)
}

my @nums = randnorm(1, 0.5) xx 1000;

# Checking
say my \$mean = @nums R/ [+] @nums;
say my \$stddev = sqrt \$mean**2 R- @nums R/ [+] @nums X** 2;
```

## Raven

```define PI
-1 acos

define rand1
9999999 choose 1 + 10000000.0 /

define randNormal
rand1 PI * 2 * cos
rand1 log -2 * sqrt
*
2 / 1 +

1000 each drop randNormal "%f\n" print```

Quick Check (on linux with code in file rand.rv)

```raven rand.rv | awk '{sum+=\$1; sumsq+=\$1*\$1;} END {print "stdev = " sqrt(sumsq/NR - (sum/NR)**2); print "mean = " sum/NR}'
stdev = 0.497773
mean = 1.01497```

## ReScript

Translation of: OCaml
```let pi = 4.0 *. atan(1.0)

let random_gaussian = () => {
1.0 +.
sqrt(-2.0 *. log(Random.float(1.0))) *.
cos(2.0 *. pi *. Random.float(1.0))
}

let a = Belt.Array.makeBy(1000, (_) => random_gaussian ())

for i in 0 to 10 {
Js.log(a[i])
}```

## REXX

The REXX language doesn't have any "higher math" functions like SQRT/SIN/COS/LN/LOG/EXP/POW/etc.,
so we hoi polloi REXX programmers have to roll our own.

Programming note:   note the range of the random numbers:   (0,1]
(that is, random numbers from   zero──►unity,   excluding zero, including unity).

```/*REXX pgm generates 1,000 normally distributed numbers:  mean=1,  standard deviation=½.*/
numeric digits 20                                /*the default decimal digit precision=9*/
parse arg n seed .                               /*allow specification of N and the seed*/
if n==''  |  n==","    then n=1000               /*N:    is the size of the array.      */
if datatype(seed,'W')  then call random ,,seed   /*SEED: for repeatable random numbers. */
newMean=1                                        /*the desired new mean (arithmetic avg)*/
sd=1/2                                           /*the desired new standard deviation.  */
do g=1  for n                             /*generate  N uniform random #'s (0,1].*/
#.g = random(1, 1e5)  /  1e5              /*REXX's RANDOM BIF generates integers.*/
end   /*g*/                               /* [↑]  random integers ──► fractions. */
say '              old mean='   mean()
say 'old standard deviation='   stdDev()
call pi;       pi2=pi * 2                        /*define   pi    and also    2 * pi.   */
say
do j=1  to n-1  by 2;    m=j+1            /*step through the iterations by two.  */
_=sd *  sqrt(ln(#.j) * -2)            /*calculate the  used-twice expression.*/
#.j=_ * cos(pi2 * #.m)  +  newMean        /*utilize the  Box─Muller method.      */
#.m=_ * sin(pi2 * #.m)  +  newMean        /*random number must be:      (0,1]    */
end   /*j*/
say '              new mean='     mean()
say 'new standard deviation='     stdDev()
exit                                             /*stick a fork in it,  we're all done. */
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/
mean:   _=0;                   do k=1  for n;  _=_ + #.k;              end;                return      _/n
stdDev: _avg=mean();  _=0;     do k=1  for n;  _=_ + (#.k - _avg)**2;  end;                return sqrt(_/n)
e:      e =2.7182818284590452353602874713526624977572470936999595749669676277240766303535; return e   /*digs overkill*/
pi:     pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862; return pi  /*  "      "   */
r2r:    return arg(1)  //  (pi() * 2)                                                                 /*normalize ang*/
sin:    procedure; parse arg x;x=r2r(x);numeric fuzz min(5,digits()-3);if abs(x)=pi then return 0;return .sincos(x,x,1)
.sincos:parse arg z,_,i; x=x*x; p=z;    do k=2 by 2; _=-_*x/(k*(k+i)); z=z+_; if z=p then leave; p=z; end;     return z
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/
ln:     procedure; parse arg x,f;   call e;   ig= x>1.5;     is=1 - 2 * (ig\==1);           ii=0;             xx=x
do while ig&xx>1.5|\ig&xx<.5;_=e;do k=-1;iz=xx*_**-is;if k>=0&(ig&iz<1|\ig&iz>.5) then leave;_=_*_;izz=iz;end
xx=izz;ii=ii+is*2**k;end;x=x*e**-ii-1;z=0;_=-1;p=z;do k=1;_=-_*x;z=z+_/k;if z=p then leave;p=z;end; return z+ii
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/
cos:    procedure; parse arg x;       x=r2r(x);        a=abs(x);               hpi=pi * .5
numeric fuzz min(6, digits() - 3);      if a=pi    then return -1
if a=hpi | a=hpi*3  then return 0;      if a=pi/3  then return .5
if a=pi * 2/3       then return -.5;                    return .sinCos(1,1,-1)
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/
sqrt:   procedure; parse arg x; if x=0  then return 0;  d=digits();  numeric digits; h=d+6
numeric form; parse value format(x,2,1,,0) 'E0'  with  g 'E' _ .;  g=g * .5'e'_ %2
m.=9;     do j=0  while h>9;       m.j=h;                h=h%2 + 1;      end /*j*/
do k=j+5  to 0  by -1;   numeric digits m.k;   g=(g+x/g)*.5;   end /*k*/
numeric digits d;     return g/1
```

output   when using the default inputs:

```              old mean= 0.5015724
old standard deviation= 0.28652466389342471402

new mean= 0.98807025356443262689
new standard deviation= 0.50002924192766720838
```

## Ring

```for i = 1 to 10
see random(i) + nl
next i```

## RPL

```≪ RAND LN NEG 2 * √
RAND 2 * π * COS *
→NUM 2 / 1 +
≫ 'RANDN' STO

≪ CL∑
1 1000 START RANDN ∑+ NEXT
MEAN PSDEV
```
Output:
```1: .990779804949
2: .487204045227
```

The collection is stored in a predefined array named `∑DAT`, which is automatically created/updated when using the `∑+` instruction and remains available until the user decides to purge it, typically by calling the `CL∑` command.

## Ruby

```Array.new(1000) { 1 + Math.sqrt(-2 * Math.log(rand)) * Math.cos(2 * Math::PI * rand) }
```

## Rust

Library: rand

Using a for-loop:

```extern crate rand;
use rand::distributions::{Normal, IndependentSample};

fn main() {
let mut rands = [0.0; 1000];
let normal = Normal::new(1.0, 0.5);
for num in rands.iter_mut() {
*num = normal.ind_sample(&mut rng);
}
}
```

Using iterators:

```extern crate rand;
use rand::distributions::{Normal, IndependentSample};

fn main() {
let rands: Vec<_> = {
let normal = Normal::new(1.0, 0.5);
(0..1000).map(|_| normal.ind_sample(&mut rng)).collect()
};
}
```

## SAS

```/* Generate 1000 random numbers with mean 1 and standard deviation 0.5.
SAS version 9.2 was used to create this code.*/

data norm1000;
call streaminit(123456);
/* Set the starting point, so we can replicate results.
If you want different results each time, comment the above line. */
do i=1 to 1000;
r=rand('normal',1,0.5);
output;
end;
run;
```

Results:

``` The MEANS Procedure

Analysis Variable : r

Mean         Std Dev
----------------------------
0.9907408       0.4844051
----------------------------
```

## Sather

```class MAIN is
main is
a:ARRAY{FLTD} := #(1000);
i:INT;

RND::seed(2010);
loop i := 1.upto!(1000) - 1;
a[i] := 1.0d + 0.5d * RND::standard_normal;
end;

-- testing the distribution
mean ::= a.reduce(bind(_.plus(_))) / a.size.fltd;
#OUT + "mean " + mean + "\n";
a.map(bind(_.minus(mean)));
a.map(bind(_.pow(2.0d)));
dev ::= (a.reduce(bind(_.plus(_))) / a.size.fltd).sqrt;
#OUT + "dev  " + dev + "\n";
end;
end;```

## Scala

### One liner

```List.fill(1000)(1.0 + 0.5 * scala.util.Random.nextGaussian)
```

```object RandomNumbers extends App {

val distribution: LazyList[Double] = {
def randomNormal: Double = 1.0 + 0.5 * scala.util.Random.nextGaussian

def normalDistribution(a: Double): LazyList[Double] = a #:: normalDistribution(randomNormal)

normalDistribution(randomNormal)
}

/*
* Let's test it
*/
def calcAvgAndStddev[T](ts: Iterable[T])(implicit num: Fractional[T]): (T, Double) = {
val mean: T =
num.div(ts.sum, num.fromInt(ts.size)) // Leaving with type of function T

// Root of mean diffs
val stdDev = Math.sqrt(ts.map { x =>
val diff = num.toDouble(num.minus(x, mean))
diff * diff
}.sum / ts.size)

(mean, stdDev)
}

println(calcAvgAndStddev(distribution.take(1000))) // e.g. (1.0061433267806525,0.5291834867560893)
}
```

## Scheme

```; linear congruential generator given in C99 section 7.20.2.1
(define ((c-rand seed)) (set! seed (remainder (+ (* 1103515245 seed) 12345) 2147483648)) (quotient seed 65536))

; uniform real numbers in open interval (0, 1)
(define (unif-rand seed) (let ((r (c-rand seed))) (lambda () (/ (+ (r) 1) 32769.0))))

; Box-Muller method to generate normal distribution
(define (normal-rand unif m s)
(let ((? #t) (! 0.0) (twopi (* 2.0 (acos -1.0))))
(lambda ()
(set! ? (not ?))
(if ? !
(let ((a (sqrt (* -2.0 (log (unif))))) (b (* twopi (unif))))
(set! ! (+ m (* s a (sin b))))
(+ m (* s a (cos b))))))))

(define rnorm (normal-rand (unif-rand 0) 1.0 0.5))

; auxiliary function to get a list of 'n random numbers from generator 'r
(define (rand-list r n) = (if (zero? n) '() (cons (r) (rand-list r (- n 1)))))

(define v (rand-list rnorm 1000))

v
#|
(-0.27965824722565835
-0.8870860825789542
0.6499618744638194
0.31336141955110863
...
0.5648743998193049
0.8282656735558756
0.6399951934564637
0.7699535302478072)
|#

; check mean and standard deviation
(define (mean-sdev v)
(let loop ((v v) (a 0) (b 0) (n 0))
(if (null? v)
(let ((mean (/ a n)))
(list mean (sqrt (/ (- b (* n mean mean)) (- n 1)))))
(let ((x (car v)))
(loop (cdr v) (+ a x) (+ b (* x x)) (+ n 1))))))

(mean-sdev v)
; (0.9562156817697293 0.5097087109575911)
```

## Seed7

```\$ include "seed7_05.s7i";
include "float.s7i";
include "math.s7i";

const func float: frand is func  # Uniform distribution, (0..1]
result
var float: frand is 0.0;
begin
repeat
frand := rand(0.0, 1.0);
until frand <> 0.0;
end func;

const func float: randomNormal is  # Normal distribution, centered on 0, std dev 1
return sqrt(-2.0 * log(frand)) * cos(2.0 * PI * frand);

const proc: main is func
local
var integer: i is 0;
var array float: rands is 1000 times 0.0;
begin
for i range 1 to length(rands) do
rands[i] := 1.0 + 0.5 * randomNormal;
end for;
end func;```

## Sidef

```var arr = 1000.of { 1 + (0.5 * sqrt(-2 * 1.rand.log) * cos(Num.tau * 1.rand)) }
arr.each { .say }
```

## Standard ML

Works with: SML/NJ

SML/NJ has two structures for random numbers:

1) Rand (a linear congruential generator). You create the generator by calling `Rand.mkRandom` with a seed (of `word` type). You can call the generator with `()` repeatedly to get a word in the range `[Rand.randMin, Rand.randMax]`. You can use the `Rand.norm` function to transform the output into a `real` from 0 to 1, or use the `Rand.range (i,j)` function to transform the output into an `int` of the given range.

```val seed = 0w42;
val gen = Rand.mkRandom seed;
fun random_gaussian () =
1.0 + Math.sqrt (~2.0 * Math.ln (Rand.norm (gen ()))) * Math.cos (2.0 * Math.pi * Rand.norm (gen ()));
val a = List.tabulate (1000, fn _ => random_gaussian ());
```

2) Random (a subtract-with-borrow generator). You create the generator by calling `Random.rand` with a seed (of a pair of `int`s). You can use the `Random.randInt` function to generate a random int over its whole range; `Random.randNat` to generate a non-negative random int; `Random.randReal` to generate a `real` between 0 and 1; or `Random.randRange (i,j)` to generate an `int` in the given range.

```val seed = (47,42);
val gen = Random.rand seed;
fun random_gaussian () =
1.0 + Math.sqrt (~2.0 * Math.ln (Random.randReal gen)) * Math.cos (2.0 * Math.pi * Random.randReal gen);
val a = List.tabulate (1000, fn _ => random_gaussian ());
```

Other implementations of Standard ML have their own random number generators. For example, Moscow ML has a `Random` structure that is different from the one from SML/NJ.

Works with: Poly/ML

The SML Basis Library does not provide a routine for uniform deviate generation, and PolyML does not have one. Using a routine from "Monte Carlo" by Fishman (Springer), in the function uniformdeviate, and avoiding the slow IntInf's:

```val urandomlist =  fn seed => fn n =>
let
val uniformdeviate = fn seed =>
let
val in31m = (Real.fromInt o Int32.toInt ) (getOpt (Int32.maxInt,0) );
val in31 = in31m +1.0;
val s1 = 41160.0;
val s2 = 950665216.0;
val v = Real.realFloor seed;
val val1 = v*s1;
val val2 = v*s2;
val next1 = Real.fromLargeInt (Real.toLargeInt IEEEReal.TO_NEGINF (val1/in31)) ;
val next2 = Real.rem(Real.realFloor(val2/in31) , in31m );
val valt = val1+val2 - (next1+next2)*in31m;
val nextt = Real.realFloor(valt/in31m);
val valt = valt - nextt*in31m;
in
(valt/in31m,valt)
end;
val store =  ref (0.0,0.0);
val rec u =  fn S => fn 0 => [] | n=> (store:=uniformdeviate S; (#1 (!store)):: (u (#2 (!store)) (n-1))) ;
in
u seed n
end;

local
open Math
in
val bmconv = fn urand => fn vrand => 1.0+0.5*(sqrt(~2.0*ln urand)*cos (2.0*pi*vrand) )
end;

val rec makeNormals = fn once => fn u::v::[] => [once u v] |
u::v::rm => (once u v )::(makeNormals once rm );

val anyrealseed=1009.0 ;
makeNormals bmconv (urandomlist anyrealseed 2000);
```

## Stata

```clear all
set obs 1000
gen x=rnormal(1,0.5)
```

### Mata

```a = rnormal(1000,1,1,0.5)
```

## Tcl

```package require Tcl 8.5
variable ::pi [expr acos(0)]
proc ::tcl::mathfunc::nrand {} {
expr {sqrt(-2*log(rand())) * cos(2*\$::pi*rand())}
}

set mean 1.0
set stddev 0.5
for {set i 0} {\$i < 1000} {incr i} {
lappend result [expr {\$mean + \$stddev*nrand()}]
}
```

## TorqueScript

```for (%i = 0; %i < 1000; %i++)
%list[%i] = 1 + mSqrt(-2 * mLog(getRandom())) * mCos(2 * \$pi * getRandom());```

## Uiua

```# Generate normal distribution with mean = 1, sd = 0.5
Gauss ← (×(∿+η××2π⚂) (√×¯2ₙe⚂))
[⍥(+1×0.5Gauss)1000] # -> mean = 1, sd = 0.5
Mean ← ÷⧻⟜/+
Sd ← √÷⊃(⋅⧻|/+×.-)Mean.
⊸⊃Sd Mean```
Output:
```[0.5549431582333864 0.6318745755493541 0.243308532933335 ...etc...]
1.026433943526625
0.4953224276300106
```

## Ursala

There are two ways of interpreting the task, either to simulate sampling a population described by the given statistics, or to construct a sample exhibiting the given statistics. Both are illustrated below. The functions parameterized by the mean and standard deviation take a sample size and return a sample of that size, represented as a list of floating point numbers. The Z library function simulates a draw from a standard normal distribution. Mean and standard deviation library functions are also used in this example.

```#import nat
#import flo

pop_stats("mu","sigma") = plus/*"mu"+ times/*"sigma"+ Z*+ iota

sample_stats("mu","sigma") = plus^*D(minus/"mu"+ mean,~&)+ vid^*D(div\"sigma"+ stdev,~&)+ Z*+ iota

#cast %eWL

test =

^(mean,stdev)* <
pop_stats(1.,0.5) 1000,
sample_stats(1.,0.5) 1000>```

The output shows the mean and standard deviation for both sample vectors, the latter being exact by construction.

```<
(1.004504e+00,4.915525e-01),
(1.000000e+00,5.000000e-01)>```

## Visual FoxPro

```LOCAL i As Integer, m As Double, n As Integer, sd As Double
py = PI()
SET TALK OFF
SET DECIMALS TO 6
CREATE CURSOR gdev (deviate B(6))
RAND(-1)
n = 1000
m = 1
sd = 0.5
CLEAR
FOR i = 1 TO n
INSERT INTO gdev VALUES (GaussDev(m, 1/sd))
ENDFOR
CALCULATE AVG(deviate), STD(deviate) TO m, sd
? "Mean", m, "Std Dev", sd
SET TALK ON
SET DECIMALS TO

FUNCTION GaussDev(mean As Double, sdev As Double) As Double
LOCAL z As Double
z = SQRT(-2*LOG(RAND()))*COS(py*RAND())
IF sdev # 0
z = mean + z/sdev
ENDIF
RETURN z
ENDFUNC
```

## V (Vlang)

```import crypto.rand

fn main() {
mut nums := []u64{}
for _ in 0..1000 {
nums << rand.int_u64(10000) or {0} // returns random unsigned 64-bit integer from real OS source of entropy
}
println(nums)
}```

## Wren

```import "random" for Random

var rand = Random.new()

var randNormal = Fn.new { (-2 * rand.float().log).sqrt * (2 * Num.pi * rand.float()).cos }

var stdDev = Fn.new { |a, m|
var c = a.count
return ((a.reduce(0) { |acc, x| acc + x*x } - m*m*c) / c).sqrt
}

var n = 1000
var numbers = List.filled(n, 0)
var mu = 1
var sigma = 0.5
var sum = 0
for (i in 0...n) {
numbers[i] = mu + sigma*randNormal.call()
sum = sum + numbers[i]
}
var mean = sum / n
System.print("Actual mean   : %(mean)")
System.print("Actual std dev: %(stdDev.call(numbers, mean))")
```
Output:

Sample run:

```Actual mean   : 1.0053988699746
Actual std dev: 0.4961645117026
```

## XPL0

Translation of: C
```define PI = 3.14159265358979323846;

func real DRand;        \Uniform distribution, [0..1]
return float(Ran(1_000_000)) / 1e6;

func real RandomNormal; \Normal distribution, centered on 0, std dev 1
return sqrt(-2.*Log(DRand)) * Cos(2.*PI*DRand);

int  I;
real Rands(1000);
for I:= 0 to 1000-1 do
Rands(I):= 1.0 + 0.5*RandomNormal```

## Yorick

Returns array of count random numbers with mean 0 and standard deviation 1.

```func random_normal(count) {
return sqrt(-2*log(random(count))) * cos(2*pi*random(count));
}```

Example of basic use:

```> nums = random_normal(1000); // create an array 1000 random numbers
> nums(avg); // show the mean
0.00901216
> nums(rms); // show the standard deviation
0.990265```

Example with a mean of 1.0 and a standard deviation of 0.5:

```> nums = random_normal(1000) * 0.5 + 1;
> nums(avg);
1.00612
> nums(rms);
0.496853```

## zkl

```fcn mkRand(mean,sd){ //normally distributed random w/mean & standard deviation
pi:=(0.0).pi;    // using the Box–Muller transform
rz1:=fcn{1.0-(0.0).random(1)}  // from [0,1) to (0,1]
return('wrap(){((-2.0*rz1().log()).sqrt() * (2.0*pi*rz1()).cos())*sd + mean })
}```

This creates a new random number generator, now to use it:

```var g=mkRand(1,0.5);
ns:=(0).pump(1000,List,g); // 1000 rands with mean==1 & sd==1/2
mean:=(ns.sum(0.0)/1000);  //-->1.00379
// calc sd of list of numbers:
(ns.reduce('wrap(p,n){p+(n-mean).pow(2)},0.0)/1000).sqrt() //-->0.494844```