Return to 3.i CVS log

Up to [Research Unix] / researchv10no / cmd / ccom / vax / tests

Annotation of researchv10no/cmd/ccom/vax/tests/3.i, revision 1.1

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.