Mathematical evolution experiment in PARI/GP

evolution.pari

\\ An Evolution Game
\\
\\ Implemented by @promovicz on 2023-11-25.
\\
\\ Inspired by Martin Escardo on Mastodon:
\\   https://mathstodon.xyz/@MartinEscardo/111467732542708261
\\

\\ Count derangements in world $x$
derangements(x) = sum(i=1,#x,x[i]!=i)
\\ Count populations in world $x$
populations(x) = my(n=#x,p=vector(n));for(i=1,n,p[x[i]]++);p
\\ List species in world $x$
species(x) = my(p=vecsort(x,,9));vector(#p,i,x[p[i]])
\\ Measure evident growth in world $x$
evidentgrowth(x) = my(p=populations(x));sum(i=1,#p,if(p[i]>1,p[i]-1,0))
\\ Determine lower bound for number of rounds experienced by world $x$
minrounds(x) = my(d=derangements(x));if(d==0,0,max(d,evidentgrowth()))

\\ Run a minimal evolution game
evolution(N=10,R=1000) = {
    my(x,a,vxi,vxj,xn);
    \\ Old world state
    x=vector(N,i,i);
    \\ Cell ages
    a=vector(N);
    \\ Pre-determined evolution history
    vxi=vector(R,i,random(N)+1);
    vxj=vector(R,i,random(N)+1);
    \\ Run up to $R$ rounds
    for(t=1,R,
        \\ Examine the world
        my(p=populations(x),s=species(x),m=minrounds(x));
        \\ Print current state
        print("t = ",t,", x = ",x,", ages = ",a,", pops = ",p,", numspecies = ",#s,", minrounds = ",m);
        \\ Terminate when one species dominates
        if(vecmax(p)==N,break());
        \\ Compute new state
        xn=vector(N,i,if(i==vxi[t],x[vxj[t]],x[i]));
        \\ Update ages
        a=vector(N,i,if(i==vxi[t],0,a[i]+1));
        \\ Begin next round
        x=xn;
        );
}