-rw-r--r-- 7469 cryptattacktester-20231020/isd2_prob.cpp raw
#include <map>
#include <cassert>
#include "transition.h"
#include "queue_prob.h"
#include "collision_prob.h"
#include "isd2_prob.h"
using namespace std;
bigfloat isd2_prob(const vector<bigint> ¶ms,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 PIJ = attackparams.at(pos++);
bigint PI = attackparams.at(pos++);
bigint L0 = attackparams.at(pos++);
bigint L1 = attackparams.at(pos++);
bigint CHECKPI = attackparams.at(pos++);
bigint CHECKSUM = attackparams.at(pos++);
bigint D = attackparams.at(pos++);
bigint Z = attackparams.at(pos++);
bigint QU0 = attackparams.at(pos++);
bigint QF0 = attackparams.at(pos++); auto PE0 = QF0*QU0;
bigint WI0 = attackparams.at(pos++);
bigint QU1 = attackparams.at(pos++);
bigint QF1 = attackparams.at(pos++); auto PE1 = QF1*QU1;
bigint WI1 = attackparams.at(pos++);
bigint FW = attackparams.at(pos++);
assert(CHECKPI || CHECKSUM);
bigint K = K_orig-FW;
bigint L = L0+L1;
bigint left = (K+L-Z)/2;
bigint right = K+L-Z-left;
bigint listsize0 = binomial(left,PIJ);
bigint listsize1 = binomial(right,PIJ);
bigint listsize = listsize0+listsize1;
WI0 = min(WI0,bigint(listsize-1));
bigint pool = (2*listsize-WI0-1)*WI0/2;
bigint queue_clears = (pool+PE0-1)/PE0;
bigint leftovers = pool%PE0;
for (bigint prec = 32;;prec *= 2) {
bigfloat::set_default_precision(prec);
map<pair<bigint,bigint>,bigfloat> prsplit;
map<pair<bigint,bigint>,bigint> representations;
for (bigint E = 0;E <= 2*PIJ;E += 2) {
if (CHECKPI) if (E != PI) continue;
for (bigint F = 0;F <= 2*PIJ;F += 2) {
if (CHECKPI) if (F != PI) continue;
pair<bigint,bigint> EF = make_pair(E,F);
// ----- basic error split
// chance of (E,F,0,W-E-F) errors in (left,right,Z,N-K-L):
// binomial(left,E)*binomial(right,F)*binomial(N-K-L,W-E-F)/binomial(N,W)
// compare to chance of having (E+F,W-E-F) errors in (K+L,N-K-L):
// binomial(K+L,E+F)*binomial(N-K-L,W-E-F)/binomial(N,W)
bigint prnumer = binomial(left,E)*binomial(right,F);
bigint prdenom = binomial(K+L,E+F);
prsplit[EF] = bigfloat(prnumer)/bigfloat(prdenom);
// _conditional_ probability of (E,F,0,W-E-F) _given_ (E+F,W-E-F)
// ----- representations
// isd2 searches for one of many representations of the E+F errors
// representation: xor(first PIJ,second PIJ) on left; xor(first PIJ,second PIJ) on right
// on left: binomial(E,E/2) ways to split E into (first part,second part)
// and also need PIJ-E/2 further bits shared between first PIJ,second PIJ
// overall binomial(left-E,PIJ-E/2)*binomial(E,E/2) left decompositions
representations[EF] = binomial(left-E,PIJ-E/2)*binomial(right-F,PIJ-F/2)*binomial(E,E/2)*binomial(F,F/2);
}
}
// ----- restricted differences
// isd2 restricts to D choices of differences Delta, with 1<=D<=2^L0
// and finds a representation only if it is suitable
// meaning first PIJ on left and first PIJ on right match Delta in L0 bits
// and (equivalently) second PIJ on left and second PIJ on right match Delta+syndrome in L0 bits
bigfloat halfL0 = exp2(bigfloat(-L0));
bigfloat reprfound = bigfloat(D)*halfL0;
// ----- window limit and queue limit for first merge and second merge
bigfloat B = bigfloat(listsize0*listsize1)*halfL0;
// expect B = listsize0*listsize1/2^L0 collisions
bigfloat C = collision_average(listsize0,listsize1,(bigint) WI0,halfL0);
// expect C collisions to survive window limit
bigfloat q = C/bigfloat(pool);
bigfloat full_queue_clears = bigfloat(pool/PE0);
bigfloat D = full_queue_clears*queue_average(PE0,q,QU0)+queue_average(leftovers,q,QU0);
// expect D collisions to survive queue-length limit
if (!D.definitely_positive())
continue;
bigfloat DB = D/B; // chance that a collision survives window limit and queue-length limit
reprfound *= DB*DB; // have to survive two merges
if (D.definitely_less(1)) {
; // model 1 collision as definitely being found
} else {
// ----- valid portion of the root list
// nominal root list size is 2*queue_clears*QU0
// but expect only first 2*D entries to be valid
bigint effectiveWI1 = WI1;
while (effectiveWI1 > 0 && (2*D).definitely_less(bigfloat(effectiveWI1)))
--effectiveWI1;
bigfloat validrootpool = (4*D-bigfloat(effectiveWI1)-1)*bigfloat(effectiveWI1)/2;
bigfloat validrootfullqueueclears = bigfloat(floor_as_bigint(validrootpool/bigfloat(PE1)));
bigfloat validrootleftovers = validrootpool-validrootfullqueueclears*bigfloat(PE1);
// ----- window limit and queue limit for root merge
bigfloat halfL1 = exp2(bigfloat(-L1));
bigfloat rootB = D*D*halfL1; // expect this many root collisions
bigfloat rootC = collision_average(D,D,(bigint) effectiveWI1,halfL1);
bigfloat rootq = rootC/validrootpool;
bigfloat rootD = validrootfullqueueclears*queue_average(PE1,rootq,QU1)+queue_average(validrootleftovers,rootq,QU1);
if (!rootB.definitely_positive()) continue;
if (!rootD.definitely_positive()) continue;
reprfound *= rootD/rootB;
}
// ----- putting everything together
vector<bool> prsuccessnonzero(4*PIJ+1,0);
vector<bigfloat> prsuccess(4*PIJ+1,0);
for (bigint E = 0;E <= 2*PIJ;E += 2) {
if (CHECKPI) if (E != PI) continue;
for (bigint F = 0;F <= 2*PIJ;F += 2) {
if (CHECKPI) if (F != PI) continue;
pair<bigint,bigint> EF = make_pair(E,F);
bigfloat pr = -expm1(bigfloat(representations[EF])*log1p(-reprfound));
pr *= prsplit[EF];
// pr is now the conditional probability
// of E,F,0 errors in left,right,Z with suitable representation
// surviving window limit and queue limit for first merge, second merge, root merge
// given E+F errors in K+L
prsuccess.at(E+F) += pr;
prsuccessnonzero.at(E+F) = 1;
}
}
// ----- markov chain
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_identity(W);
for (bigint EplusF = 0;EplusF <= 4*PIJ;EplusF += 2)
if (prsuccessnonzero.at(EplusF))
T2 = transition_multiply(W,T2,transition_found(W,EplusF,prsuccess.at(EplusF)));
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;
}
}