Random numbers
The goal of this task is to generate a collection filled with 1000 normally distributed random (or pseudorandom) numbers with a mean of 1.0 and a standard deviation of 0.5
You are encouraged to solve this task according to the task description, using any language you may know.
Many libraries only generate uniformly distributed random numbers. If so, use this formula to convert them to a normal distribution.
Ada
<lang ada>with Ada.Numerics; use Ada.Numerics; with Ada.Numerics.Float_Random; use Ada.Numerics.Float_Random; with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;
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;</lang>
ALGOL 68
<lang algol68>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]))
)</lang> Output sample:
( 0.693461, 0.948424, 0.482261, 1.045939, 0.890818, 1.467935, 0.604153, 0.804811, 0.690227, 0.83462 ... )
AutoHotkey
contributed by Laszlo on the ahk forum <lang AutoHotkey>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
}</lang>
AWK
<lang awk>$ 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}'</lang>
BASIC
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
BBC BASIC
<lang bbcbasic> 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</lang>
Output:
Mean = 1.01848064 Standard deviation = 0.503551814
C
<lang 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;
}</lang>
C#
This is the Javascript translation into C#. <lang csharp> private static double randomNormal() { return Math.Cos(2 * Math.PI * tRand.NextDouble()) * Math.Sqrt(-2 * Math.Log(tRand.NextDouble())); } </lang>
Then the methods in Metafont are used to calculate the average and the Standard Deviation: <lang csharp> 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());
Console.ReadLine(); } </lang>
An example result:
Average: 1,00510073053613 Standard Deviation: 0,502540443430955
C++
<lang cpp>#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() // 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:
double mu, sigma;
};
int main() {
double array[1000]; std::generate_n(array, 1000, normal_distribution(1.0, 0.5)); return 0;
}</lang>
This example used Mersenne Twister generator. It can be changed by changing the typedef.
<lang cpp>
- 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;
} </lang>
C++0x
The new C++ standard looks very similar to the Boost library here.
<lang cpp>#include <random>
- include <functional>
- include <vector>
- include <algorithm>
using namespace std;
int main() {
random_device seed; mt19937 engine(seed()); normal_distribution<> dist(1.0, 0.5); auto rnd = bind(dist, engine);
vector<double> v(1000); generate(v.begin(), v.end(), rnd);
}</lang>
Clojure
<lang lisp>(import '(java.util Random)) (def normals
(let [r (Random.)] (take 1000 (repeatedly #(-> r .nextGaussian (* 0.5) (+ 1.0))))))</lang>
Common Lisp
<lang lisp>(loop for i from 1 to 1000
collect (1+ (* (sqrt (* -2 (log (random 1.0)))) (cos (* 2 pi (random 1.0))) 0.5)))</lang>
D
<lang d>import std.stdio, std.random, std.math;
struct NormalRandom {
double mean, stdDev;
this(double mean_, double stdDev_) { this.mean = mean_; this.stdDev = stdDev_; }
double opCall() { double r1 = uniform(0.0, 1.0); double r2 = uniform(0.0, 1.0); return mean + stdDev * sqrt(-2 * log(r1)) * cos(2 * PI * r2); }
}
void main() {
double[1000] array; auto nrnd = NormalRandom(1.0, 0.5); foreach (ref x; array) x = nrnd();
}</lang>
<lang d>import tango.math.random.Random; double [1000]list; auto r = new Random(); foreach(ref l; list) {
r.normalSource!(double)()(l); l = 1.0 + 0.5 * l;
}</lang>
Delphi
Delphi has RandG function which generates random numbers with normal distribution using Marsaglia-Bray algorithm:
<lang Delphi>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); Readln;
end.</lang> Output:
Mean = 1.0098 Std deviation = 0.5016
E
<lang e>accum [] for _ in 1..1000 { _.with(entropy.nextGaussian()) }</lang>
Euphoria
<lang euphoria>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</lang>
Factor
<lang factor>1000 [ 1.0 0.5 normal-random-float ] replicate</lang>
Falcon
<lang 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</lang>
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:
<lang fantom> 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) } }
} </lang>
The second calls out to Java's Gaussian random-number generator:
<lang fantom> 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) } }
} </lang>
Forth
<lang forth>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</lang>
Fortran
<lang fortran>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</lang>
Output
Mean = 0.995112 Standard Deviation = 0.503373
F#
<lang fsharp>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() ]</lang>
Go
<lang go>package main import "rand" import "time"
func main() {
var list [1000]float64 rand.Seed(time.Nanoseconds()) for i := range list { list[i] = 1.0 + 0.5 * rand.NormFloat64() }
}</lang>
Groovy
<lang groovy>rnd = new Random() result = (1..1000).inject([]) { r, i -> r << rnd.nextGaussian() }</lang>
Haskell
<lang haskell>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</lang>
HicEst
<lang 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 </lang>
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. <lang icon> 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 </lang>
IDL
<lang idl>result = 1.0 + 0.5*randomn(seed,1000)</lang>
J
Solution: <lang j>urand=: ?@$ 0: zrand=: (2 o. 2p1 * urand) * [: %: _2 * [: ^. urand
1 + 0.5 * zrand 100</lang>
Alternative Solution:
Using the normal script from the stats/distribs addon.
<lang j> require 'stats/distribs/normal'
1 0.5 rnorm 1000
1.44868803 1.21548637 0.812460657 1.54295452 1.2470606 ...</lang>
Java
<lang 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();
}</lang>
JavaScript
<lang 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
}</lang>
LabVIEW
Liberty BASIC
<lang lb>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</lang>
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. <lang logo>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 ?] []</lang>
lua
<lang 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 </lang>
Mathematica
Built-in function RandomReal with built-in distribution NormalDistribution as an argument: <lang Mathematica>RandomReal[NormalDistribution[1, 1/2], 1000]</lang>
MATLAB
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.
<lang MATLAB>function randNum = randNorm(mu0,chi2)
radiusSquared = +Inf;
while (radiusSquared >= 1) u = ( 2 * rand(1) ) - 1; v = ( 2 * rand(1) ) - 1;
radiusSquared = u^2 + v^2; end
scaleFactor = sqrt( ( -2*log(radiusSquared) )/ radiusSquared ); randNum = (v * scaleFactor * chi2) + mu0;
end</lang>
Output: <lang MATLAB>>> randNorm(1,.5)
ans =
0.693984121077029</lang>
MAXScript
<lang 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
)</lang>
Metafont
Metafont has normaldeviate
which produces pseudorandom normal distributed numbers with mean 0 and variance one. So the following complete the task:
<lang metafont>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</lang>
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
Mirah
<lang 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 </lang>
Modula-3
<lang modula3>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.</lang>
Objeck
<lang 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(); } }
} </lang>
OCaml
<lang 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 ());;</lang>
Octave
<lang octave>p = normrnd(1.0, 0.5, 1000, 1); disp(mean(p)); disp(sqrt(sum((p - mean(p)).^2)/numel(p));</lang>
Output:
1.0209 0.51048
PARI/GP
<lang parigp>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*u1) \\ Could easily be extended with a second normal at very little cost. }; vector(1000,unused,rnormal()/2+1)</lang>
Pascal
See Delphi
Perl
<lang perl>my $PI = 2 * atan2 1, 0;
my @nums = map {
1 + 0.5 * sqrt(-2 * log rand) * cos(2 * $PI * rand)
} 1..1000;</lang>
Perl 6
<lang perl6>sub randnorm ($mean, $stddev) {
$mean + $stddev * sqrt(-2 * log rand) * cos(2 * pi * rand)
}
my @nums = map { randnorm 1, 0.5 }, ^1000;</lang>
PHP
<lang 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);
}</lang>
PicoLisp
<lang PicoLisp>(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) " ") ) )</lang>
Output:
1.500334 1.212931 1.095283 0.433122 0.459116 1.302446 0.402477
PL/I
<lang 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; </lang>
Pop11
<lang pop11>;;; Choose radians as arguments to trigonometic functions true -> popradians;
- 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;</lang>
PureBasic
<lang 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.f(n-1) For i = 0 To n-1
a(i) = 1 + 0.5 * RandomNormal()
Next</lang>
Python
<lang python>import random values = [random.gauss(1, .5) for i in range(1000)]
- or [ random.normalvariate(1, 0.5) for i in range(1000)]</lang>
Note that the random module in the Python standard library supports a number of statistical distribution methods.
Quick check <lang python>>>> mean = sum(values)/1000 >>> sdeviation = (sum((i - mean)**2 for i in values)/1000)**0.5 >>> mean, sdeviation (1.0127861555468178, 0.5006682783828207) </lang>
R
<lang r>result <- rnorm(1000, mean=1, sd=0.5)</lang>
Ruby
<lang ruby>Array.new(1000) { 1 + Math.sqrt(-2 * Math.log(rand)) * Math.cos(2 * Math::PI * rand) }</lang>
Sather
<lang 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;</lang>
Scala
<lang scala>List.fill(1000)(1.0 + 0.5 * scala.util.Random.nextGaussian)</lang>
Seed7
<lang 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;</lang>
Standard ML
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.
<lang sml>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 ());</lang>
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.
<lang sml>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 ());</lang>
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.
Tcl
<lang 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()}]
}</lang>
TI-83 BASIC
Calculator symbol translations:
"STO" arrow: →
Square root sign: √
ClrList L1 Radian For(A,1,1000) √(-2*ln(rand))*cos(2*π*A)→L1(A) End
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. <lang Ursala>#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></lang>
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)>
Yorick
Returns array of count random numbers with mean 0 and standard deviation 1. <lang yorick>func random_normal(count) {
return sqrt(-2*log(random(count))) * cos(2*pi*random(count));
}</lang>
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
ZX Spectrum Basic
Here we have converted the QBasic code to suit the ZX Spectrum:
<lang zxbasic>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</lang>