1 
  2 /* bitreversal performance comparison
  3    Bracewell-Buneman and Yu-Loh-Miller
  4    John Bryan 2018
  5    -lm -lrt -O3 
  6    gcc 4.6.3
  7 */
  8 
  9 
 10 #include <stdlib.h>
 11 #include <stdio.h>
 12 #include <time.h>
 13 #include <complex.h>
 14 #include <math.h>
 15 #include <stdint.h>
 16 #include <inttypes.h>
 17 #include <assert.h>
 18 #include <float.h>
 19 
 20 
 21 #define BILLION 1000000000L
 22 // rows is the number of lengths
 23 #define rows 9    
 24 // columns is the number of iterations averaged
 25 #define columns 20000  
 26 
 27 
 28 void bracewell_buneman(complex double *x, int m, int N)
 29 {
 30      int M=1, i,  T, k, r[N], muplus, q, p, upper_range, nu;
 31      int nprime, rprime, rprimeprime;
 32      complex double temp;
 33      muplus = (m+1)>>1;
 34      for (i=0; i < N; i++)
 35          {r[i] = 0;}
 36      upper_range = muplus + 1;
 37      for (nu=1; nu < upper_range; nu++)
 38      {
 39          for (k=0; k<M; k++)
 40          {
 41              T = r[k] << 1;
 42              r[k] = T;
 43              r[k+M] = T + 1;
 44           }
 45           M = M << 1;
 46      }
 47      if (m & 0x01) {M = M >> 1;};
 48      for (q=1; q < M; q++)
 49      {
 50           nprime = q - M;
 51           rprimeprime = r[q] * M;
 52           for (p=0; p < r[q]; p++)
 53           {
 54               nprime = nprime + M;
 55               rprime = rprimeprime + r[p];
 56               temp = x[nprime];
 57               x[nprime] = x[rprime];
 58               x[rprime] = temp;
 59           }
 60       }
 61       return;
 62 }
 63 
 64 
 65 void yu_loh_miller(complex double *x, int p, int n)
 66 {
 67     int n1, i, j, k, w, u, v, y;
 68     int r[n];
 69     complex double temp;
 70     if (!(p & 0x01))
 71         n1 = (int)sqrt(n);       //seed table size 
 72     else
 73         n1 = (int)sqrt(n>>1);
 74     // algorithm 2, compute seed table  
 75     for (i=0; i<n; i++)
 76        {r[i] = 0;}
 77     r[1] = n>>1;
 78     for (i=1; i<(n1>>1); i++)
 79     {
 80         r[i<<1] = r[i]>>1;
 81         r[(i<<1)+1] = r[(i<<1)] + r[1];
 82     }
 83     //algorithm 1
 84     for (i=0; i < (n1-1); i++)
 85     {
 86         for (j = i+1; j < n1; j++)
 87         {
 88             u = i + r[j];
 89             v = j + r[i];
 90             temp = x[u];
 91             x[u] = x[v];
 92             x[v] = temp;
 93             if (p & 0x01)
 94             {
 95                 w = u + n1;
 96                 y = v + n1;
 97                 temp = x[w];
 98                 x[w] = x[y];
 99                 x[y] = temp;
100             }
101         }
102     }
103     return;
104 }
105 
106 
107 double randf(void)
108 {
109     //return random value between -1 and 1
110     double low = -1.0, high = 1.0;
111     return (rand()/(double)(RAND_MAX))*abs(low-high)+low;
112 }
113 
114 
115 
116 uint64_t dAvg(uint64_t sum){
117      return sum/columns;
118 }
119 
120 
121 int main()
122 {
123    int j, i, N, m, k = 0, b = 2;
124    uint64_t  bb_mean[rows], ylm_mean[rows];
125    uint64_t bb_sum,ylm_sum;
126    complex double *x, *y;
127    uint64_t ylm_diff,bb_diff;
128    struct timespec start,end;
129    FILE *f = fopen("file2.txt", "w");
130    for (j = 0; j < rows; j++)
131    {
132       m = j+5;
133       printf("%d more sequence length(s) to compute performance... \n",rows-j);
134       N = (int)powf(2.0, (float)(m));
135       x = malloc(sizeof(complex double)*N);
136       y = malloc(sizeof(complex double)*N);
137       ylm_sum = 0;
138       bb_sum = 0;
139       for (k = 0; k < columns; k++)
140       {
141           srand(k);
142           for (i = 0; i < N; i++)
143               {x[i] = randf()+(randf()*I);}
144           srand(k);
145           for (i = 0; i < N; i++)
146               {y[i] = randf()+(randf()*I);}
147           clock_gettime(CLOCK_MONOTONIC, &start);
148           bracewell_buneman(x, m, N);
149           clock_gettime(CLOCK_MONOTONIC, &end);
150           bb_diff = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
151           bb_sum = bb_sum + bb_diff;
152           clock_gettime(CLOCK_MONOTONIC, &start);
153           yu_loh_miller(x, m, N);
154           clock_gettime(CLOCK_MONOTONIC, &end);
155           ylm_diff = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
156           ylm_sum = ylm_sum + ylm_diff;
157           // involution check: bitreversal(bitreversal(x))=x.
158           for (i = 0; i < N; i++)
159               if (x[i] != y[i])
160                   assert(0);
161       }
162       free(x);
163       free(y);
164       ylm_mean[j]=dAvg(ylm_sum);
165       bb_mean[j]=dAvg(bb_sum);
166    }
167    for ( i = 0; i < rows; i++ )
168        fprintf (f, "  %3d   "  "%"PRIu64 "\n",(int)powf(2.0, (float)(i+5)), ylm_mean[i]);
169    for ( i = 0; i < rows; i++ )
170        fprintf (f, "  %3d   "  "%"PRIu64 "\n",(int)powf(2.0, (float)(i+5)), bb_mean[i]);
171    fclose(f);
172    system("python plot_br.py");
173    exit(0);
174 }