-rw-r--r-- 3315 cryptattacktester-20230614/isd1_prob.cpp raw
#include <cassert>
#include "transition.h"
#include "queue_prob.h"
#include "collision_prob.h"
#include "isd1_prob.h"

using namespace std;

bigfloat isd1_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 PI = 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 WINDOW = attackparams.at(pos++);
  bigint FW = attackparams.at(pos++);

  bigint K = K_orig-FW;

  bigint left = (K+L-Z)/2;
  bigint right = K+L-Z-left;

  // simple upper bound on isd1 success probability for one iteration:
  // binomial(left,PI)*binomial(right,PI)*binomial(N-K-L,W-2*PI)/binomial(N,W)

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

  bigint listsize0 = binomial(left,PI);
  bigint listsize1 = binomial(right,PI);
  bigint listsize = listsize0+listsize1;
  bigint prdenom = binomial(K+L,2*PI);
  bigint WINDOW1 = min(WINDOW,bigint(listsize-1));
  bigint pool = (2*listsize-WINDOW1-1)*WINDOW1/2;

  for (bigint prec = 32;;prec *= 2) {
    bigfloat::set_default_precision(prec);
  
    bigfloat prnumer = bigfloat(listsize0*listsize1);
    bigfloat pr = prnumer/bigfloat(prdenom);

    // pr is the _conditional_ probability of the right error distribution
    // _given_ 2*PI errors in K+L positions
  
    bigfloat p = exp2(bigfloat(-L));
    bigfloat B = p*prnumer;
    // expect B = listsize0*listsize1/2^L collisions
    bigfloat C = collision_average(listsize0,listsize1,(bigint) WINDOW,p);
    // expect C collisions to survive window limit
    bigfloat q = C/bigfloat(pool);
    bigfloat full_queue_clears = bigfloat(pool/PERIOD);
    bigint leftovers = pool%PERIOD;
    bigfloat D = full_queue_clears*queue_average(PERIOD,q,QUEUE_SIZE)+queue_average(leftovers,q,QUEUE_SIZE);
    // expect D collisions to survive queue-length limit
    pr *= D/B;

    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,2*PI,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);

    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;
  }
}