Seven-sided dice from five-sided dice: Difference between revisions
(Update) |
(Updated D entry) |
||
| Line 420: | Line 420: | ||
/// Generates a random number in [1, 5]. |
/// Generates a random number in [1, 5]. |
||
int dice5() /*pure nothrow*/ { |
int dice5() /*pure nothrow*/ @safe { |
||
return uniform(1, 6); |
return uniform(1, 6); |
||
} |
} |
||
/// Naive, generates a random number in [1, 7] using dice5. |
/// Naive, generates a random number in [1, 7] using dice5. |
||
int fiveToSevenNaive() /*pure nothrow*/ { |
int fiveToSevenNaive() /*pure nothrow*/ @safe { |
||
immutable int r = dice5() + dice5() * 5 - 6; |
immutable int r = dice5() + dice5() * 5 - 6; |
||
return (r < 21) ? (r % 7) + 1 : fiveToSevenNaive(); |
return (r < 21) ? (r % 7) + 1 : fiveToSevenNaive(); |
||
| Line 434: | Line 434: | ||
minimizing calls to dice5. |
minimizing calls to dice5. |
||
*/ |
*/ |
||
int fiveToSevenSmart() { |
int fiveToSevenSmart() @safe { |
||
static int rem = 0, max = 1; |
static int rem = 0, max = 1; |
||
| Line 457: | Line 457: | ||
} |
} |
||
void main() { |
void main() /*@safe*/ { |
||
enum int N = 400_000; |
enum int N = 400_000; |
||
distCheck(&dice5, N, 1); |
distCheck(&dice5, N, 1); |
||
Revision as of 15:48, 5 January 2015
You are encouraged to solve this task according to the task description, using any language you may know.
Given an equal-probability generator of one of the integers 1 to 5
as dice5; create dice7 that generates a pseudo-random integer from
1 to 7 in equal probability using only dice5 as a source of random
numbers, and check the distribution for at least 1000000 calls using the function created in Simple Random Distribution Checker.
Implementation suggestion:
dice7 might call dice5 twice, re-call if four of the 25
combinations are given, otherwise split the other 21 combinations
into 7 groups of three, and return the group index from the rolls.
(Task adapted from an answer here)
Ada
The specification of a package Random_57: <lang Ada>package Random_57 is
type Mod_7 is mod 7;
function Random7 return Mod_7;
-- a "fast" implementation, minimazing the calls to the Random5 generator
function Simple_Random7 return Mod_7;
-- a simple implementation
end Random_57;</lang> Implementation of Random_57: <lang Ada> with Ada.Numerics.Discrete_Random;
package body Random_57 is
type M5 is mod 5;
package Rand_5 is new Ada.Numerics.Discrete_Random(M5);
Gen: Rand_5.Generator;
function Random7 return Mod_7 is
N: Natural;
begin
loop
N := Integer(Rand_5.Random(Gen))* 5 + Integer(Rand_5.Random(Gen));
-- N is uniformly distributed in 0 .. 24
if N < 21 then
return Mod_7(N/3);
else -- (N-21) is in 0 .. 3
N := (N-21) * 5 + Integer(Rand_5.Random(Gen)); -- N is in 0 .. 19
if N < 14 then
return Mod_7(N / 2);
else -- (N-14) is in 0 .. 5
N := (N-14) * 5 + Integer(Rand_5.Random(Gen)); -- N is in 0 .. 29
if N < 28 then
return Mod_7(N/4);
else -- (N-28) is in 0 .. 1
N := (N-28) * 5 + Integer(Rand_5.Random(Gen)); -- 0 .. 9
if N < 7 then
return Mod_7(N);
else -- (N-7) is in 0, 1, 2
N := (N-7)* 5 + Integer(Rand_5.Random(Gen)); -- 0 .. 14
if N < 14 then
return Mod_7(N/2);
else -- (N-14) is 0. This is not useful for us!
null;
end if;
end if;
end if;
end if;
end if;
end loop;
end Random7;
function Simple_Random7 return Mod_7 is
N: Natural :=
Integer(Rand_5.Random(Gen))* 5 + Integer(Rand_5.Random(Gen));
-- N is uniformly distributed in 0 .. 24
begin
while N > 20 loop
N := Integer(Rand_5.Random(Gen))* 5 + Integer(Rand_5.Random(Gen));
end loop; -- Now I <= 20
return Mod_7(N / 3);
end Simple_Random7;
begin
Rand_5.Reset(Gen);
end Random_57;</lang> A main program, using the Random_57 package: <lang Ada>with Ada.Text_IO, Random_57;
procedure R57 is
use Random_57;
type Fun is access function return Mod_7;
function Rand return Mod_7 renames Random_57.Random7; -- change this to "... renames Random_57.Simple_Random;" if you like
procedure Test(Sample_Size: Positive; Rand: Fun; Precision: Float := 0.3) is
Counter: array(Mod_7) of Natural := (others => 0);
Expected: Natural := Sample_Size/7;
Small: Mod_7 := Mod_7'First;
Large: Mod_7 := Mod_7'First;
Result: Mod_7;
begin
Ada.Text_IO.New_Line;
Ada.Text_IO.Put_Line("Sample Size: " & Integer'Image(Sample_Size));
Ada.Text_IO.Put( " Bins:");
for I in 1 .. Sample_Size loop
Result := Rand.all;
Counter(Result) := Counter(Result) + 1;
end loop;
for J in Mod_7 loop
Ada.Text_IO.Put(Integer'Image(Counter(J)));
if Counter(J) < Counter(Small) then Small := J; end if;
if Counter(J) > Counter(Large) then Large := J; end if;
end loop;
Ada.Text_IO.New_Line;
Ada.Text_IO.Put_Line(" Small Bin:" & Integer'Image(Counter(Small)));
Ada.Text_IO.Put_Line(" Large Bin: " & Integer'Image(Counter(Large)));
if Float(Counter(Small)*7) * (1.0+Precision) < Float(Sample_Size) then
Ada.Text_IO.Put_Line("Failed! Small too small!");
elsif Float(Counter(Large)*7) * (1.0-Precision) > Float(Sample_Size) then
Ada.Text_IO.Put_Line("Failed! Large too large!");
else
Ada.Text_IO.Put_Line("Passed");
end if;
end Test;
begin
Test( 10_000, Rand'Access, 0.08); Test( 100_000, Rand'Access, 0.04); Test( 1_000_000, Rand'Access, 0.02); Test(10_000_000, Rand'Access, 0.01);
end R57;</lang> A sample output:
Sample Size: 10000 Bins: 1368 1404 1435 1491 1483 1440 1379 Small Bin: 1368 Large Bin: 1491 Passed Sample Size: 100000 Bins: 14385 14110 14362 14404 14362 14206 14171 Small Bin: 14110 Large Bin: 14404 Passed Sample Size: 1000000 Bins: 143765 142384 142958 142684 142799 142956 142454 Small Bin: 142384 Large Bin: 143765 Passed Sample Size: 10000000 Bins: 1429266 1428214 1428753 1427032 1428418 1428699 1429618 Small Bin: 1427032 Large Bin: 1429618 Passed
ALGOL 68
- note: This specimen retains the original C coding style.
C's version using no multiplications, divisions, or mod operators: <lang algol68>PROC dice5 = INT:
1 + ENTIER (5*random);
PROC mulby5 = (INT n)INT:
ABS (BIN n SHL 2) + n;
PROC dice7 = INT: (
INT d55 := 0; INT m := 1; WHILE m := ABS ((2r1 AND BIN m) SHL 2) + ABS (BIN m SHR 1); # repeats 4 - 2 - 1 # d55 := mulby5(mulby5(d55)) + mulby5(dice5) + dice5 - 6;
- WHILE # d55 < m DO SKIP OD;
m := 1; WHILE d55>0 DO d55 +:= m; m := ABS (BIN d55 AND 2r111); # modulas by 8 # d55 := ABS (BIN d55 SHR 3) # divide by 8 # OD; m
);
PROC distcheck = (PROC INT dice, INT count, upb)VOID: (
[upb]INT sum; FOR i TO UPB sum DO sum[i] := 0 OD;
FOR i TO count DO sum[dice]+:=1 OD;
FOR i TO UPB sum WHILE print(whole(sum[i],0)); i /= UPB sum DO print(", ") OD;
print(new line)
);
main: (
distcheck(dice5, 1000000, 5); distcheck(dice7, 1000000, 7)
)</lang> Sample output:
200598, 199852, 199939, 200602, 199009 143529, 142688, 142816, 142747, 142958, 142802, 142460
AutoHotkey
<lang AutoHotkey>dice5() { Random, v, 1, 5
Return, v
}
dice7() { Loop
{ v := 5 * dice5() + dice5() - 6
IfLess v, 21, Return, (v // 3) + 1
}
}</lang>
Distribution check: Total elements = 10000 Margin = 3% --> Lbound = 1386, Ubound = 1471 Bucket 1 contains 1450 elements. Bucket 2 contains 1374 elements. Skewed. Bucket 3 contains 1412 elements. Bucket 4 contains 1465 elements. Bucket 5 contains 1370 elements. Skewed. Bucket 6 contains 1485 elements. Skewed. Bucket 7 contains 1444 elements.
BBC BASIC
<lang bbcbasic> MAXRND = 7
FOR r% = 2 TO 5
check% = FNdistcheck(FNdice7, 10^r%, 0.1)
PRINT "Over "; 10^r% " runs dice7 ";
IF check% THEN
PRINT "failed distribution check with "; check% " bin(s) out of range"
ELSE
PRINT "passed distribution check"
ENDIF
NEXT
END
DEF FNdice7
LOCAL x% : x% = FNdice5 + 5*FNdice5
IF x%>26 THEN = FNdice7 ELSE = (x%+1) MOD 7 + 1
DEF FNdice5 = RND(5)
DEF FNdistcheck(RETURN func%, repet%, delta)
LOCAL i%, m%, r%, s%, bins%()
DIM bins%(MAXRND)
FOR i% = 1 TO repet%
r% = FN(^func%)
bins%(r%) += 1
IF r%>m% m% = r%
NEXT
FOR i% = 1 TO m%
IF bins%(i%)/(repet%/m%) > 1+delta s% += 1
IF bins%(i%)/(repet%/m%) < 1-delta s% += 1
NEXT
= s%</lang>
Output:
Over 100 runs dice7 failed distribution check with 4 bin(s) out of range Over 1000 runs dice7 failed distribution check with 2 bin(s) out of range Over 10000 runs dice7 passed distribution check Over 100000 runs dice7 passed distribution check
C
<lang c>int rand5() { int r, rand_max = RAND_MAX - (RAND_MAX % 5); while ((r = rand()) >= rand_max); return r / (rand_max / 5) + 1; }
int rand5_7() { int r; while ((r = rand5() * 5 + rand5()) >= 27); return r / 3 - 1; }
int main() { printf(check(rand5, 5, 1000000, .05) ? "flat\n" : "not flat\n"); printf(check(rand7, 7, 1000000, .05) ? "flat\n" : "not flat\n"); return 0; }</lang>output
flat flat
C++
This solution tries to minimize calls to the underlying d5 by reusing information from earlier calls. <lang cpp>template<typename F> class fivetoseven { public:
fivetoseven(F f): d5(f), rem(0), max(1) {}
int operator()();
private:
F d5; int rem, max;
};
template<typename F>
int fivetoseven<F>::operator()()
{
while (rem/7 == max/7)
{
while (max < 7)
{
int rand5 = d5()-1;
max *= 5;
rem = 5*rem + rand5;
}
int groups = max / 7;
if (rem >= 7*groups)
{
rem -= 7*groups;
max -= 7*groups;
}
}
int result = rem % 7; rem /= 7; max /= 7; return result+1;
}
int d5() {
return 5.0*std::rand()/(RAND_MAX + 1.0) + 1;
}
fivetoseven<int(*)()> d7(d5);
int main() {
srand(time(0)); test_distribution(d5, 1000000, 0.001); test_distribution(d7, 1000000, 0.001);
}</lang>
Clojure
Uses the verify function defined in Verify distribution uniformity/Naive#Clojure <lang Clojure>(def dice5 #(rand-int 5))
(defn dice7 []
(quot (->> dice5 ; do the following to dice5
(repeatedly 2) ; call it twice
(apply #(+ %1 (* 5 %2))) ; d1 + 5*d2 => 0..24
#() ; wrap that up in a function
repeatedly ; make infinite sequence of the above
(drop-while #(> % 20)) ; throw away anything > 20
first) ; grab first acceptable element
3)) ; divide by three rounding down
(doseq [n [100 1000 10000] [num count okay?] (verify dice7 n)]
(println "Saw" num count "times:"
(if okay? "that's" " not") "acceptable"))</lang>
Saw 0 10 times: not acceptable Saw 1 19 times: not acceptable Saw 2 12 times: not acceptable Saw 3 15 times: that's acceptable Saw 4 11 times: not acceptable Saw 5 11 times: not acceptable Saw 6 22 times: not acceptable Saw 0 142 times: that's acceptable Saw 1 158 times: not acceptable Saw 2 151 times: that's acceptable Saw 3 153 times: that's acceptable Saw 4 118 times: not acceptable Saw 5 139 times: that's acceptable Saw 6 139 times: that's acceptable Saw 0 1498 times: that's acceptable Saw 1 1411 times: that's acceptable Saw 2 1436 times: that's acceptable Saw 3 1434 times: that's acceptable Saw 4 1414 times: that's acceptable Saw 5 1408 times: that's acceptable Saw 6 1399 times: that's acceptable
Common Lisp
<lang lisp>(defun d5 ()
(1+ (random 5)))
(defun d7 ()
(loop for d55 = (+ (* 5 (d5)) (d5) -6)
until (< d55 21)
finally (return (1+ (mod d55 7)))))</lang>
> (check-distribution 'd7 1000)
Distribution potentially skewed for 1: expected around 1000/7 got 153.
Distribution potentially skewed for 2: expected around 1000/7 got 119.
Distribution potentially skewed for 3: expected around 1000/7 got 125.
Distribution potentially skewed for 7: expected around 1000/7 got 156.
T
#<EQL Hash Table{7} 200B5A53>
> (check-distribution 'd7 10000)
NIL
#<EQL Hash Table{7} 200CB5BB>
D
<lang d>import std.random; import verify_distribution_uniformity_naive: distCheck;
/// Generates a random number in [1, 5]. int dice5() /*pure nothrow*/ @safe {
return uniform(1, 6);
}
/// Naive, generates a random number in [1, 7] using dice5. int fiveToSevenNaive() /*pure nothrow*/ @safe {
immutable int r = dice5() + dice5() * 5 - 6; return (r < 21) ? (r % 7) + 1 : fiveToSevenNaive();
}
/** Generates a random number in [1, 7] using dice5, minimizing calls to dice5.
- /
int fiveToSevenSmart() @safe {
static int rem = 0, max = 1;
while (rem / 7 == max / 7) {
while (max < 7) {
immutable int rand5 = dice5() - 1;
max *= 5;
rem = 5 * rem + rand5;
}
immutable int groups = max / 7;
if (rem >= 7 * groups) {
rem -= 7 * groups;
max -= 7 * groups;
}
}
immutable int result = rem % 7; rem /= 7; max /= 7; return result + 1;
}
void main() /*@safe*/ {
enum int N = 400_000; distCheck(&dice5, N, 1); distCheck(&fiveToSevenNaive, N, 1); distCheck(&fiveToSevenSmart, N, 1);
}</lang>
- Output:
1 80365 2 79941 3 80065 4 79784 5 79845 1 57186 2 57201 3 57180 4 57231 5 57124 6 56832 7 57246 1 57367 2 56869 3 57644 4 57111 5 57157 6 56809 7 57043
E
<lang e>def dice5() {
return entropy.nextInt(5) + 1
}
def dice7() {
var d55 := null
while ((d55 := 5 * dice5() + dice5() - 6) >= 21) {}
return d55 %% 7 + 1
}</lang> <lang e>def bins := ([0] * 7).diverge() for x in 1..1000 {
bins[dice7() - 1] += 1
} println(bins.snapshot())</lang>
Erlang
<lang Erlang> -module( dice ).
-export( [dice5/0, dice7/0, task/0] ).
dice5() -> random:uniform( 5 ).
dice7() -> dice7_small_enough( dice5() * 5 + dice5() - 6 ). % 0 - 24
task() ->
verify_distribution_uniformity:naive( fun dice7/0, 1000000, 1 ).
dice7_small_enough( N ) when N < 21 -> N div 3 + 1; dice7_small_enough( _N ) -> dice7(). </lang>
- Output:
76> dice:task(). ok
Fortran
<lang fortran>module rand_mod
implicit none
contains
function rand5()
integer :: rand5 real :: r
call random_number(r) rand5 = 5*r + 1
end function
function rand7()
integer :: rand7
do
rand7 = 5*rand5() + rand5() - 6
if (rand7 < 21) then
rand7 = rand7 / 3 + 1
return
end if
end do
end function end module
program Randtest
use rand_mod implicit none integer, parameter :: samples = 1000000 call distcheck(rand7, samples, 0.005) write(*,*) call distcheck(rand7, samples, 0.001)
end program</lang> Output
Distribution Uniform Distribution potentially skewed for bucket 1 Expected: 142857 Actual: 143142 Distribution potentially skewed for bucket 2 Expected: 142857 Actual: 143454 Distribution potentially skewed for bucket 3 Expected: 142857 Actual: 143540 Distribution potentially skewed for bucket 4 Expected: 142857 Actual: 142677 Distribution potentially skewed for bucket 5 Expected: 142857 Actual: 142511 Distribution potentially skewed for bucket 6 Expected: 142857 Actual: 142163 Distribution potentially skewed for bucket 7 Expected: 142857 Actual: 142513
Go
<lang go>package main
import (
"fmt" "math" "math/rand" "time"
)
// "given" func dice5() int {
return rand.Intn(5) + 1
}
// function specified by task "Seven-sided dice from five-sided dice" func dice7() (i int) {
for {
i = 5*dice5() + dice5()
if i < 27 {
break
}
}
return (i / 3) - 1
}
// function specified by task "Verify distribution uniformity/Naive" // // Parameter "f" is expected to return a random integer in the range 1..n. // (Values out of range will cause an unceremonious crash.) // "Max" is returned as an "indication of distribution achieved." // It is the maximum delta observed from the count representing a perfectly // uniform distribution. // Also returned is a boolean, true if "max" is less than threshold // parameter "delta." func distCheck(f func() int, n int,
repeats int, delta float64) (max float64, flatEnough bool) {
count := make([]int, n)
for i := 0; i < repeats; i++ {
count[f()-1]++
}
expected := float64(repeats) / float64(n)
for _, c := range count {
max = math.Max(max, math.Abs(float64(c)-expected))
}
return max, max < delta
}
// Driver, produces output satisfying both tasks. func main() {
rand.Seed(time.Now().UnixNano())
const calls = 1000000
max, flatEnough := distCheck(dice7, 7, calls, 500)
fmt.Println("Max delta:", max, "Flat enough:", flatEnough)
max, flatEnough = distCheck(dice7, 7, calls, 500)
fmt.Println("Max delta:", max, "Flat enough:", flatEnough)
}</lang> Output:
Max delta: 356.1428571428696 Flat enough: true Max delta: 787.8571428571304 Flat enough: false
Groovy
<lang groovy>random = new Random()
int rand5() {
random.nextInt(5) + 1
}
int rand7From5() {
def raw = 25
while (raw > 21) {
raw = 5*(rand5() - 1) + rand5()
}
(raw % 7) + 1
}</lang> Test: <lang groovy>def test = {
(1..6). each {
def counts = [0g, 0g, 0g, 0g, 0g, 0g, 0g]
def target = 10g**it
def popSize = 7*target
(0..<(popSize)).each {
def i = rand7From5() - 1
counts[i] = counts[i] + 1g
}
BigDecimal stdDev = (counts.collect { it - target}.collect { it * it }.sum() / popSize) ** 0.5g
def countMap = (0..<counts.size()).inject([:]) { map, index -> map + [(index+1):counts[index]] }
println """\
counts: ${countMap}
population size: ${popSize}
std dev: ${stdDev.round(new java.math.MathContext(3))}
"""
}
}
4.times {
println """
TRIAL #${it+1} =============="""
test(it)
}</lang> Output:
TRIAL #1
==============
counts: [1:16, 2:10, 3:9, 4:7, 5:12, 6:8, 7:8]
population size: 70
std dev: 0.910
counts: [1:85, 2:97, 3:108, 4:110, 5:95, 6:105, 7:100]
population size: 700
std dev: 0.800
counts: [1:990, 2:1008, 3:992, 4:1060, 5:1008, 6:997, 7:945]
population size: 7000
std dev: 0.995
counts: [1:9976, 2:10007, 3:10009, 4:9858, 5:10109, 6:9988, 7:10053]
population size: 70000
std dev: 0.714
counts: [1:100310, 2:99783, 3:99843, 4:100353, 5:99804, 6:99553, 7:100354]
population size: 700000
std dev: 0.968
counts: [1:999320, 2:1000942, 3:1000201, 4:1000878, 5:999181, 6:999632, 7:999846]
population size: 7000000
std dev: 0.654
TRIAL #2
==============
counts: [1:10, 2:8, 3:9, 4:9, 5:14, 6:7, 7:13]
population size: 70
std dev: 0.756
counts: [1:104, 2:101, 3:97, 4:108, 5:100, 6:87, 7:103]
population size: 700
std dev: 0.619
counts: [1:995, 2:970, 3:1001, 4:953, 5:1006, 6:1081, 7:994]
population size: 7000
std dev: 1.18
counts: [1:10013, 2:10063, 3:9843, 4:9984, 5:9986, 6:10059, 7:10052]
population size: 70000
std dev: 0.711
counts: [1:100048, 2:99647, 3:100240, 4:100683, 5:99813, 6:100320, 7:99249]
population size: 700000
std dev: 1.39
counts: [1:1000579, 2:1000541, 3:999497, 4:1000805, 5:999708, 6:999161, 7:999709]
population size: 7000000
std dev: 0.586
TRIAL #3
==============
counts: [1:9, 2:8, 3:11, 4:14, 5:10, 6:11, 7:7]
population size: 70
std dev: 0.676
counts: [1:100, 2:92, 3:105, 4:107, 5:111, 6:91, 7:94]
population size: 700
std dev: 0.733
counts: [1:1010, 2:1053, 3:967, 4:981, 5:1027, 6:959, 7:1003]
population size: 7000
std dev: 0.984
counts: [1:9857, 2:10037, 3:9992, 4:10231, 5:9828, 6:10140, 7:9915]
population size: 70000
std dev: 1.37
counts: [1:99650, 2:99580, 3:99848, 4:100507, 5:99916, 6:100212, 7:100287]
population size: 700000
std dev: 1.01
counts: [1:1001710, 2:999667, 3:1000685, 4:1000411, 5:999369, 6:998469, 7:999689]
population size: 7000000
std dev: 0.965
TRIAL #4
==============
counts: [1:12, 2:7, 3:11, 4:12, 5:7, 6:9, 7:12]
population size: 70
std dev: 0.676
counts: [1:97, 2:96, 3:101, 4:93, 5:96, 6:124, 7:93]
population size: 700
std dev: 1.01
counts: [1:985, 2:1023, 3:1018, 4:1023, 5:995, 6:973, 7:983]
population size: 7000
std dev: 0.615
counts: [1:9948, 2:9968, 3:10131, 4:10050, 5:9990, 6:10039, 7:9874]
population size: 70000
std dev: 0.764
counts: [1:100125, 2:99616, 3:99912, 4:100286, 5:99674, 6:100190, 7:100197]
population size: 700000
std dev: 0.787
counts: [1:1001267, 2:999911, 3:1000602, 4:999483, 5:1000549, 6:998725, 7:999463]
population size: 7000000
std dev: 0.798
Haskell
<lang haskell>import System.Random import Data.List
sevenFrom5Dice = do
d51 <- randomRIO(1,5) :: IO Int
d52 <- randomRIO(1,5) :: IO Int
let d7 = 5*d51+d52-6
if d7 > 20 then sevenFrom5Dice
else return $ 1 + d7 `mod` 7</lang>
Output: <lang haskell>*Main> replicateM 10 sevenFrom5Dice [2,3,1,1,6,2,5,6,5,3]</lang> Test: <lang haskell>*Main> mapM_ print .sort =<< distribCheck sevenFrom5Dice 1000000 3 (1,(142759,True)) (2,(143078,True)) (3,(142706,True)) (4,(142403,True)) (5,(142896,True)) (6,(143028,True)) (7,(143130,True))</lang>
Icon and Unicon
Uses verify_uniform from here.
<lang Icon>
$include "distribution-checker.icn"
- return a uniformly distributed number from 1 to 7,
- but only using a random number in range 1 to 5.
procedure die_7 ()
while rnd := 5*?5 + ?5 - 6 do {
if rnd < 21 then suspend rnd % 7 + 1
}
end
procedure main ()
if verify_uniform (create (|die_7()), 1000000, 0.01)
then write ("uniform")
else write ("skewed")
end </lang>
Output:
5 142870 2 142812 7 142901 4 142960 1 143113 6 142706 3 142638 uniform
J
The first step is to create 7-sided dice rolls from 5-sided dice rolls (rollD5):
<lang j>rollD5=: [: >: ] ?@$ 5: NB. makes a y shape array of 5s, "rolls" the array and increments.
roll2xD5=: [: rollD5 2 ,~ */ NB. rolls D5 twice for each desired D7 roll (y rows, 2 cols)
toBase10=: 5 #. <: NB. decrements and converts rows from base 5 to 10
keepGood=: #~ 21&> NB. compress out values not less than 21
groupin3s=: [: >. >: % 3: NB. increments, divides by 3 and takes ceiling
getD7=: groupin3s@keepGood@toBase10@roll2xD5</lang> Here are a couple of variations on the theme that achieve the same result: <lang j>getD7b=: 0 8 -.~ 3 >.@%~ 5 #. [: <:@rollD5 2 ,~ ] getD7c=: [: (#~ 7&>:) 3 >.@%~ [: 5&#.&.:<:@rollD5 ] , 2:</lang> The trouble is that we probably don't have enough D7 rolls yet because we compressed out any double D5 rolls that evaluated to 21 or more. So we need to accumulate some more D7 rolls until we have enough. J has two types of verb definition - tacit (arguments not referenced) and explicit (more conventional function definitions) illustrated below:
Here's an explicit definition that accumulates rolls from getD7:
<lang j>rollD7x=: monad define
n=. */y NB. product of vector y is total number of D7 rolls required rolls=. NB. initialize empty noun rolls while. n > #res do. NB. checks if if enough D7 rolls accumulated rolls=. rolls, getD7 >. 0.75 * n NB. calcs 3/4 of required rolls and accumulates getD7 rolls end. y $ rolls NB. shape the result according to the vector y
)</lang> Here's a tacit definition that does the same thing: <lang j>getNumRolls=: [: >. 0.75 * */@[ NB. calc approx 3/4 of the required rolls accumD7Rolls=: ] , getD7@getNumRolls NB. accumulates getD7 rolls isNotEnough=: */@[ > #@] NB. checks if enough D7 rolls accumulated
rollD7t=: ] $ (accumD7Rolls ^: isNotEnough ^:_)&</lang>
The verb1 ^: verb2 ^:_ construct repeats x verb1 y while x verb2 y is true. It is like saying "Repeat accumD7Rolls while isNotEnough".
Example usage: <lang j> rollD7t 10 NB. 10 rolls of D7 6 4 5 1 4 2 4 5 2 5
rollD7t 2 5 NB. 2 by 5 array of D7 rolls
5 1 5 1 3 3 4 3 5 6
rollD7t 2 3 5 NB. 2 by 3 by 5 array of D7 rolls
4 7 7 5 7 3 7 1 4 5 5 4 5 7 6
1 1 7 6 3 4 4 1 4 4 1 1 1 6 5
NB. check results from rollD7x and rollD7t have same shape
($@rollD7x -: $@rollD7t) 10
1
($@rollD7x -: $@rollD7t) 2 3 5
1</lang>
Java
<lang Java>import java.util.Random; public class SevenSidedDice { private static final Random rnd = new Random(); public static void main(String[] args) { SevenSidedDice now=new SevenSidedDice(); System.out.println("Random number from 1 to 7: "+now.seven()); } int seven() { int v=21; while(v>20) v=five()+five()*5-6; return 1+v%7; } int five() { return 1+rnd.nextInt(5); } }</lang>
JavaScript
<lang javascript>function dice5() {
return 1 + Math.floor(5 * Math.random());
}
function dice7() {
while (true)
{
var dice55 = 5 * dice5() + dice5() - 6;
if (dice55 < 21)
return dice55 % 7 + 1;
}
}
distcheck(dice5, 1000000); print(); distcheck(dice7, 1000000);</lang> output
1 199792 2 200425 3 199243 4 200407 5 200133 1 143617 2 142209 3 143023 4 142990 5 142894 6 142648 7 142619
Julia
<lang Julia>dice5() = rand(1:5)
function dice7()
r = 5*dice5() + dice5() - 6 r < 21 ? (r%7 + 1) : dice7()
end</lang> Distribution check:
julia> hist([dice5() for i=1:10^6]) (0:1:5,[199932,200431,199969,199925,199743]) julia> hist([dice7() for i=1:10^6]) (0:1:7,[142390,143032,142837,142999,142800,142642,143300])
Lua
<lang lua>dice5 = function() return math.random(5) end
function dice7()
x = dice5() * 5 + dice5() - 6 if x > 20 then return dice7() end return x%7 + 1
end</lang>
Mathematica
<lang Mathematica>sevenFrom5Dice := (tmp$ = 5*RandomInteger[{1, 5}] + RandomInteger[{1, 5}] - 6;
If [tmp$ < 21, 1 + Mod[tmp$ , 7], sevenFrom5Dice])</lang>
CheckDistribution[sevenFrom5Dice, 1000000, 5]
->Expected: 142857., Generated :{142206,142590,142650,142693,142730,143475,143656}
->"Flat"
OCaml
<lang ocaml>let dice5() = 1 + Random.int 5 ;;
let dice7 =
let rolls2answer = Hashtbl.create 25 in
let n = ref 0 in
for roll1 = 1 to 5 do
for roll2 = 1 to 5 do
Hashtbl.add rolls2answer (roll1,roll2) (!n / 3 +1);
incr n
done;
done;
let rec aux() =
let trial = Hashtbl.find rolls2answer (dice5(),dice5()) in
if trial <= 7 then trial else aux()
in
aux
- </lang>
PARI/GP
<lang parigp>dice5()=random(5)+1;
dice7()={
my(t); while((t=dice5()*5+dice5()) > 21,); (t+2)\3
};</lang>
Perl
Using dice5 twice to generate numbers in the range 0 to 24. If we consider these modulo 8 and re-call if we get zero, we have eliminated 4 cases and created the necessary number in the range from 1 to 7. <lang perl>sub dice5 { 1+int rand(5) }
sub dice7 {
while(1) {
my $d7 = (5*dice5()+dice5()-6) % 8;
return $d7 if $d7;
}
}
my %count7; my $n = 1000000; $count7{dice7()}++ for 1..$n; printf "%s: %5.2f%%\n", $_, 100*($count7{$_}/$n*7-1) for sort keys %count7; </lang>
- Output:
1: 0.05% 2: 0.16% 3: -0.43% 4: 0.11% 5: 0.01% 6: -0.15% 7: 0.24%
Perl 6
Since rakudo is still pretty slow, we've done some interesting bits of optimization. We factor out the range object construction so that it doesn't have to be recreated each time, and we sneakily subtract the 1's from the 5's, which takes us back to 0 based without having to subtract 6. <lang perl6>my $d5 = 1..5; sub d5() { $d5.roll; } # 1d5
sub d7() {
my $flat = 21; $flat = 5 * d5() - d5() until $flat < 21; $flat % 7 + 1;
}</lang>
Here's the test. We use a C-style for loop, except it's named loop, because it's currently faster than the other loops--and, er, doesn't segfault the GC on a million iterations...
<lang perl6>my @dist;
my $n = 1_000_000;
my $expect = $n / 7;
loop ($_ = $n; $n; --$n) { @dist[d7()]++; }
say "Expect\t",$expect.fmt("%.3f"); for @dist.kv -> $i, $v {
say "$i\t$v\t" ~ (($v - $expect)/$expect*100).fmt("%+.2f%%") if $v;
}</lang> And the output:
Expect 142857.143 1 142835 -0.02% 2 143021 +0.11% 3 142771 -0.06% 4 142706 -0.11% 5 143258 +0.28% 6 142485 -0.26% 7 142924 +0.05%
PicoLisp
<lang PicoLisp>(de dice5 ()
(rand 1 5) )
(de dice7 ()
(use R
(until (> 21 (setq R (+ (* 5 (dice5)) (dice5) -6))))
(inc (% R 7)) ) )</lang>
Output:
: (let R NIL (do 1000000 (accu 'R (dice7) 1)) (sort R) ) -> ((1 . 142295) (2 . 142491) (3 . 143448) (4 . 143129) (5 . 142701) (6 . 143142) (7 . 142794))
PureBasic
<lang PureBasic>Procedure dice5()
ProcedureReturn Random(4) + 1
EndProcedure
Procedure dice7()
Protected x x = dice5() * 5 + dice5() - 6 If x > 20 ProcedureReturn dice7() EndIf ProcedureReturn x % 7 + 1
EndProcedure</lang>
Python
<lang python>from random import randint
def dice5():
return randint(1, 5)
def dice7():
r = dice5() + dice5() * 5 - 6 return (r % 7) + 1 if r < 21 else dice7()</lang>
Distribution check using Simple Random Distribution Checker:
>>> distcheck(dice5, 1000000, 1)
{1: 200244, 2: 199831, 3: 199548, 4: 199853, 5: 200524}
>>> distcheck(dice7, 1000000, 1)
{1: 142853, 2: 142576, 3: 143067, 4: 142149, 5: 143189, 6: 143285, 7: 142881}
Racket
<lang Racket>
- lang racket
(define (dice5) (add1 (random 5)))
(define (dice7)
(define res (+ (* 5 (dice5)) (dice5) -6)) (if (< res 21) (+ 1 (modulo res 7)) (dice7)))
</lang>
Checking the uniformity using math library:
<lang racket> -> (require math/statistics) -> (samples->hash (for/list ([i 700000]) (dice7))) '#hash((7 . 100392)
(6 . 100285)
(5 . 99774)
(4 . 100000)
(3 . 100000)
(2 . 99927)
(1 . 99622))
</lang>
R
5-sided die. <lang r>dice5 <- function(n=1) sample(5, n, replace=TRUE)</lang> Simple but slow 7-sided die, using a while loop. <lang r>dice7.while <- function(n=1) {
score <- numeric()
while(length(score) < n)
{
total <- sum(c(5,1) * dice5(2)) - 3
if(total < 24) score <- c(score, total %/% 3)
}
score
} system.time(dice7.while(1e6)) # longer than 4 minutes</lang> More complex, but much faster vectorised version. <lang r>dice7.vec <- function(n=1, checkLength=TRUE) {
morethan2n <- 3 * n + 10 + (n %% 2) #need more than 2*n samples, because some are discarded
twoDfive <- matrix(dice5(morethan2n), nrow=2)
total <- colSums(c(5, 1) * twoDfive) - 3
score <- ifelse(total < 24, total %/% 3, NA)
score <- score[!is.na(score)]
#If length is less than n (very unlikely), add some more samples
if(checkLength)
{
while(length(score) < n)
{
score <- c(score, dice7(n, FALSE))
}
score[1:n]
} else score
} system.time(dice7.vec(1e6)) # ~1 sec</lang>
REXX
<lang rexx>/*REXX program to simulate 7-sided die based on a 5-sided throw. */ parse arg trials sample . /*get arguments from command line*/ if trials== then trials=1 /*Not specified? Then use default*/ if sample== then sample=1000000 /* " " " " " */
do t=1 for trials /*performe the number of trials. */
die.=0; k=0
do until k==sample; r=5*random(1,5)+random(1,5)-6
if r>20 then iterate
k=k+1; r=r//7+1; die.r=die.r+1
end /*until*/
expect=sample%7
say
say center('trial:'right(t,4) ' ' sample 'samples, expect='expect,79,'─')
say
do j=1 for 7
say ' side' j "had " die.j ' occurences',
' difference from expected:'right(die.j-expect,length(sample))
end /*j*/
end /*t*/
/*stick a fork in it, we're done.*/</lang>
output when the input is specified as 11:
─────────────────trial: 1 1000000 samples, expect=142857─────────────────
side 1 had 142879 occurences difference from expected: 22
side 2 had 142980 occurences difference from expected: 123
side 3 had 142366 occurences difference from expected: -491
side 4 had 143129 occurences difference from expected: 272
side 5 had 143316 occurences difference from expected: 459
side 6 had 142742 occurences difference from expected: -115
side 7 had 142588 occurences difference from expected: -269
─────────────────trial: 2 1000000 samples, expect=142857─────────────────
side 1 had 142583 occurences difference from expected: -274
side 2 had 142797 occurences difference from expected: -60
side 3 had 143093 occurences difference from expected: 236
side 4 had 142168 occurences difference from expected: -689
side 5 had 143203 occurences difference from expected: 346
side 6 had 143260 occurences difference from expected: 403
side 7 had 142896 occurences difference from expected: 39
─────────────────trial: 3 1000000 samples, expect=142857─────────────────
side 1 had 143089 occurences difference from expected: 232
side 2 had 142815 occurences difference from expected: -42
side 3 had 142645 occurences difference from expected: -212
side 4 had 142959 occurences difference from expected: 102
side 5 had 142764 occurences difference from expected: -93
side 6 had 143205 occurences difference from expected: 348
side 7 had 142523 occurences difference from expected: -334
─────────────────trial: 4 1000000 samples, expect=142857─────────────────
side 1 had 142535 occurences difference from expected: -322
side 2 had 142316 occurences difference from expected: -541
side 3 had 143101 occurences difference from expected: 244
side 4 had 143159 occurences difference from expected: 302
side 5 had 143130 occurences difference from expected: 273
side 6 had 142551 occurences difference from expected: -306
side 7 had 143208 occurences difference from expected: 351
─────────────────trial: 5 1000000 samples, expect=142857─────────────────
side 1 had 143335 occurences difference from expected: 478
side 2 had 142296 occurences difference from expected: -561
side 3 had 142500 occurences difference from expected: -357
side 4 had 142775 occurences difference from expected: -82
side 5 had 142581 occurences difference from expected: -276
side 6 had 143363 occurences difference from expected: 506
side 7 had 143150 occurences difference from expected: 293
─────────────────trial: 6 1000000 samples, expect=142857─────────────────
side 1 had 142593 occurences difference from expected: -264
side 2 had 143012 occurences difference from expected: 155
side 3 had 142703 occurences difference from expected: -154
side 4 had 143702 occurences difference from expected: 845
side 5 had 142867 occurences difference from expected: 10
side 6 had 142158 occurences difference from expected: -699
side 7 had 142965 occurences difference from expected: 108
─────────────────trial: 7 1000000 samples, expect=142857─────────────────
side 1 had 143065 occurences difference from expected: 208
side 2 had 142892 occurences difference from expected: 35
side 3 had 142681 occurences difference from expected: -176
side 4 had 142637 occurences difference from expected: -220
side 5 had 142563 occurences difference from expected: -294
side 6 had 142934 occurences difference from expected: 77
side 7 had 143228 occurences difference from expected: 371
─────────────────trial: 8 1000000 samples, expect=142857─────────────────
side 1 had 142745 occurences difference from expected: -112
side 2 had 142906 occurences difference from expected: 49
side 3 had 142894 occurences difference from expected: 37
side 4 had 142976 occurences difference from expected: 119
side 5 had 142265 occurences difference from expected: -592
side 6 had 142909 occurences difference from expected: 52
side 7 had 143305 occurences difference from expected: 448
─────────────────trial: 9 1000000 samples, expect=142857─────────────────
side 1 had 142320 occurences difference from expected: -537
side 2 had 142510 occurences difference from expected: -347
side 3 had 142867 occurences difference from expected: 10
side 4 had 143711 occurences difference from expected: 854
side 5 had 142732 occurences difference from expected: -125
side 6 had 143115 occurences difference from expected: 258
side 7 had 142745 occurences difference from expected: -112
─────────────────trial: 10 1000000 samples, expect=142857─────────────────
side 1 had 142725 occurences difference from expected: -132
side 2 had 142865 occurences difference from expected: 8
side 3 had 143076 occurences difference from expected: 219
side 4 had 142759 occurences difference from expected: -98
side 5 had 142593 occurences difference from expected: -264
side 6 had 142750 occurences difference from expected: -107
side 7 had 143232 occurences difference from expected: 375
─────────────────trial: 11 1000000 samples, expect=142857─────────────────
side 1 had 142825 occurences difference from expected: -32
side 2 had 142870 occurences difference from expected: 13
side 3 had 142919 occurences difference from expected: 62
side 4 had 143098 occurences difference from expected: 241
side 5 had 142627 occurences difference from expected: -230
side 6 had 143163 occurences difference from expected: 306
side 7 had 142498 occurences difference from expected: -359
Ruby
Uses distcheck from here.
<lang ruby>require './distcheck.rb'
def d5
1 + rand(5)
end
def d7
loop do d55 = 5*d5 + d5 - 6 return (d55 % 7 + 1) if d55 < 21 end
end
distcheck(1_000_000) {d5} distcheck(1_000_000) {d7}</lang>
- Output:
1 200227 2 200264 3 199777 4 199387 5 200345 1 143175 2 143031 3 142731 4 142716 5 142931 6 142605 7 142811
Tcl
Any old D&D hand will know these as a D5 and a D7... <lang tcl>proc D5 {} {expr {1 + int(5 * rand())}}
proc D7 {} {
while 1 {
set d55 [expr {5 * [D5] + [D5] - 6}]
if {$d55 < 21} {
return [expr {$d55 % 7 + 1}]
}
}
}</lang> Checking:
% distcheck D5 1000000 1 199893 2 200162 3 200075 4 199630 5 200240 % distcheck D7 1000000 1 143121 2 142383 3 143353 4 142811 5 142172 6 143291 7 142869
VBScript
<lang vb>Option Explicit
function dice5 dice5 = int(rnd*5) + 1 end function
function dice7 dim j do j = 5 * dice5 + dice5 - 6 loop until j < 21 dice7 = j mod 7 + 1 end function</lang>
zkl
<lang zkl>var die5=(1).random.fp(6); // [1..5] fcn die7{ while((r:=5*die5() + die5())>=27){} r/3-1 }
fcn rtest(N){ //test spread over [0..9]
dist:=L(0,0,0,0,0,0,0,0,0,0);
do(N){ dist[die7()]+=1 }
sum:=dist.sum();
dist=dist.apply('wrap(n){ "%.2f%%".fmt(n.toFloat()/sum*100) }).println();
}
println("Looking for ",100.0/7,"%"); rtest(0d1_000_000);</lang>
- Output:
Looking for 14.2857%
L("0.00%","14.28%","14.36%","14.22%","14.26%","14.34%","14.33%","14.21%","0.00%","0.00%")