-rw-r--r-- 2816 cryptattacktester-20230614/isd0_prob.cpp raw
#include "transition.h"
#include "queue_prob.h"
#include "isd0_prob.h"

using namespace std;

bigfloat isd0_prob(const vector<bigint> &params,const vector<bigint> &attackparams)
{
  bigint N = params.at(0);
  bigint K_orig = params.at(1);
  bigint W = params.at(2);
  bigint pos = 0;
  bigint ITERS = attackparams.at(pos++);
  bigint RESET = attackparams.at(pos++);
  bigint X = attackparams.at(pos++);
  bigint YX = attackparams.at(pos++); auto Y = X+YX;
  bigint P = attackparams.at(pos++);
  bigint L = attackparams.at(pos++);
  bigint Z = attackparams.at(pos++);
  bigint QUEUE_SIZE = attackparams.at(pos++);
  bigint QF = attackparams.at(pos++); auto PERIOD = QF*QUEUE_SIZE;
  bigint FW = attackparams.at(pos++);

  bigint K = K_orig-FW;

  // simple upper bound on isd0 success probability for one iteration:
  // binomial(K+L-Z,P)*binomial(N-K-L,W-P)/binomial(N,W)

  // extra complications below are for three reasons:
  // * account for non-randomness from X
  // * account for losses from Y
  // * account for losses from QUEUE_SIZE vs. PERIOD

  // first compute pr: the _conditional_ probability of isd0 working
  // _given_ P errors in K+L positions
  // isd0 works exactly when the P errors are in K+L-Z positions
  bigint prnumer = binomial(K+L-Z,P);
  bigint prdenom = binomial(K+L,P);

  for (bigint prec = 32;;prec *= 2) {
    bigfloat::set_default_precision(prec);
  
    bigfloat pr = bigfloat(prnumer)/bigfloat(prdenom);
  
    if (P > 0 && L > 0 && PERIOD > QUEUE_SIZE) {
      bigfloat p = exp2(bigfloat(-L));
      bigint pool = binomial(K+L-Z,P);
      bigfloat full_queue_clears = bigfloat(pool/PERIOD);
      bigint leftovers = pool%PERIOD;
      bigfloat average = full_queue_clears*queue_average(PERIOD,p,QUEUE_SIZE)+queue_average(leftovers,p,QUEUE_SIZE);
      pr *= average/(p*bigfloat(pool));
    }
    vector<vector<bigfloat>> T0 = transition_weightstep(N,K+L,W,X);
    vector<vector<bigfloat>> T1 = transition_checksystematic(W,X,Y);
    vector<vector<bigfloat>> T2 = transition_found(W,P,pr);
    vector<vector<bigfloat>> Titer = transition_multiply(W,T0,transition_multiply(W,T1,T2));
    vector<vector<bigfloat>> Tmanyiters = transition_power(W,Titer,RESET-1);
    Tmanyiters = transition_multiply(W,T2,Tmanyiters); // first iteration skips weightstep and checksystematic

    bigfloat result = 0;
    for (bigint u = 0;u <= W;++u) {
      bigint prstartnumer = binomial(K+L,u)*binomial(N-K-L,W-u);
      bigint prstartdenom = binomial(N,W);
      bigfloat prstart = bigfloat(prstartnumer)/bigfloat(prstartdenom);
      result += prstart*Tmanyiters.at(u).at(W+1);
    }

    result = -expm1(bigfloat(ITERS/RESET)*log1p(-result));

    if (FW)
      result *= bigfloat(1)-exp2(bigfloat(-K_orig));

    if (result.definitely_positive_and_reasonably_precise()) return result;
  }
}