![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to poly_atan.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [MW Coherent from dump] / coherent / b / kernel / emulator|
Return to poly_atan.c CVS log | Up to [MW Coherent from dump] / coherent / b / kernel / emulator |
1.1 root 1: /*---------------------------------------------------------------------------+
2: | p_atan.c |
3: | |
4: | Compute the tan of a FPU_REG, using a polynomial approximation. |
5: | |
6: | Copyright (C) 1992 W. Metzenthen, 22 Parker St, Ormond, Vic 3163, |
7: | Australia. E-mail [email protected] |
8: | |
9: | |
10: +---------------------------------------------------------------------------*/
11:
12: #include "fpu_system.h"
13: #include "exception.h"
14: #include "reg_constant.h"
15: #include "fpu_emu.h"
16:
17:
18: #define HIPOWERon 6 /* odd poly, negative terms */
19: static unsigned oddnegterms[HIPOWERon][2] =
20: {
21: { 0x00000000, 0x00000000 }, /* for + 1.0 */
22: { 0x763b6f3d, 0x1adc4428 },
23: { 0x20f0630b, 0x0502909d },
24: { 0x4e825578, 0x0198ce38 },
25: { 0x22b7cb87, 0x008da6e3 },
26: { 0x9b30ca03, 0x00239c79 }
27: } ;
28:
29: #define HIPOWERop 6 /* odd poly, positive terms */
30: static unsigned oddplterms[HIPOWERop][2] =
31: {
32: { 0xa6f67cb8, 0x94d910bd },
33: { 0xa02ffab4, 0x0a43cb45 },
34: { 0x04265e6b, 0x02bf5655 },
35: { 0x0a728914, 0x00f280f7 },
36: { 0x6d640e01, 0x004d6556 },
37: { 0xf1dd2dbf, 0x000a530a }
38: };
39:
40:
41: static unsigned denomterm[2] =
42: { 0xfc4bd208, 0xea2e6612 };
43:
44:
45:
46: /*--- poly_atan() -----------------------------------------------------------+
47: | |
48: +---------------------------------------------------------------------------*/
49: void poly_atan(FPU_REG *arg)
50: {
51: char recursions = 0;
52: short exponent;
53: FPU_REG odd_poly, even_poly, pos_poly, neg_poly;
54: FPU_REG argSq;
55: long long arg_signif, argSqSq;
56:
57:
58: #ifdef PARANOID
59: if ( arg->sign != 0 ) /* Can't hack a number < 0.0 */
60: { arith_invalid(arg); return; }
61: #endif PARANOID
62:
63: exponent = arg->exp - EXP_BIAS;
64:
65: if ( arg->tag == TW_Zero )
66: {
67: /* Return 0.0 */
68: reg_move(&CONST_Z, arg);
69: return;
70: }
71:
72: if ( exponent >= -2 )
73: {
74: /* argument is in the range [0.25 .. 1.0] */
75: if ( exponent >= 0 )
76: {
77: #ifdef PARANOID
78: if ( (exponent == 0) &&
79: (arg->sigl == 0) && (arg->sigh == 0x80000000) )
80: #endif PARANOID
81: {
82: reg_move(&CONST_PI4, arg);
83: return;
84: }
85: #ifdef PARANOID
86: EXCEPTION(EX_INTERNAL|0x104); /* There must be a logic error */
87: #endif PARANOID
88: }
89:
90: /* If the argument is greater than sqrt(2)-1 (=0.414213562...) */
91: /* convert the argument by an identity for atan */
92: if ( (exponent >= -1) || (arg->sigh > 0xd413ccd0) )
93: {
94: FPU_REG numerator, denom;
95:
96: recursions++;
97:
98: arg_signif = *(long long *)&(arg->sigl);
99: if ( exponent < -1 )
100: {
101: if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
102: arg_signif++; /* round up */
103: }
104: *(long long *)&(numerator.sigl) = -arg_signif;
105: numerator.exp = EXP_BIAS - 1;
106: normalize(&numerator); /* 1 - arg */
107:
108: arg_signif = *(long long *)&(arg->sigl);
109: if ( shrx(&arg_signif, -exponent) >= 0x80000000U )
110: arg_signif++; /* round up */
111: *(long long *)&(denom.sigl) = arg_signif;
112: denom.sigh |= 0x80000000; /* 1 + arg */
113:
114: arg->exp = numerator.exp;
115: reg_u_div((long long *)&(numerator.sigl),
116: (long long *)&(denom.sigl), arg);
117:
118: exponent = arg->exp - EXP_BIAS;
119: }
120: }
121:
122: *(long long *)&arg_signif = *(long long *)&(arg->sigl);
123:
124: #ifdef PARANOID
125: /* This must always be true */
126: if ( exponent >= -1 )
127: {
128: EXCEPTION(EX_INTERNAL|0x120); /* There must be a logic error */
129: }
130: #endif PARANOID
131:
132: /* shift the argument right by the required places */
133: if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
134: arg_signif++; /* round up */
135:
136: /* Now have arg_signif with binary point at the left
137: .1xxxxxxxx */
138: mul64(&arg_signif, &arg_signif, (long long *)(&argSq.sigl));
139: mul64((long long *)(&argSq.sigl), (long long *)(&argSq.sigl), &argSqSq);
140:
141: /* will be a valid positive nr with expon = 0 */
142: *(short *)&(pos_poly.sign) = 0;
143: pos_poly.exp = EXP_BIAS;
144:
145: /* Do the basic fixed point polynomial evaluation */
146: polynomial(&pos_poly.sigl, (unsigned *)&argSqSq,
147: (unsigned short (*)[4])oddplterms, HIPOWERop-1);
148: mul64((long long *)(&argSq.sigl), (long long *)(&pos_poly.sigl),
149: (long long *)(&pos_poly.sigl));
150:
151: /* will be a valid positive nr with expon = 0 */
152: *(short *)&(neg_poly.sign) = 0;
153: neg_poly.exp = EXP_BIAS;
154:
155: /* Do the basic fixed point polynomial evaluation */
156: polynomial(&neg_poly.sigl, (unsigned *)&argSqSq,
157: (unsigned short (*)[4])oddnegterms, HIPOWERon-1);
158:
159: /* Subtract the mantissas */
160: *((long long *)(&pos_poly.sigl)) -= *((long long *)(&neg_poly.sigl));
161:
162: reg_move(&pos_poly, &odd_poly);
163: poly_add_1(&odd_poly);
164:
165: reg_u_mul(&odd_poly, arg, &odd_poly); /* The complete odd polynomial */
166: odd_poly.exp -= EXP_BIAS - 1;
167:
168: /* will be a valid positive nr with expon = 0 */
169: *(short *)&(even_poly.sign) = 0;
170:
171: mul64((long long *)(&argSq.sigl),
172: (long long *)(&denomterm), (long long *)(&even_poly.sigl));
173:
174: poly_add_1(&even_poly);
175:
176: reg_div(&odd_poly, &even_poly, arg);
177:
178: if ( recursions )
179: reg_sub(&CONST_PI4, arg, arg);
180: }
181:
182:
183: /* The argument to this function must be polynomial() compatible,
184: i.e. have an exponent (not checked) of EXP_BIAS-1 but need not
185: be normalized.
186: This function adds 1.0 to the (assumed positive) argument. */
187: void poly_add_1(FPU_REG *src)
188: {
189: /* Rounding in a consistent direction produces better results
190: for the use of this function in poly_atan. Simple truncation
191: is used here instead of round-to-nearest. */
192:
193: #ifdef OBSOLETE
194: char round = (src->sigl & 3) == 3;
195: #endif OBSOLETE
196:
197: shrx(&src->sigl, 1);
198:
199: #ifdef OBSOLETE
200: if ( round ) (*(long long *)&src->sigl)++; /* Round to even */
201: #endif OBSOLETE
202:
203: src->sigh |= 0x80000000;
204:
205: src->exp = EXP_BIAS;
206:
207: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.