Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Z rotation approximation algorithm based on number theoretic methods …
…added
  • Loading branch information
vadym-kl committed Jun 15, 2013
1 parent 354dd07 commit 8f8be24
Show file tree
Hide file tree
Showing 21 changed files with 2,209 additions and 72 deletions.
26 changes: 26 additions & 0 deletions CMakeLists.txt
Expand Up @@ -36,8 +36,34 @@ add_executable(sqct
main.cpp
)

add_executable(sqct2
hprhelpers.cpp
resring.cpp
rint.cpp
matrix2x2.cpp
vector2.cpp
output.cpp
gatelibrary.cpp
symbolic_angle.cpp

es/numbersgen.cpp
es/optsequencegenerator.cpp
es/seqlookupcliff.cpp
es/exactdecomposer.cpp

appr/zrot_cache.cpp
appr/toptzrot2.cpp
appr/normsolver.cpp
appr/approxlist.cpp

mainA.cpp
)

set(BOOST_LIBS program_options timer chrono system)
find_package(Boost COMPONENTS ${BOOST_LIBS} REQUIRED)

target_link_libraries(sqct ${Boost_LIBRARIES})
target_link_libraries(sqct gomp pthread mpfr gmpxx gmp rt)

target_link_libraries(sqct2 ${Boost_LIBRARIES})
target_link_libraries(sqct2 gomp pthread pari mpfr gmpxx gmp rt )
8 changes: 7 additions & 1 deletion README
Expand Up @@ -19,4 +19,10 @@ The program code based on results of http://arxiv.org/abs/1206.5236. It also imp
the version of Solovay Kitaev algorithm described in http://arxiv.org/abs/quant-ph/0505030.
In addition to Boost, The GNU Multiple Precision Arithmetic Library, The GNU MPFR Library the library
mpfr::real by Christian Schneider <software(at)chschneider(dot)eu> is used for high precision
floating point arithmetic.
floating point arithmetic.

DIRECTORY STRUCTURE
sk -- implementation of the Solovay-Kitaev algorithm
es -- exact synthesis algorithm
theory -- numerical proof of result from arXiv:1206.5236, tests of exact synthesis algorithm
appr -- optimal round off of unitaries
95 changes: 95 additions & 0 deletions appr/approxlist.cpp
@@ -0,0 +1,95 @@
// Copyright (c) 2013 Vadym Kliuchnikov sqct(dot)software(at)gmail(dot)com
//
// 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/>.
//
// Based on "Practical approximation of single-qubit unitaries by single-qubit quantum Clifford and T circuits"
// by Vadym Kliuchnikov, Dmitri Maslov, Michele Mosca
// http://arxiv.org/abs/1212.6964v1, new version of the paper describing this modification: [to be published]

#include "approxlist.h"

#include <stdexcept>
#include <gmpxx.h>

using namespace std;

typedef __int128 ul;

static long calc_b_min( const approxlist_params& in)
{
hprr tmp = max( -pow2(in.n),
( pow2(in.n) * (in.s - hprr(1.)) - hprr(0.5) ) / sqrt( hprr(2.) ) );

if( !mpfr_fits_slong_p(tmp._x,MPFR_RNDD) )
throw logic_error("too big search space");

return mpfr_get_si(tmp._x,MPFR_RNDD);
}

static long calc_b_max( const approxlist_params& in )
{
hprr tmp = min( pow2(in.n),
( pow2(in.n) * (in.s + hprr(1.)) + hprr(0.5) ) / sqrt( hprr(2.) ) );

if( !mpfr_fits_slong_p(tmp._x,MPFR_RNDU) )
throw logic_error("too big search space");

return mpfr_get_si(tmp._x ,MPFR_RNDU);
}

bool approx_list_builder::norm_condition( long a, long b )
{
ul aa = a;
aa *= a;
ul b2b = b;
b2b *= b;
b2b <<= 1;
return (b2b + aa) <= pow4n;
}

approx_list_builder::approx_list_builder(const approxlist_params &_in, approx_list &_out) :
in(_in),out(_out)
{
if( in.delta > 0.5 || in.delta < 0.0 )
{
cout << "delta:" << in.delta << endl;
throw logic_error("Unsupported value of delta");
}

pow4n = 1; pow4n <<= (in.n * 2);

out.clear();
long b_min = calc_b_min(in);
long b_max = calc_b_max(in);

hprr sqrt2 ( sqrt(hprr(2)) );
hprr sqrt2st ( sqrt(hprr(2)) * hprr(in.step) );
hprr s2n ( in.s * pow2(in.n) );
hprr t ( s2n - hprr(b_min + in.init) * sqrt2 );
long double sld = to_ld(in.s);

hprr r = 0.;

int st = in.step;
for( long b = (b_min + in.init) ; b <= b_max ; b+=st )
{
long a = mpfr_get_si(t._x, MPFR_RNDN);
r = t - a;
if( abs(r) < in.delta && norm_condition(a,b) )
out.push_back( make_pair(to_ld(r) *sld,b) );
t -= sqrt2st;
}
}
76 changes: 76 additions & 0 deletions appr/approxlist.h
@@ -0,0 +1,76 @@
// Copyright (c) 2013 Vadym Kliuchnikov sqct(dot)software(at)gmail(dot)com
//
// 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/>.
//
// Based on "Practical approximation of single-qubit unitaries by single-qubit quantum Clifford and T circuits"
// by Vadym Kliuchnikov, Dmitri Maslov, Michele Mosca
// http://arxiv.org/abs/1212.6964v1, new version of the paper describing this modification: [to be published]

#ifndef APPROXLIST_H
#define APPROXLIST_H

#include "hprhelpers.h"
#include <vector>

/// \brief Specifies input for approx_list_builder struct.
/// See approx_list_builder struct for the meaning of each member.
struct approxlist_params
{
int n;
hprr s;
hprr delta;
int step;
int init;

approxlist_params( int _n , hprr _s, hprr _delta, int _step = 1, int _init = 0 ) :
n( _n ), s(_s), delta(_delta), step(_step), init(_init)
{}

};

typedef std::vector< std::pair< long double, long > > approx_list;

/// \brief Finds an approximations of a real number <in.s>
/// using numbers of the from ( a + \sqrt{2} b ) / 2^<in.n>;
/// b is chosen from the set ( b_min + <in.init> + N * <in.step> ), where N is a non-negative integer;
/// b_min calculated based on <in.n> and <in.s>;
/// a is an integer;
///
/// <in.init> and <in.step> members we introduced to run approximation using multiple threads
///
/// For any approximation found it is true that:
/// | a + \sqrt{2} b - 2^<in.n> * s | <= <in.delta>
/// and a^2 + b^2 <= 4^<in.n>
///
/// Implemented algorithm supports only <in.delta> < 0.5. In this case a uniquely defined by b
///
/// The result of the execution is a vector of pairs:
/// ( ( 2^<in.n> * <in.s> - ( a + \sqrt{2} b ) ) * <in.s> , b )
///
/// in.<name> refers to the member named <name> of approxlist_params struct
struct approx_list_builder
{
approx_list_builder( const approxlist_params& in,
approx_list& out );
bool norm_condition( long a, long b );

const approxlist_params& in;
approx_list& out;
__int128 pow4n;

};

#endif // APPROXLIST_H
162 changes: 162 additions & 0 deletions appr/normsolver.cpp
@@ -0,0 +1,162 @@
// Copyright (c) 2013 Vadym Kliuchnikov sqct(dot)software(at)gmail(dot)com
//
// 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/>.
//
// Based on "Practical approximation of single-qubit unitaries by single-qubit quantum Clifford and T circuits"
// by Vadym Kliuchnikov, Dmitri Maslov, Michele Mosca
// http://arxiv.org/abs/1212.6964v1, new version of the paper describing this modification: [to be published]

#include "normsolver.h"
#include <iostream>

using namespace std;

normSolver::normSolver()
{
pari_init(8000000l * 128l, 500509);
rnf = gp_read_str("rnfisnorminit(z^2-2,x^2+1)");
}

static void genToMpz( GEN gen, mpz_class& out )
{
if( gen != 0 )
{
char* gen_str = GENtostr(gen);
out = mpz_class(gen_str);
pari_free(gen_str);
}
else
out = 0;
}

static void genToMpz( GEN gen, mpz_class& out, bool& div )
{
div = false;
if( gen != 0 )
{
if( typ(gen) == t_INT )
{
char* gen_str = GENtostr(gen);
out = mpz_class(gen_str);
pari_free(gen_str);
}
else if( typ(gen) == t_FRAC )
{
char* gen_str = GENtostr(gel(gen,1));
out = mpz_class(gen_str);
pari_free(gen_str);
div = true;
if( ! gequalgs(gel(gen,2),2) )
throw std::logic_error("unexpected denominator");
}
else
throw std::logic_error("wrong data type");
}
else
out = 0;
}


bool normSolver::solve(const ring_int_real<mpz_class>& rhs, ring_int<mpz_class> &res) const
{
stringstream ss;
ss << "(" << rhs[0] << ")+(" << rhs[1] << ")*z";
GEN in = gp_read_str(ss.str().c_str());
GEN solution = rnfisnorm(rnf,in,0);

bool ok = gequal1(gel(solution,2));
if( !ok ) //there is no solution
return false;

GEN re_a = 0, re_b = 0, im_a = 0, im_b = 0;
GEN sln = gel(gel(solution,1),2);

if( degree(sln) >= 0 )
{
GEN re = gel(gel(sln,2),2);
if( degree(re) >= 0 )
re_a = gel(re,2);
if( degree(re) == 1 )
re_b = gel(re,3);
}

if( degree(sln) == 1 )
{
GEN im = gel(gel(sln,3),2);
if( degree(im) >= 0 )
im_a = gel(im,2);
if( degree(im) == 1 )
im_b = gel(im,3);
}

genToMpz( re_a, res[0] );
genToMpz( im_a, res[2] );

bool re_b_fr, im_b_fr;
mpz_class re_bz, im_bz;
genToMpz( re_b, re_bz, re_b_fr );
genToMpz( im_b, im_bz, im_b_fr );

if( re_b_fr && im_b_fr )
{
res[1] = re_bz + im_bz;
res[1] /= 2;
res[3] = im_bz - re_bz;
res[3] /= 2;
}
else if( !im_b_fr && !re_b_fr )
{
res[1] = re_bz + im_bz;
res[3] = im_bz - re_bz;
}
else
{
cout << ss.str().c_str() << endl;
pari_printf("(%Ps)+Sqrt[2]*(%Ps)+I*((%Ps)+Sqrt[2]*(%Ps))\n",re_a,re_b,im_a,im_b);
throw std::logic_error("unexpected solution");
}


return true;
}

bool normSolver::solve( const ring_int<mpz_class>& u00, int denompower, m& matr ) const
{
mpz_class pow2 = 1;
pow2 <<= denompower * 2;
ring_int_real<mpz_class> norm(u00.abs2());
if( norm[0] > pow2 )
return false;
norm[0] = pow2 - norm[0];
norm[1] = -norm[1];
norm[3] = -norm[3];

ring_int<mpz_class> u10;
if( !solve(norm,u10) )
return false;

matr.set(u00,-u10.conjugate(),
u10, u00.conjugate(), denompower * 2);

matr.reduce();
return true;
}

const normSolver &normSolver::instance()
{
static normSolver inst;
return inst;
}

0 comments on commit 8f8be24

Please sign in to comment.