fft.cpp

/**
 * Copyright (c) 2016-2018 Koninklijke Philips N.V. All rights reserved. A
 * copyright license for redistribution and use in source and binary forms,
 * with or without modification, is hereby granted for non-commercial,
 * experimental and research purposes, provided that the following conditions
 * are met:
 * - Redistributions of source code must retain the above copyright notice,
 *   this list of conditions and the following disclaimers.
 * - Redistributions in binary form must reproduce the above copyright notice,
 *   this list of conditions and the following disclaimers in the
 *   documentation and/or other materials provided with the distribution. If
 *   you wish to use this software commercially, kindly contact
 *   info.licensing@philips.com to obtain a commercial license.
 *
 * This license extends only to copyright and does not include or grant any
 * patent license or other license whatsoever.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

#include <gmpxx.h>
#define MIE_ATE_USE_GMP
#include <bn.h>

using namespace bn;
using namespace std;

#include "fft.h"
#include "modp.h"

template<typename T> void fft(
  const vector<T>& roots,
  vector<T>& X, int X0,
  const vector<T>& x, int N, int x0, int s, int rs) {

  int Ndiv2 = N/2;
  int NR = roots.size();

  if (N==1) { X[X0] = x[x0]; return; }

  fft(roots, X, X0, x, Ndiv2, x0, 2*s, 2*rs);
  fft(roots, X, X0+Ndiv2, x, Ndiv2, x0+s, 2*s, 2*rs);

  for (int k = 0; k < Ndiv2; k++) {
    int rix = (((k*rs) % NR) + NR) % NR;
    T t = X[X0+k], delta = roots[rix]*X[X0+Ndiv2+k];
    X[X0+k] = t + delta;
    X[X0+Ndiv2+k] = t - delta;
  }
}

template<typename T> void fft(const vector<T>& roots, vector<T>& ret, const vector<T>& coeff) {
  fft(roots, ret, 0, coeff, coeff.size(), 0, 1, 1);
}

template<typename T> void fftinv(const vector<T>& roots, vector<T>& ret, const vector<T>& coeff) {
  fft(roots, ret, 0, coeff, coeff.size(), 0, 1, -1);

  T n(ret.size()), ninv = n.inv();
  for (unsigned int i = 0; i < ret.size(); i++)
    ret[i] = ninv*ret[i];
}

template void fft<modp>(const vector<modp>& roots, std::vector<modp>& ret, const std::vector<modp>& coeff);
template void fftinv<modp>(const vector<modp>& roots, std::vector<modp>& ret, const std::vector<modp>& coeff);