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
|
#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
#define NUM_STEPS 200
typedef enum {
NONE = 0,
COS = 1,
SIN = 2
} function_t;
static clock_t DoFPUTransIteration(double *abase,
double *bbase,
unsigned long arraysize);
static inline double TrapezoidIntegrate(double x0,
double x1,
int n,
function_t 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 = 64;
if (is_adjusted == false) {
is_adjusted = true;
do {
array_size += 64;
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 NUM_STEPS fourier coefficients of the function (x+1)^x on
** the interval 0,2.
*/
static clock_t
DoFPUTransIteration(double *abase, double *bbase, unsigned long arraysize)
{
clock_t start, stop;
int i;
start = clock();
/*
** Calculate the fourier series. Begin by
** calculating A[0], B[0]
*/
abase[0] = TrapezoidIntegrate(0.0, 2.0, 0, NONE) / 2.0;
bbase[0] = TrapezoidIntegrate(0.0, 2.0, 0, NONE) / 2.0;
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, i, COS);
/*
** Calculate the B[i] terms.
*/
bbase[i] = TrapezoidIntegrate(0.0, 2.0, i, SIN);
}
stop = clock();
return stop - start;
}
/***********************
** TrapezoidIntegrate **
************************
** Perform a simple trapezoid integration on the
** function (x+1)**x.
** double x0 - lower bound
** double x1 - upper bound
** int n - series number
** int select - select functions FIXME: this is dumb
*/
static inline double
TrapezoidIntegrate(double x0, double x1, int n, function_t select)
{
double dx = (x1 - x0) / (double)NUM_STEPS; /* Stepsize */
double rvalue = pow(x0 + 1.0, x0);
int num_steps = NUM_STEPS - 2; /* Already done 1 step */
switch (select) {
case NONE:
while(num_steps) {
x0 += dx;
rvalue += pow(x0 + 1.0, x0);
num_steps--;
}
rvalue += pow(x1 + 1.0, x1);
break;
case COS:
rvalue *= cos(M_PI * n * x0);
while(num_steps) {
x0 += dx;
rvalue += pow(x0 + 1.0, x0) * cos(M_PI * n * x0);
num_steps--;
}
rvalue += pow(x1 + 1.0, x1) * cos(M_PI * n * x1);
break;
case SIN:
rvalue *= sin(M_PI * n * x0);
while(num_steps) {
x0 += dx;
rvalue += pow(x0 + 1.0, x0) * sin(M_PI * n * x0);
num_steps--;
}
rvalue += pow(x1 + 1.0, x1) * sin(M_PI * n * x1);
break;
}
return rvalue / 2.0 * dx;
}
|
