matrix2x2.cpp
// Copyright (c) 2012 Vadym Kliuchnikov sqct(dot)software(at)gmail(dot)com, Dmitri Maslov, Michele Mosca
//
// This file is part of SQCT.
//
// SQCT is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// SQCT is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with SQCT. If not, see <http://www.gnu.org/licenses/>.
//
#include "matrix2x2.h"
#include "resring.h"
#include "hprhelpers.h"
#include <iostream>
#include <limits>
#include <algorithm>
using namespace std;
typedef ring_int<int>::mpclass mpclass;
template < class TInt>
matrix2x2<TInt> matrix2x2<TInt>::operator*( const matrix2x2<TInt>& b ) const
{
return matrix2x2<TInt>(
b.d[0][0]* d[0][0]+b.d[1][0]* d[0][1], b.d[0][1]* d[0][0]+b.d[1][1]* d[0][1],
b.d[0][0]* d[1][0]+b.d[1][0]* d[1][1], b.d[0][1]* d[1][0]+b.d[1][1]* d[1][1],
b.de + de);
}
template < class TInt>
matrix2x2<TInt> matrix2x2<TInt>::operator*( const ring_int<TInt>& b ) const
{
return matrix2x2<TInt>( b * d[0][0], b*d[0][1], b*d[1][0], b*d[1][1], de);
}
template < class TInt>
void matrix2x2<TInt>::mul_eq_w( int k )
{
d[0][0].mul_eq_w(k);
d[0][1].mul_eq_w(k);
d[1][0].mul_eq_w(k);
d[1][1].mul_eq_w(k);
}
template < class TInt>
matrix2x2<TInt>& matrix2x2<TInt>::div_eq_sqrt2_exp(int a)
{
d[0][0].div_eq_sqrt2(a) ; d[0][1].div_eq_sqrt2(a);
d[1][0].div_eq_sqrt2(a) ; d[1][1].div_eq_sqrt2(a);
return *this;
}
template < class TInt>
void matrix2x2<TInt>::set( const ring_int<TInt>& z, const ring_int<TInt>& u, const ring_int<TInt>& w, const ring_int<TInt>& v, int denom_exp )
{
d[0][0] = z; d[0][1] = u;
d[1][0] = w; d[1][1] = v;
de = denom_exp;
}
template < class TInt>
matrix2x2<TInt>::matrix2x2( const ring_int<TInt>& z, const ring_int<TInt>& u, const ring_int<TInt>& w, const ring_int<TInt>& v, int denom_exp )
{
set(z,u,
w,v,
denom_exp);
}
template < class TInt>
matrix2x2<TInt>::matrix2x2()
{
set( ring_int<TInt>(),ring_int<TInt>(),ring_int<TInt>(),ring_int<TInt>(),0 );
}
template < class TInt>
int matrix2x2<TInt>::min_gde() const
{
int min_gde = std::numeric_limits<int>::max();
for( int i = 0; i < 2; ++i )
for( int j = 0; j < 2; ++j )
min_gde = std::min( d[i][j].gde(), min_gde );
return min_gde;
}
template < class TInt>
int matrix2x2<TInt>::min_gde_abs2() const
{
int min_gde = std::numeric_limits<int>::max();
for( int i = 0; i < 2; ++i )
for( int j = 0; j < 2; ++j )
min_gde = std::min( d[i][j].abs2().gde(), min_gde );
return min_gde;
}
template < class TInt>
int matrix2x2<TInt>::reduce()
{
matrix2x2<TInt>& A = *this;
int gde = min_gde();
if( gde != std::numeric_limits<int>::max() )
{
int red_power = std::min ( gde, de );
A.div_eq_sqrt2_exp(red_power);
de -= red_power;
}
return de;
}
template < class TInt>
int matrix2x2<TInt>::reduce(int red_power)
{
div_eq_sqrt2_exp(red_power);
de -= red_power;
return de;
}
template < class TInt >
int matrix2x2<TInt>::max_sde_abs2() const
{
return sde( 2 * de, min_gde_abs2() );
}
template < class TInt >
void matrix2x2<TInt>::mul_TkH(int k, matrix2x2<TInt> &out)
{
ring_int<TInt> d01wk;
ring_int<TInt> d11wk;
d[0][1].mul_w(k,d01wk);
d[1][1].mul_w(k,d11wk);
out.d[0][0] = d[0][0];
out.d[0][0] += d01wk;
out.d[1][0] = d[1][0];
out.d[1][0] += d11wk;
out.d[0][1] = d[0][0];
out.d[0][1] -= d01wk;
out.d[1][1] = d[1][0];
out.d[1][1] -= d11wk;
out.de = de + 1;
}
template < class TInt >
matrix2x2<TInt> matrix2x2<TInt>::T( int power )
{
return matrix2x2<TInt>( ring_int<TInt>(1), ring_int<TInt>(0),
ring_int<TInt>(0), ring_int<TInt>::omega(power) );
}
template < class TInt>
matrix2x2<TInt> matrix2x2<TInt>::H()
{
return matrix2x2<TInt>( ring_int<TInt>(1), ring_int<TInt>( 1),
ring_int<TInt>(1), ring_int<TInt>(-1),
1 );
}
template < class TInt>
matrix2x2<TInt> matrix2x2<TInt>::X()
{
return matrix2x2<TInt>( ring_int<TInt>(0), ring_int<TInt>( 1),
ring_int<TInt>(1), ring_int<TInt>(0),
0 );
}
template < class TInt>
matrix2x2<TInt> matrix2x2<TInt>::Z()
{
return matrix2x2<TInt>( ring_int<TInt>(1), ring_int<TInt>( 0),
ring_int<TInt>(0), ring_int<TInt>(-1),
0 );
}
template < class TInt>
matrix2x2<TInt> matrix2x2<TInt>::Y()
{
return matrix2x2<TInt>( ring_int<TInt>(0), -ring_int<TInt>::omega(2),
ring_int<TInt>::omega(2), ring_int<TInt>(0),
0 );
}
template < class TInt>
matrix2x2<TInt> matrix2x2<TInt>::Id( int power )
{
return matrix2x2<TInt>( ring_int<TInt>::omega(power), ring_int<TInt>(0),
ring_int<TInt>(0), ring_int<TInt>::omega(power),
0 );
}
template < class TInt>
matrix2x2<TInt> matrix2x2<TInt>::P()
{
return matrix2x2<TInt>( ring_int<TInt>(1), ring_int<TInt>(0),
ring_int<TInt>(0), ring_int<TInt>::omega(2),
0 );
}
template < class TInt>
matrix2x2<TInt> matrix2x2<TInt>::conjugateTranspose()
{
return matrix2x2<TInt>( d[0][0].conjugate(), d[1][0].conjugate(),
d[0][1].conjugate(), d[1][1].conjugate(), de );
}
template < class TInt>
bool matrix2x2<TInt>::is_unitary() const
{
auto r1 = d[0][0].abs2() + d[0][1].abs2();
auto r2 = d[1][0].abs2() + d[1][1].abs2();
auto r3 = d[0][0].abs2() + d[1][0].abs2();
auto r4 = d[0][1].abs2() + d[1][1].abs2();
decltype(r1) sum( 1 );
sum <<= de;
return (r1 == sum) && (r2 == sum) && (r3 == sum) && (r4 == sum);
}
template < class TInt>
matrix2x2<TInt> &matrix2x2<TInt>::operator =(const vector2<TInt> &v)
{
d[0][0] = v[0]; d[0][1] = - v[1].conjugate();
d[1][0] = v[1]; d[1][1] = v[0].conjugate();
de = v.de;
}
template < class TInt>
bool matrix2x2<TInt>::operator <( const matrix2x2<TInt>& B ) const
{
if ( de > B.de )
return false;
else if( de < B.de )
return true;
for( int i = 0; i < 2; ++i )
for( int j = 0; j < 2; ++j )
{
int le = d[i][j].le(B.d[i][j]); // +1 -- less, 0 -- equal, -1 greater
if( le < 0 )
return false;
else if( le > 0 )
return true;
}
return false;
}
template < class TInt >
template < class Tcl >
matrix2x2<TInt>::matrix2x2( const matrix2x2<Tcl>& val )
{
typedef ring_int<TInt> r;
set( (r) val.d[0][0], (r) val.d[0][1],(r) val.d[1][0],(r) val.d[1][1], val.de );
}
template < class TInt>
bool matrix2x2<TInt>::operator ==( const matrix2x2<TInt>& B ) const
{
for( int i = 0; i < 2; ++i )
for( int j = 0; j < 2; ++j )
{
if ( d[i][j].le(B.d[i][j]) != 0 ) // +1 -- less, 0 -- equal, -1 greater
return false;
}
return true;
}
template < class TInt >
matrix2x2hpr::matrix2x2hpr(const matrix2x2<TInt> &val)
{
for( int i = 0; i < 2; ++i )
for( int j = 0; j < 2; ++j )
{
d[i][j] = val.d[i][j].toHprComplex( val.de );
}
}
matrix2x2hpr::matrix2x2hpr( scalar a, scalar b, scalar c, scalar _d)
{
d[0][0] = a;
d[0][1] = b;
d[1][0] = c;
d[1][1] = _d;
}
matrix2x2hpr::matrix2x2hpr( cd a, cd b, cd c, cd _d)
{
d[0][0] = a;
d[0][1] = b;
d[1][0] = c;
d[1][1] = _d;
}
matrix2x2hpr::matrix2x2hpr( const matrix2x2cd& val )
{
d[0][0] = val[0][0];
d[0][1] = val[0][1];
d[1][0] = val[1][0];
d[1][1] = val[1][1];
}
matrix2x2hpr::matrix2x2hpr() :
d( {
{scalar(mpclass(0),mpclass(0)),scalar(mpclass(0),mpclass(0))},
{scalar(mpclass(0),mpclass(0)),scalar(mpclass(0),mpclass(0))}}
)
{
}
matrix2x2hpr::matrix2x2hpr( const matrix2x2hpr& val )
{
*this = val;
}
matrix2x2hpr::scalar matrix2x2hpr::trace() const
{
return d[0][0] + d[1][1];
}
matrix2x2hpr matrix2x2hpr::operator-( const matrix2x2hpr& b ) const
{
matrix2x2hpr res;
for( int i = 0; i < 2; ++i )
for( int j = 0; j < 2; ++j )
{
res.d[i][j] = d[i][j] - b.d[i][j];
}
return res;
}
double matrix2x2hpr::dist(const matrix2x2hpr& matr) const
{
mpclass dist(0);
for( int i = 0 ; i < 2; ++i )
for( int j = 0; j < 2; ++j )
{
mpclass d2 = std::norm( matr.d[i][j] - d[i][j] );
dist += d2;
}
return hprHelpers::toMachine( sqrt(dist) );
}
const matrix2x2hpr& matrix2x2hpr::Id()
{
static matrix2x2hpr Idm(1,0,
0,1);
return Idm;
}
const matrix2x2hpr& matrix2x2hpr::X()
{
static matrix2x2hpr Xm(0,1,
1,0);
return Xm;
}
const matrix2x2hpr &matrix2x2hpr::Y()
{
static matrix2x2hpr Ym(cd(0,0),cd(0,-1),
cd(0,1),cd(0, 0));
return Ym;
}
const matrix2x2hpr &matrix2x2hpr::Z()
{
static matrix2x2hpr Zm(1,0,
0,-1);
return Zm;
}
matrix2x2hpr matrix2x2hpr::operator*( const matrix2x2hpr& b ) const
{
return matrix2x2hpr(
b.d[0][0]* d[0][0]+b.d[1][0]* d[0][1], b.d[0][1]* d[0][0]+b.d[1][1]* d[0][1],
b.d[0][0]* d[1][0]+b.d[1][0]* d[1][1], b.d[0][1]* d[1][0]+b.d[1][1]* d[1][1]);
}
matrix2x2hpr matrix2x2hpr::operator+( const matrix2x2hpr& b ) const
{
matrix2x2hpr res;
for( int i = 0; i < 2; ++i )
for( int j = 0; j < 2; ++j )
{
res.d[i][j] = d[i][j] + b.d[i][j];
}
return res;
}
matrix2x2hpr matrix2x2hpr::operator*( const scalar& val ) const
{
matrix2x2hpr res;
for( int i = 0; i < 2; ++i )
for( int j = 0; j < 2; ++j )
{
res.d[i][j] = d[i][j] * val;
}
return res;
}
matrix2x2hpr matrix2x2hpr::operator*( const mpclass& val ) const
{
matrix2x2hpr res;
for( int i = 0; i < 2; ++i )
for( int j = 0; j < 2; ++j )
{
res.d[i][j] = d[i][j] * val;
}
return res;
}
matrix2x2hpr operator*( const mpclass& val, const matrix2x2hpr& rhs )
{
return rhs * val;
}
matrix2x2hpr operator*( const matrix2x2hpr::scalar& val, const matrix2x2hpr& rhs )
{
return rhs * val;
}
matrix2x2hpr& matrix2x2hpr::operator= ( const matrix2x2hpr& v )
{
for( int i = 0; i < 2; ++i )
for( int j = 0; j < 2; ++j )
{
d[i][j]=v.d[i][j];
}
return *this;
}
matrix2x2hpr matrix2x2hpr::adjoint() const
{
return matrix2x2hpr( std::conj( d[0][0] ), std::conj( d[1][0]),
std::conj( d[0][1]), std::conj( d[1][1]) );
}
//////// template compilation requests //////////////////
template class matrix2x2<int>;
template class matrix2x2<long int>;
template class matrix2x2<mpz_class>;
template class matrix2x2< resring<8> >;
template matrix2x2hpr::matrix2x2hpr( const matrix2x2<int> &val );
template matrix2x2hpr::matrix2x2hpr( const matrix2x2<long int> &val );
template matrix2x2hpr::matrix2x2hpr( const matrix2x2<mpz_class> &val );
template matrix2x2<int>::matrix2x2( const matrix2x2<mpz_class>& );
template matrix2x2<long>::matrix2x2( const matrix2x2<mpz_class>& );
template matrix2x2<mpz_class>::matrix2x2( const matrix2x2<int>& );
template matrix2x2<mpz_class>::matrix2x2( const matrix2x2<long>& );
template matrix2x2< resring<8> >::matrix2x2( const matrix2x2<mpz_class>& );
//////////////// Helper functions //////////////////////
double trace_dist( const matrix2x2hpr& a, const matrix2x2hpr& b )
{
static const mpclass& half = hprHelpers::half();
static const mpclass& one = hprHelpers::one();
// One could use Tailor expansion of the whole formula when a and b close to each other.
// This would give better results in terms of precision.
// We, intstead, take advantage of high precision arithmetic and use straightforward calculation
return hprHelpers::toMachine( sqrt(one - abs( ( a * b.adjoint() ).trace() ) * half) );
}
void convert(const matrix2x2hpr& in, matrix2x2cd &out)
{
for( int i =0; i < 2; ++i )
for( int j =0; j < 2; ++j )
hprHelpers::convert( in.d[i][j],out[i][j] );
}