/* amoeba.tc - tct - 9/14/12 */

/* We use the pseudorandom number generator seen below instead of the one */
/* in the standard library which seems to be OK only for the DOS */
/* (16-bit integer) interpreter. */

int gseed

randint [
 int k1,ix;ix=gseed
 k1=ix/127773
 ix=16807*(ix-k1*127773)-k1*2836
 if (ix<0) ix=ix+2147483647
 gseed=ix
 ]

amoeba [
 pl "amoeba.tc - tct - 9/14/12";pl ""
 pl "An amoeba either dies, lives, splits into two or splits into three,"
 pl "each with probability 1/4. It can be proved that out of 100"
 pl "simulations, about 100*(sqrt(2)-1), which equals about 41,"
 pl "populations will die out. It can also be proved that 1.5^n"
 pl "is the expected number of amoebas after n generations. If n=10,"
 pl "58 amoebas are expected.";pl ""
 pl "Ref: http://primepuzzle.com/tc/amoeba.html"
 int c,p,pp,psp,s,gens,n,ptot
 pl "";pl "seed? ";gseed=gn
 gens=10
 for (s=1;s<=100;++s) [
  p=1;n=gens
  while (n*p) [
   --n;psp=0
   for (pp=1;pp<=p;++pp) [
    randint;psp=psp+gseed%4
    ]
   p=psp;if (!p) ++c
   ]
  ptot=ptot+p;pl "simulation ";pn s;ps ": population";pn p
  ]
 pl "";pl "";pn ptot/100;ps " amoebas, on average, were living after"
 pn gens;ps " generations.";pl ""
 pn c;ps " population(s) died out.";pl ""
 ]