![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to rhmath.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd / PDP11 / fpp / aux|
Return to rhmath.c CVS log | Up to [Research Unix] / researchv10no / cmd / PDP11 / fpp / aux |
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.