Randil's math problem revisited (C++ code)

//this version updates by bayesian probability instead of using maximum likelihood
#include <iostream>
#include <iomanip>
#include <ctime>
#include <random>
#include <cmath>

#define TRIES 5000

using namespace std;

int main()
{
    mt19937 mersenne(static_cast<unsigned int>(time(0)));

    for(int lo=1; lo<=18; lo++)
        for(int hi=lo+1;hi<=19;hi++){

            double pr=0.05*hi, pw=0.05*lo;
            double qr=1-pr, qw=1-pw;

            double bgreen=pr/pw; //should be greater than 1
            double bred=qr/qw; //should be less than 1

            int countsum=0;
            for(int temp=1;temp<=TRIES;temp++){

                unsigned int thenum=mersenne()%10000;
                unsigned int num[]={thenum/1000, (thenum/100)%10, (thenum/10)%10, thenum%10};

                int flags[10000]{};
                double digitprobs[40]{};

                double bayesprobs[10000]{};
                int countsgreen[40]{};
                int countsred[40]{};

                //first guess

                unsigned int theguess=mersenne()%10000;
                unsigned int guess[]={theguess/1000, (theguess/100)%10, (theguess/10)%10, theguess%10};


                int count=1;
                while(true){
                    //check if guess matches num
                    int check=1;
                    for(int i=0;i<4;i++)
                        if(guess[i]!=num[i]) {check=0; break;}
                    if(check==1) //match
                        break;

                    //no match

                    //add occurrences
                    for(int i=0;i<4;i++){
                        if(guess[i]==num[i]){
                            //green with prob pr, red with prob qr
                            if(((double) mersenne())/mersenne.max() < pr) //green
                                countsgreen[10*i+guess[i]]++;
                            else
                                countsred[10*i+guess[i]]++;
                        }
                        else{
                            //green with prob pw, red with prob qw
                            if(((double) mersenne())/mersenne.max() < pw) //green
                                countsgreen[10*i+guess[i]]++;
                            else
                                countsred[10*i+guess[i]]++;
                        }
                    }

                    //set flag so as to not guess it again
                    flags[theguess]=1;

                    //calculate normalized likelihood for each digit
                    for(int i=0;i<40;i++)
                        digitprobs[i]=pow(bgreen,countsgreen[i])*pow(bred,countsred[i]);


                    //calculate normalized likelihood for each 4-digit number
                    //bayesprobs will store the cumulative sum (likelihood function is proportional to bayesian probability function)
                    double cumulative=0;
                    for(int i1=0;i1<10;i1++)
                        for(int i2=0;i2<10;i2++)
                            for(int i3=0;i3<10;i3++)
                                for(int i4=0;i4<10;i4++){
                                    double s=0;
                                    int fourdigitval=1000*i1+100*i2+10*i3+i4;
                                    if(!flags[fourdigitval]){
                                        s=digitprobs[i1]*digitprobs[10+i2]*digitprobs[20+i3]*digitprobs[30+i4];
                                    }
                                    cumulative+=s;
                                    bayesprobs[fourdigitval]=cumulative;
                                }

                    //now take a new random number (uniform distribution from 0 to bayesprobs[9999]) and use bayesprobs to determine which number to guess
                    //use halving strategy to find the first number larger than rndnum or equal to bayesprobs[9999]
                    double rndnum=(((double) mersenne())/mersenne.max())*bayesprobs[9999];
                    int min=-1, max=9999;
                    while(max>min+1){


                        int mid=floor((max+min)/2);
                        if(bayesprobs[mid]>rndnum || bayesprobs[mid]==bayesprobs[9999])
                            max=mid;
                        else
                            min=mid;
                    }

                    //set guess to value of max (first number larger than rndnum or equal to bayesprobs[9999])

                    theguess=max;
                    guess[0]=max/1000;
                    guess[1]=(max/100)%10;
                    guess[2]=(max/10)%10;
                    guess[3]=max%10;

                    count++;

                    if(count>=10002){
                        cout << "failed";
                        exit(1);
                    }


                }

                countsum+=count;
            }

            cout << fixed << setprecision(2) << "Average tries for pr=" << pr << " , pw=" << pw << ": " << ((double) countsum)/TRIES << endl;
        }
    return 0;
}