From c471cfe8ed04c68bd3ba96de578160018676966f Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@fftw.org>
Date: Sat, 1 Mar 2003 19:15:18 -0500
Subject: [PATCH] use K for constants

---
 dft/zero.c                |  6 +++---
 kernel/trig1.c            | 10 +++++-----
 rdft/generic.c            | 16 ++++++++--------
 rdft/problem.c            |  6 +++---
 rdft/rdft-dht.c           |  6 +++---
 rdft/rdft2-radix2.c       | 22 +++++++++++-----------
 reodft/redft00e-r2hc.c    |  6 +++---
 reodft/reodft010e-r2hc.c  | 22 +++++++++++-----------
 reodft/reodft11e-r2hc.c   | 14 +++++++-------
 reodft/reodft11e-radix2.c | 22 +++++++++++-----------
 reodft/rodft00e-r2hc.c    |  6 +++---
 11 files changed, 68 insertions(+), 68 deletions(-)

diff --git a/dft/zero.c b/dft/zero.c
index 07dd3ee8d..5e9512c6f 100644
--- a/dft/zero.c
+++ b/dft/zero.c
@@ -18,7 +18,7 @@
  *
  */
 
-/* $Id: zero.c,v 1.3 2003-01-15 02:10:25 athena Exp $ */
+/* $Id: zero.c,v 1.4 2003-03-02 00:15:18 stevenj Exp $ */
 
 #include "dft.h"
 
@@ -28,7 +28,7 @@ static void recur(const iodim *dims, int rnk, R *ri, R *ii)
      if (rnk == RNK_MINFTY)
           return;
      else if (rnk == 0)
-          ri[0] = ii[0] = 0.0;
+          ri[0] = ii[0] = K(0.0);
      else if (rnk > 0) {
           int i, n = dims[0].n;
           int is = dims[0].is;
@@ -36,7 +36,7 @@ static void recur(const iodim *dims, int rnk, R *ri, R *ii)
 	  if (rnk == 1) {
 	       /* this case is redundant but faster */
 	       for (i = 0; i < n; ++i)
-		    ri[i * is] = ii[i * is] = 0.0;
+		    ri[i * is] = ii[i * is] = K(0.0);
 	  } else {
 	       for (i = 0; i < n; ++i)
 		    recur(dims + 1, rnk - 1, ri + i * is, ii + i * is);
diff --git a/kernel/trig1.c b/kernel/trig1.c
index e73a55fc7..e155cfc6d 100644
--- a/kernel/trig1.c
+++ b/kernel/trig1.c
@@ -18,7 +18,7 @@
  *
  */
 
-/* $Id: trig1.c,v 1.1 2002-09-23 22:49:10 athena Exp $ */
+/* $Id: trig1.c,v 1.2 2003-03-02 00:15:18 stevenj Exp $ */
 
 /* trigonometric functions */
 #include "ifftw.h"
@@ -48,10 +48,10 @@ static const trigreal K2PI =
 trigreal X(sincos)(trigreal m, trigreal n, int sinp)
 {
      /* waiting for C to get tail recursion... */
-     trigreal half_n = n * 0.5;
-     trigreal quarter_n = half_n * 0.5;
-     trigreal eighth_n = quarter_n * 0.5;
-     trigreal sgn = 1.0;
+     trigreal half_n = n * KTRIG(0.5);
+     trigreal quarter_n = half_n * KTRIG(0.5);
+     trigreal eighth_n = quarter_n * KTRIG(0.5);
+     trigreal sgn = KTRIG(1.0);
 
      if (sinp) goto sin;
  cos:
diff --git a/rdft/generic.c b/rdft/generic.c
index 1e4e0b291..da34b67b6 100644
--- a/rdft/generic.c
+++ b/rdft/generic.c
@@ -67,8 +67,8 @@ static void apply_dit(const plan *ego_, R *I, R *O)
 
      /* compute the transform of the r 0th elements (which are real) */
      for (i = 0; i + i < r; ++i) {
-	  rsum = 0.0;
-	  isum = 0.0;
+	  rsum = K(0.0);
+	  isum = K(0.0);
 	  wincr = m * i;
 	  for (j = 0, wp = 0; j < r; ++j) {
 	       E tw_r = W[2*wp];
@@ -98,8 +98,8 @@ static void apply_dit(const plan *ego_, R *I, R *O)
      /* compute the transform of the middle elements (which are complex) */
      for (k = 1; k + k < m; ++k, X += os, YI -= os, YO -= os) {
 	  for (i = 0; i < r; ++i) {
-	       rsum = 0.0;
-	       isum = 0.0;
+	       rsum = K(0.0);
+	       isum = K(0.0);
 	       wincr = k + m * i;
 	       for (j = 0, wp = 0; j < r; ++j) {
 		    E tw_r = W[2*wp];
@@ -168,7 +168,7 @@ static void apply_dif(const plan *ego_, R *I, R *O)
      }
 
      for (i = 0; i < r; ++i) {
-	  rsum = 0.0;
+	  rsum = K(0.0);
 	  wincr = m * i;
 	  for (j = 1, wp = wincr; j + j < r; ++j) {
 	       E tw_r = W[2*wp];
@@ -180,7 +180,7 @@ static void apply_dif(const plan *ego_, R *I, R *O)
 	       if (wp >= n)
 		    wp -= n;
 	  }
-	  X[i * ism] = 2.0 * rsum + buf[0];
+	  X[i * ism] = K(2.0) * rsum + buf[0];
      }
 
      X += is;
@@ -200,8 +200,8 @@ static void apply_dif(const plan *ego_, R *I, R *O)
 	  }
 
 	  for (i = 0; i < r; ++i) {
-	       rsum = 0.0;
-	       isum = 0.0;
+	       rsum = K(0.0);
+	       isum = K(0.0);
 	       wincr = m * i;
 	       for (j = 0, wp = k * i; j < r; ++j) {
 		    E tw_r = W[2*wp];
diff --git a/rdft/problem.c b/rdft/problem.c
index 767df6a26..dcd692885 100644
--- a/rdft/problem.c
+++ b/rdft/problem.c
@@ -18,7 +18,7 @@
  *
  */
 
-/* $Id: problem.c,v 1.40 2003-02-27 18:44:19 stevenj Exp $ */
+/* $Id: problem.c,v 1.41 2003-03-02 00:15:18 stevenj Exp $ */
 
 #include "rdft.h"
 #include <stddef.h>
@@ -57,7 +57,7 @@ static void recur(const iodim *dims, int rnk, R *I)
      if (rnk == RNK_MINFTY)
           return;
      else if (rnk == 0)
-          I[0] = 0.0;
+          I[0] = K(0.0);
      else if (rnk > 0) {
           int i, n = dims[0].n;
           int is = dims[0].is;
@@ -65,7 +65,7 @@ static void recur(const iodim *dims, int rnk, R *I)
 	  if (rnk == 1) {
 	       /* this case is redundant but faster */
 	       for (i = 0; i < n; ++i)
-		    I[i * is] = 0.0;
+		    I[i * is] = K(0.0);
 	  } else {
 	       for (i = 0; i < n; ++i)
 		    recur(dims + 1, rnk - 1, I + i * is);
diff --git a/rdft/rdft-dht.c b/rdft/rdft-dht.c
index 125a43284..3db4e04b2 100644
--- a/rdft/rdft-dht.c
+++ b/rdft/rdft-dht.c
@@ -18,7 +18,7 @@
  *
  */
 
-/* $Id: rdft-dht.c,v 1.15 2003-02-28 23:28:58 stevenj Exp $ */
+/* $Id: rdft-dht.c,v 1.16 2003-03-02 00:15:18 stevenj Exp $ */
 
 /* Solve an R2HC/HC2R problem via post/pre processing of a DHT.  This
    is mainly useful because we can use Rader to compute DHTs of prime
@@ -53,8 +53,8 @@ static void apply_r2hc(const plan *ego_, R *I, R *O)
      os = ego->os;
      for (i = 1; i < n - i; ++i) {
 	  E a, b;
-	  a = 0.5 * O[os * i];
-	  b = 0.5 * O[os * (n - i)];
+	  a = K(0.5) * O[os * i];
+	  b = K(0.5) * O[os * (n - i)];
 	  O[os * i] = a + b;
 #if FFT_SIGN == -1
 	  O[os * (n - i)] = b - a;
diff --git a/rdft/rdft2-radix2.c b/rdft/rdft2-radix2.c
index b12932f2d..ca9e8ef78 100644
--- a/rdft/rdft2-radix2.c
+++ b/rdft/rdft2-radix2.c
@@ -18,7 +18,7 @@
  *
  */
 
-/* $Id: rdft2-radix2.c,v 1.22 2003-02-28 23:28:58 stevenj Exp $ */
+/* $Id: rdft2-radix2.c,v 1.23 2003-03-02 00:15:18 stevenj Exp $ */
 
 /*
   Compute RDFT2 of even size via either a DFT or a vector RDFT of
@@ -129,8 +129,8 @@ static void k_f_dft(R *rio, R *iio, const R *W, int n, int dist)
           E rop = pp[0], iop = pp[im];
           pp[0] = rop + iop;
           pm[0] = rop - iop;
-          pp[im] = 0.0;
-          pm[im] = 0.0;
+          pp[im] = K(0.0);
+          pm[im] = K(0.0);
 	  pp += dist; pm -= dist;
      }
 
@@ -144,10 +144,10 @@ static void k_f_dft(R *rio, R *iio, const R *W, int n, int dist)
           E id = iop + iom;
           E tr = rd * wr - id * wi;
           E ti = id * wr + rd * wi;
-          pp[0] = 0.5 * (re + ti);
-          pp[im] = 0.5 * (ie + tr);
-          pm[0] = 0.5 * (re - ti);
-          pm[im] = 0.5 * (tr - ie);
+          pp[0] = K(0.5) * (re + ti);
+          pp[im] = K(0.5) * (ie + tr);
+          pm[0] = K(0.5) * (re - ti);
+          pm[im] = K(0.5) * (tr - ie);
 	  pp += dist; pm -= dist;
      }
 
@@ -201,8 +201,8 @@ static void k_f_rdft(R *rio, R *iio, const R *W, int n, int dist)
           E rop = pp[0], iop = pp[im];
           pp[0] = rop + iop;
           pm[0] = rop - iop;
-          pp[im] = 0.0;
-          pm[im] = 0.0;
+          pp[im] = K(0.0);
+          pm[im] = K(0.0);
 	  pp += dist; pm -= dist;
      }
 
@@ -291,7 +291,7 @@ static void k_b_dft(R *rio, R *iio, const R *W, int n, int dist)
      }
 
      /* i = n/2 when n is even */
-     if (!(n & 1)) { pp[0] *= 2.0; pp[im] *= -2.0; }
+     if (!(n & 1)) { pp[0] *= K(2.0); pp[im] *= -K(2.0); }
 }
 
 static void apply_b_dft(const plan *ego_, R *r, R *rio, R *iio)
@@ -357,7 +357,7 @@ static void k_b_rdft(R *rio, R *iio, const R *W, int n, int dist)
      }
 
      /* i = n/2 when n is even */
-     if (!(n & 1)) { pp[0] *= 2.0; pp[im] *= -2.0; }
+     if (!(n & 1)) { pp[0] *= K(2.0); pp[im] *= -K(2.0); }
 }
 
 static void apply_b_rdft(const plan *ego_, R *r, R *rio, R *iio)
diff --git a/reodft/redft00e-r2hc.c b/reodft/redft00e-r2hc.c
index 5a97ee804..89db5f4bb 100644
--- a/reodft/redft00e-r2hc.c
+++ b/reodft/redft00e-r2hc.c
@@ -18,7 +18,7 @@
  *
  */
 
-/* $Id: redft00e-r2hc.c,v 1.21 2003-02-28 23:28:58 stevenj Exp $ */
+/* $Id: redft00e-r2hc.c,v 1.22 2003-03-02 00:15:18 stevenj Exp $ */
 
 /* Do a REDFT00 problem via an R2HC problem, with some pre/post-processing. */
 
@@ -61,14 +61,14 @@ static void apply(const plan *ego_, R *I, R *O)
 	       E a, b, apb, amb;
 	       a = I[is * i];
 	       b = I[is * (n - i)];
-	       csum += W[2*i] * (amb = 2.0*(a - b));
+	       csum += W[2*i] * (amb = K(2.0)*(a - b));
 	       amb = W[2*i+1] * amb;
 	       apb = (a + b);
 	       buf[i] = apb - amb;
 	       buf[n - i] = apb + amb;
 	  }
 	  if (i == n - i) {
-	       buf[i] = 2.0 * I[is * i];
+	       buf[i] = K(2.0) * I[is * i];
 	  }
 	  
 	  {
diff --git a/reodft/reodft010e-r2hc.c b/reodft/reodft010e-r2hc.c
index cbe145c22..9d8eb6f53 100644
--- a/reodft/reodft010e-r2hc.c
+++ b/reodft/reodft010e-r2hc.c
@@ -18,7 +18,7 @@
  *
  */
 
-/* $Id: reodft010e-r2hc.c,v 1.24 2003-02-28 23:28:58 stevenj Exp $ */
+/* $Id: reodft010e-r2hc.c,v 1.25 2003-03-02 00:15:18 stevenj Exp $ */
 
 /* Do an R{E,O}DFT{01,10} problem via an R2HC problem, with some
    pre/post-processing ala FFTPACK. */
@@ -100,7 +100,7 @@ static void apply_re01(const plan *ego_, R *I, R *O)
 	       buf[n - i] = wa * apb - wb * amb; 
 	  }
 	  if (i == n - i) {
-	       buf[i] = 2.0 * I[is * i] * W[2*i];
+	       buf[i] = K(2.0) * I[is * i] * W[2*i];
 	  }
 	  
 	  {
@@ -154,7 +154,7 @@ static void apply_ro01(const plan *ego_, R *I, R *O)
 	       buf[n - i] = wa * apb - wb * amb; 
 	  }
 	  if (i == n - i) {
-	       buf[i] = 2.0 * I[is * (i - 1)] * W[2*i];
+	       buf[i] = K(2.0) * I[is * (i - 1)] * W[2*i];
 	  }
 	  
 	  {
@@ -211,18 +211,18 @@ static void apply_re10(const plan *ego_, R *I, R *O)
 	       cld->apply((plan *) cld, buf, buf);
 	  }
 	  
-	  O[0] = 2.0 * buf[0];
+	  O[0] = K(2.0) * buf[0];
 	  for (i = 1; i < n - i; ++i) {
 	       E a, b, wa, wb;
-	       a = 2.0 * buf[i];
-	       b = 2.0 * buf[n - i];
+	       a = K(2.0) * buf[i];
+	       b = K(2.0) * buf[n - i];
 	       wa = W[2*i];
 	       wb = W[2*i + 1];
 	       O[os * i] = wa * a + wb * b;
 	       O[os * (n - i)] = wb * a - wa * b;
 	  }
 	  if (i == n - i) {
-	       O[os * i] = 2.0 * buf[i] * W[2*i];
+	       O[os * i] = K(2.0) * buf[i] * W[2*i];
 	  }
      }
 
@@ -262,18 +262,18 @@ static void apply_ro10(const plan *ego_, R *I, R *O)
 	       cld->apply((plan *) cld, buf, buf);
 	  }
 	  
-	  O[os * (n - 1)] = 2.0 * buf[0];
+	  O[os * (n - 1)] = K(2.0) * buf[0];
 	  for (i = 1; i < n - i; ++i) {
 	       E a, b, wa, wb;
-	       a = 2.0 * buf[i];
-	       b = 2.0 * buf[n - i];
+	       a = K(2.0) * buf[i];
+	       b = K(2.0) * buf[n - i];
 	       wa = W[2*i];
 	       wb = W[2*i + 1];
 	       O[os * (n - 1 - i)] = wa * a + wb * b;
 	       O[os * (i - 1)] = wb * a - wa * b;
 	  }
 	  if (i == n - i) {
-	       O[os * (i - 1)] = 2.0 * buf[i] * W[2*i];
+	       O[os * (i - 1)] = K(2.0) * buf[i] * W[2*i];
 	  }
      }
 
diff --git a/reodft/reodft11e-r2hc.c b/reodft/reodft11e-r2hc.c
index 25bbf1791..d1dc3d9a1 100644
--- a/reodft/reodft11e-r2hc.c
+++ b/reodft/reodft11e-r2hc.c
@@ -18,7 +18,7 @@
  *
  */
 
-/* $Id: reodft11e-r2hc.c,v 1.23 2003-02-28 23:28:58 stevenj Exp $ */
+/* $Id: reodft11e-r2hc.c,v 1.24 2003-03-02 00:15:18 stevenj Exp $ */
 
 /* Do an R{E,O}DFT11 problem via an R2HC problem, with some
    pre/post-processing ala FFTPACK.  Use a trick from: 
@@ -67,10 +67,10 @@ static void apply_re11(const plan *ego_, R *I, R *O)
      for (iv = 0; iv < vl; ++iv, I += ivs, O += ovs) {
 	  /* I wish that this didn't require an extra pass. */
 	  /* FIXME: use recursive/cascade summation for better stability? */
-	  buf[n - 1] = cur = 2.0 * I[is * (n - 1)];
+	  buf[n - 1] = cur = K(2.0) * I[is * (n - 1)];
 	  for (i = n - 1; i > 0; --i) {
 	       E curnew;
-	       buf[(i - 1)] = curnew = 2.0 * I[is * (i - 1)] - cur;
+	       buf[(i - 1)] = curnew = K(2.0) * I[is * (i - 1)] - cur;
 	       cur = curnew;
 	  }
 	  
@@ -87,7 +87,7 @@ static void apply_re11(const plan *ego_, R *I, R *O)
 	       buf[n - i] = wa * apb - wb * amb; 
 	  }
 	  if (i == n - i) {
-	       buf[i] = 2.0 * buf[i] * W[2*i];
+	       buf[i] = K(2.0) * buf[i] * W[2*i];
 	  }
 	  
 	  {
@@ -132,10 +132,10 @@ static void apply_ro11(const plan *ego_, R *I, R *O)
      for (iv = 0; iv < vl; ++iv, I += ivs, O += ovs) {
 	  /* I wish that this didn't require an extra pass. */
 	  /* FIXME: use recursive/cascade summation for better stability? */
-	  buf[n - 1] = cur = 2.0 * I[0];
+	  buf[n - 1] = cur = K(2.0) * I[0];
 	  for (i = n - 1; i > 0; --i) {
 	       E curnew;
-	       buf[(i - 1)] = curnew = 2.0 * I[is * (n - i)] - cur;
+	       buf[(i - 1)] = curnew = K(2.0) * I[is * (n - i)] - cur;
 	       cur = curnew;
 	  }
 	  
@@ -152,7 +152,7 @@ static void apply_ro11(const plan *ego_, R *I, R *O)
 	       buf[n - i] = wa * apb - wb * amb; 
 	  }
 	  if (i == n - i) {
-	       buf[i] = 2.0 * buf[i] * W[2*i];
+	       buf[i] = K(2.0) * buf[i] * W[2*i];
 	  }
 	  
 	  {
diff --git a/reodft/reodft11e-radix2.c b/reodft/reodft11e-radix2.c
index 7876995bb..707fced20 100644
--- a/reodft/reodft11e-radix2.c
+++ b/reodft/reodft11e-radix2.c
@@ -18,7 +18,7 @@
  *
  */
 
-/* $Id: reodft11e-radix2.c,v 1.5 2003-03-01 00:01:20 stevenj Exp $ */
+/* $Id: reodft11e-radix2.c,v 1.6 2003-03-02 00:15:18 stevenj Exp $ */
 
 /* Do an R{E,O}DFT11 problem of *even* size by a pair of R2HC problems
    of half the size, plus some pre/post-processing.  Use a trick from:
@@ -71,8 +71,8 @@ static void apply_re11(const plan *ego_, R *I, R *O)
      buf = (R *) MALLOC(sizeof(R) * n, BUFFERS);
 
      for (iv = 0; iv < vl; ++iv, I += ivs, O += ovs) {
-	  buf[0] = 2.0 * I[0];
-	  buf[n2] = 2.0 * I[is * (n - 1)];
+	  buf[0] = K(2.0) * I[0];
+	  buf[n2] = K(2.0) * I[is * (n - 1)];
 	  for (i = 1; i + i < n2; ++i) {
 	       int k = i + i;
 	       E a, b, a2, b2;
@@ -114,8 +114,8 @@ static void apply_re11(const plan *ego_, R *I, R *O)
 	       E u, v;
 	       u = I[is * (n2 - 1)];
 	       v = I[is * n2];
-	       buf[i] = 2.0 * (u + v) * W[2*i];
-	       buf[n - i] = 2.0 * (u - v) * W[2*i];
+	       buf[i] = K(2.0) * (u + v) * W[2*i];
+	       buf[n - i] = K(2.0) * (u - v) * W[2*i];
 	  }
 
 
@@ -204,8 +204,8 @@ static void apply_re11(const plan *ego_, R *I, R *O)
      buf = (R *) MALLOC(sizeof(R) * n, BUFFERS);
 
      for (iv = 0; iv < vl; ++iv, I += ivs, O += ovs) {
-	  buf[0] = 2.0 * I[0];
-	  buf[n/2] = 2.0 * I[is * (n - 1)];
+	  buf[0] = K(2.0) * I[0];
+	  buf[n/2] = K(2.0) * I[is * (n - 1)];
 	  for (i = 1; i + i < n; ++i) {
 	       int k = i + i;
 	       E a, b;
@@ -279,8 +279,8 @@ static void apply_ro11(const plan *ego_, R *I, R *O)
      buf = (R *) MALLOC(sizeof(R) * n, BUFFERS);
 
      for (iv = 0; iv < vl; ++iv, I += ivs, O += ovs) {
-	  buf[0] = 2.0 * I[is * (n - 1)];
-	  buf[n2] = 2.0 * I[0];
+	  buf[0] = K(2.0) * I[is * (n - 1)];
+	  buf[n2] = K(2.0) * I[0];
 	  for (i = 1; i + i < n2; ++i) {
 	       int k = i + i;
 	       E a, b, a2, b2;
@@ -322,8 +322,8 @@ static void apply_ro11(const plan *ego_, R *I, R *O)
 	       E u, v;
 	       u = I[is * n2];
 	       v = I[is * (n2 - 1)];
-	       buf[i] = 2.0 * (u + v) * W[2*i];
-	       buf[n - i] = 2.0 * (u - v) * W[2*i];
+	       buf[i] = K(2.0) * (u + v) * W[2*i];
+	       buf[n - i] = K(2.0) * (u - v) * W[2*i];
 	  }
 
 
diff --git a/reodft/rodft00e-r2hc.c b/reodft/rodft00e-r2hc.c
index 4e09d3734..5bf2859b3 100644
--- a/reodft/rodft00e-r2hc.c
+++ b/reodft/rodft00e-r2hc.c
@@ -18,7 +18,7 @@
  *
  */
 
-/* $Id: rodft00e-r2hc.c,v 1.22 2003-02-28 23:28:58 stevenj Exp $ */
+/* $Id: rodft00e-r2hc.c,v 1.23 2003-03-02 00:15:18 stevenj Exp $ */
 
 /* Do a RODFT00 problem via an R2HC problem, with some pre/post-processing. */
 
@@ -59,13 +59,13 @@ static void apply(const plan *ego_, R *I, R *O)
 	       E a, b, apb, amb;
 	       a = I[is * (i - 1)];
 	       b = I[is * ((n - i) - 1)];
-	       apb =  2.0 * W[i] * (a + b);
+	       apb =  K(2.0) * W[i] * (a + b);
 	       amb = (a - b);
 	       buf[i] = apb + amb;
 	       buf[n - i] = apb - amb;
 	  }
 	  if (i == n - i) {
-	       buf[i] = 4.0 * I[is * (i - 1)];
+	       buf[i] = K(4.0) * I[is * (i - 1)];
 	  }
 	  
 	  {