-rw-r--r-- 4362 cryptattacktester-20230614/transition.cpp raw
#include <map>
#include <cassert>
#include "transition.h"

using namespace std;

// consider the following markov chain:
// state 0: search has not finished, and target has 0 errors within KK positions
// state 1: search has not finished, and target has 1 error within KK positions
// ...
// state W: search has not finished, and target has W errors within KK positions
// state W+1: have finished by finding target
// state W+2: have finished with systematic-form failure

// initial state is 0...W
// chance of state e: binomial(KK,e)*binomial(N-KK,W-e)/binomial(N,W)

vector<vector<bigfloat>> transition_zero(bigint W)
{
  vector<vector<bigfloat>> result;
  vector<bigfloat> zerorow;

  for (bigint v = 0;v <= W+2;++v) zerorow.push_back(0);
  for (bigint u = 0;u <= W+2;++u) result.push_back(zerorow);
  return result;
}

vector<vector<bigfloat>> transition_identity(bigint W)
{
  vector<vector<bigfloat>> result = transition_zero(W);
  for (bigint u = 0;u <= W+2;++u) result.at(u).at(u) = 1;
  return result;
}

static map<tuple<bigint,bigint,bigint,bigint,bigint>,vector<vector<bigfloat>>> transition_weightstep_cache;

// output[u][v] is chance that a type-3 exchange of X columns moves state u to state v
vector<vector<bigfloat>> transition_weightstep(bigint N,bigint KK,bigint W,bigint X)
{
  tuple<bigint,bigint,bigint,bigint,bigint> key = make_tuple(N,KK,W,X,bigfloat::get_default_precision());

  if (transition_weightstep_cache.count(key))
    return transition_weightstep_cache[key];

  vector<vector<bigfloat>> result = transition_identity(W);

  bigfloat denom = bigfloat(binomial(N-KK,X)*binomial(KK,X));

  for (bigint u = 0;u <= W;++u) {
    for (bigint v = 0;v <= W;++v) {
      bigint resultuv = 0;
      bigint d = u-v;
      if (d > X && d < -X) continue;
      for (bigint i = 0;i <= W-u;++i) {
        if (X-i < 0) continue;
        if (d+i < 0) continue;
        if (X-d-i < 0) continue;
        if (X-i > N-KK-W+u) continue;
        if (d+i > u) continue;
        if (X-d-i > KK-u) continue;
        resultuv += binomial(W-u,i)*binomial(N-KK-W+u,X-i)*binomial(u,d+i)*binomial(KK-u,X-d-i);
      }
      result.at(u).at(v) = bigfloat(resultuv)/denom;
    }
  }

  transition_weightstep_cache[key] = result;

  return result;
}

vector<vector<bigfloat>> transition_checksystematic(bigint W,bigint X,bigint Y)
{
  vector<vector<bigfloat>> result = transition_identity(W);
  bigfloat success;

  assert(X <= Y);
  // what's the chance that an XxY matrix has rank X?
  // 1/2^Y chance of first row failing
  // if it's ok: 1/2^(Y-1) chance of second row failing
  // ...
  // so overall chance (1-1/2^Y)(1-1/2^(Y-1))...(1-1/2^(Y-X+1))

  bigfloat log05 = -log1p(bigfloat(1));
  success = 1;
  for (bigint pos = Y;pos > Y-X;--pos)
    success *= -expm1(pos*log05);

  for (bigint u = 0;u <= W;++u) {
    result.at(u).at(u) = success;
    result.at(u).at(W+2) = 1-success;
  }
  return result;
}

// pr: conditional probability of finding target given P errors
vector<vector<bigfloat>> transition_found(bigint W,bigint P,bigfloat pr)
{
  vector<vector<bigfloat>> result = transition_identity(W);
  assert(P <= W);
  result.at(P).at(P) = 1-pr;
  result.at(P).at(W+1) = pr;
  return result;
}

vector<vector<bigfloat>> transition_multiply(bigint W,const vector<vector<bigfloat>> &T0,const vector<vector<bigfloat>> &T1)
{
  vector<vector<bigfloat>> result = transition_zero(W);
  for (bigint t = 0;t <= W+2;++t)
    for (bigint u = 0;u <= W+2;++u) {
      bigfloat T0tu = T0.at(t).at(u);
      if (T0tu.definitely_zero()) continue;
      for (bigint v = 0;v <= W+2;++v)
        result.at(t).at(v) += T0tu*T1.at(u).at(v);
    }

#ifdef TRANSITION_RENORMALIZE
  // renormalize for better floating-point stability
  for (bigint t = 0;t <= W+2;++t) {
    bigfloat x = 0;
    for (bigint u = 0;u <= W+2;++u)
      x += result.at(t).at(u);
    if (x <= 0) continue;
    x = 1/x;
    for (bigint u = 0;u <= W+2;++u)
      result.at(t).at(u) *= x;
  }
#endif
  return result;
}

vector<vector<bigfloat>> transition_power(bigint W,const vector<vector<bigfloat>> &T,bigint e)
{
  assert(e >= 0);
  if (e == 0) return transition_identity(W);
  if (e == 1) return T;
  vector<vector<bigfloat>> result = transition_power(W,T,e/2);
  result = transition_multiply(W,result,result);
  if ((e&1) != 0) result = transition_multiply(W,result,T);
  return result;
}