! ! 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(-2LOG(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



```Euler Math Toolbox

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

Euphoria

{{trans|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

=={{header|F_Sharp|F#}}==

let n = MathNet.Numerics.Distributions.Normal(1.0,0.5)
List.init 1000 (fun _->n.Sample())

{{out}}

  [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; ...]

Factor

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

=={{header|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|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.

Fortran

{{works with|Fortran|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

{{out}}

 Mean = 0.995112
 Standard Deviation = 0.503373

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: {{out}}

Mean is               1.000763573902885
Standard Deviation is 0.500653063426955

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;

=={{header|F_Sharp|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() ]

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

{{out}}

mean 0.995, stdv 0.503
**
****
******
********
************
************
*************
************
**********
********
*****
***
**

FutureBasic

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

include "ConsoleWindow"

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

local fn RandomGaussian as double
dim as double r

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

dim as long i
dim as 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

Output:

           Average: 1.0258434498
Standard Deviation: 0.2771047023

Groovy

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

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

Or using Data.Random from random-fu package:

replicateM 1000 $ normal 1 0.5

To print them:

import  Data.Random
import Control.Monad

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

=={{header|Icon}} and {{header|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 [[j:Addons/stats/distribs|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|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

{{out}} $ 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: {{out}}

Mean is 1.0071784073168768
S.D. is 0.48567118114896807

LabVIEW

{{works with|LabVIEW|8.6}} [[File:LV_array_of_randoms_with_given_mean_and_stdev.png]]

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

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 ?] []

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

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

Mathematica

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)

    radiusSquared = +Inf;

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

        radiusSquared = u.^2 + v.^2;
    end

    scaleFactor = sqrt( ( -2*log(radiusSquared) )./ radiusSquared );
    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

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

=={{header|MK-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:

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

=={{header|Modula-3}}== {{trans|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.

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()
  numberList.add(BigDecimal(rn, MathContext.DECIMAL128))
  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

{{out}}

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 math, strutils

const  precisn = 5
var rs: TRunningStat

proc normGauss: float {.inline.} = 1 + 0.76 * cos(2*PI*random(1.0)) * sqrt(-2*log10(random(1.0)))

randomize()

for j in 0..5:
   for i in 0..1000:
      rs.push(normGauss())
   echo("mean: ", $formatFloat(rs.mean,ffDecimal,precisn),
        " stdDev: ", $formatFloat(rs.standardDeviation(),ffDecimal,precisn))

{{out}}

mean: 1.01703 stdDev: 0.50324
mean: 1.01187 stdDev: 0.50060
mean: 1.00216 stdDev: 0.49969
mean: 1.00335 stdDev: 0.50184
mean: 1.00120 stdDev: 0.49830
mean: 1.00217 stdDev: 0.49911

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

{{out}}

1.0209
0.51048

ooRexx

{{trans|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

{{out}}

              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*u1)
	\\ 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);
end;

[[#Delphi | Delphi]] and [[#Free Pascal|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;

Perl 6

{{works with|Rakudo|#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;

Phix

{{Trans|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

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

}

PicoLisp

{{trans|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) " ") ) )

{{out}}

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;

{{out}}

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

PowerShell

Equation adapted from Liberty BASIC

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

{{out}}

Count             : 1000
Average           : 1.01206560135809
StandardDeviation : 0.489099623426272

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

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

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

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

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

Ruby

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

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

Rust

{{libheader|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);
    let mut rng = rand::thread_rng();
    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);
        let mut rng = rand::thread_rng();
        (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)

Academic

val distrubution = {
  def randomNormal = 1.0 + 0.5 * scala.util.Random.nextGaussian

  def normalDistribution(a: Double): Stream[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 = sqrt(ts.map { x =>
      val diff = num.toDouble(num.minus(x, mean))
      diff * diff
    }.sum / ts.size)

    (mean, stdDev)
  }

println(calcAvgAndStddev(distrubution.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 ints). 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|PolyML}} 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()}]
}

=={{header|TI-83 BASIC}}== Builtin function: randNorm() randNorm(1,.5)

Or by a program:

Calculator symbol translations:

"STO" arrow: →

Square root sign: √

ClrList L₁ Radian For(A,1,1000) √(-2ln(rand))cos(2πA)→L₁(A) End

TorqueScript

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

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

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

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