sleef-sys 0.1.2

//          Copyright Naoki Shibata 2010 - 2018.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE.txt or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <time.h>

#if defined(POWER64_UNDEF_USE_EXTERN_INLINES)
// This is a workaround required to cross compile for PPC64 binaries
#include <features.h>
#ifdef __USE_EXTERN_INLINES
#undef __USE_EXTERN_INLINES
#endif
#endif

#include <math.h>
#include <complex.h>

#include "sleef.h"
#include "sleefdft.h"

#include <fftw3.h>

#ifndef MODE
#define MODE SLEEF_MODE_DEBUG
#endif

#if BASETYPEID == 1
#define THRES 1e-30
#define SleefDFT_init1d SleefDFT_double_init1d
#define SleefDFT_execute SleefDFT_double_execute
typedef double real;
#elif BASETYPEID == 2
#define THRES 1e-13
#define SleefDFT_init1d SleefDFT_float_init1d
#define SleefDFT_execute SleefDFT_float_execute
typedef float real;
#else
#error BASETYPEID not set
#endif

static double squ(double x) { return x * x; }

// complex forward
double check_cf(int n) {
  fftw_complex *in  = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
  fftw_complex *out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
  fftw_plan w = fftw_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

  real *sx = (real *)Sleef_malloc(n*2*sizeof(real));
  real *sy = (real *)Sleef_malloc(n*2*sizeof(real));
  struct SleefDFT *p = SleefDFT_init1d(n, sx, sy, MODE);

  for(int i=0;i<n;i++) {
    real re = (2.0 * random() - 1) / (real)RAND_MAX;
    real im = (2.0 * random() - 1) / (real)RAND_MAX;
    sx[(i*2+0)] = re;
    sx[(i*2+1)] = im;
    in[i] = re + im * _Complex_I;
  }

  SleefDFT_execute(p, NULL, NULL);
  fftw_execute(w);

  double rmsn = 0, rmsd = 0;

  for(int i=0;i<n;i++) {
    rmsn += squ(sy[i*2+0] - creal(out[i])) + squ(sy[i*2+1] - cimag(out[i]));
    rmsd += squ(            creal(out[i])) + squ(            cimag(out[i]));
  }

  fftw_destroy_plan(w);
  fftw_free(in);
  fftw_free(out);

  Sleef_free(sx);
  Sleef_free(sy);
  SleefDFT_dispose(p);

  return rmsn / rmsd;
}

// complex backward
double check_cb(int n) {
  fftw_complex *in  = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
  fftw_complex *out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
  fftw_plan w = fftw_plan_dft_1d(n, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);

  real *sx = (real *)Sleef_malloc(n*2*sizeof(real));
  real *sy = (real *)Sleef_malloc(n*2*sizeof(real));
  struct SleefDFT *p = SleefDFT_init1d(n, sx, sy, SLEEF_MODE_BACKWARD | MODE);

  for(int i=0;i<n;i++) {
    real re = (2.0 * random() - 1) / (real)RAND_MAX;
    real im = (2.0 * random() - 1) / (real)RAND_MAX;
    sx[(i*2+0)] = re;
    sx[(i*2+1)] = im;
    in[i] = re + im * _Complex_I;
  }

  SleefDFT_execute(p, NULL, NULL);
  fftw_execute(w);

  double rmsn = 0, rmsd = 0;

  for(int i=0;i<n;i++) {
    rmsn += squ(sy[i*2+0] - creal(out[i])) + squ(sy[i*2+1] - cimag(out[i]));
    rmsd += squ(            creal(out[i])) + squ(            cimag(out[i]));
  }

  fftw_destroy_plan(w);
  fftw_free(in);
  fftw_free(out);

  Sleef_free(sx);
  Sleef_free(sy);
  SleefDFT_dispose(p);

  return rmsn / rmsd;
}

// real forward
double check_rf(int n) {
  double       *in  = (double *)      fftw_malloc(sizeof(double) * n);
  fftw_complex *out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (n/2+1));
  fftw_plan w       = fftw_plan_dft_r2c_1d(n, in, out, FFTW_ESTIMATE);

  real *sx = (real *)Sleef_malloc(n*sizeof(real));
  real *sy = (real *)Sleef_malloc((n/2+1)*sizeof(real)*2);
  struct SleefDFT *p = SleefDFT_init1d(n, sx, sy, SLEEF_MODE_REAL | MODE);

  for(int i=0;i<n;i++) {
    real re = (2.0 * random() - 1) / (real)RAND_MAX;
    sx[i] = re;
    in[i] = re;
  }

  SleefDFT_execute(p, NULL, NULL);
  fftw_execute(w);

  double rmsn = 0, rmsd = 0;

  for(int i=0;i<n/2+1;i++) {
    rmsn += squ(sy[i*2+0] - creal(out[i])) + squ(sy[i*2+1] - cimag(out[i]));
    rmsd += squ(            creal(out[i])) + squ(            cimag(out[i]));
  }

  fftw_destroy_plan(w);
  fftw_free(in);
  fftw_free(out);

  Sleef_free(sx);
  Sleef_free(sy);
  SleefDFT_dispose(p);

  return rmsn / rmsd;
}

// real backward
double check_rb(int n) {
  fftw_complex *in  = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (n/2+1));
  double       *out = (double *)      fftw_malloc(sizeof(double) * n);
  fftw_plan w = fftw_plan_dft_c2r_1d(n, in, out, FFTW_ESTIMATE);

  real *sx = (real *)Sleef_malloc((n/2+1) * sizeof(real)*2);
  real *sy = (real *)Sleef_malloc(sizeof(real)*n);
  struct SleefDFT *p = SleefDFT_init1d(n, sx, sy, SLEEF_MODE_REAL | SLEEF_MODE_BACKWARD | MODE);

  for(int i=0;i<n/2;i++) {
    if (i == 0) {
      in[0  ] = (2.0 * (rand() / (real)RAND_MAX) - 1);
      in[n/2] = (2.0 * (rand() / (real)RAND_MAX) - 1);
    } else {
      in[i  ] = (2.0 * (rand() / (real)RAND_MAX) - 1) + (2.0 * (rand() / (real)RAND_MAX) - 1) * _Complex_I;
    }
  }

  for(int i=0;i<n/2+1;i++) {
    sx[2*i+0] = creal(in[i]);
    sx[2*i+1] = cimag(in[i]);
  }
  
  SleefDFT_execute(p, NULL, NULL);
  fftw_execute(w);

  double rmsn = 0, rmsd = 0;

  for(int i=0;i<n;i++) {
    rmsn += squ(sy[i] - out[i]);
    rmsd += squ(        out[i]);
  }
  
  fftw_destroy_plan(w);
  fftw_free(in);
  fftw_free(out);

  Sleef_free(sx);
  Sleef_free(sy);
  SleefDFT_dispose(p);

  return rmsn / rmsd;
}

int main(int argc, char **argv) {
  if (argc != 2) {
    fprintf(stderr, "%s <log2n>\n", argv[0]);
    exit(-1);
  }

  const int n = 1 << atoi(argv[1]);

  srand(time(NULL));

  SleefDFT_setPlanFilePath(NULL, NULL, SLEEF_PLAN_RESET | SLEEF_PLAN_READONLY);

  //

  int success = 1;
  double e;

  e = check_cf(n);
  success = success && e < THRES;
  printf("complex forward   : %s (%g)\n", e < THRES ? "OK" : "NG", e);
  e = check_cb(n);
  success = success && e < THRES;
  printf("complex backward  : %s (%g)\n", e < THRES ? "OK" : "NG", e);
  e = check_rf(n);
  success = success && e < THRES;
  printf("real forward      : %s (%g)\n", e < THRES ? "OK" : "NG", e);
  e = check_rb(n);
  success = success && e < THRES;
  printf("real backward     : %s (%g)\n", e < THRES ? "OK" : "NG", e);

  exit(success ? 0 : -1);
}