1 #include <math.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <string.h>
5 #include <assert.h>
6
7 #include "sphinxbase/f2c.h"
8
9 #ifdef _MSC_VER
10 #pragma warning (disable: 4244)
11 #endif
12
13
14 extern void
s_wsfe(cilist * f)15 s_wsfe(cilist * f)
16 {;
17 }
18 extern void
e_wsfe(void)19 e_wsfe(void)
20 {;
21 }
22 extern void
do_fio(integer * c,char * s,ftnlen l)23 do_fio(integer * c, char *s, ftnlen l)
24 {;
25 }
26
27 /* You'll want this if you redo the *_lite.c files with the -C option
28 * to f2c for checking array subscripts. (It's not suggested you do that
29 * for production use, of course.) */
30 extern int
s_rnge(char * var,int index,char * routine,int lineno)31 s_rnge(char *var, int index, char *routine, int lineno)
32 {
33 fprintf(stderr,
34 "array index out-of-bounds for %s[%d] in routine %s:%d\n", var,
35 index, routine, lineno);
36 fflush(stderr);
37 assert(2+2 == 5);
38 return 0;
39 }
40
41
42 #ifdef KR_headers
43 extern double sqrt();
44 float
f__cabs(real,imag)45 f__cabs(real, imag)
46 float real, imag;
47 #else
48 #undef abs
49
50 float
51 f__cabs(float real, float imag)
52 #endif
53 {
54 float temp;
55
56 if (real < 0)
57 real = -real;
58 if (imag < 0)
59 imag = -imag;
60 if (imag > real) {
61 temp = real;
62 real = imag;
63 imag = temp;
64 }
65 if ((imag + real) == real)
66 return ((float) real);
67
68 temp = imag / real;
69 temp = real * sqrt(1.0 + temp * temp); /*overflow!! */
70 return (temp);
71 }
72
73
74 VOID
75 #ifdef KR_headers
s_cnjg(r,z)76 s_cnjg(r, z)
77 complex *r, *z;
78 #else
79 s_cnjg(complex * r, complex * z)
80 #endif
81 {
82 r->r = z->r;
83 r->i = -z->i;
84 }
85
86
87 #ifdef KR_headers
88 float
r_imag(z)89 r_imag(z)
90 complex *z;
91 #else
92 float
93 r_imag(complex * z)
94 #endif
95 {
96 return (z->i);
97 }
98
99
100 #define log10e 0.43429448190325182765
101
102 #ifdef KR_headers
103 double log();
104 float
r_lg10(x)105 r_lg10(x)
106 real *x;
107 #else
108 #undef abs
109
110 float
111 r_lg10(real * x)
112 #endif
113 {
114 return (log10e * log(*x));
115 }
116
117
118 #ifdef KR_headers
119 float
r_sign(a,b)120 r_sign(a, b)
121 real *a, *b;
122 #else
123 float
124 r_sign(real * a, real * b)
125 #endif
126 {
127 float x;
128 x = (*a >= 0 ? *a : -*a);
129 return (*b >= 0 ? x : -x);
130 }
131
132
133 #ifdef KR_headers
134 double floor();
135 integer
i_dnnt(x)136 i_dnnt(x)
137 real *x;
138 #else
139 #undef abs
140
141 integer
142 i_dnnt(real * x)
143 #endif
144 {
145 return ((*x) >= 0 ? floor(*x + .5) : -floor(.5 - *x));
146 }
147
148
149 #ifdef KR_headers
150 double pow();
151 double
pow_dd(ap,bp)152 pow_dd(ap, bp)
153 doublereal *ap, *bp;
154 #else
155 #undef abs
156
157 double
158 pow_dd(doublereal * ap, doublereal * bp)
159 #endif
160 {
161 return (pow(*ap, *bp));
162 }
163
164
165 #ifdef KR_headers
166 float
pow_ri(ap,bp)167 pow_ri(ap, bp)
168 real *ap;
169 integer *bp;
170 #else
171 float
172 pow_ri(real * ap, integer * bp)
173 #endif
174 {
175 float pow, x;
176 integer n;
177 unsigned long u;
178
179 pow = 1;
180 x = *ap;
181 n = *bp;
182
183 if (n != 0) {
184 if (n < 0) {
185 n = -n;
186 x = 1 / x;
187 }
188 for (u = n;;) {
189 if (u & 01)
190 pow *= x;
191 if (u >>= 1)
192 x *= x;
193 else
194 break;
195 }
196 }
197 return (pow);
198 }
199
200 /* Unless compiled with -DNO_OVERWRITE, this variant of s_cat allows the
201 * target of a concatenation to appear on its right-hand side (contrary
202 * to the Fortran 77 Standard, but in accordance with Fortran 90).
203 */
204 #define NO_OVERWRITE
205
206
207 #ifndef NO_OVERWRITE
208
209 #undef abs
210 #ifdef KR_headers
211 extern char *F77_aloc();
212 extern void free();
213 extern void exit_();
214 #else
215
216 extern char *F77_aloc(ftnlen, char *);
217 #endif
218
219 #endif /* NO_OVERWRITE */
220
221 VOID
222 #ifdef KR_headers
s_cat(lp,rpp,rnp,np,ll)223 s_cat(lp, rpp, rnp, np, ll)
224 char *lp, *rpp[];
225 ftnlen rnp[], *np, ll;
226 #else
227 s_cat(char *lp, char *rpp[], ftnlen rnp[], ftnlen * np, ftnlen ll)
228 #endif
229 {
230 ftnlen i, nc;
231 char *rp;
232 ftnlen n = *np;
233 #ifndef NO_OVERWRITE
234 ftnlen L, m;
235 char *lp0, *lp1;
236
237 lp0 = 0;
238 lp1 = lp;
239 L = ll;
240 i = 0;
241 while (i < n) {
242 rp = rpp[i];
243 m = rnp[i++];
244 if (rp >= lp1 || rp + m <= lp) {
245 if ((L -= m) <= 0) {
246 n = i;
247 break;
248 }
249 lp1 += m;
250 continue;
251 }
252 lp0 = lp;
253 lp = lp1 = F77_aloc(L = ll, "s_cat");
254 break;
255 }
256 lp1 = lp;
257 #endif /* NO_OVERWRITE */
258 for (i = 0; i < n; ++i) {
259 nc = ll;
260 if (rnp[i] < nc)
261 nc = rnp[i];
262 ll -= nc;
263 rp = rpp[i];
264 while (--nc >= 0)
265 *lp++ = *rp++;
266 }
267 while (--ll >= 0)
268 *lp++ = ' ';
269 #ifndef NO_OVERWRITE
270 if (lp0) {
271 memmove(lp0, lp1, L);
272 free(lp1);
273 }
274 #endif
275 }
276
277
278 /* compare two strings */
279
280 #ifdef KR_headers
281 integer
s_cmp(a0,b0,la,lb)282 s_cmp(a0, b0, la, lb)
283 char *a0, *b0;
284 ftnlen la, lb;
285 #else
286 integer
287 s_cmp(char *a0, char *b0, ftnlen la, ftnlen lb)
288 #endif
289 {
290 register unsigned char *a, *aend, *b, *bend;
291 a = (unsigned char *) a0;
292 b = (unsigned char *) b0;
293 aend = a + la;
294 bend = b + lb;
295
296 if (la <= lb) {
297 while (a < aend)
298 if (*a != *b)
299 return (*a - *b);
300 else {
301 ++a;
302 ++b;
303 }
304
305 while (b < bend)
306 if (*b != ' ')
307 return (' ' - *b);
308 else
309 ++b;
310 }
311
312 else {
313 while (b < bend)
314 if (*a == *b) {
315 ++a;
316 ++b;
317 }
318 else
319 return (*a - *b);
320 while (a < aend)
321 if (*a != ' ')
322 return (*a - ' ');
323 else
324 ++a;
325 }
326 return (0);
327 }
328
329 /* Unless compiled with -DNO_OVERWRITE, this variant of s_copy allows the
330 * target of an assignment to appear on its right-hand side (contrary
331 * to the Fortran 77 Standard, but in accordance with Fortran 90),
332 * as in a(2:5) = a(4:7) .
333 */
334
335
336
337 /* assign strings: a = b */
338
339 #ifdef KR_headers
340 VOID
s_copy(a,b,la,lb)341 s_copy(a, b, la, lb)
342 register char *a, *b;
343 ftnlen la, lb;
344 #else
345 void
346 s_copy(register char *a, register char *b, ftnlen la, ftnlen lb)
347 #endif
348 {
349 register char *aend, *bend;
350
351 aend = a + la;
352
353 if (la <= lb)
354 #ifndef NO_OVERWRITE
355 if (a <= b || a >= b + la)
356 #endif
357 while (a < aend)
358 *a++ = *b++;
359 #ifndef NO_OVERWRITE
360 else
361 for (b += la; a < aend;)
362 *--aend = *--b;
363 #endif
364
365 else {
366 bend = b + lb;
367 #ifndef NO_OVERWRITE
368 if (a <= b || a >= bend)
369 #endif
370 while (b < bend)
371 *a++ = *b++;
372 #ifndef NO_OVERWRITE
373 else {
374 a += lb;
375 while (b < bend)
376 *--a = *--bend;
377 a += lb;
378 }
379 #endif
380 while (a < aend)
381 *a++ = ' ';
382 }
383 }
384
385
386 #ifdef KR_headers
387 float f__cabs();
388 float
389 z_abs(z)
390 complex *z;
391 #else
392 float f__cabs(float, float);
393 float
z_abs(complex * z)394 z_abs(complex * z)
395 #endif
396 {
397 return (f__cabs(z->r, z->i));
398 }
399
400
401 #ifdef KR_headers
402 extern void sig_die();
403 VOID
z_div(c,a,b)404 z_div(c, a, b)
405 complex *a, *b, *c;
406 #else
407 extern void sig_die(char *, int);
408 void
409 z_div(complex * c, complex * a, complex * b)
410 #endif
411 {
412 float ratio, den;
413 float abr, abi;
414
415 if ((abr = b->r) < 0.)
416 abr = -abr;
417 if ((abi = b->i) < 0.)
418 abi = -abi;
419 if (abr <= abi) {
420 /*Let IEEE Infinties handle this ;( */
421 /*if(abi == 0)
422 sig_die("complex division by zero", 1); */
423 ratio = b->r / b->i;
424 den = b->i * (1 + ratio * ratio);
425 c->r = (a->r * ratio + a->i) / den;
426 c->i = (a->i * ratio - a->r) / den;
427 }
428
429 else {
430 ratio = b->i / b->r;
431 den = b->r * (1 + ratio * ratio);
432 c->r = (a->r + a->i * ratio) / den;
433 c->i = (a->i - a->r * ratio) / den;
434 }
435
436 }
437
438
439 #ifdef KR_headers
440 double sqrt();
441 double f__cabs();
442 VOID
z_sqrt(r,z)443 z_sqrt(r, z)
444 complex *r, *z;
445 #else
446 #undef abs
447
448 extern float f__cabs(float, float);
449 void
450 z_sqrt(complex * r, complex * z)
451 #endif
452 {
453 float mag;
454
455 if ((mag = f__cabs(z->r, z->i)) == 0.)
456 r->r = r->i = 0.;
457 else if (z->r > 0) {
458 r->r = sqrt(0.5 * (mag + z->r));
459 r->i = z->i / r->r / 2;
460 }
461 else {
462 r->i = sqrt(0.5 * (mag - z->r));
463 if (z->i < 0)
464 r->i = -r->i;
465 r->r = z->i / r->i / 2;
466 }
467 }
468
469 #ifdef __cplusplus
470 extern "C" {
471 #endif
472
473 #ifdef KR_headers
pow_ii(ap,bp)474 integer pow_ii(ap, bp) integer *ap, *bp;
475 #else
476 integer pow_ii(integer * ap, integer * bp)
477 #endif
478 {
479 integer pow, x, n;
480 unsigned long u;
481
482 x = *ap;
483 n = *bp;
484
485 if (n <= 0) {
486 if (n == 0 || x == 1)
487 return 1;
488 if (x != -1)
489 return x != 0 ? 1 / x : 0;
490 n = -n;
491 } u = n;
492 for (pow = 1;;) {
493 if (u & 01)
494 pow *= x;
495 if (u >>= 1)
496 x *= x;
497 else
498 break;
499 }
500 return (pow);
501 }
502 #ifdef __cplusplus
503 }
504 #endif
505
506 #ifdef KR_headers
507 extern void f_exit();
508 VOID
s_stop(s,n)509 s_stop(s, n)
510 char *s;
511 ftnlen n;
512 #else
513 #undef abs
514 #undef min
515 #undef max
516 #ifdef __cplusplus
517 extern "C" {
518 #endif
519 #ifdef __cplusplus
520 extern "C" {
521 #endif
522 void f_exit(void);
523
524 int s_stop(char *s, ftnlen n)
525 #endif
526 {
527 int i;
528
529 if (n > 0) {
530 fprintf(stderr, "STOP ");
531 for (i = 0; i < n; ++i)
532 putc(*s++, stderr);
533 fprintf(stderr, " statement executed\n");
534 }
535 #ifdef NO_ONEXIT
536 f_exit();
537 #endif
538 exit(0);
539
540 /* We cannot avoid (useless) compiler diagnostics here: */
541 /* some compilers complain if there is no return statement, */
542 /* and others complain that this one cannot be reached. */
543
544 return 0; /* NOT REACHED */
545 }
546 #ifdef __cplusplus
547 }
548 #endif
549 #ifdef __cplusplus
550 }
551 #endif
552