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
|
// gcc tutorial.c -lsleef -lsleefdft -lm
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <complex.h>
#include "sleef.h"
#include "sleefdft.h"
#define THRES 1e-4
typedef double complex cmpl;
cmpl omega(double n, double kn) {
return cexp((-2 * M_PI * _Complex_I / n) * kn);
}
void forward(cmpl *ts, cmpl *fs, int len) {
for(int k=0;k<len;k++) {
fs[k] = 0;
for(int n=0;n<len;n++) fs[k] += ts[n] * omega(len, n*k);
}
}
int main(int argc, char **argv) {
int n = 256;
if (argc == 2) n = 1 << atoi(argv[1]);
SleefDFT_setPlanFilePath("plan.txt", NULL, SLEEF_PLAN_AUTOMATIC);
double *sx = (double *)Sleef_malloc(n*2 * sizeof(double));
double *sy = (double *)Sleef_malloc(n*2 * sizeof(double));
struct SleefDFT *p = SleefDFT_double_init1d(n, sx, sy, SLEEF_MODE_FORWARD);
if (p == NULL) {
printf("SleefDFT initialization failed\n");
exit(-1);
}
cmpl *ts = (cmpl *)malloc(sizeof(cmpl)*n);
cmpl *fs = (cmpl *)malloc(sizeof(cmpl)*n);
for(int i=0;i<n;i++) {
ts[i] =
(2.0 * (rand() / (double)RAND_MAX) - 1) * 1.0 +
(2.0 * (rand() / (double)RAND_MAX) - 1) * _Complex_I;
sx[(i*2+0)] = creal(ts[i]);
sx[(i*2+1)] = cimag(ts[i]);
}
forward(ts, fs, n);
SleefDFT_double_execute(p, NULL, NULL);
int success = 1;
for(int i=0;i<n;i++) {
if ((fabs(sy[(i*2+0)] - creal(fs[i])) > THRES) ||
(fabs(sy[(i*2+1)] - cimag(fs[i])) > THRES)) {
success = 0;
}
}
printf("%s\n", success ? "OK" : "NG");
free(fs); free(ts);
Sleef_free(sy); Sleef_free(sx);
SleefDFT_dispose(p);
exit(success);
}
|
