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
|
/* test of fftsg3d.c */
#include <math.h>
#include <stdio.h>
#include "alloc.h"
#define MAX(x,y) ((x) > (y) ? (x) : (y))
/* random number generator, 0 <= RND < 1 */
#define RND(p) ((*(p) = (*(p) * 7141 + 54773) % 259200) * (1.0 / 259200))
int main()
{
void cdft3d(int, int, int, int, double ***, double *, int *, double *);
void rdft3d(int, int, int, int, double ***, double *, int *, double *);
void ddct3d(int, int, int, int, double ***, double *, int *, double *);
void ddst3d(int, int, int, int, double ***, double *, int *, double *);
void putdata3d(int n1, int n2, int n3, double ***a);
double errorcheck3d(int n1, int n2, int n3, double scale, double ***a);
int *ip, n1, n2, n3, n, nt, i, j;
double ***a, *w, err;
printf("data length n1=? (n1 = power of 2) \n");
scanf("%d", &n1);
printf("data length n2=? (n2 = power of 2) \n");
scanf("%d", &n2);
printf("data length n3=? (n3 = power of 2) \n");
scanf("%d", &n3);
a = alloc_3d_double(n1, n2, n3);
nt = MAX(n1, n2);
n = MAX(nt, n3 / 2);
ip = alloc_1d_int(2 + (int) sqrt(n + 0.5));
n = MAX(nt, n3) * 3 / 2;
w = alloc_1d_double(n);
ip[0] = 0;
/* check of CDFT */
putdata3d(n1, n2, n3, a);
cdft3d(n1, n2, n3, 1, a, NULL, ip, w);
cdft3d(n1, n2, n3, -1, a, NULL, ip, w);
err = errorcheck3d(n1, n2, n3, 2.0 / n1 / n2 / n3, a);
printf("cdft3d err= %g \n", err);
/* check of RDFT */
putdata3d(n1, n2, n3, a);
rdft3d(n1, n2, n3, 1, a, NULL, ip, w);
rdft3d(n1, n2, n3, -1, a, NULL, ip, w);
err = errorcheck3d(n1, n2, n3, 2.0 / n1 / n2 / n3, a);
printf("rdft3d err= %g \n", err);
/* check of DDCT */
putdata3d(n1, n2, n3, a);
ddct3d(n1, n2, n3, 1, a, NULL, ip, w);
ddct3d(n1, n2, n3, -1, a, NULL, ip, w);
for (i = 0; i <= n1 - 1; i++) {
for (j = 0; j <= n2 - 1; j++) {
a[i][j][0] *= 0.5;
}
for (j = 0; j <= n3 - 1; j++) {
a[i][0][j] *= 0.5;
}
}
for (i = 0; i <= n2 - 1; i++) {
for (j = 0; j <= n3 - 1; j++) {
a[0][i][j] *= 0.5;
}
}
err = errorcheck3d(n1, n2, n3, 8.0 / n1 / n2 / n3, a);
printf("ddct3d err= %g \n", err);
/* check of DDST */
putdata3d(n1, n2, n3, a);
ddst3d(n1, n2, n3, 1, a, NULL, ip, w);
ddst3d(n1, n2, n3, -1, a, NULL, ip, w);
for (i = 0; i <= n1 - 1; i++) {
for (j = 0; j <= n2 - 1; j++) {
a[i][j][0] *= 0.5;
}
for (j = 0; j <= n3 - 1; j++) {
a[i][0][j] *= 0.5;
}
}
for (i = 0; i <= n2 - 1; i++) {
for (j = 0; j <= n3 - 1; j++) {
a[0][i][j] *= 0.5;
}
}
err = errorcheck3d(n1, n2, n3, 8.0 / n1 / n2 / n3, a);
printf("ddst3d err= %g \n", err);
free_1d_double(w);
free_1d_int(ip);
free_3d_double(a);
return 0;
}
void putdata3d(int n1, int n2, int n3, double ***a)
{
int j1, j2, j3, seed = 0;
for (j1 = 0; j1 <= n1 - 1; j1++) {
for (j2 = 0; j2 <= n2 - 1; j2++) {
for (j3 = 0; j3 <= n3 - 1; j3++) {
a[j1][j2][j3] = RND(&seed);
}
}
}
}
double errorcheck3d(int n1, int n2, int n3, double scale, double ***a)
{
int j1, j2, j3, seed = 0;
double err = 0, e;
for (j1 = 0; j1 <= n1 - 1; j1++) {
for (j2 = 0; j2 <= n2 - 1; j2++) {
for (j3 = 0; j3 <= n3 - 1; j3++) {
e = RND(&seed) - a[j1][j2][j3] * scale;
err = MAX(err, fabs(e));
}
}
}
return err;
}
|
