/* sox.tc - tct - 3/20/9

dpf int n,d,p [
  /* display proper fraction n/d as decimal number
  /* with p digit precision 
  int w
  w=10*n/d
  if (p) [
    MC w,14
    dpf 10*n%d,d,p-1 /* recursive call
    ]
  ]

int last

rand [ /* returns random number in 1..2^31-1
int m;m=2147483647
int a;a=16807	
last=last*a%m
if (last<0) last=-last
return last
]	

/* this was motivated by DeeDee's post
/* http://tobee-interpres.blogspot.com/2009/03/sox.html

sox [
pl "";ps "sox.tc - 3/20/9 - tct";pl "";pl ""
ps "seed (e.g. 314) ";last=gn
pl "number of feet";pl ""
pl "2 (huminds), 4 (kitehs), 6 (flize),"
pl "8 (spidees), 100 (centapeeds) etc. "
int pair;pair=gn
pl "number of simulations (e.g. 1, 100, 200, 2000, 15000 etc.) "
int ul;ul=gn
pl ""
int i,j,xbar
int w,b,p,m /* white, black, "pairs," max
for (j=1;j<=ul;++j) [
 p=w=b=0
 for (i=1;i<=3*pair-1;++i) [
  if (rand()>1073741823) [
   ++w;if (w==pair) [++p;w=0]
   ]
  else [
   ++b;if (b==pair) [++p;b=0]
   ]
  if (p==2) [ /* got 2 "pair"
   if (i>m) [m=i;pn m] /* update, show new max
   /* update xbar (xbar/ul will equal average number of pulls required to get 2 "pair")
   xbar=xbar+i
   break
   ]
  ]
 ]
pl ""
pl "total number of sox pulled to get 2 "
putchar '"';ps "pair";putchar '"'
pl "in ";MC ul,14;ps " simulations was";pn xbar;ps ",";pl ""
ps "giving an expected value of"
pn xbar/ul;ps ".";dpf xbar%ul,ul,4;ps " pulls.";pl ""
]