attack.h

/*
attack.h
Written by Walter O. Krawec, Michael Nelson, and Eric Geiss
Copyright (c) 2017
 
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
 
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
 
 
#pragma once
#include <vector>
#include "Quantum/Algebra.h"
 
namespace QKD
{
struct Basis
{
  std::vector < std::vector < Complex > > B;
  void orthogonalize();
 
  Complex innerProduct(std::vector < Complex >& a, std::vector < Complex >& b)
  {
    Complex sum = 0;
 
    for(int i=0; i<a.size(); ++i)
      sum = sum + a[i] * conj(b[i]);
 
    return sum;
  }
 
  double norm(std::vector <Complex> & a)
  {
    double sum = 0;
    Complex b;
    for(int i=0; i<a.size(); ++i){
      b = a[i] * conj(a[i]);
      sum += creal(b);
    }
 
    return sqrt(sum);
  }
};
 
struct Stats
{
  double p00, p01, p10, p11;
  double Qx, Qy; // Qx = 1-paa; Qy = 1-pbb
  double p0a, p1a, pa0;
  double p0b, p1b, pb0;
 
  double BB84Rate();
  double SafeLog(double x)
  {
    if(x < 0.0000001)
      return 0;
    return log(x)/log(2.0);
  }
  double entropy(double x)
  {
    if(x > 1.00001){
      std::cout << "Entropy error: x = " << x << "\n";
      return 0;
      //throw std::string("Entropy Error.");
    }
 
    return -x*SafeLog(x) - (1-x)*SafeLog(1-x);
  }
  double HAB(double key00, double key01, double key10, double key11)
  {
    double first = -key00*SafeLog(key00) - key01*SafeLog(key01) - key10*SafeLog(key10) - key11*SafeLog(key11);
 
    if(key00+key10 > 1.00001)
      std::cout << "HAB ERROR!\n";
    double second = entropy(key00+key10);
 
    return first-second;
  }
 
 
  void cnot();
  void symmetric(double Q);
  void testChannel1()
  {
    p00 = 0.752;
    p01 = 1-p00;
    p10 = .055;
    p11 = 1-p10;
 
    Qx = .163;
    Qy = .115;
 
    p0a = .435;
    p1a = .592;
    pa0 = .507;
 
    p0b = .660;
    p1b = .357;
    pb0 = .322;
  }
  void testChannel2()
  {
    p00 = 0.790;
    p01 = 1-p00;
    p10 = .099;
    p11 = 1-p10;
 
    Qx = .121;
    Qy = .283;
 
    p0a = .581;
    p1a = .381;
    pa0 = .373;
 
    p0b = .337;
    p1b = .453;
    pb0 = .499;
  }
 
  void printChannel(std::ostream& f)
  {
    f << "CHANNEL STATS:\n";
    f << "p00 = " << p00 << "\n";
    f << "p01 = " << p01 << "\n";
    f << "p10 = " << p10 << "\n";
    f << "p11 = " << p11 << "\n\n";
 
    f << "Qx = " << Qx << "\n";
    f << "Qy = " << Qy << "\n\n";
 
    f << "p0+ = " << p0a << "\n";
    f << "p1+ = " << p1a << "\n";
    f << "p+0 = " << pa0 << "\n\n";
 
    f << "p0Y = " << p0b << "\n";
    f << "p1Y = " << p1b << "\n";
    f << "pY0 = " << pb0 << "\n\n";
  }
 
  double randFloat(double small, double large)
  {
    int r = rand()%RAND_MAX;
    double r2 = r / (double)RAND_MAX;
    
    return (large-small)*r2 + small;
  }
 
  void generateRandom(double p01, double p10);
  void chooseRandomVector(double minScale, double maxScale, std::vector <Complex>& v, int startIndex, int endIndex);
};
 
class Attack
{
  static const double precision = .01;
 public:
  Attack(Stats& _stats) : stats(_stats)
  {
    setLocalVars();
 
    v0.create(4,1);
    v1.create(4,1);
    v2.create(4,1);
    v3.create(4,1);
 
    v0.zero();
    v1.zero();
    v2.zero();
    v3.zero();
 
    v0(0,0) = 1;
    v1(1,0) = 1;
    v2(2,0) = 1;
    v3(3,0) = 1;
  }
 
  // sets E03 and E12 to their start values
  void begin()
  {
    double im03lb = -SafeSqrt(E00*E33 - ReE03*ReE03);
    double im12lb = -SafeSqrt(E11*E22 - ReE12*ReE12);
 
    E03 = creal(E03) + I*im03lb;
    E12 = creal(E12) + I*im12lb;
 
    debugStart = true;
  }
 
  // gets the next possible attack; returns TRUE when this is the last legal attack; FALSE o/w
  bool getNext(algebra::mat & e0, algebra::mat& e1, algebra::mat& e2, algebra::mat& e3)
  {
    bool done = false;
    do{
      //E12 = quantum::Complex(-.0361068, .00601651);
      //E03 = quantum::Complex(.73823, -.230657);
 
      done = buildAttackVec(e0, e1, e2, e3);
 
      // add to E12.imag:
      E12 = creal(E12) + I*(cimag(E12)+precision);
      if(creal(E12*conj(E12)) > E11*E22){
    double im12lb = -SafeSqrt(E11*E22 - ReE12*ReE12);
    E12 = creal(E12) + I*im12lb;
 
    E03 = creal(E03) + I*(cimag(E03) + precision);  
    if(creal(E03*conj(E03)) > E00*E33){
      if(debugStart)
        throw(std::string("NO ATTACKS"));
      return true; // no more attacks to consider
    }
      }
    }while(!done);
 
    debugStart= false;
    return false;
  }
 
    bool debugStart;
 
  double SafeSqrt(double x)
  {
    if(x < -0.00001){
      std::cout << "\n\nSafeSQRT Error: x = " << x << "\n\n";
      throw(std::string("SafeSQRT Error; X negative."));
    }
    if(x < .00000001)
      return 0;
 
    return sqrt(x);
  }
 
  bool buildAttackVec(algebra::mat & e0, algebra::mat& e1, algebra::mat& e2, algebra::mat& e3)
  {
    // e0:
    Complex a0 = SafeSqrt(stats.p00);
 
    // e3:
    Complex a3 = E03 / creal(a0);
 
    if(stats.p11 - creal(a3)*creal(a3) - cimag(a3)*cimag(a3) < 0)
      return false;
    Complex b3 = SafeSqrt(stats.p11 - creal(a3)*creal(a3) - cimag(a3)*cimag(a3));
 
 
    // e1:
    Complex a1 = E01 / creal(a0);
    Complex b1 = E13 - conj(a1)*a3;
    b1 = b1 * (1.0/creal(b3));
    b1 = conj(b1);
    if(creal(b3) < .0000001){
      b1 = 0;
    }
 
    if(stats.p01 - creal(a1*conj(a1) + b1*conj(b1)) < 0)
      return false;
 
    Complex c1 = SafeSqrt(stats.p01 - creal(a1*conj(a1) + b1*conj(b1)));
 
 
    // e2:
    Complex a2 = E02/creal(a0);
    Complex b2 = E23 - conj(a2)*a3;
    b2 = b2 * (1.0/creal(b3));
    b2 = conj(b2);
 
    if(creal(b3) < .0000001){
      // check to make sure this is a valid solution:
      Complex check = conj(a1)*a3;
      if(fabs(creal(check) - creal(E13)) > .0000001 ||
     fabs(cimag(check) - cimag(E13)) > .0000001){
    //std::cout << "INVALID\n";
    return false;
      }
      b2 = 0;
    }
 
    Complex c2 = E12 - conj(a1)*a2 - conj(b1)*b2;
 
 
    c2 = c2 * (1.0/creal(c1));
 
    if(creal(c1) < .0000001){
      // check:
      Complex check = conj(a1)*a2 + conj(b1)*b2;
      if(fabs(creal(check) - creal(E12)) > .0000001 ||
     fabs(cimag(check) - cimag(E12)) > .0000001){
    //std::cout << "INVALID12\n";
    return false;
      }
      c2 = 0;
    }
 
    
 
    double test = stats.p10 - creal( a2*conj(a2) + b2*conj(b2) + c2*conj(c2) );
    if(test < 0)
      return false;
    Complex d2 = SafeSqrt(test);
 
 
    /*
    // a test:
    std::cout << "\n--test--\n";
    // e0e2:
    quantum::Complex t = a0.conj()*a2;
    std::cout << "E02: ";
    t.print();
    //e1e3:
    t = a1.conj()*a3 + b1.conj()*b3;
    std::cout << "\nE13: ";
    t.print();
    std::cout << "--end test--\n\n";
    */
 
    e0.zero();
    e1.zero();
    e2.zero();
    e3.zero();
 
    e0.add(a0, v0);
 
    e3.add(a3, v0);
    e3.add(b3, v1);
 
    e1.add(a1, v0);
    e1.add(b1, v1);
    e1.add(c1, v2);
 
    e2.add(a2, v0);
    e2.add(b2, v1);
    e2.add(c2, v2);
    e2.add(d2, v3);
 
    return true;
  }
 
 private:
  void setLocalVars()
  {
    E00 = stats.p00;
    E11 = stats.p01;
    E22 = stats.p10;
    E33 = stats.p11;
 
    E01 = stats.p0a - .5 + I*(stats.p0b-.5);
    E23 = stats.p1a - .5 + I*(stats.p1b-.5);
    E02 = stats.pa0 - .5*(E00+E22) + I*(.5*(E00+E22)-stats.pb0);
    E13 = E02*(-1);
 
    ReE03 = 1 - stats.Qx - stats.Qy - .5*(creal(E01) + cimag(E01) + creal(E23) + cimag(E23));
    ReE12 = stats.Qy - stats.Qx + .5*(cimag(E01) - creal(E01) + cimag(E23) - creal(E23));
 
    E03 = ReE03;
    E12 = ReE12;
  }
 
  Stats stats;
  algebra::mat v0, v1, v2, v3; // standard basis vectors
 
  double E00, E11, E22, E33;
  Complex E01, E23, E02, E13;
 
  double ReE03, ReE12;
  Complex E03, E12;
};
}