Return to rhmath.c CVS log

Up to [Research Unix] / researchv10no / cmd / PDP11 / fpp / aux

Annotation of researchv10no/cmd/PDP11/fpp/aux/rhmath.c, revision 1.1

1.1     ! root        1: /* These are unguaranteed math routines lifted from NORGEN "K" code */
        !             2: 
        !             3: static float alogtab[] = {
        !             4:        -.00645354,
        !             5:        .03608849,
        !             6:        -.09532938,
        !             7:        .16765407,
        !             8:        -.24073380,
        !             9:        .33179902,
        !            10:        -.49987412,
        !            11:        .999964239,
        !            12: };
        !            13: 
        !            14: extern long _fsign2;
        !            15: float exp(exparg)
        !            16: float exparg;
        !            17: {
        !            18:        float ffr(); 
        !            19:        long fex();
        !            20:        float expd,expz;
        !            21:        float two();
        !            22:        long expc;
        !            23:        expz = exparg*1.44269504;
        !            24:        /* 1/LN(2) */
        !            25:        expc=expz;
        !            26:        expd=expc-expz;
        !            27:        if(expz>=0.0){
        !            28:                expd = expd+1.0;
        !            29:                expc = expc+1;
        !            30:        }
        !            31:        expz=expd*expd;
        !            32:        return(two(expc)*(1.-(2.*expd)/(expd+.034657359*(expz)+9.9545948-
        !            33:            617.97227/(expz+87.417497))));
        !            34: }
        !            35: float atan(atanx)
        !            36: float atanx;
        !            37: {
        !            38:        if(atanx<0.0)   return(-atan(-atanx));
        !            39:        if(atanx>1.0)   return(1.57079632-atan(1./atanx));
        !            40:        if(atanx>2.67949192e-1) 
        !            41:                return(5.23598775e-1+atan((7.32050808e-1*atanx-1.+atanx)/
        !            42:                            (1.732050808+atanx)));
        !            43:        return(atanx*(.60310579-.05160454*atanx*atanx+.55913709/(atanx*atanx+
        !            44:                            1.4087812)));
        !            45: }
        !            46: float atan2(atan2x1,atan2x2)
        !            47: float atan2x1,atan2x2;
        !            48: {
        !            49:        float atan2z;
        !            50:        if(atan2x1<0.0)         return(-atan2(-atan2x1,atan2x2));
        !            51:        if(!atan2x2)            return(1.5707963);
        !            52:        atan2z = (atan2x1/atan2x2);
        !            53:        if(atan2z<0.)           atan2z= -atan2z;
        !            54:        if(atan2z>1.67777e7)    return(1.5707963);
        !            55:        if(atan2x2>0.0)         return(atan(atan2z));
        !            56:        return(3.14159265-atan(atan2z));
        !            57: }
        !            58: float log(alogxx)
        !            59: float alogxx;
        !            60: {
        !            61:        float ffr(); 
        !            62:        long fex();
        !            63:        float p();
        !            64:        float alogx;
        !            65:        float temp;
        !            66:        alogx=ffr(alogxx)*2.-1.;
        !            67:        /* "c" won't take this lovely expression! 
        !            68:        return(6.9314718e-1*(fex(alogxx)-1.)+alogx*(.999964239+
        !            69:                alogx*(-.49987412+alogx*(.33179902+alogx*(-.24073380+alogx*(.16765407+
        !            70:                alogx*(-.09532938+alogx*(.03608849+alogx*(-.00645354)))))))));
        !            71:        */
        !            72: 
        !            73:        return(6.9314718e-1*(fex(alogxx)-1)+p(alogx,alogtab,8));
        !            74: }
        !            75: float sin(sinx)
        !            76: float sinx;
        !            77: {
        !            78:        float psin(),pcos();
        !            79:        long sinq,sinqz;
        !            80:        float sinr,sinz;
        !            81:        sinz=1.2732395*(sinx>=0.0?sinx:-sinx);
        !            82:        sinq=sinz;
        !            83:        sinr=sinz-sinq;
        !            84:        if(sinx<0.0)sinq = sinq+4;
        !            85:        sinqz=sinq&07;
        !            86:        if(!sinqz)      return(psin(sinr));
        !            87:        if(sinqz==1)    return(pcos(1.0-sinr));
        !            88:        if(sinqz==2)    return(pcos(sinr));
        !            89:        if(sinqz==3)    return(psin(1.0-sinr));
        !            90:        if(sinqz==4)    return(-psin(sinr));
        !            91:        if(sinqz==5)    return(-pcos(1.0-sinr));
        !            92:        if(sinqz==6)    return(-pcos(sinr));
        !            93:        return(-psin(1.0-sinr));
        !            94: }
        !            95: float cos(cosx)
        !            96: float cosx;
        !            97: {
        !            98:        float sin();
        !            99:        return(sin(cosx+1.5707963));
        !           100: }
        !           101: float psin(psinss)
        !           102: float psinss;
        !           103: {
        !           104:        float psinst,psinst2;
        !           105:        psinst2=.78539816*psinss;
        !           106:        psinst=psinst2*psinst2;
        !           107:        return(psinst2*(1.-psinst*(1.66666667e-1-
        !           108:            psinst*(8.33333333e-3-psinst*(1.98412698e-4-psinst*
        !           109:            (2.75573192e-6))))));
        !           110: }
        !           111: float pcos(pcoscs)
        !           112: float pcoscs;
        !           113: {
        !           114:        float pcosct;
        !           115:        pcosct=.616850275*pcoscs*pcoscs;
        !           116:        return(1.-pcosct*(.5-pcosct*(4.16666667e-2-pcosct*
        !           117:            (1.38888888e-3-pcosct*(2.48015873e-5-pcosct*(2.75573192e-7))))));
        !           118: }
        !           119: float sqrt(sqrtarg)
        !           120: float sqrtarg;
        !           121: {
        !           122:        float ffr(); 
        !           123:        long fex();
        !           124:        float sqrtiarg;
        !           125:        float two();
        !           126:        sqrtiarg = two((long)(fex(sqrtarg )*.5-1.))*(1.+ffr(sqrtarg));
        !           127:        sqrtiarg = .5*(sqrtiarg+sqrtarg/sqrtiarg);
        !           128:        sqrtiarg = .5*(sqrtiarg+sqrtarg/sqrtiarg);
        !           129:        sqrtiarg = .5*(sqrtiarg+sqrtarg/sqrtiarg);
        !           130:        return(.5*(sqrtiarg+sqrtarg/sqrtiarg));
        !           131: }
        !           132: 
        !           133: static long fex(arg)
        !           134: long arg;
        !           135: {
        !           136:        if(_fsign2==0100000L)   arg >>= 7; 
        !           137:        else                    arg >>= 23;
        !           138:        return((arg&0377)-0200);
        !           139: }
        !           140: 
        !           141: static long ffr(arg)
        !           142: long arg;
        !           143: {
        !           144:        if(arg==0)      return(0);
        !           145:        if(_fsign2==0100000L) {
        !           146:                arg = arg & (~077600L);
        !           147:                return(arg|040000L);
        !           148:        }
        !           149:        arg = arg & (~(077600L<<16));
        !           150:        return(arg|(040000L<<16));
        !           151: }
        !           152: 
        !           153: static long two(arg)
        !           154: long arg;
        !           155: {
        !           156:        if(_fsign2==0100000L)   return((arg+0201L)<<7);
        !           157:        return((arg+0201L)<<23);
        !           158: }
        !           159: static float p(x,t,n)
        !           160: float x,t[];
        !           161: int n;
        !           162: {
        !           163:        int i;
        !           164:        float ans;
        !           165:        ans=t[0]; 
        !           166:        for(i=1; i<n; i++)      ans=ans*x+t[i];
        !           167:        return(x*ans);
        !           168: }

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.