C++ c_re函数代码示例

/* ---------------------------------------------------------------------------- */
COMPLEX
Cadd(COMPLEX x, COMPLEX y) {
  COMPLEX z;
  c_re(z) = c_re(x) + c_re(y);
  c_im(z) = c_im(x) + c_im(y);
  return z;
}

/* ---------------------------------------------------------------------------- */
COMPLEX
Cscl(COMPLEX x, REAL a) {
  COMPLEX z;
  c_re(z) = c_re(x) * a;
  c_im(z) = c_im(x) * a;
  return z;
}

/* ---------------------------------------------------------------------------- */
COMPLEX
Csub(COMPLEX x, COMPLEX y) {
  COMPLEX z;
  c_re(z) = c_re(x) - c_re(y);
  c_im(z) = c_im(x) - c_im(y);
  return z;
}

/* ---------------------------------------------------------------------------- */
COMPLEX
Cmul(COMPLEX x, COMPLEX y) {
  COMPLEX z;
  c_re(z) = c_re(x) * c_re(y) - c_im(x) * c_im(y);
  c_im(z) = c_im(x) * c_re(y) + c_re(x) * c_im(y);
  return z;
}

void testnd_out_of_place(int rank, int *n, fftw_direction dir,
			 fftwnd_plan validated_plan)
{
     int istride, ostride;
     int N, dim, i;
     fftw_complex *in1, *in2, *out1, *out2;
     fftwnd_plan p;
     int flags = measure_flag | wisdom_flag;

     if (coinflip())
	  flags |= FFTW_THREADSAFE;

     N = 1;
     for (dim = 0; dim < rank; ++dim)
	  N *= n[dim];

     in1 = (fftw_complex *) fftw_malloc(N * MAX_STRIDE * sizeof(fftw_complex));
     out1 = (fftw_complex *) fftw_malloc(N * MAX_STRIDE * sizeof(fftw_complex));
     in2 = (fftw_complex *) fftw_malloc(N * sizeof(fftw_complex));
     out2 = (fftw_complex *) fftw_malloc(N * sizeof(fftw_complex));

     p = fftwnd_create_plan(rank, n, dir, flags);

     for (istride = 1; istride <= MAX_STRIDE; ++istride) {
	  /* generate random inputs */
	  for (i = 0; i < N; ++i) {
	       int j;
	       c_re(in2[i]) = DRAND();
	       c_im(in2[i]) = DRAND();
	       for (j = 0; j < istride; ++j) {
		    c_re(in1[i * istride + j]) = c_re(in2[i]);
		    c_im(in1[i * istride + j]) = c_im(in2[i]);
	       }
	  }

	  for (ostride = 1; ostride <= MAX_STRIDE; ++ostride) {
	       int howmany = (istride < ostride) ? istride : ostride;

	       if (howmany != 1 || istride != 1 || ostride != 1 || coinflip())
		    fftwnd_threads(nthreads, p, howmany, in1, istride, 1,
				   out1, ostride, 1);
	       else
		    fftwnd_threads_one(nthreads, p, in1, out1);

	       fftwnd(validated_plan, 1, in2, 1, 1, out2, 1, 1);

	       for (i = 0; i < howmany; ++i)
		    CHECK(compute_error_complex(out1 + i, ostride, out2, 1, N)
			  < TOLERANCE,
			  "testnd_out_of_place: wrong answer");
	  }
     }

     fftwnd_destroy_plan(p);

     fftw_free(out2);
     fftw_free(in2);
     fftw_free(out1);
     fftw_free(in1);
}

double linear(dofft_closure *k, int realp,
	      int n, C *inA, C *inB, C *inC, C *outA,
	      C *outB, C *outC, C *tmp, int rounds, double tol)
{
     int j;
     double e = 0.0;

     for (j = 0; j < rounds; ++j) {
	  C alpha, beta;
	  c_re(alpha) = mydrand();
	  c_im(alpha) = realp ? 0.0 : mydrand();
	  c_re(beta) = mydrand();
	  c_im(beta) = realp ? 0.0 : mydrand();
	  arand(inA, n);
	  arand(inB, n);
	  k->apply(k, inA, outA);
	  k->apply(k, inB, outB);

	  ascale(outA, alpha, n);
	  ascale(outB, beta, n);
	  aadd(tmp, outA, outB, n);
	  ascale(inA, alpha, n);
	  ascale(inB, beta, n);
	  aadd(inC, inA, inB, n);
	  k->apply(k, inC, outC);

	  e = dmax(e, acmp(outC, tmp, n, "linear", tol));
     }
     return e;
}

void realrft (
  double *in,
  unsigned n,
  double *out
)
{
	COMPLEX *c_in, *c_out;
	unsigned i;

	if (n == 0 ||
	    (c_in = (COMPLEX *) malloc (n * sizeof (COMPLEX))) == 0 ||
	    (c_out = (COMPLEX *) malloc (n * sizeof (COMPLEX))) == 0)
		return;
	
	c_re (c_in [0]) = in [0];		/* dc */
	c_im (c_in [0]) = 0;
	for (i = 1; i < (n + 1) / 2; i++) {	/* geconj. symm. harmonischen */
		c_re (c_in [i]) = in [2 * i - 1] / 2;
		c_im (c_in [i]) = in [2 * i] / -2;
		c_re (c_in [n - i]) = in [2 * i - 1] / 2;
		c_im (c_in [n - i]) = in [2 * i] / 2;
	}
	if (n % 2 == 0) {			/* Nyquist */
		c_re (c_in [n / 2]) = in [n - 1];
		c_im (c_in [n / 2]) = 0;
	}

	rft (c_in, n, c_out);

	for (i = 0; i < n; i++)
		out [i] = c_re (c_out [i]);

	free ((char *) c_in);
	free ((char *) c_out);
}

void fftwi_twiddle_2(fftw_complex *A, const fftw_complex *W, int iostride,
		     int m, int dist)
{
     int i;
     fftw_complex *inout;
     inout = A;
     for (i = m; i > 0; i = i - 1, inout = inout + dist, W = W + 1) {
	  fftw_real tmp1;
	  fftw_real tmp8;
	  fftw_real tmp6;
	  fftw_real tmp7;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp1 = c_re(inout[0]);
	  tmp8 = c_im(inout[0]);
	  {
	       fftw_real tmp3;
	       fftw_real tmp5;
	       fftw_real tmp2;
	       fftw_real tmp4;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp3 = c_re(inout[iostride]);
	       tmp5 = c_im(inout[iostride]);
	       tmp2 = c_re(W[0]);
	       tmp4 = c_im(W[0]);
	       tmp6 = (tmp2 * tmp3) + (tmp4 * tmp5);
	       tmp7 = (tmp2 * tmp5) - (tmp4 * tmp3);
	  }
	  c_re(inout[iostride]) = tmp1 - tmp6;
	  c_re(inout[0]) = tmp1 + tmp6;
	  c_im(inout[0]) = tmp7 + tmp8;
	  c_im(inout[iostride]) = tmp8 - tmp7;
     }
}

void realfft (
  double *in,
  unsigned n,
  double *out
)
{
	COMPLEX *c_in, *c_out;
	unsigned i;

	if (n == 0 ||
	    (c_in = (COMPLEX *) malloc (n * sizeof (COMPLEX))) == 0 ||
	    (c_out = (COMPLEX *) malloc (n * sizeof (COMPLEX))) == 0)
		return;
	
	for (i = 0; i < n; i++) {
		c_re (c_in [i]) = in [i];
		c_im (c_in [i]) = 0;
	}

	fft (c_in, n, c_out);

	out [0] = c_re (c_out [0]);		/* cos van dc */
	for (i = 1; i < (n + 1) / 2; i++) {	/* cos/sin i-de harmonische */
		out [2 * i - 1] = c_re (c_out [i]) * 2;
		out [2 * i] = c_im (c_out [i]) * -2;
	}
	if (n % 2 == 0)				/* cos van Nyquist */
		out [n - 1] = c_re (c_out [n / 2]);

	free ((char *) c_in);
	free ((char *) c_out);
}

static void mkre01(C *a, int n) /* n should be be multiple of 4 */
{
     R a0;
     a0 = c_re(a[0]);
     mko00(a, n/2, 0);
     c_re(a[n/2]) = -(c_re(a[0]) = a0);
     mkre00(a, n);
}

/* C = A - B */
void asub(C *c, C *a, C *b, int n)
{
     int i;

     for (i = 0; i < n; ++i) {
	  c_re(c[i]) = c_re(a[i]) - c_re(b[i]);
	  c_im(c[i]) = c_im(a[i]) - c_im(b[i]);
     }
}

/* C = A */
void acopy(C *c, C *a, int n)
{
     int i;

     for (i = 0; i < n; ++i) {
	  c_re(c[i]) = c_re(a[i]);
	  c_im(c[i]) = c_im(a[i]);
     }
}

/* C = A - B */
void array_sub(fftw_complex *C, fftw_complex *A, fftw_complex *B, int n)
{
     int i;

     for (i = 0; i < n; ++i) {
	  c_re(C[i]) = c_re(A[i]) - c_re(B[i]);
	  c_im(C[i]) = c_im(A[i]) - c_im(B[i]);
     }
}

void accuracy_test(dofft_closure *k, aconstrain constrain,
		   int sign, int n, C *a, C *b, int rounds, int impulse_rounds,
		   double t[6])
{
     int r, i;
     int ntests = 0;
     bench_complex czero = {0, 0};

     for (i = 0; i < 6; ++i) t[i] = 0.0;

     for (r = 0; r < rounds; ++r) {
	  arand(a, n);
	  if (one_accuracy_test(k, constrain, sign, n, a, b, t))
	       ++ntests;
     }

     /* impulses at beginning of array */
     for (r = 0; r < impulse_rounds; ++r) {
	  if (r > n - r - 1)
	       continue;
	  
	  caset(a, n, czero);
	  c_re(a[r]) = c_im(a[r]) = 1.0;
	  
	  if (one_accuracy_test(k, constrain, sign, n, a, b, t))
	       ++ntests;
     }
     
     /* impulses at end of array */
     for (r = 0; r < impulse_rounds; ++r) {
	  if (r <= n - r - 1)
	       continue;
	  
	  caset(a, n, czero);
	  c_re(a[n - r - 1]) = c_im(a[n - r - 1]) = 1.0;
	  
	  if (one_accuracy_test(k, constrain, sign, n, a, b, t))
	       ++ntests;
     }
     
     /* randomly-located impulses */
     for (r = 0; r < impulse_rounds; ++r) {
	  caset(a, n, czero);
	  i = rand() % n;
	  c_re(a[i]) = c_im(a[i]) = 1.0;
	  
	  if (one_accuracy_test(k, constrain, sign, n, a, b, t))
	       ++ntests;
     }

     t[0] /= ntests;
     t[1] = sqrt(t[1] / ntests);
     t[3] /= ntests;
     t[4] = sqrt(t[4] / ntests);

     fftaccuracy_done();
}

void fftw_no_twiddle_1(const fftw_complex *input, fftw_complex *output, int istride, int ostride)
{
     fftw_real tmp1;
     fftw_real tmp2;
     ASSERT_ALIGNED_DOUBLE;
     tmp1 = c_re(input[0]);
     c_re(output[0]) = tmp1;
     tmp2 = c_im(input[0]);
     c_im(output[0]) = tmp2;
}