E
August 13, 2010
We begin with a function f to count the random numbers required to reach a sum greater than one:
(define (f)
(do ((i 0 (+ i 1))
(s 0.0 (+ s (rand))))
((< 1 s) i)))
Now we run the simulation n times:
(define (e n)
(do ((i 0 (+ i 1))
(s 0.0 (+ s (f))))
((= i n) (/ s i))))
And the answer is:
> (e #e1e6)
2.718844
Run f enough times, and the surprising average is the transcendental number e. A good explanation of the math involved is given by John Walker. Mathematically, this puzzle is known as the uniform sum distribution.
We used rand from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/AWzS1d7U.
Programming Praxis 
[…] Praxis – E By Remco Niemeijer In today’s Programming Praxis exercise our task is to determine how many random numbers between 0 and 1 we […]
My Haskell solution (see http://bonsaicode.wordpress.com/2010/08/13/programming-praxis-e/ for a version with comments):
#include <stdio.h> #include <stdlib.h> main() { int trials = 1000000 ; int ttrials = 0; int i, j ; double sum ; for (i=0; i<trials; i++) for (sum=0., ttrials++; (sum += drand48()) < 1.; ttrials++) ; printf("%f trials on average.\n", (double) ttrials / trials) ; }A couple of improvements could be made. I left the j variable in accidently, it is unused. In C++, the scope of the sum value could be confined to the for loop in which it appears. A do { } while loop would mean that I didnt have two different increments of the ttrial variable.
Here’s my solution in Clojure:
(defn to-one [] (loop [times 0 sum 0] (cond (> sum 1) times :else (recur (inc times) (+ sum (rand)))))) (defn to-one-avg [n] (/ (reduce + (repeatedly n to-one)) n 1.0)) user=> (to-one-avg 10000000) 2.7184935PLT “racket”:
#! /bin/sh
#| Hey Emacs, this is -*-scheme-*- code!
#$Id$
exec mzscheme -l errortrace –require “$0″ –main — ${1+”$@”}
|#
#lang scheme
;; https://programmingpraxis.com/2010/08/13/e/
(define (trial)
(let loop ([number-of-numbers 0]
[sum 0])
(if (inexact
(/ (apply + outcomes)
(length outcomes)))))
(provide main)
/me tips his hat to Mark VandeWettering, Oregon ’86 or thereabouts
Gaah. Stupid wordpad.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
gistfile1
hosted with ❤ by GitHub
Here’s a Python implementation:
import random
import sys
random.seed()
numIterations = int(sys.argv[1])
count = 0.0
sum = 0.0
for i in range(numIterations):
sum += random.random()
if sum > 1.0:
count += 1.0
sum = 0.0
print numIterations/count
A naive ruby implementation
count_sum = 0 iterations = 1000 (1..iterations).each do sum = 0 count = 0 while sum < 1.0 do sum += rand count += 1 end count_sum += count end puts "#{count_sum.to_f / iterations }"Here’s a Java solution:
package progpraxis; import java.util.ArrayList; import java.util.List; import java.util.Random; public class E { public static void main(String[] args) { Random generator = new Random(); List<Integer> totalAttempts = new ArrayList<Integer>(); for (int i = 1; i <= 100000; i++) { double total = generator.nextDouble(); int attempts = 1; while (total <= 1.0) { total += generator.nextDouble(); attempts++; } totalAttempts.add(attempts); } double average = 0d; for (Integer number : totalAttempts) { average += number; } average = average / totalAttempts.size(); System.out.println("Average before surpassing 1: " + average); } }The output, as expected, is 2.72197.
(* OCaml *)
let test () =
let acc = ref 0. and cpt = ref 0 in
while !acc < 1. do
incr cpt;
acc := !acc +. (Random.float 1.);
done;
!cpt;;
let mean n =
let rec aux n =
if n = 0
then 0
else test () + aux (n-1)
in
(float_of_int (aux n)) /. (float_of_int n);;
print_string ((string_of_float (mean 10000))^”\n”);;
(* Purely functional looks even better (well, as in “no references”…) *)
let rec test acc cpt =
if acc < 1.
then test (acc +. (Random.float 1.)) (cpt + 1)
else cpt
let mean n =
let rec aux n =
if n = 0
then 0
else test 0. 0 + aux (n-1)
in (float_of_int (aux n)) /. (float_of_int n);;
print_string ((string_of_float (mean 10000))^”\n”);;
[…] about computing Euler’s Number, e, using Clojure. It wasn’t until I saw the idea on Programming Praxis, though, that I decided to just do […]
//quick and dirty javascript
var getRandomCount = function(ctr,buf){
while(buf < 1){
var nbr = Math.random();
ctr = ctr + 1;
buf = buf + nbr;
}
return ctr;
}
var getAverageNumberOfTimes = function(iter){
var total = 0;
for(i=0;i<iter;i++){
total = total + getRandomCount(0,0);
}
return total;
}
alert(getAverageNumberOfTimes(1000)/1000);
C#
Random randNum = new Random();
int trials = 1000000;
long totalTrials = 0;
double i = 0;
for (int runs = 0; runs < trials; i = 0, runs++)
for (totalTrials++ ; (i += randNum.NextDouble()) < 1.0; totalTrials++) ;
Console.WriteLine((float)totalTrials / trials);
Console.ReadLine();
Here is a FORTH version, without floating point.
A random 32 bit number is treated as binary random between 0-1. When overflow occurs (requires 33 bits), it is interpreted as exceeding 1.0.
\ multiply with carry RNG \ quickie MWC-0 (no lag, 2^63 period is good enough) create seed counter , create carry 48313 , : random ( -- r ) seed @ 4164903690 um* carry @ 0 d+ carry ! dup seed ! ; : trial ( -- n ) 0 0. ( stack: count rng ) BEGIN random 0 D+ rot 1+ -rot dup 0> UNTIL 2drop ; : trials ( n -- ) locals| n | 0 n 0 DO trial + LOOP 10000 n */ s>d <# # # # # [char] . hold #s #> type space ;And here is session with SwiftForth :-