1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
|
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include <limits.h>
#include <time.h>
#include "cleanbench.h"
#include "randnum.h"
/*************************
** FOURIER COEFFICIENTS **
*************************/
/* M_PI isn't defined if compiled with -ansi */
#ifndef M_PI
#define M_PI 3.14159265358979323846 /* pi */
#endif
static clock_t DoFPUTransIteration(double *abase,
double *bbase,
unsigned long arraysize);
static double TrapezoidIntegrate(double x0,
double x1,
int nsteps,
double omega_n,
int select);
static double thefunction(double x,
double omega_n,
int select);
/**************
** DoFourier **
***************
** Perform the transcendental/trigonometric portion of the
** benchmark. This benchmark calculates the first n
** fourier coefficients of the function (x+1)^x defined
** on the interval 0,2.
*/
double
DoFourier(void)
{
double* abase = NULL;
double* bbase = NULL;
clock_t total_time = 0;
int iterations = 0;
static bool is_adjusted = false;
static int array_size = 50;
if (is_adjusted == false) {
is_adjusted = true;
do {
array_size += 50;
abase = realloc(abase, array_size * sizeof(double));
bbase = realloc(bbase, array_size * sizeof(double));
/*
** Do an iteration of the tests. If the elapsed time is
** less than or equal to the permitted minimum, re-allocate
** larger arrays and try again.
*/
} while (DoFPUTransIteration(abase,bbase, array_size) <= MINIMUM_TICKS);
} else {
/*
** Don't need self-adjustment. Just allocate the
** arrays, and go.
*/
abase = malloc(array_size * sizeof(double));
bbase = malloc(array_size * sizeof(double));
}
do {
total_time += DoFPUTransIteration(abase,bbase,array_size);
iterations += array_size * 2 - 1;
} while (total_time < MINIMUM_SECONDS * CLOCKS_PER_SEC);
free(abase);
free(bbase);
return (double)(iterations * CLOCKS_PER_SEC) / (double)total_time;
}
/************************
** DoFPUTransIteration **
*************************
** Perform an iteration of the FPU Transcendental/trigonometric
** benchmark. Here, an iteration consists of calculating the
** first n fourier coefficients of the function (x+1)^x on
** the interval 0,2. n is given by arraysize.
** NOTE: The # of integration steps is fixed at
** 200.
*/
static clock_t
DoFPUTransIteration(double *abase, double *bbase, unsigned long arraysize)
{
clock_t start, stop;
double omega; /* Fundamental frequency */
unsigned long i; /* Index */
start = clock();
/*
** Calculate the fourier series. Begin by
** calculating A[0].
*/
*abase = TrapezoidIntegrate(0.0, 2.0, 200, 0.0, 0) / 2.0;
/*
** Calculate the fundamental frequency.
** ( 2 * pi ) / period...and since the period
** is 2, omega is simply pi.
*/
omega = M_PI;
for(i = 1; i < arraysize; i++) {
/*
** Calculate A[i] terms. Note, once again, that we
** can ignore the 2/period term outside the integral
** since the period is 2 and the term cancels itself
** out.
*/
*(abase + i) = TrapezoidIntegrate(0.0, 2.0, 200, omega * (double)i, 1);
/*
** Calculate the B[i] terms.
*/
*(bbase + i) = TrapezoidIntegrate(0.0, 2.0, 200, omega * (double)i, 2);
}
stop = clock();
return stop - start;
}
/***********************
** TrapezoidIntegrate **
************************
** Perform a simple trapezoid integration on the
** function (x+1)**x.
** x0,x1 set the lower and upper bounds of the
** integration.
** nsteps indicates # of trapezoidal sections
** omega_n is the fundamental frequency times
** the series member #
** select = 0 for the A[0] term, 1 for cosine terms, and
** 2 for sine terms.
** Returns the value.
** double x0 - lower bound
** double x1 - upper bound
** int nsteps - number of steps
** double omega_n - omega * n
** int select - select functions FIXME: this is dumb
*/
static double
TrapezoidIntegrate( double x0, double x1, int nsteps, double omega_n, int select)
{
double x = x0; /* Independent variable */
double dx = (x1 - x0) / (double)nsteps; /* Stepsize */
double rvalue = thefunction(x0, omega_n, select) / 2.0;
/*
** Compute the other terms of the integral.
*/
if(nsteps!=1)
{ --nsteps; /* Already done 1 step */
while(--nsteps )
{
x+=dx;
rvalue+=thefunction(x,omega_n,select);
}
}
/*
** Finish computation
*/
rvalue=(rvalue+thefunction(x1,omega_n,select)/(double)2.0)*dx;
return(rvalue);
}
/****************
** thefunction **
*****************
** This routine selects the function to be used
** in the Trapezoid integration.
** x is the independent variable
** omega_n is omega * n
** select chooses which of the sine/cosine functions
** are used. note the special case for select=0.
*/
static double thefunction(double x, /* Independent variable */
double omega_n, /* Omega * term */
int select) /* Choose term */
{
/*
** Use select to pick which function we call.
*/
switch(select)
{
case 0: return(pow(x+(double)1.0,x));
case 1: return(pow(x+(double)1.0,x) * cos(omega_n * x));
case 2: return(pow(x+(double)1.0,x) * sin(omega_n * x));
}
/*
** We should never reach this point, but the following
** keeps compilers from issuing a warning message.
*/
return(0.0);
}
|
