![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to 3.i CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd / ccom / vax / tests|
Return to 3.i CVS log | Up to [Research Unix] / researchv10no / cmd / ccom / vax / tests |
1.1 root 1: extern double fabs(), floor(), ceil(), fmod();
2: extern double frexp(), ldexp(), modf();
3: extern double sqrt(), hypot(), atof(), strtod();
4: extern double sin(), cos(), tan(), asin(), acos(), atan(), atan2();
5: extern double exp(), log(), log10(), pow();
6: extern double sinh(), cosh(), tanh();
7: extern double gamma();
8: extern double j0(), j1(), jn(), y0(), y1(), yn();
9: extern double frand();
10: extern double erf(), erfc();
11: extern struct _iobuf {
12: int _cnt;
13: unsigned char *_ptr;
14: unsigned char *_base;
15: short _flag;
16: char _file;
17: } _iob[120];
18: struct _iobuf *fopen();
19: struct _iobuf *fdopen();
20: struct _iobuf *freopen();
21: struct _iobuf *popen();
22: long ftell();
23: char *fgets();
24: typedef double real;
25: realfft ( log2n, x, y, sign )
26: int sign ;
27: real log2n, x[], y[] ;
28: {
29: static int prevn = 0 ;
30: static double *cost, *sint ;
31: int s[13], u[13] ;
32: int j1, j2, j3, j4, jt, n, m4 ;
33: int a, b, c, d, e, f, g, h, i, j, k, l, m ;
34: double arg, c1, c2, c3, s1, s2, s3, tpim4 ;
35: double r1, r2, r3, r4, r5, r6, r7, r8 ;
36: real *w, *z, t ;
37: w = &x[0] ;
38: z = &x[1] ;
39: n = pow(2.,log2n -1.) + 0.2 ;
40:
41: if ( n / 4 != prevn ) {
42: if ( prevn != 0 ) {
43: free ((char *)cost ) ;
44: free ((char *)sint ) ;
45: }
46: prevn = n / 4 ;
47:
48: if ((( cost = (double *)calloc((unsigned)prevn,sizeof(*cost)))
49: == (double *)0 ) ||
50: (( sint = (double *)calloc((unsigned)prevn,sizeof(*sint)))
51: == (double *)0 ) ) {
52: fprintf((&_iob[2]), " alloc error in realfft\n");
53: exit(1);
54: }
55: tpim4 = 6.283185307179586476925287 / ( 2 * n ) ;
56: for ( j = 1 ; j < prevn ; j++ ) {
57: arg = tpim4 * j ;
58: cost[j] = cos ( arg ) ;
59: sint[j] = sin ( arg ) ;
60: }
61: cost[0] = 1. ;
62: sint[0] = 0. ;
63: }
64: if (log2n > 2.1 ) {
65: for (m = n/4; m >= 1; m /= 4) {
66: m4 = 4 * m ;
67: e = n / m4 ;
68: for (j = 0; j < m; j++ ) {
69: f = e * j ;
70: c1 = cost[f] * cost[f] * 2. - 1. ;
71: s1 = cost[f] * sint[f] * 2. ;
72: c2 = c1 * c1 * 2. - 1. ;
73: s2 = c1 * s1 * 2. ;
74: c3 = c2 * c1 - s2 * s1 ;
75: s3 = c2 * s1 + s2 * c1 ;
76: for (i = m4; i <= n; i += m4) {
77: j1 = (i + j - m4) * 2 ;
78: a = m * 2 ;
79: j2 = j1 + a ;
80: j3 = j2 + a ;
81: j4 = j3 + a ;
82: r1 = w[j1] + w[j3] ;
83: r2 = w[j1] - w[j3] ;
84: r3 = z[j1] + z[j3] ;
85: r4 = z[j1] - z[j3] ;
86: r5 = w[j2] + w[j4] ;
87: r6 = w[j2] - w[j4] ;
88: r7 = z[j2] + z[j4] ;
89: r8 = z[j2] - z[j4] ;
90: w[j1] = r1 + r5 ;
91: z[j1] = r3 + r7 ;
92: if (arg != 0. ) {
93: w[j3] = (r2+r8)*c1 + (r4-r6)*s1 ;
94: z[j3] = (r4-r6)*c1 - (r2+r8)*s1 ;
95: w[j2] = (r1-r5)*c2 + (r3-r7)*s2 ;
96: z[j2] = (r3-r7)*c2 - (r1-r5)*s2 ;
97: w[j4] = (r2-r8)*c3 + (r4+r6)*s3 ;
98: z[j4] = (r4+r6)*c3 - (r2-r8)*s3 ;
99: }
100: else {
101: w[j3] = r2+r8 ;
102: z[j3] = r4-r6 ;
103: w[j2] = r1-r5 ;
104: z[j2] = r3-r7 ;
105: w[j4] = r2-r8 ;
106: z[j4] = r4+r6 ;
107: }
108: }
109: }
110: }
111: }
112:
113: if (((int)(log2n+.2)-1) % 2) {
114: m = n * 2 ;
115: for (j=2, i = 0; i < m ; i += 4, j += 4 ) {
116: r1 = w[i]+w[j] ;
117: r2 = w[i]-w[j] ;
118: r3 = z[i]+z[j] ;
119: r4 = z[i]-z[j] ;
120: w[i] = r1 ;
121: z[i] = r3 ;
122: w[j] = r2 ;
123: z[j] = r4 ;
124: }
125: }
126:
127: s[12] = n ;
128: u[12] = n * 2 ;
129: for (j = 11; j > 0; --j) {
130: s[j] = (s[j+1] > 2) ? s[j+1] / 2 : 2 ;
131: u[j] = s[j+1] ;
132: }
133: u[0] = s[1] ;
134: jt = -2 ;
135: for (a = 0; a < u[0]; a += 2 )
136: for (b = a; b < u[1]; b += s[1])
137: for (c = b; c < u[2]; c += s[2])
138: for (d = c; d < u[3]; d += s[3])
139: for (e = d; e < u[4]; e += s[4])
140: for (f = e; f < u[5]; f += s[5])
141: for (g = f; g < u[6]; g += s[6])
142: for (h = g; h < u[7]; h += s[7])
143: for (i = h; i < u[8]; i += s[8])
144: for (j = i; j < u[9]; j += s[9])
145: for (k = j; k < u[10]; k += s[10])
146: for (l = k; l < u[11]; l += s[11])
147: for (m = l; m < u[12]; m += s[12]) {
148: jt += 2 ;
149: if ( jt > m) {
150: t = w[jt] ;
151: w[jt] = w[m] ;
152: w[m] = t ;
153: t = z[jt] ;
154: z[jt] = z[m] ;
155: z[m] = t ;
156: }
157: }
158:
159: n *= 2 ;
160: tpim4 = w[0] - z[0] ;
161: x[0] = w[0] + z[0] ;
162: y[0] = 0. ;
163: m = n / 8 ;
164: for (l=n/2-1, i=n-2, j=2, g=n/4-1, f=n/2-2, e=n/2+2, d=n/4+1, k=1;
165: k < m ;
166: --l, i -= 2, j += 2, --g, f -= 2, e += 2, d++, k++ ) {
167:
168: r1 = z[j] + z[i] ;
169: r2 = w[i] - w[j] ;
170: r3 = cost[ k ] ;
171: r4 = sint[ k ] ;
172: r5 = r3 * r2 - r4 * r1 ;
173: y[k] = .5 * ((z[j] - z[i]) + r5) ;
174: y[l] = .5 * ((z[i] - z[j]) + r5) ;
175: r5 = r3 * r1 + r4 * r2 ;
176: r6 = w[j] + w[i] ;
177: x[k] = .5 * (r6 + r5) ;
178: x[i] = .5 * (r6 - r5) ;
179:
180: r1 = z[f] + z[e] ;
181: r2 = w[e] - w[f] ;
182: r5 = r4 * r2 - r3 * r1 ;
183: y[g] = .5 * ((z[f] - z[e]) + r5) ;
184: y[d] = .5 * ((z[e] - z[f]) + r5) ;
185: r5 = r4 * r1 + r3 * r2 ;
186: r6 = w[f] + w[e] ;
187: x[f] = .5 * (r6 + r5) ;
188: x[e] = .5 * (r6 - r5) ;
189: }
190:
191: r1 = ( z[j] + z[i] ) * 0.70710678118654752440 ;
192: r2 = ( w[i] - w[j] ) * 0.70710678118654752440 ;
193: r5 = r2 - r1 ;
194: y[k] = .5 * ((z[j] - z[i]) + r5) ;
195: y[l] = .5 * ((z[i] - z[j]) + r5) ;
196: r5 = r1 + r2 ;
197: r6 = w[j] + w[i] ;
198: x[k] = .5 * (r6 + r5) ;
199: x[i] = .5 * (r6 - r5) ;
200:
201: m = n / 2 ;
202: for (i=n/4+2, l=n/8+1; l< m ; i += 2, l++)
203: x[l] = x[i] ;
204:
205: y[n/4] = -z[n/2] ;
206: x[n/2] = tpim4 ;
207:
208: if ( sign <= 0 ) {
209: r1 = 1. / (double) n ;
210: r2 = - r1 ;
211: m = n / 2 ;
212: for (i=0; i<m ; i++) {
213: x[i] *= r1 ;
214: y[i] *= r2 ;
215: }
216: x[n/2] *= r1 ;
217: }
218: return ;
219: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.