From 8f8be241b51e34296e4e908bbd0fbacb69d952ee Mon Sep 17 00:00:00 2001 From: Vadym Kliuchnikov Date: Fri, 14 Jun 2013 20:39:15 -0400 Subject: [PATCH] Z rotation approximation algorithm based on number theoretic methods added --- CMakeLists.txt | 26 ++ README | 8 +- appr/approxlist.cpp | 95 ++++++ appr/approxlist.h | 76 +++++ appr/normsolver.cpp | 162 ++++++++++ appr/normsolver.h | 52 +++ appr/toptzrot2.cpp | 697 +++++++++++++++++++++++++++++++++++++++++ appr/toptzrot2.h | 188 +++++++++++ appr/zrot_cache.cpp | 364 +++++++++++++++++++++ appr/zrot_cache.h | 86 +++++ es/exactdecomposer.cpp | 23 +- es/exactdecomposer.h | 14 +- gatelibrary.cpp | 57 ++-- gatelibrary.h | 17 +- hprhelpers.cpp | 45 ++- hprhelpers.h | 22 +- mainA.cpp | 163 ++++++++++ matrix2x2.cpp | 117 ++++++- matrix2x2.h | 28 +- symbolic_angle.cpp | 24 ++ symbolic_angle.h | 17 + 21 files changed, 2209 insertions(+), 72 deletions(-) create mode 100644 appr/approxlist.cpp create mode 100644 appr/approxlist.h create mode 100644 appr/normsolver.cpp create mode 100644 appr/normsolver.h create mode 100644 appr/toptzrot2.cpp create mode 100644 appr/toptzrot2.h create mode 100644 appr/zrot_cache.cpp create mode 100644 appr/zrot_cache.h create mode 100644 mainA.cpp create mode 100644 symbolic_angle.cpp create mode 100644 symbolic_angle.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 2d90582..ca0dbaf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 ) diff --git a/README b/README index 1da96a2..f04b715 100644 --- a/README +++ b/README @@ -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 is used for high precision -floating point arithmetic. \ No newline at end of file +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 \ No newline at end of file diff --git a/appr/approxlist.cpp b/appr/approxlist.cpp new file mode 100644 index 0000000..d26b23c --- /dev/null +++ b/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 . +// +// 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 +#include + +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; + } +} diff --git a/appr/approxlist.h b/appr/approxlist.h new file mode 100644 index 0000000..f0364fc --- /dev/null +++ b/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 . +// +// 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 + +/// \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 +/// using numbers of the from ( a + \sqrt{2} b ) / 2^; +/// b is chosen from the set ( b_min + + N * ), where N is a non-negative integer; +/// b_min calculated based on and ; +/// a is an integer; +/// +/// and members we introduced to run approximation using multiple threads +/// +/// For any approximation found it is true that: +/// | a + \sqrt{2} b - 2^ * s | <= +/// and a^2 + b^2 <= 4^ +/// +/// Implemented algorithm supports only < 0.5. In this case a uniquely defined by b +/// +/// The result of the execution is a vector of pairs: +/// ( ( 2^ * - ( a + \sqrt{2} b ) ) * , b ) +/// +/// in. refers to the member named 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 diff --git a/appr/normsolver.cpp b/appr/normsolver.cpp new file mode 100644 index 0000000..5eb1625 --- /dev/null +++ b/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 . +// +// 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 + +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& rhs, ring_int &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& u00, int denompower, m& matr ) const +{ + mpz_class pow2 = 1; + pow2 <<= denompower * 2; + ring_int_real norm(u00.abs2()); + if( norm[0] > pow2 ) + return false; + norm[0] = pow2 - norm[0]; + norm[1] = -norm[1]; + norm[3] = -norm[3]; + + ring_int 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; +} diff --git a/appr/normsolver.h b/appr/normsolver.h new file mode 100644 index 0000000..759191d --- /dev/null +++ b/appr/normsolver.h @@ -0,0 +1,52 @@ +// 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 . +// +// 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 NORMSOLVER_H +#define NORMSOLVER_H + +#include "gatelibrary.h" +#include "rint.h" + +#include + +/// \brief Solves norm equation |z|^2 = A + \sqrt{2} B, where z is from Z[exp(i pi/4)]. +/// PARI (http://en.wikipedia.org/wiki/PARI/GP) C interface is used for this. +class normSolver +{ +public: + typedef matrix2x2 m; + normSolver(); + + /// \brief rsh specifies A + \sqrt{2}B and res contains one of possible solutions to the norm equation + bool solve( const ring_int_real& rhs, ring_int& res ) const; + + /// \brief Uses norm equation solver to find an entry y of the unitary + /// 1/2^denompower * {{u00,-y^{\dagger}},{y,u00^{\dagger}}}; + /// \return true if such a unitary exists and writes it down into matr + bool solve( const ring_int& u00, int denompower, m& matr ) const; + + /// \brief Instance of normSolver to be used for all calls; allows to avoid multiple initialization of PARI + static const normSolver& instance(); +private: + GEN rnf; ///< PARI object for the extension Z[exp(i pi/4)] / Z[sqrt{2}] +}; + +#endif // NORMSOLVER_H diff --git a/appr/toptzrot2.cpp b/appr/toptzrot2.cpp new file mode 100644 index 0000000..c0da30d --- /dev/null +++ b/appr/toptzrot2.cpp @@ -0,0 +1,697 @@ +// 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 . +// +// 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 "toptzrot2.h" +#include +#include +#include +#include "output.h" +#include + +#include +#include + +#include + +using namespace std; +using namespace std::chrono; +typedef matrix2x2 matr; + +////////////////////////////////////////////////////////////// + +unitary::unitary(const matr &m) : + x(m.d[0][0]),y(m.d[1][0]),n(m.de), k(0) +{ + matr sp(m); + sp.d[1][1] = m.d[0][0].conjugate(); + sp.d[0][1] = - m.d[1][0].conjugate(); + + for( int i = 0; i < 8; ++i ) + { + if( sp * m.T(i) == m ) + { + k = i; + break; + } + } +} + +unitary::operator matr() +{ + matr tmp(x, - y.conjugate() * y.omega(k) , + y, x.conjugate() * y.omega(k) ,n); + tmp.reduce(); + return tmp; +} + +std::ostream& operator<<(std::ostream& out, const unitary& u ) +{ + out << u.n << ","; + for( int i = 0; i < 4; ++i ) + out << u.x[i] << ","; + for( int i = 0; i < 4; ++i ) + out << u.y[i] << ","; + out << u.k ; + return out; +} + +////////////////////////////////////////////////////////////// + +static const int maximum_sde2 = 250; + +struct time_logger +{ + time_logger( const string& event_name ) : name(event_name) + { + t1 = high_resolution_clock::now(); + } + + ~time_logger() + { + high_resolution_clock::time_point t2 = high_resolution_clock::now(); + ofstream ofs("time.log",ios_base::app ); + duration time_span = duration_cast>(t2 - t1); + //ofs << "[" ; ofs.width(40); + ofs << name << ","; + //ofs.width(20); + ofs << time_span.count() << endl; + } + + high_resolution_clock::time_point t1; + string name; +}; + +//////////////////////////////////// + +static long integer_part( const hprr& num, long b ) +{ + hprr tmp = ( num - hprHelpers::sqrt2() * hprr(b) ); + if( ! mpfr_fits_slong_p( tmp._x, MPFR_RNDN ) ) + throw std::logic_error("overflow"); + return mpfr_get_si( tmp._x, MPFR_RNDN ); +} + +static ring_int_real sqrt2_conjugate( const ring_int_real& val ) +{ + return ring_int_real(val[0],-val[1]); +} + +//////////////////////////////////// + +template< class T> +static inline void sort( vector< T >& vec ) +{ + sort(begin(vec),end(vec)); +} + +static pair vec_values_range( const vector< pair< long double, long> >& vec ) +{ + return make_pair( vec.front().first, (vec.end() - 1)->first ); +} + +static inline auto lower_bound( const vector< pair< long double, long> >& vec, long double val ) -> decltype(vec.begin()) +{ + return lower_bound(begin(vec),end(vec), make_pair(val,0l)); +} + +static inline auto upper_bound( const vector< pair< long double, long> >& vec, long double val ) -> decltype(vec.begin()) +{ + return upper_bound(begin(vec),end(vec), make_pair(val,0l)); +} + +static long double vector_threshold( const vector< pair< long double, long> >& vec ) +{ + pair zp(0,0); + auto r1 = lower_bound( begin(vec), end(vec), zp ); + auto r2 = upper_bound( begin(vec), end(vec), zp ); + if( r1 == vec.end() || r2 == vec.end() ) + throw std::logic_error("initial threshold estimate failed"); + return max( abs(r1->first), abs(r2->first) ); +} + +/////////////////////////// + +toptzrot2::toptzrot2(const approx_params& _in, approx_result& _out) : + in(_in), out(_out), precisions(maximum_sde2,2.0), matrices(maximum_sde2), + solver( normSolver::instance() ) +{ + theta[0] = - in.phi / hprHelpers::two(); + theta[1] = - in.phi / hprHelpers::two() + hprHelpers::pi() / hprr(8.0); + n = ceil( (in.sde2 + 2.) / 4. ); + delta = in.epsilon * pow2(n) * sqrt(hprr(2.0)); + + for( int i = 0; i < 2 ;++i ) + { + pow2nSinTh[i] = pow2(n) * sin ( theta[i] ); + pow2nCosTh[i] = pow2(n) * cos ( theta[i] ); + } + + sde2min = 4 * n - 5; + sde2max = 4 * n - 2; + + std::cout << "sde2 in [" << sde2min << "," << sde2max << "]" << endl; + + RzTh(theta[0],target); + + search(); +} + +void toptzrot2::search() +{ + +//before parallelization: +// { +// time_logger l("approx_list_builder," + to_string(n)); +// for( int i = 0; i < 2; ++i ) +// { +// approx_list_builder lb0({n, cos(theta[i]),delta},cosTheta[i]); +// approx_list_builder lb1({n, sin(theta[i]),delta},sinTheta[i]); +// } +// } + + + { + time_logger l("approx_list_builder," + to_string(n)); + int threads = omp_get_max_threads(); + cout << "Using " << threads << " threads" << endl; + + vector< approx_list > cl(threads * 2); + vector< approx_list > sl(threads * 2); + + #pragma omp parallel + //for( int thid = 0; thid < threads; ++thid ) + { + int thid = omp_get_thread_num(); + for( int i = 0; i < 2; ++i ) + { + approx_list_builder lb0({n, cos(theta[i]),delta,threads,thid},cl[ 2 * thid + i]); + approx_list_builder lb1({n, sin(theta[i]),delta,threads,thid},sl[ 2 * thid + i]); + } + } + + size_t ssz[2] = {0,0}; + size_t csz[2] = {0,0}; + + for( int th = 0; th < threads ; ++th ) + { + for( int i = 0; i < 2; ++i ) + { + csz[i] += cl[ 2 * th + i].size(); + ssz[i] += sl[ 2 * th + i].size(); + } + } + + approx_list::iterator ci[2]; + approx_list::iterator si[2]; + + for( int i = 0; i < 2; ++i ) + { + cosTheta[i].resize(csz[i]); + ci[i] = cosTheta[i].begin(); + + sinTheta[i].resize(ssz[i]); + si[i] = sinTheta[i].begin(); + } + + for( int th = 0; th < threads ; ++th ) + { + for( int i = 0; i < 2; ++i ) + { + ci[i] = copy(cl[ 2 * th + i].begin(),cl[ 2 * th + i].end(),ci[i]); + si[i] = copy(sl[ 2 * th + i].begin(),sl[ 2 * th + i].end(),si[i]); + } + } + + } + + { + time_logger l("sort," + to_string(n)); + for( int i = 0; i < 2; ++i ) + { + sort(cosTheta[i]); + sort(sinTheta[i]); + } + } + + { + time_logger l("combine_lists," + to_string(n)); + combine_lists(); + } +} + +void toptzrot2::combine_lists() +{ + bool not_done = true; + long double search_threshold = to_ld( pow2(n) * in.epsilon * in.epsilon ); + + long double eps1 = 0; + long double eps2 = min( search_initial_threshold(), search_threshold ); + + while( not_done && eps1 <= search_threshold ) + { + candidates.clear(); + for( int i = 0; i < 2; ++i ) + add_candidates_in_range( i, eps1, eps2 ); + + sort( candidates ); + for( const auto& t : candidates ) + { + if( try_solution( get<1>(t), get<2>(t) ) ) + not_done = ! done(); + } + + eps1 = eps2; + eps2 *= 2.0; + } + + //if( !not_done ) + { + cout << "epsilon = " << in.epsilon << endl; + for( int k = sde2min; k <= sde2max; ++k ) + { + cout << k << "," << precisions[k] ; + auto cpr = precisions[k]; + if( cpr > in.epsilon ) + { + if( cpr < 1. ) + { + out.set_hint(k,precisions[k],matrices[k]); + cout << "( > epsilon )" ; + } + out.set_certified_epsilon(k,to_ld(in.epsilon)); + } + else + out.set_result(k,precisions[k],matrices[k]); + + cout << endl; + } + + for( size_t k = sde2max + 1; k < precisions.size(); ++k ) + { + if( precisions[k] < 1. ) + { + out.set_hint(k,precisions[k],matrices[k]); + cout << "hint:" << k << "," << precisions[k] << endl; + } + } + //return; + } + + if( not_done ) + cout << "search threshold reached" << endl; +} + +long double toptzrot2::search_initial_threshold() +{ + long double initial_threshold[2]; + for( int i = 0; i < 2; ++i ) + { + initial_threshold[i] = max( vector_threshold(sinTheta[i]), + vector_threshold(cosTheta[i])); + } + return max( initial_threshold[0], initial_threshold[1] ) * 16.; +} + +bool toptzrot2::try_solution( long b, long d ) +{ + int theta_id = b & 1; + //note: we use last bit of b and d to store which angle they approximate : Theta1 or Theta2 + hprr pow2nCT = pow2nCosTh[ theta_id ]; + hprr pow2nST = pow2nSinTh[ theta_id ]; + b >>= 1; + d >>= 1; + long a = integer_part( pow2nCT, b ); + long c = integer_part( pow2nST, d ); + + mpz_class max = 1; + max <<= 2 * n; + + ring_int rr(a,d+b,c,d-b); + auto rrabs2 = rr.abs2(); + + if( rrabs2[0] > max ) + return false; + + if( abs(rrabs2.toComplex(0)) > pow2(2*n) ) + return false; + + if( abs(sqrt2_conjugate(rrabs2).toComplex(0)) > pow2(2*n) ) + return false; + + matr m; + int sde2 = 4 * n - rrabs2.gde(); + + if( sde2 < sde2min ) + return false; + + if(precisions[sde2] > 1. ) + { + if( solver.solve(rr,n,m) ) + { + if( theta_id == 0 ) + matrices[sde2] = m; + else + matrices[sde2] = m * m.T(1); + + + precisions[sde2] = trace_dist(target,matrices[sde2]); + return true; + } + } + return false; +} + +bool toptzrot2::done() +{ + for( int k = sde2min; k <= sde2max; ++k ) + if( precisions[k] > 1. ) + return false; + return true; +} + +void toptzrot2::add_candidates_in_range(int i, long double eps1, long double eps2) +{ + auto im_range = vec_values_range( sinTheta[i] ); + // auto re_range = vec_values_range( cosTheta[i] ); + // cout << im_range.first << ":" << im_range.second << endl; + // cout << re_range.first << ":" << re_range.second << endl; + long double re_min = eps1 - im_range.second; + long double re_max = eps2 - im_range.first; + auto re_beg = lower_bound(cosTheta[i], re_min ); + auto re_end = upper_bound(cosTheta[i], re_max ); + for( auto k = re_beg; k < re_end; ++k ) + { + long double im_min = eps1 - k->first; + long double im_max = eps2 - k->first; + auto im_beg = lower_bound( sinTheta[i], im_min ); + auto im_end = upper_bound( sinTheta[i], im_max ); + for( auto j = im_beg; j < im_end; ++j ) + candidates.push_back( make_tuple( k->first + j->first, + (k->second << 1) + i, + (j->second << 1) + i ) ); + } +} + + + +approx_result_block::approx_result_block() : + precisions(maximum_sde2,2.), + matrices(maximum_sde2) +{ +} + + +approx_result::approx_result() : + certified_epsilon(maximum_sde2,2.) +{ +} + +void approx_result::set_hint(int sde2, long double precision, const matr &m) +{ + long double cpr = hints.precisions[sde2]; + if( cpr < precision ) + { + hints.precisions[sde2] = precision; + hints.matrices[sde2] = m; + } +} + +void approx_result::set_certified_epsilon(int sde2, long double eps) +{ + if( certified_epsilon[sde2] > eps ) + certified_epsilon[sde2] = eps; +} + +void approx_result::set_result(int sde2, long double precision, const matr &m) +{ + precisions[sde2] = precision; + matrices[sde2] = m; +} + + +optzrot_driver::optzrot_driver(const optzrot_driver_params &_in, approx_result &_out) : + in(_in), out(_out) +{ + search(); +} + +static void print_vec( const std::vector& v) +{ + for( auto i : v ) + cout << i << "," ; + cout << endl; +} + +static long double update( long double epsilon, int n ) +{ + long double delta = to_ld( hprr(epsilon) * pow2(n) * hprHelpers::sqrt2() ); + if( delta > 0.5 ) + { + epsilon = to_ld( hprr(0.4999999) / ( pow2(n) * hprHelpers::sqrt2() ) ); + cout << "epsilon was updated to:" << epsilon << endl; + } + return epsilon; +} + +void optzrot_driver::search() +{ + out.precisions[0] = trace_dist(matrix2x2hpr::Id(), RzTh(in.phi) ); + out.matrices[0] = matr::Id(); + + int n = ceil( ( in.sde2initial + 2.0 ) / 4. ); + int n_max = ceil( ( in.sde2max + 2.0 ) / 4. ); + + for( ; n <= n_max; ++n ) + { + long double epsilon = *min_element(out.precisions.begin(), + out.precisions.begin() + 4*n - 1 ); + cout << "epsilon:" << epsilon << endl; + epsilon = update( epsilon, n ); //in the case if algorithm is not applicable for given precision + + approx_params inp{in.phi,epsilon, 4 * n - 2 }; + try + { + + toptzrot2 tr(inp,out); + } + catch(const std::exception& e ) + { + cerr << in.phi << "," << n << "," << epsilon << "," << ":" << e.what() << endl; + } + //here we need to decide which epsilon to feed next + + cout << "epsilon:" << epsilon << endl; + } + +} + + +unitary::unitary() : + x(1,0,0,0), y(0,0,0,0), n(0), k(0) +{ + +} + + + +toptimizer::toptimizer(const approx_result &_in, topt_result& _out ) : + in(_in), out(_out) +{ + for (size_t i = 0; i < out.size(); ++i) + { + auto& o = out[i]; + + auto eps = min_element(out.begin(), out.begin() + i, + [] ( const result_entry& a, const result_entry& b ) + { + return a.precision < b.precision; + } + )->precision; + + if( in.precisions[i] < eps ) + { + o.precision = in.precisions[i]; + o.u = optimize_t_gates(in.matrices[i]); + matr m(o.u); + o.t_gates = m.t(); + o.delta_t = m.t() - m.h() + 1; + } + else if( in.precisions[i] < 1.0 ) + { + o.cert_precision = eps; + o.cert_gap = false; + } else { + if( in.certified_epsilon[i] < 1.0 ) + { + o.cert_precision = in.certified_epsilon[i]; + o.cert_gap = in.certified_epsilon[i] < eps ; + if( o.cert_gap ) + o.required_precision = eps; + } + } + } +} + +unitary toptimizer::optimize_t_gates( const matr& m ) +{ + int min_i = 0; + bool is_conjugate = false; + unitary u(m); + + int t_min = m.t(); + + for( int i = 0; i < 8; ++i ) + { + u.y.mul_eq_w(); + int tc = matr(u).t(); + if( t_min > tc ) + { + + t_min = tc; + min_i = i + 1; + } + } + + u.y.conjugate_eq(); + + for( int i = 0; i < 8; ++i ) + { + u.y.mul_eq_w(); + int tc = matr(u).t(); + if( t_min > tc ) + { + t_min = tc; + min_i = i + 1; + is_conjugate = true; + } + } + + u = (unitary)m; + if( is_conjugate ) + u.y.conjugate_eq(); + u.y.mul_eq_w(min_i); + + return u; +} + + +topt_result::topt_result() : + vector(maximum_sde2) +{ +} + +void topt_result::write(const string &filename) +{ + size_t max_items = 0; + for( size_t i = 0; i < size(); ++i ) + if( at(i).precision < 1.0 ) + max_items = i + 1; + + ofstream fout(filename + ".data", ios_base::app ); + ofstream fc(filename , ios_base::app ); + ofstream ftc(filename + ".tc.csv", ios_base::app ); + + for( size_t i = 0; i < max_items; ++i ) + { + if( at(i).precision < 1.0 || at(i).cert_gap ) + { + fout << angle << "," << i << "," << at(i) << endl; + if( at(i).precision < 1.0 ) + { + fc << angle.numerator << "," << angle.denominator << "," << + at(i).precision << "," << exactDecomposer::decompose(at(i).u) << endl; + } + + if( at(i).precision < 1.0 && at(i).cert_gap ) + { + ftc << at(i).precision << "," << at(i).t_gates << endl; + } + } + } +} + + + + + +result_entry::result_entry() : + precision(2.0), + cert_precision(2.0), + cert_gap(false), + required_precision(2.0), + t_gates(0), + delta_t(0) +{ +} + +std::ostream& operator<<(std::ostream& out, const result_entry& re ) +{ + out << re.precision << "," << re.t_gates << "," << re.delta_t << "," + << re.cert_precision << "," << re.cert_gap << "," << re.required_precision << "," + << re.u; + return out; +} + + + + + +void topt_app::run() +{ + ifstream ifs; + ifs.exceptions ( std::ifstream::badbit ); + ifs.open(params.in_filename); + vector angles; + while( ifs ) + { + string line; + getline(ifs,line); + if(line.size() == 0 ) + continue; + if( line[0] == '#' ) + continue; + + istringstream ss(line); + symbolic_angle a{0,1,true}; + //char delim; + ss >> a.numerator >> a.denominator >> a.pi; + angles.push_back(a); + } + + cout << "scheduled to find approximations for " << angles.size() << " unitaries"; + + for( const symbolic_angle& a : angles ) + { + optzrot_driver_params ap; + + ap.sde2initial = 15; + ap.sde2max = params.max_sde; + ap.phi = a; + + approx_result ares; + optzrot_driver drv(ap,ares); + + topt_result out; + out.angle = a; + toptimizer tpz(ares,out); + out.write(params.out_filename); + } +} diff --git a/appr/toptzrot2.h b/appr/toptzrot2.h new file mode 100644 index 0000000..729464a --- /dev/null +++ b/appr/toptzrot2.h @@ -0,0 +1,188 @@ +// 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 . +// +// 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 TOPTZROT2_H +#define TOPTZROT2_H + +#include "approxlist.h" +#include "rint.h" +#include "matrix2x2.h" +#include "normsolver.h" +#include "symbolic_angle.h" +#include +#include + + +typedef matrix2x2 matr; + +/// \brief Represents a unitary matrix +/// 1 / sqrt{2}^n {{x,-y* w^k},{y, x* w^k}}, for w = Exp[i pi/4] +struct unitary +{ + unitary(); + explicit unitary( const matr& m ); + operator matr(); + ring_int x; + ring_int y; + int n; + int k; +}; + +std::ostream& operator<<(std::ostream& out, const unitary& re ); + +///////////////////////////////////////////////////////////////////// + +struct approx_params +{ + hprr phi; + hprr epsilon; + int sde2; +}; + +struct approx_result_block +{ + approx_result_block(); + std::vector precisions; + std::vector< matr > matrices; +}; + +struct approx_result : public approx_result_block +{ + approx_result(); + void set_hint( int sde2, long double precision, const matr& m ); + void set_certified_epsilon( int sde2, long double eps ); + void set_result( int sde2, long double precision, const matr& m ); + + approx_result_block hints; // ocasionally found results or results produced by other methods + std::vector certified_epsilon; //if certified_epsilon[k] > 0 there is no approximation better than epsilon below certified_epsilon[k] for sde2 = k +}; + +/// \brief Finds an approximation of R_z() rotation by unitary over the ring Z[i,1/sqrt{2}]. +/// +/// Searches for unitaries with sde in the range [,]. +/// Considers unitaries of the form 1/2^n {{x,-y* w^k},{y, x* w^k}}, for w = Exp[i pi/4] and x,y from Z[w]. +/// Uses global phase invariant distance to determine the quality of the approximation. +/// +/// This results in two approximation subproblems corresponding to +/// theta[0] = - / 2 and +/// theta[1] = - / 2 + pi / 8 +/// +/// Each subproblem is stated as following: +/// Find number x over the ring Z[i,1/sqrt{2}] such that: +/// +struct toptzrot2 +{ + toptzrot2(const approx_params& in, approx_result& out); + void search(); + void combine_lists(); + long double search_initial_threshold(); + bool try_solution( long b, long d ); + bool done(); + void add_candidates_in_range( int i,long double eps_min, long double eps_max ); + + const approx_params& in; + approx_result& out; + + hprr theta[2]; + approx_list cosTheta[2];///< cos(theta[i]) + approx_list sinTheta[2];///< sin(theta[i]) + hprr delta; + int n; + int sde2min; + int sde2max; + + hprr pow2nCosTh[2]; + hprr pow2nSinTh[2]; + + std::vector< std::tuple > candidates; + + std::vector precisions; + std::vector< matr > matrices; + matrix2x2hpr target; + const normSolver& solver; +}; + +///////////////////////////////////////////////////////////////////// + +struct optzrot_driver_params +{ + int sde2initial; + int sde2max; + hprr phi; +}; + +struct optzrot_driver +{ + optzrot_driver( const optzrot_driver_params& in, approx_result& out ); + + const optzrot_driver_params& in; + approx_result& out; + void search(); +}; + +struct result_entry +{ + result_entry(); + long double precision; //precision achieved by unitary + long double cert_precision; //if there is no unitary, shows which precision was assumed + bool cert_gap; + long double required_precision; + unitary u; //U[x,y,k] + int t_gates; // number of gates we off from minimum for given x in U[x,y,k] + int delta_t; +private: + result_entry( const result_entry& ); +}; + +std::ostream& operator<<(std::ostream& out, const result_entry& re ); + +struct topt_result : public std::vector< result_entry > +{ + topt_result(); + symbolic_angle angle; + void write( const std::string& filename ); +}; + +/// \brief tweaks results of toptzrot2 to achieve sde2 - 2 or sde2 - 1 T gates +struct toptimizer +{ + toptimizer( const approx_result& in, topt_result& out ); + unitary optimize_t_gates( const matr& m ); + const approx_result& in; + topt_result& out; +}; + +struct topt_app_params +{ + int max_sde; + std::string in_filename; + std::string out_filename; +}; + +struct topt_app +{ + topt_app( const topt_app_params& in ) : + params(in) {} + void run(); + const topt_app_params& params; +}; + +#endif // TOPTZROT2_H diff --git a/appr/zrot_cache.cpp b/appr/zrot_cache.cpp new file mode 100644 index 0000000..41d4017 --- /dev/null +++ b/appr/zrot_cache.cpp @@ -0,0 +1,364 @@ +// 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 . +// +// 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 "zrot_cache.h" + +#include +#include +#include +#include "output.h" + +using namespace std; + + +static void normalize( step_function &sf1 ) +{ + step_function tmp; + vector tc(sf1.size()); + for( size_t k = 0; k < sf1.size(); ++k ) + tc[k] = ((matr)sf1[k].second).t(); + + tmp.push_back(sf1[0]); + for( size_t k = 1; k < sf1.size(); ++k ) + { + if( tc[k - 1] > tc[k] ) + tmp.push_back(sf1[k]); + } + + swap(tmp,sf1); +} + +zrot_cache::zrot_cache() : recalculate(false) +{} + +zrot_cache::zrot_cache(const string &filename, bool rec ) : recalculate(rec) +{ + read_file(filename); +} + +void zrot_cache::read_file(const std::string &filename) +{ + ifstream ifs; + ifs.exceptions ( std::ifstream::badbit ); + ifs.open(filename); + while( ifs ) + { + string line; + getline(ifs,line); + + if(line.size() == 0 ) + continue; + if( line[0] == '#' ) + continue; + + //cout << line << endl; + + istringstream ss(line); + symbolic_angle a{0,1,false}; + circuit c; + double file_precision; + char delim; + ss >> a.numerator >> delim >> a.denominator >> delim >> file_precision >> delim; + c.fromStream(ss,false); + if(recalculate) + (*this)[a].push_back({trace_dist(Rz(a),c),c}); + else + (*this)[a].push_back({file_precision,c}); + } + + for ( auto& sf : *this ) + { + sort(sf.second.begin(),sf.second.end()); + normalize(sf.second); + } +} + +circuit zrot_cache::lookup(const symbolic_angle &angle, double precision) const throw ( not_found_exception ) +{ + auto it = find(angle); + if ( it == end() ) + throw not_found_exception(angle,1.); + + const auto& sf = it->second; + step_function::value_type v{precision,{}}; + auto lb = std::lower_bound( std::begin(sf), std::end(sf), v ); + if( lb == std::begin(sf) ) + throw not_found_exception(angle,precision); + + return (lb - 1)->second; +} + +std::ostream& operator<< ( std::ostream& out , const zrot_cache& zc ) +{ + out << "Precision, T count, circuit" << endl; + for( const auto& v : zc ) + { + out << v.first << ":" << endl; + + for( const auto& r : v.second ) + { + int TC = gateLibrary::toCliffordT(r.second.count())[gateLibrary::T]; + //out << " " << r.first << "," << TC << "," << r.second << endl; + out.setf ( ios_base::fixed ); + out.precision( 100 ); + out << "{" << r.first << "," << TC << "},"<< endl; + } + } + + return out; +} + + + + +void DiffApplication::run() +{ + if( params.filenames.size() != 2 ) + throw logic_error("wrong number of inputs"); + zrot_cache zc1(params.filenames[0],true); + zrot_cache zc2(params.filenames[1],true); + + cout << "File contents:" << endl; + + cout << "file:" << params.filenames[0] << endl; + cout << zc1 << endl; + cout << "**************************" << endl; + + cout << "file:" << params.filenames[1] << endl; + cout << zc2 << endl; + cout << "**************************" << endl; + + // find the angles that are common for both files + + vector angles; + + for( const auto& rec : zc1 ) + { + if(zc2.count(rec.first) > 0 ) + angles.push_back(rec.first); + } + + cout << "Angles in common:" << angles.size() << endl; + for( const auto& angle : angles ) + cout << " " << angle << endl; + + for( const auto& angle : angles ) + compare(zc1[angle],zc2[angle],angle); + +} + +const long double low_k = 0.99999999; +const long double high_k = 1.00000001; + + + +static bool belongs( const step_function &sf1, long double val) +{ + auto en = upper_bound(sf1.begin(), sf1.end(), pair{val * low_k,{}} ); + if( en == sf1.end() ) return false; + return en->first < val * high_k; +} + +static const circuit& circ( const step_function &sf1, long double val ) +{ + auto en = upper_bound(sf1.begin(), sf1.end(), pair{val * low_k,{}} ); + return en->second; +} + + +static bool compare_circuits( const circuit& a, const circuit& b ) +{ + matr ma(a); + matr mb(b); + ma.reduce(); + mb.reduce(); + + if( ma.h() != mb.h() ) + { + cout << "totally wrong: different sde|.|^2" << endl; + throw std::logic_error(""); + return false; + } + + bool found = false; + int glob_ph = 0; + for( int i = 0; i < 8; ++i ) + { + if ( ma.d[0][0] == mb.d[0][0] * mb.d[0][0].omega(i) ) + { + found = true; + glob_ph = i; + } + } + + if( !found ) + { + cout << "unexpected: different x" << endl; + throw std::logic_error(""); + return false; + } + + mb.mul_eq_w(glob_ph); + + if( mb.t() != ma.t() ) + { + cout << endl << "different number of t gates:" << ma.t() << " vs " << mb.t() << endl; + if(abs( mb.t() - ma.t() ) > 1 ) + throw logic_error("too big difference in the number of T gates"); + return false; + } + + return true; +} + +template< class T> +ostream& operator<< ( ostream& out, const vector& vec ) +{ + for( const auto& a : vec ) + out << a << " ,"; + return out; +} + + +template< class T> +ostream& operator<< ( ostream& out, const vector< pair >& vec ) +{ + for( const auto& a : vec ) + out << a.first << " ,"; + return out; +} + + +void DiffApplication::compare(const step_function &sf1, const step_function &sf2, const symbolic_angle &angle) +{ + + double dist_id = trace_dist(Rz(angle),matrix2x2hpr::Id()) * low_k; + cout << "comparing angle:" << angle << " " ; + pair vid {dist_id,{}}; + long double max1 = ( lower_bound( sf1.begin(), sf1.end(), vid ) - 1 )->first ; + long double max2 = ( lower_bound( sf2.begin(), sf2.end(), vid ) - 1 )->first ; + + long double minl = max( sf1.begin()->first, sf2.begin()->first ) * low_k; + long double maxl = min(max1,max2) * high_k ; + + auto bg = lower_bound(sf1.begin(), sf1.end(), pair{minl,{}} ); + auto en = upper_bound(sf1.begin(), sf1.end(), pair{maxl,{}} ); + + std::vector mfs; + std::vector diff_circs; + int records = 0; + for( auto i = bg; i < en; ++i ) + { + if( belongs(sf2,i->first) ) + { + records++; + if( !compare_circuits(i->second,circ(sf2,i->first)) ) + diff_circs.push_back(i->first); + } + else + { + mfs.push_back(i->first); + } + } + + bg = lower_bound(sf2.begin(), sf2.end(), pair{minl,{}} ); + en = upper_bound(sf2.begin(), sf2.end(), pair{maxl,{}} ); + + std::vector mff; + for( auto i = bg; i < en; ++i ) + { + if( !belongs(sf1,i->first) ) + { + mff.push_back(i->first); + } + } + + cout << " total:" << records << " "; + cout << "different in t gates: " << diff_circs.size() ; + cout << " missing: from 1st: " << mff.size() << " ,from 2nd: " << mfs.size() << endl; + + if( mff.size() > 0 ) + { + cout << " missed from 1st: " << mff << endl; + } + + if( mfs.size() > 0 ) + { + cout << " missed from 2nd: " << mfs << endl; + } + + + + } + + + void sqct_light_app::run() + { + ifstream ifs; + ifs.exceptions ( std::ifstream::badbit ); + ifs.open(par.in_filename); + vector angles; + vector precisions; + vector names; + while( ifs ) + { + string line; + getline(ifs,line); + if(line.size() == 0 ) + continue; + if( line[0] == '#' ) + continue; + + istringstream ss(line); + symbolic_angle a{0,1,true}; + //char delim; + double precision = 0.; + string name; + ss >> a.numerator >> a.denominator >> a.pi >> precision >> name; + angles.push_back(a); + precisions.push_back(precision); + names.push_back(name); + } + + cout << "Loading cahce ... "; + zrot_cache zc(par.cache_filename); + + cout << "Processing " << angles.size() << " requests" << endl; + + ofstream ofs; + ifs.exceptions ( std::ifstream::badbit ); + ofs.open(par.out_filename); + + + for( size_t i = 0; i < angles.size(); ++i ) + { + try + { + circuit c = zc.lookup(angles[i],precisions[i]); + } + catch(not_found_exception& e) + { + cout << "Not found:" << e.angle << ":" << e.precision << endl; + } + } + + //write footer + } diff --git a/appr/zrot_cache.h b/appr/zrot_cache.h new file mode 100644 index 0000000..88e0b44 --- /dev/null +++ b/appr/zrot_cache.h @@ -0,0 +1,86 @@ +// 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 . +// +// 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 ZROT_CACHE_H +#define ZROT_CACHE_H + +#include "toptzrot2.h" +#include "gatelibrary.h" +#include +#include +#include + +struct not_found_exception : public std::logic_error +{ + not_found_exception( const symbolic_angle& _angle, double _precision ) + : logic_error("approximation not found"), + angle(_angle), + precision(_precision) {}; + + symbolic_angle angle; + double precision; +}; + +typedef std::vector > step_function; + +struct zrot_cache : public std::map< symbolic_angle, step_function > +{ + zrot_cache(); + zrot_cache( const std::string& filename, bool recalculate = false ); + // loads cached data from file + void read_file( const std::string& filename ); + // looks up approximation + circuit lookup(const symbolic_angle& angle, double precision ) const throw ( not_found_exception ) ; + bool recalculate; +}; + +std::ostream& operator<< ( std::ostream& out , const zrot_cache& zc ); + +struct DiffApplicationParams +{ + std::vector< std::string > filenames; +}; + +struct DiffApplication +{ + DiffApplication( const DiffApplicationParams& in ) : params(in) {} + void run(); + void compare(const step_function& sf1, const step_function& sf2, const symbolic_angle& angle ); + + const DiffApplicationParams& params; +}; + +struct sqct_light_params +{ + std::string in_filename; + std::string out_filename; + std::string cache_filename; +}; + +struct sqct_light_app +{ + sqct_light_app( const sqct_light_params& in ) : + par(in) {} + void run(); + const sqct_light_params& par; +}; + +#endif // ZROT_CACHE_H diff --git a/es/exactdecomposer.cpp b/es/exactdecomposer.cpp index 35df64e..a900803 100644 --- a/es/exactdecomposer.cpp +++ b/es/exactdecomposer.cpp @@ -1,20 +1,20 @@ // 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 . -// +// #include "exactdecomposer.h" #include "matrix2x2.h" @@ -35,7 +35,7 @@ exactDecomposer::exactDecomposer() } } -void exactDecomposer::decompose( const matrix2x2& matr, circuit& c) +void exactDecomposer::decompose( const matrix2x2& matr, circuit& c) const { typedef ring_int< resring<8> > rr8; typedef matrix2x2< resring<8> > mrr8; @@ -82,3 +82,16 @@ void exactDecomposer::decompose( const matrix2x2& matr, circuit& c) slC.find(current,r); c.push_back(r); } + +circuit exactDecomposer::decompose(const matrix2x2 &matr) +{ + circuit tmp; + instance().decompose(matr,tmp); + return tmp; +} + +const exactDecomposer &exactDecomposer::instance() +{ + static exactDecomposer ed; + return ed; +} diff --git a/es/exactdecomposer.h b/es/exactdecomposer.h index 80fd60f..bc1f498 100644 --- a/es/exactdecomposer.h +++ b/es/exactdecomposer.h @@ -1,20 +1,20 @@ // 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 . -// +// #ifndef EXACTDECOMPOSER_H #define EXACTDECOMPOSER_H @@ -30,7 +30,9 @@ class exactDecomposer /// \brief Initializes generators exactDecomposer(); /// \brief Decomposes matr into circuit c - void decompose( const matrix2x2 &matr, circuit& c ); + void decompose( const matrix2x2 &matr, circuit& c ) const; + static circuit decompose( const matrix2x2 &matr ); + static const exactDecomposer& instance(); private: /// \brief Matrix over the ring type used during decomposition typedef matrix2x2 M; @@ -43,7 +45,7 @@ class exactDecomposer /// \brief Inverses of generators used for decomposition M generators_ctr[gen_count]; /// \brief Object of class to lookup T-optimal sequences for unitaries with \f$ sde(|\cdot|^2) < 4 \f$ - seqLookupCliff slC; + mutable seqLookupCliff slC; }; #endif // EXACTDECOMPOSER_H diff --git a/gatelibrary.cpp b/gatelibrary.cpp index 470272a..e0c7198 100644 --- a/gatelibrary.cpp +++ b/gatelibrary.cpp @@ -1,20 +1,20 @@ // 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 . -// +// #include "gatelibrary.h" #include "output.h" @@ -24,18 +24,18 @@ using namespace std; gateLibrary::gateLibrary() : - name(GLw7 + 1), name_qc(GLw7 + 1), - matrix_str( GLw7 + 1), matrix( GLw7 + 1), - inverse(GLw7 + 1), cost( GLw7 + 1,0 ) + name(GLw7 + 10), name_qc(GLw7 + 10), + matrix_str( GLw7 + 10), matrix( GLw7 + 10), + inverse(GLw7 + 10), cost( GLw7 + 10,0 ) { - name[Id] = "Id"; + name[Id] = ""; name[T] = "T"; name[P] = "P"; - name[TP] = "Td Z"; + name[TP] = "T*Z"; name[Z] = "Z"; - name[TZ] = "T Z"; - name[Pd] = "Pd"; - name[Td] = "Td"; + name[TZ] = "TZ"; + name[Pd] = "P*"; + name[Td] = "T*"; name[H] = "H"; name[X] = "X"; name[Y] = "Y"; @@ -63,11 +63,11 @@ gateLibrary::gateLibrary() : matrix_str[Id] = "Id"; matrix_str[T] = "T"; matrix_str[P] = "P"; - matrix_str[TP] = "Td.Z"; + matrix_str[TP] = "T*.Z"; matrix_str[Z] = "Z"; matrix_str[TZ] = "T.Z"; - matrix_str[Pd] = "Pd"; - matrix_str[Td] = "Td"; + matrix_str[Pd] = "P*"; + matrix_str[Td] = "T*"; matrix_str[H] = "H"; matrix_str[X] = "X"; matrix_str[Y] = "Y"; @@ -130,6 +130,9 @@ gateLibrary::gateLibrary() : symbols['H'] = H; symbols['h'] = H; + symbols['I'] = Id; + symbols['i'] = Id; + cost[Td] = cost[T] = 1000; cost[H] = 10; cost[P] = cost[Pd] = 40; @@ -170,7 +173,7 @@ void circuit::convert(circuit::m& res) const const circuit& c = *this; static const gateLibrary& gl = gateLibrary::instance(); res = gl.matrix[gl.Id]; - for( int i = 0; i < size(); ++i ) + for( size_t i = 0; i < size(); ++i ) { m tmp = gl.matrix[ c[i] ] * res; res = tmp; @@ -196,7 +199,7 @@ void circuit::toStream(ostream &out) const return; } - for( int i = 0; i < size(); ++i ) + for( size_t i = 0; i < size(); ++i ) out << gl.name_qc[ c[i] ] << endl; } @@ -211,8 +214,8 @@ void circuit::toStreamSym(ostream &out) const return; } - for( int i = 0; i < size() - 1; ++i ) - out << gl.name[ c[i] ] << " "; + for( size_t i = 0; i < size() - 1; ++i ) + out << gl.name[ c[i] ] ; out << gl.name[ c.back() ]; } @@ -254,7 +257,7 @@ void circuit::push_back( const circuit& c ) push_back(i); } -vector circuit::count() +vector circuit::count() const { static const gateLibrary& gl = gateLibrary::instance(); vector res( gl.GLw7 + 1,0 ); @@ -298,8 +301,20 @@ int circuit::cost() { static const gateLibrary& gl = gateLibrary::instance(); int sum = 0; - for( int i = 0; i < size(); ++i ) + for( size_t i = 0; i < size(); ++i ) sum += gl.cost[ at(i) ]; return sum; } + + +circuit::operator matrix2x2hpr() const +{ + return matrix2x2hpr(matrix2x2(*this) ); +} + +std::ostream& operator<< ( std::ostream& out , const circuit& c) +{ + c.toStreamSym(out); + return out; +} diff --git a/gatelibrary.h b/gatelibrary.h index d8f56ba..2260b09 100644 --- a/gatelibrary.h +++ b/gatelibrary.h @@ -1,20 +1,20 @@ // 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 . -// +// #ifndef GATELIBRARY_H #define GATELIBRARY_H @@ -88,6 +88,8 @@ class circuit : public std::deque typedef matrix2x2 m; /// \brief Conversion operator to unitary operator matrix2x2 () const; + /// \brief Conversion operator to unitary + operator matrix2x2hpr () const; /// \brief Converts circuit to unitary and writes result into res void convert( m& res ) const; /// \brief Outputs circuit to stream in Dot QC friendly way @@ -109,12 +111,17 @@ class circuit : public std::deque /// For resulting vector \b res, the number \b res[gateId] gives a number of gates with /// specific id in the circuit. /// \see gateLibrary class for the list of ids - std::vector count(); + std::vector count() const; /// \brief Reads circuit from stream in line by line mode. Lines beginning with # /// are skipped. All charachters that are not gate symbols or '*','d' are ignored. /// If gate symbol is followed by * or d it repplaced by its inverse. void fromStream( std::istream& in, bool reverse = false ); /// \brief Returns total cost of the circuit defined by gateLibrary::cost vector int cost(); + }; + +std::ostream& operator<< ( std::ostream& out , const circuit& c); + + #endif // GATELIBRARY_H diff --git a/hprhelpers.cpp b/hprhelpers.cpp index d449f73..854abec 100644 --- a/hprhelpers.cpp +++ b/hprhelpers.cpp @@ -1,23 +1,24 @@ // 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 . -// +// #define MPFR_REAL_DATA_PUBLIC #include "hprhelpers.h" +#include void hprHelpers::convert(const hprHelpers::hpr_complex &in, std::complex &out) { @@ -46,14 +47,14 @@ struct piData piData() { pi = std::string( "3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798" - "2148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786" - "78316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882" - "04665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733" - "62440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637" - "1787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317" - "328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083" - "8142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642" - "0198938095257201065485863279" ); + "2148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786" + "78316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882" + "04665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733" + "62440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637" + "1787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317" + "328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083" + "8142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642" + "0198938095257201065485863279" ); } }; @@ -92,3 +93,23 @@ const hprHelpers::hpr_real &hprHelpers::sqrt2ov2() static hpr_real m_s( sqrt(two()) / two() ); return m_s; } + +const hprHelpers::hpr_real &hprHelpers::sqrt2() +{ + static hpr_real m_s( sqrt(two()) ); + return m_s; +} + +hprr pow2(int n) +{ + mpz_class pow = 1; + pow <<= n; + hprr powf = 1; + mpfr_set_z(powf._x,pow.get_mpz_t(),MPFR_RNDN); //powf = power + return powf; +} + +long double to_ld( const hprr& a ) +{ + return mpfr_get_ld( a._x , MPFR_RNDN ); +} diff --git a/hprhelpers.h b/hprhelpers.h index 22db933..7ea5253 100644 --- a/hprhelpers.h +++ b/hprhelpers.h @@ -1,33 +1,35 @@ // 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 . -// +// #ifndef HPRHELPERS_H #define HPRHELPERS_H +#define MPFR_REAL_DATA_PUBLIC #include "real.hpp" #include +typedef mpfr::real<256> hprr; + /// \brief Helper functions and data for high precision arithmetic -class hprHelpers +struct hprHelpers { -public: + typedef mpfr::real<256> hpr_real; /// \brief High precision real type used in the project - typedef mpfr::real<4096> hpr_real; /// \brief High precision complex type used in the project typedef std::complex hpr_complex; /// \brief Transforms high precision complex number into machine complex @@ -51,6 +53,12 @@ class hprHelpers static const hpr_real& mhalf(); /// \brief High precision \f$ \frac{\sqrt{2}}{2} \f$ static const hpr_real& sqrt2ov2(); + /// \brief High precision \f$ \frac{\sqrt{2}}{2} \f$ + static const hpr_real& sqrt2(); }; +hprr pow2( int n ); +long double to_ld( const hprr& a ); + + #endif // HPRHELPERS_H diff --git a/mainA.cpp b/mainA.cpp new file mode 100644 index 0000000..1ed052d --- /dev/null +++ b/mainA.cpp @@ -0,0 +1,163 @@ +// Copyright (c) 2012 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 . +// + +#include "appr/toptzrot2.h" +#include "output.h" + +#include "appr/zrot_cache.h" + +#include +#include +#include +#include + +#include +#include +#include + +// command line parcing code based on http://www.boost.org/doc/libs/1_49_0/doc/html/program_options/tutorial.html#id2499896 +// and BOOST_ROOT/libs/program_options/example/first.cpp +namespace po = boost::program_options; +namespace btm = boost::timer; + +using namespace std; + +static bool print_help( const string& topic ) +{ + map hi; + hi["in"]; + if( hi.count(topic) ) + { + cout << hi[topic] << endl; + return true; + } + + return false; +} + +//////////////////////////////////////////////////////////////////// + +static void print_about_message() +{ +cout << "Copyright (c) 2013 Vadym Kliuchnikov sqct(dot)software(at)gmail(dot)com" << endl << endl; +cout << "SQCT is free software: you can redistribute it and/or modify" << endl; +cout << "it under the terms of the GNU Lesser General Public License as published by" << endl; +cout << "the Free Software Foundation, either version 3 of the License, or" << endl; +cout << "(at your option) any later version." << endl; +cout << "" << endl; +cout << "SQCT is distributed in the hope that it will be useful," << endl; +cout << "but WITHOUT ANY WARRANTY; without even the implied warranty of" << endl; +cout << "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the" << endl; +cout << "GNU Lesser General Public License for more details." << endl; +cout << "" << endl; +cout << "You should have received a copy of the GNU Lesser General Public License" << endl; +cout << "along with SQCT. If not, see ." << endl; +cout << "" << endl; +} + +//////////////////////////////////////////////////////////////////// + +int main(int ac, char* av[]) +{ + + string help_topic; + DiffApplicationParams diff_params; + topt_app_params topt{80,"in","rz"}; + sqct_light_params sp; + + try { + + po::options_description desc("Allowed options"); + desc.add_options() + + ("help,H", po::value< string >(&help_topic)->implicit_value(""), + "Produce help message, see help