Skip to content

Commit

Permalink
Merge pull request #279 from scality/fix/fft_add
Browse files Browse the repository at this point in the history 
Fix bugs in fft::Additive
  • Loading branch information
@vrancurel
vrancurel authored Apr 4, 2019
2 parents 6211b25 + 136d2a7 commit 621f3f5
Show file tree
Hide file tree
Showing 3 changed files with 291 additions and 244 deletions.
186 changes: 98 additions & 88 deletions src/fft_add.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,8 +53,10 @@ class Additive : public FourierTransform<T> {
using FourierTransform<T>::ifft;
using FourierTransform<T>::fft_inv;

Additive(const gf::Field<T>& gf, T m, vec::Vector<T>* betas = nullptr);
~Additive();
Additive(
const gf::Field<T>& gf,
T m,
std::shared_ptr<vec::Vector<T>> betas = nullptr);
void compute_basis();
void compute_beta_m_powers();
void compute_G();
Expand All @@ -63,7 +65,7 @@ class Additive : public FourierTransform<T> {
void fft(vec::Vector<T>& output, vec::Vector<T>& input) override;
void ifft(vec::Vector<T>& output, vec::Vector<T>& input) override;
void fft_inv(vec::Vector<T>& output, vec::Vector<T>& input) override;
void taylor_expand_t2(vec::Vector<T>& input, int n, bool do_copy = false);
void taylor_expand_t2(vec::Vector<T>& input, int n);
void
taylor_expand(vec::Vector<T>& output, vec::Vector<T>& input, int n, int t);
void inv_taylor_expand_t2(vec::Vector<T>& output);
Expand All @@ -83,21 +85,24 @@ class Additive : public FourierTransform<T> {
T m_k, deg0, deg1, deg2;
T beta_1, inv_beta_1;
T beta_m, inv_beta_m;
vec::Vector<T>* betas = nullptr;
vec::Vector<T>* gammas = nullptr;
vec::Vector<T>* deltas = nullptr;
vec::Vector<T>* beta_m_powers = nullptr;
vec::Vector<T>* G = nullptr;
vec::Vector<T>* g0 = nullptr;
vec::Vector<T>* g1 = nullptr;
vec::Vector<T>* u = nullptr;
vec::Vector<T>* v = nullptr;
vec::Vector<T>* mem = nullptr;
Additive<T>* fft_add = nullptr;
std::shared_ptr<vec::Vector<T>> betas = nullptr;
std::unique_ptr<vec::Vector<T>> gammas = nullptr;
std::shared_ptr<vec::Vector<T>> deltas = nullptr;
std::unique_ptr<vec::Vector<T>> beta_m_powers = nullptr;
std::unique_ptr<vec::Vector<T>> G = nullptr;
std::unique_ptr<vec::Vector<T>> g0 = nullptr;
std::unique_ptr<vec::Vector<T>> g1 = nullptr;
std::unique_ptr<vec::Vector<T>> u = nullptr;
std::unique_ptr<vec::Vector<T>> v = nullptr;
std::unique_ptr<vec::Vector<T>> mem = nullptr;
std::unique_ptr<Additive<T>> fft_add = nullptr;
};

template <typename T>
Additive<T>::Additive(const gf::Field<T>& gf, T m, vec::Vector<T>* betas)
Additive<T>::Additive(
const gf::Field<T>& gf,
T m,
std::shared_ptr<vec::Vector<T>> betas)
: FourierTransform<T>(gf, arith::exp2<T>(m), true)
{
assert(m >= 1);
Expand All @@ -108,7 +113,7 @@ Additive<T>::Additive(const gf::Field<T>& gf, T m, vec::Vector<T>* betas)
assert(gf.get_sub_field().card() == 2);

create_betas = true;
betas = new vec::Vector<T>(gf, m);
betas = std::make_shared<vec::Vector<T>>(gf, m);
T beta = this->gf->get_primitive_root();
betas->set(0, beta);
for (T i = 1; i < m - 1; i++) {
Expand All @@ -122,54 +127,28 @@ Additive<T>::Additive(const gf::Field<T>& gf, T m, vec::Vector<T>* betas)
this->beta_m = betas->get(m - 1);
this->inv_beta_m = this->gf->inv(this->beta_m);

this->beta_m_powers = new vec::Vector<T>(gf, this->n);
this->beta_m_powers = std::make_unique<vec::Vector<T>>(gf, this->n);
this->compute_beta_m_powers();
if (m > 1) {
this->m_k = arith::exp2<T>(m - 1);

this->m_k = arith::exp2<T>(m - 1);
this->g0 = std::make_unique<vec::Vector<T>>(gf, this->m_k);
this->g1 = std::make_unique<vec::Vector<T>>(gf, this->m_k);

if (m > 1) {
// compute gammas and deltas
this->gammas = new vec::Vector<T>(gf, m - 1);
this->deltas = new vec::Vector<T>(gf, m - 1);
this->G = new vec::Vector<T>(gf, this->m_k);
this->gammas = std::make_unique<vec::Vector<T>>(gf, m - 1);
this->deltas = std::make_shared<vec::Vector<T>>(gf, m - 1);
this->G = std::make_unique<vec::Vector<T>>(gf, this->m_k);
this->compute_basis();

this->u = new vec::Vector<T>(gf, this->m_k);
this->v = new vec::Vector<T>(gf, this->m_k);
this->g0 = new vec::Vector<T>(gf, this->m_k);
this->g1 = new vec::Vector<T>(gf, this->m_k);
this->u = std::make_unique<vec::Vector<T>>(gf, this->m_k);
this->v = std::make_unique<vec::Vector<T>>(gf, this->m_k);

this->mem = new vec::Vector<T>(gf, this->n);
this->fft_add = new Additive(gf, m - 1, this->deltas);
this->mem = std::make_unique<vec::Vector<T>>(gf, this->n);
this->fft_add = std::make_unique<Additive>(gf, m - 1, this->deltas);
}
}

template <typename T>
Additive<T>::~Additive()
{
if (create_betas && betas)
delete betas;
if (gammas)
delete gammas;
if (deltas)
delete deltas;
if (beta_m_powers)
delete beta_m_powers;
if (G)
delete G;
if (u)
delete u;
if (v)
delete v;
if (g0)
delete g0;
if (g1)
delete g1;
if (mem)
delete mem;
if (fft_add)
delete fft_add;
}

template <typename T>
void Additive<T>::compute_beta_m_powers()
{
Expand Down Expand Up @@ -261,8 +240,9 @@ template <typename T>
void Additive<T>::_fft(vec::Vector<T>& output, vec::Vector<T>& input)
{
mem->copy(&input, this->n);
if (beta_m > 1)
mem->hadamard_mul(beta_m_powers);
if (beta_m > 1) {
mem->hadamard_mul(beta_m_powers.get());
}

// compute taylor expansion of g(x) at (x^2 - x)
// outputs are this->g0 and this->g1
Expand All @@ -272,15 +252,15 @@ void Additive<T>::_fft(vec::Vector<T>& output, vec::Vector<T>& input)
this->fft_add->fft(*v, *g1);

// copy output = (undefined, v)
output.copy(v, m_k, m_k);
output.copy(v.get(), m_k, m_k);
// perform G * v hadamard multiplication
v->hadamard_mul(G);
v->hadamard_mul(G.get());
// v += u
v->add(u);
v->add(u.get());
// output = (u + G*v, v)
output.copy(v, m_k);
output.copy(v.get(), m_k);
// output = (u + G*v, (u + G*v) + v)
output.add(v, m_k);
output.add(v.get(), m_k);
}

template <typename T>
Expand All @@ -297,9 +277,9 @@ void Additive<T>::_ifft(vec::Vector<T>& output, vec::Vector<T>& input)
// v = w0 + w1
v->add(&w0);
// u = v
u->copy(v);
u->copy(v.get());
// u = G * v
u->hadamard_mul(G);
u->hadamard_mul(G.get());
// u = w0 + G * v;
u->add(&w0);

Expand Down Expand Up @@ -418,37 +398,37 @@ void Additive<T>::ifft(vec::Vector<T>& output, vec::Vector<T>& input)
*
* This function is used for FFT over n=2^m, hence n must be power of 2
*
* @param input input vector that would be modified by the function
* @param n a power of 2
* @param do_copy flag to do a copy of input or not since this function will
* modify input vector
*/
template <typename T>
void Additive<T>::taylor_expand_t2(vec::Vector<T>& input, int n, bool do_copy)
void Additive<T>::taylor_expand_t2(vec::Vector<T>& input, int n)
{
assert(g0 != nullptr);
assert(n >= 1);
assert(input.get_n() <= n);
assert(arith::is_power_of_2<int>(n));

if (n <= 2) {
g0->set(0, input.get(0));

if (n == 2) {
g1->set(0, input.get(1));
}

return;
}

// find k s.t. t2^k < n <= 2 *t2^k
int k = find_k(n, 2);

vec::Vector<T>* _input;
if (do_copy) {
_input = new vec::Vector<T>(*(this->gf), input.get_n());
_input->copy(&input, input.get_n());
} else
_input = &input;

_taylor_expand_t2(*_input, n, k, 0);
_taylor_expand_t2(input, n, k, 0);

// get g0, g1 from mem
for (int i = 0; i < this->n; i += 2) {
g0->set(i / 2, _input->get(i));
g1->set(i / 2, _input->get(i + 1));
for (int i = 0; i < n; i += 2) {
g0->set(i / 2, input.get(i));
g1->set(i / 2, input.get(i + 1));
}

if (do_copy)
delete _input;
}

/**
Expand Down Expand Up @@ -508,9 +488,23 @@ void Additive<T>::inv_taylor_expand_t2(vec::Vector<T>& output)
{
assert(g0 != nullptr);
assert(g0->get_n() == g1->get_n());

const int o_len = output.get_n();
assert(arith::is_power_of_2<int>(o_len));

if (o_len <= 2) {
output.set(0, g0->get(0));

if (o_len == 2) {
output.set(1, g1->get(0));
}

return;
}

output.zero_fill();

int i = output.get_n() / 2 - 1;
int i = o_len / 2 - 1;
output.set(0, g0->get(i));
output.set(1, g1->get(i));
while (--i >= 0) {
Expand Down Expand Up @@ -548,6 +542,8 @@ void Additive<T>::mul_xt_x(vec::Vector<T>& vec, int t)
template <typename T>
inline int Additive<T>::find_k(int n, int t)
{
assert(t < n);

int k;
int t2k = t; // init for t*2^k with k = 0
for (k = 0; k < n; k++) {
Expand Down Expand Up @@ -586,7 +582,10 @@ void Additive<T>::taylor_expand(

// set output as input
output.copy(&input);
_taylor_expand(output, n, t);

if (n > t) {
_taylor_expand(output, n, t);
}
}

template <typename T>
Expand Down Expand Up @@ -640,18 +639,29 @@ void Additive<T>::inv_taylor_expand(
vec::Vector<T>& input,
int t)
{
output.zero_fill();
const int i_len = input.get_n();
const int o_len = output.get_n();

int i = input.get_n() - t;
assert(o_len <= this->n);
assert(i_len >= o_len);

// A temporary vector storing result
vec::Vector<T> tmp(input.get_gf(), i_len);
tmp.zero_fill();

int i = i_len - t;
vec::Slice<T> hi(&input, t, i);
output.add_mutual(&hi);
tmp.add_mutual(&hi);
while (i >= t) {
i -= t;
// multiply output to (x^t - x)
mul_xt_x(output, t);
mul_xt_x(tmp, t);
hi.set_map(i);
output.add_mutual(&hi);
tmp.add_mutual(&hi);
}

// copy only first `n` coefficients to output
output.copy(&tmp, o_len);
}

} // namespace fft
Expand Down
2 changes: 1 addition & 1 deletion src/vec_vector.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -167,7 +167,7 @@ Vector<T>::Vector(Vector const& other)
{
if (new_mem) {
this->mem = this->allocator.allocate(other.mem_len);
std::copy_n(this->mem, other.mem_len, other.mem);
std::copy_n(other.mem, other.mem_len, this->mem);
} else {
this->mem = other.mem;
}
Expand Down
Loading

0 comments on commit 621f3f5

Please sign in to comment.