SphinxBase 0.6
src/libsphinxbase/util/f2c_lite.c
00001 #include <math.h>
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include <string.h>
00005 #include <assert.h>
00006 
00007 #include "sphinxbase/f2c.h"
00008 
00009 #ifdef _MSC_VER
00010 #pragma warning (disable: 4244)
00011 #endif
00012 
00013 
00014 extern void
00015 s_wsfe(cilist * f)
00016 {;
00017 }
00018 extern void
00019 e_wsfe(void)
00020 {;
00021 }
00022 extern void
00023 do_fio(integer * c, char *s, ftnlen l)
00024 {;
00025 }
00026 
00027 /* You'll want this if you redo the *_lite.c files with the -C option
00028  * to f2c for checking array subscripts. (It's not suggested you do that
00029  * for production use, of course.) */
00030 extern int
00031 s_rnge(char *var, int index, char *routine, int lineno)
00032 {
00033     fprintf(stderr,
00034             "array index out-of-bounds for %s[%d] in routine %s:%d\n", var,
00035             index, routine, lineno);
00036     fflush(stderr);
00037         assert(2+2 == 5);
00038         return 0;
00039 }
00040 
00041 
00042 #ifdef KR_headers
00043 extern double sqrt();
00044 float
00045 f__cabs(real, imag)
00046 float real, imag;
00047 #else
00048 #undef abs
00049 
00050 float
00051 f__cabs(float real, float imag)
00052 #endif
00053 {
00054     float temp;
00055 
00056     if (real < 0)
00057         real = -real;
00058     if (imag < 0)
00059         imag = -imag;
00060     if (imag > real) {
00061         temp = real;
00062         real = imag;
00063         imag = temp;
00064     }
00065     if ((imag + real) == real)
00066         return ((float) real);
00067 
00068     temp = imag / real;
00069     temp = real * sqrt(1.0 + temp * temp);      /*overflow!! */
00070     return (temp);
00071 }
00072 
00073 
00074 VOID
00075 #ifdef KR_headers
00076 s_cnjg(r, z)
00077 complex *r, *z;
00078 #else
00079 s_cnjg(complex * r, complex * z)
00080 #endif
00081 {
00082     r->r = z->r;
00083     r->i = -z->i;
00084 }
00085 
00086 
00087 #ifdef KR_headers
00088 float
00089 r_imag(z)
00090 complex *z;
00091 #else
00092 float
00093 r_imag(complex * z)
00094 #endif
00095 {
00096     return (z->i);
00097 }
00098 
00099 
00100 #define log10e 0.43429448190325182765
00101 
00102 #ifdef KR_headers
00103 double log();
00104 float
00105 r_lg10(x)
00106 real *x;
00107 #else
00108 #undef abs
00109 
00110 float
00111 r_lg10(real * x)
00112 #endif
00113 {
00114     return (log10e * log(*x));
00115 }
00116 
00117 
00118 #ifdef KR_headers
00119 float
00120 r_sign(a, b)
00121 real *a, *b;
00122 #else
00123 float
00124 r_sign(real * a, real * b)
00125 #endif
00126 {
00127     float x;
00128     x = (*a >= 0 ? *a : -*a);
00129     return (*b >= 0 ? x : -x);
00130 }
00131 
00132 
00133 #ifdef KR_headers
00134 double floor();
00135 integer
00136 i_dnnt(x)
00137 real *x;
00138 #else
00139 #undef abs
00140 
00141 integer
00142 i_dnnt(real * x)
00143 #endif
00144 {
00145     return ((*x) >= 0 ? floor(*x + .5) : -floor(.5 - *x));
00146 }
00147 
00148 
00149 #ifdef KR_headers
00150 double pow();
00151 double
00152 pow_dd(ap, bp)
00153 doublereal *ap, *bp;
00154 #else
00155 #undef abs
00156 
00157 double
00158 pow_dd(doublereal * ap, doublereal * bp)
00159 #endif
00160 {
00161     return (pow(*ap, *bp));
00162 }
00163 
00164 
00165 #ifdef KR_headers
00166 float
00167 pow_ri(ap, bp)
00168 real *ap;
00169 integer *bp;
00170 #else
00171 float
00172 pow_ri(real * ap, integer * bp)
00173 #endif
00174 {
00175     float pow, x;
00176     integer n;
00177     unsigned long u;
00178 
00179     pow = 1;
00180     x = *ap;
00181     n = *bp;
00182 
00183     if (n != 0) {
00184         if (n < 0) {
00185             n = -n;
00186             x = 1 / x;
00187         }
00188         for (u = n;;) {
00189             if (u & 01)
00190                 pow *= x;
00191             if (u >>= 1)
00192                 x *= x;
00193             else
00194                 break;
00195         }
00196     }
00197     return (pow);
00198 }
00199 
00200 /* Unless compiled with -DNO_OVERWRITE, this variant of s_cat allows the
00201  * target of a concatenation to appear on its right-hand side (contrary
00202  * to the Fortran 77 Standard, but in accordance with Fortran 90).
00203  */
00204 #define NO_OVERWRITE
00205 
00206 
00207 #ifndef NO_OVERWRITE
00208 
00209 #undef abs
00210 #ifdef KR_headers
00211 extern char *F77_aloc();
00212 extern void free();
00213 extern void exit_();
00214 #else
00215 
00216 extern char *F77_aloc(ftnlen, char *);
00217 #endif
00218 
00219 #endif                          /* NO_OVERWRITE */
00220 
00221 VOID
00222 #ifdef KR_headers
00223 s_cat(lp, rpp, rnp, np, ll)
00224 char *lp, *rpp[];
00225 ftnlen rnp[], *np, ll;
00226 #else
00227 s_cat(char *lp, char *rpp[], ftnlen rnp[], ftnlen * np, ftnlen ll)
00228 #endif
00229 {
00230     ftnlen i, nc;
00231     char *rp;
00232     ftnlen n = *np;
00233 #ifndef NO_OVERWRITE
00234     ftnlen L, m;
00235     char *lp0, *lp1;
00236 
00237     lp0 = 0;
00238     lp1 = lp;
00239     L = ll;
00240     i = 0;
00241     while (i < n) {
00242         rp = rpp[i];
00243         m = rnp[i++];
00244         if (rp >= lp1 || rp + m <= lp) {
00245             if ((L -= m) <= 0) {
00246                 n = i;
00247                 break;
00248             }
00249             lp1 += m;
00250             continue;
00251         }
00252         lp0 = lp;
00253         lp = lp1 = F77_aloc(L = ll, "s_cat");
00254         break;
00255     }
00256     lp1 = lp;
00257 #endif                          /* NO_OVERWRITE */
00258     for (i = 0; i < n; ++i) {
00259         nc = ll;
00260         if (rnp[i] < nc)
00261             nc = rnp[i];
00262         ll -= nc;
00263         rp = rpp[i];
00264         while (--nc >= 0)
00265             *lp++ = *rp++;
00266     }
00267     while (--ll >= 0)
00268         *lp++ = ' ';
00269 #ifndef NO_OVERWRITE
00270     if (lp0) {
00271         memmove(lp0, lp1, L);
00272         free(lp1);
00273     }
00274 #endif
00275 }
00276 
00277 
00278 /* compare two strings */
00279 
00280 #ifdef KR_headers
00281 integer
00282 s_cmp(a0, b0, la, lb)
00283 char *a0, *b0;
00284 ftnlen la, lb;
00285 #else
00286 integer
00287 s_cmp(char *a0, char *b0, ftnlen la, ftnlen lb)
00288 #endif
00289 {
00290     register unsigned char *a, *aend, *b, *bend;
00291     a = (unsigned char *) a0;
00292     b = (unsigned char *) b0;
00293     aend = a + la;
00294     bend = b + lb;
00295 
00296     if (la <= lb) {
00297         while (a < aend)
00298             if (*a != *b)
00299                 return (*a - *b);
00300             else {
00301                 ++a;
00302                 ++b;
00303             }
00304 
00305         while (b < bend)
00306             if (*b != ' ')
00307                 return (' ' - *b);
00308             else
00309                 ++b;
00310     }
00311 
00312     else {
00313         while (b < bend)
00314             if (*a == *b) {
00315                 ++a;
00316                 ++b;
00317             }
00318             else
00319                 return (*a - *b);
00320         while (a < aend)
00321             if (*a != ' ')
00322                 return (*a - ' ');
00323             else
00324                 ++a;
00325     }
00326     return (0);
00327 }
00328 
00329 /* Unless compiled with -DNO_OVERWRITE, this variant of s_copy allows the
00330  * target of an assignment to appear on its right-hand side (contrary
00331  * to the Fortran 77 Standard, but in accordance with Fortran 90),
00332  * as in  a(2:5) = a(4:7) .
00333  */
00334 
00335 
00336 
00337 /* assign strings:  a = b */
00338 
00339 #ifdef KR_headers
00340 VOID
00341 s_copy(a, b, la, lb)
00342 register char *a, *b;
00343 ftnlen la, lb;
00344 #else
00345 void
00346 s_copy(register char *a, register char *b, ftnlen la, ftnlen lb)
00347 #endif
00348 {
00349     register char *aend, *bend;
00350 
00351     aend = a + la;
00352 
00353     if (la <= lb)
00354 #ifndef NO_OVERWRITE
00355         if (a <= b || a >= b + la)
00356 #endif
00357             while (a < aend)
00358                 *a++ = *b++;
00359 #ifndef NO_OVERWRITE
00360         else
00361             for (b += la; a < aend;)
00362                 *--aend = *--b;
00363 #endif
00364 
00365     else {
00366         bend = b + lb;
00367 #ifndef NO_OVERWRITE
00368         if (a <= b || a >= bend)
00369 #endif
00370             while (b < bend)
00371                 *a++ = *b++;
00372 #ifndef NO_OVERWRITE
00373         else {
00374             a += lb;
00375             while (b < bend)
00376                 *--a = *--bend;
00377             a += lb;
00378         }
00379 #endif
00380         while (a < aend)
00381             *a++ = ' ';
00382     }
00383 }
00384 
00385 
00386 #ifdef KR_headers
00387 float f__cabs();
00388 float
00389 z_abs(z)
00390 complex *z;
00391 #else
00392 float f__cabs(float, float);
00393 float
00394 z_abs(complex * z)
00395 #endif
00396 {
00397     return (f__cabs(z->r, z->i));
00398 }
00399 
00400 
00401 #ifdef KR_headers
00402 extern void sig_die();
00403 VOID
00404 z_div(c, a, b)
00405 complex *a, *b, *c;
00406 #else
00407 extern void sig_die(char *, int);
00408 void
00409 z_div(complex * c, complex * a, complex * b)
00410 #endif
00411 {
00412     float ratio, den;
00413     float abr, abi;
00414 
00415     if ((abr = b->r) < 0.)
00416         abr = -abr;
00417     if ((abi = b->i) < 0.)
00418         abi = -abi;
00419     if (abr <= abi) {
00420         /*Let IEEE Infinties handle this ;( */
00421         /*if(abi == 0)
00422            sig_die("complex division by zero", 1); */
00423         ratio = b->r / b->i;
00424         den = b->i * (1 + ratio * ratio);
00425         c->r = (a->r * ratio + a->i) / den;
00426         c->i = (a->i * ratio - a->r) / den;
00427     }
00428 
00429     else {
00430         ratio = b->i / b->r;
00431         den = b->r * (1 + ratio * ratio);
00432         c->r = (a->r + a->i * ratio) / den;
00433         c->i = (a->i - a->r * ratio) / den;
00434     }
00435 
00436 }
00437 
00438 
00439 #ifdef KR_headers
00440 double sqrt();
00441 double f__cabs();
00442 VOID
00443 z_sqrt(r, z)
00444 complex *r, *z;
00445 #else
00446 #undef abs
00447 
00448 extern float f__cabs(float, float);
00449 void
00450 z_sqrt(complex * r, complex * z)
00451 #endif
00452 {
00453     float mag;
00454 
00455     if ((mag = f__cabs(z->r, z->i)) == 0.)
00456         r->r = r->i = 0.;
00457     else if (z->r > 0) {
00458         r->r = sqrt(0.5 * (mag + z->r));
00459         r->i = z->i / r->r / 2;
00460     }
00461     else {
00462         r->i = sqrt(0.5 * (mag - z->r));
00463         if (z->i < 0)
00464             r->i = -r->i;
00465         r->r = z->i / r->i / 2;
00466     }
00467 }
00468 
00469 #ifdef __cplusplus
00470 extern "C" {
00471 #endif
00472 
00473 #ifdef KR_headers
00474     integer pow_ii(ap, bp) integer *ap, *bp;
00475 #else
00476     integer pow_ii(integer * ap, integer * bp)
00477 #endif
00478     {
00479         integer pow, x, n;
00480         unsigned long u;
00481 
00482          x = *ap;
00483          n = *bp;
00484 
00485         if (n <= 0) {
00486             if (n == 0 || x == 1)
00487                 return 1;
00488             if (x != -1)
00489                 return x != 0 ? 1 / x : 0;
00490             n = -n;
00491         } u = n;
00492         for (pow = 1;;) {
00493             if (u & 01)
00494                 pow *= x;
00495             if (u >>= 1)
00496                 x *= x;
00497             else
00498                 break;
00499         }
00500         return (pow);
00501     }
00502 #ifdef __cplusplus
00503 }
00504 #endif
00505 
00506 #ifdef KR_headers
00507 extern void f_exit();
00508 VOID
00509 s_stop(s, n)
00510 char *s;
00511 ftnlen n;
00512 #else
00513 #undef abs
00514 #undef min
00515 #undef max
00516 #ifdef __cplusplus
00517 extern "C" {
00518 #endif
00519 #ifdef __cplusplus
00520     extern "C" {
00521 #endif
00522         void f_exit(void);
00523 
00524         int s_stop(char *s, ftnlen n)
00525 #endif
00526         {
00527             int i;
00528 
00529             if (n > 0) {
00530                 fprintf(stderr, "STOP ");
00531                 for (i = 0; i < n; ++i)
00532                     putc(*s++, stderr);
00533                 fprintf(stderr, " statement executed\n");
00534             }
00535 #ifdef NO_ONEXIT
00536             f_exit();
00537 #endif
00538              exit(0);
00539 
00540 /* We cannot avoid (useless) compiler diagnostics here:         */
00541 /* some compilers complain if there is no return statement,     */
00542 /* and others complain that this one cannot be reached.         */
00543 
00544              return 0;          /* NOT REACHED */
00545         }
00546 #ifdef __cplusplus
00547     }
00548 #endif
00549 #ifdef __cplusplus
00550 }
00551 #endif