https://github.com/epiqc/ScaffCC
Tip revision: fb0341d7eb6d3ae899a20f913e9f550f738d1bea authored by ah744 on 22 December 2016, 07:02:43 UTC
afree patch
afree patch
Tip revision: fb0341d
netgenerator.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 "netgenerator.h"
#include "numbersgen.h"
#include "epsilonnet.h"
#include <memory>
#include <deque>
#include <vector>
#include <algorithm>
#include <sstream>
#include <functional>
#include <cassert>
using namespace std;
void netGenerator::generateInitial()
{
typedef ring_int<int> ri;
unique_ptr< numbersGenerator<10> > png( new numbersGenerator<10>);
auto& ng = *png.get();
ng.generate_all_numbers();
epsilonnet enets[ ng.max_sde + 1 ];
for( int i = 0; i < ng.ip_range; ++i )
for( int j = 0; j < ng.max_val2 + 1; ++j )
{
const auto& vec = ng.numbers(i,j);
const auto& vec_compl = ng.numbers( ng.max_val2 - i, -j);
if( (!vec.empty() ) && (!vec_compl.empty() ))
{
int s = std::min ( ng.getSde(i,j), ng.getSde(ng.max_val2 - i, -j) );
int de = ng.max_sde / 2 - (s + 1 ) / 2;
vector<ri> v;
vector<ri> vc;
// normalize denominator and pick canonical representatives
for( auto ii : vec )
{
v.push_back( ii.canonical() );
v.back().div_eq_sqrt2( de );
}
for( auto ii : vec_compl )
{
vc.push_back( ii.canonical() );
vc.back().div_eq_sqrt2( de );
}
// sort to leave only unique elements
sort( v.begin(), v.end() );
sort( vc.begin(), vc.end() );
nodeRanges nr = {v.begin(),v.end(),vc.begin(), vc.end() };
assert( ( (i >> de) << de ) == i );
assert( ( (j >> de) << de ) == j );
// in the case when ipQxx = 0 we only add numbers such that ipxx <= ng.max_val2 / 2
// the number with ipxx >= ng.max_val2 / 2 will be in complementary part automatically
if( ( j != 0 ) || ( (j == 0) && (i <= ng.max_val2 / 2) ) )
enets[s].addNode(i >> de ,j >> de ,nr);
}
}
for( int i = 0; i <= ng.max_sde; ++i )
{
enets[i].saveToFile( fileName(i).c_str() );
}
}
string netGenerator::fileName(int sde)
{
stringstream s;
char const *folder = getenv("TMPDIR");
if (folder == 0)
folder = "/tmp";
s << folder << "/epsilon-net." << sde << ".bin";
return s.str();
}
///////////////// internal implementation of generator ///////////////////////////////////
typedef netGenerator::ri ri;
/// \brief Node of generated data
struct gdata
{
/// \brief ip.first is P(x) of num and ip.second is Q(x) of num \see ring_int::ipxx() for defintions of P(x) and Q(x)
pair<ri::pr_type,ri::pr_type> ip;
/// \brief Found number
ri num;
/// \brief Found complementary number \see enetNode::compl_offset
ri num_compl;
};
/// \brief Order on pointers to gdata based on values of P(x) and Q(x) \see ring_int::ipxx() for defintions of P(x) and Q(x)
struct gdata_comparator
{
bool operator () ( const gdata* a, const gdata* b) const
{
return a->ip < b->ip;
}
};
/// \brief Order on pointers based on data that they point to
template< class T >
struct ptr_comp : public binary_function< const T*, const T*, bool >
{
const bool operator() ( const T* a, const T* b)
{
return *a < *b;
}
};
/// \brief Internal structure for epsilon net generation
struct generation_data
{
/// \brief Sets sde end denominator exponent used during generation
generation_data( int sde, int de2 )
{
isGde1 = (sde % 2 == 1);
pow2n = 1 << (isGde1 ? de2 : de2 + 1 );
pow2nm1 = pow2n >> 1;
}
/// \brief Found values of P(x) and Q(x)
deque< gdata > ips;
/// \brief Sum of P(x) and P(y) for two complementary numbers \see epsilonnet::compl_offset
ri::pr_type pow2n;
/// \brief If P(x) is less than pow2nm1 we store (x,y), and (y,x) otherwise
ri::pr_type pow2nm1;
/// \brief Pointers to computed numbers
vector< const gdata* > order;
/// \brief When sde is even then \f$gde(|\cdot|^2)\f$ 0 and 1 otherwise. See properties of gde and sde in
/// http://arxiv.org/abs/1206.5236
bool isGde1;
/// \brief Adds all numbers to epsilon net that we can get from unit column based on (omega_k_y,x)
void add_for_all_k ( ri& omega_k_y, const ri& x )
{
for( int k = 0; k < 8; ++k )
{
omega_k_y.mul_eq_w(); //goes to next element in eq. class
ri z = x + omega_k_y;
ri w = x - omega_k_y;
auto ipxx = z.ipxx();
auto ipQxx = z.ipQxx();
if( isGde1 )
{
// in this case \f$gde(|\cdot|^2)=2\f$ and we need further reduction
if( ri::isGde2(ipxx,ipQxx) )
{
ipxx >>= 1;
ipQxx >>= 1;
z.div_eq_sqrt2();
w.div_eq_sqrt2();
push_back(ipxx,ipQxx,z,w);
}
}
else
{
if( ri::isGde1(ipxx,ipQxx) )
push_back(ipxx,ipQxx,z,w);
}
}
}
/// \brief Adds number equivalence class to epsilon net. By calling ring_int::canonical
/// we mod out multiplication by power of \f$ \omega \f$ and conjugation.
/// By comparison of ipxx with pow2nm1 we mod out permutation of column entries
void push_back( ri::pr_type ipxx, ri::pr_type ipQxx, const ri& a, const ri& ac )
{
if( ipQxx > 0 )
ips.push_back( { {ipxx, ipQxx}, a.canonical(),ac.canonical()} );
else if( ipQxx == 0 )
{
if( ipxx < pow2nm1 )
ips.push_back( {{ipxx, 0},a.canonical(),ac.canonical()} );
else
ips.push_back( {{pow2n - ipxx,0},ac.canonical(),a.canonical()} );
} else if ( ipQxx < 0 )
ips.push_back( {{pow2n - ipxx,-ipQxx},ac.canonical(),a.canonical()} );
}
/// \brief Builds epsilon net from computed data
epsilonnet* build_net()
{
epsilonnet* enet = new epsilonnet;
order.resize( ips.size() );
int sz = order.size();
for( int i = 0; i < sz; ++i )
order[i] = & ips[i];
gdata_comparator gcomp;
/// \todo This sort can be avoided with a bit more clever algorithm
sort( order.begin(), order.end(), gcomp );
int j = 0;
int j_end = 0;
for( int i = 0; i < sz - 1; ++i )
{
if( gcomp( order[i], order[i+1] ) )
{
j_end = i + 1;
add_to_net(j,j_end,enet);
j = j_end;
}
}
j_end = sz;
add_to_net(j,j_end,enet);
return enet;
}
/// \brief Adds range from \b order to epsilon net
void add_to_net( int j, int j_end ,epsilonnet* enet )
{
vector< const ri* > nums,nums_compl;
for( int i = j ; i < j_end; ++i )
{
nums.push_back( &(order[i]->num) );
nums_compl.push_back( &( order[i]->num_compl ) );
}
ptr_comp<ri> ri_comp;
sort( nums.begin(), nums.end(), ri_comp );
sort( nums_compl.begin(), nums_compl.end(), ri_comp );
nodeRangesPtr nr = { nums.begin(), nums.end(),
nums_compl.begin(), nums_compl.end() };
enet->addNode( order[j]->ip, nr);
}
};
epsilonnet* netGenerator::generate(const epsilonnet& enet)
{
generation_data gdata( enet.sde(), enet.denominatorExponent2() );
for( int i = 0; i < enet.nodes.size() - 1; ++i )
{
const enetNode& node = enet.nodes[i];
nodeRanges r;
enet.getNode( i, r );
for( auto i1 = r.nums_begin; i1 != r.nums_end; ++i1 )
{
const ri& x = *i1;
ri x_c = x.conjugate();
for( auto i2 = r.nums_compl_begin; i2 != r.nums_compl_end; ++i2 )
{
// addes all vectors that can be generated from equivalence class
ri omega_k_y(*i2);
gdata.add_for_all_k( omega_k_y, x );
ri omega_k_y_c( i2->conjugate() );
gdata.add_for_all_k( omega_k_y_c, x );
}
}
}
return gdata.build_net();
}