This source file includes following definitions.
- poly_l2
- poly_l2p1
1
2
3
4
5
6
7
8
9
10
11
12
13 #include "exception.h"
14 #include "reg_constant.h"
15 #include "fpu_emu.h"
16
17
18
19 #define HIPOWER 9
20 static unsigned short lterms[HIPOWER][4] =
21 {
22
23
24 { 0xe177, 0xb82f, 0x7652, 0x7154 },
25 { 0xee0f, 0xe80f, 0x2770, 0x7b1c },
26 { 0x0fc0, 0xbe87, 0xb143, 0x49dd },
27 { 0x78b9, 0xdadd, 0xec54, 0x34c2 },
28 { 0x003a, 0x5de9, 0x628b, 0x2909 },
29 { 0x5588, 0xed16, 0x4abf, 0x2193 },
30 { 0xb461, 0x85f7, 0x347a, 0x1c6a },
31 { 0x0975, 0x87b3, 0xd5bf, 0x1876 },
32 { 0xe85c, 0xcec9, 0x84e7, 0x187d }
33 };
34
35
36
37
38
39
40
41 void poly_l2(REG *arg, REG *result)
42 {
43 short exponent;
44 char zero;
45 unsigned short bits, shift;
46 long long Xsq;
47 REG accum, denom, num, Xx;
48
49
50 exponent = arg->exp - EXP_BIAS;
51
52 accum.tag = TW_Valid;
53
54 if ( arg->sigh > (unsigned)0xb504f334 )
55 {
56
57
58
59
60 exponent++;
61 accum.sign = 1;
62 num.exp = EXP_BIAS;
63 reg_u_div((long long *)&(CONST_1.sigl), (long long *)&(arg->sigl), &num);
64 }
65 else
66 {
67 accum.sign = 0;
68 num.sigl = arg->sigl;
69 num.sigh = arg->sigh;
70 }
71
72
73 num.sigh <<= 1;
74 if ( num.sigl & 0x80000000 ) num.sigh |= 1;
75 num.sigl <<= 1;
76
77 denom.sigl = num.sigl;
78 denom.sigh = num.sigh;
79 poly_div4((long long *)&(denom.sigl));
80 denom.sigh += 0x80000000;
81 Xx.exp = EXP_BIAS;
82 reg_u_div((long long *)&num.sigl, (long long *)&(denom.sigl), &Xx);
83
84 zero = !(Xx.sigh | Xx.sigl);
85
86 mul64((long long *)&Xx.sigl, (long long *)&Xx.sigl, &Xsq);
87 poly_div16(&Xsq);
88
89 accum.exp = -1;
90
91
92 polynomial((unsigned *)&accum.sigl, (unsigned *)&Xsq, lterms, HIPOWER-1);
93
94 if ( !exponent )
95 {
96
97
98
99 REG lXx;
100 char sign;
101
102 sign = accum.sign;
103 accum.sign = 0;
104
105
106 accum.exp = EXP_BIAS + accum.exp;
107 normalize(&accum);
108
109 if ( zero )
110 {
111 reg_move(&CONST_Z, result);
112 }
113 else
114 {
115
116 num.tag = TW_Valid;
117 num.sign = 0;
118 num.exp = EXP_BIAS - 1;
119 if ( sign )
120 {
121
122
123 *((long long *)&num.sigl) = - *((long long *)&(arg->sigl));
124 normalize(&num);
125 reg_div(&num, arg, &num);
126 }
127 else
128 {
129 normalize(&num);
130 }
131
132 denom.tag = TW_Valid;
133 denom.sign = SIGN_POS;
134 denom.exp = EXP_BIAS;
135
136 reg_div(&num, &denom, &lXx);
137
138 reg_u_mul(&lXx, &accum, &accum);
139 accum.exp += - EXP_BIAS + 1;
140
141
142 reg_u_add(&lXx, &accum, result);
143
144 normalize(result);
145 }
146 result->sign = sign;
147 return;
148 }
149
150 mul64((long long *)&accum.sigl,
151 (long long *)&Xx.sigl, (long long *)&accum.sigl);
152
153 *((long long *)(&accum.sigl)) += *((long long *)(&Xx.sigl));
154
155 if ( Xx.sigh > accum.sigh )
156 {
157
158
159 poly_div2((long long *)&accum.sigl);
160 accum.sigh |= 0x80000000;
161 accum.exp++;
162 }
163
164
165
166
167 if ( exponent && ((exponent < 0) ^ (accum.sign)) )
168 {
169
170
171 accum.sign = !accum.sign;
172
173
174 if ( accum.sigl | accum.sigh )
175 {
176
177
178 if ( accum.exp < 0 )
179 {
180 poly_div2((long long *)&accum.sigl);
181 accum.exp++;
182 }
183
184 *((long long *)&(accum.sigl)) = - *((long long *)&(accum.sigl));
185 if ( exponent < 0 )
186 exponent++;
187 else
188 exponent--;
189 }
190 }
191
192 shift = exponent >= 0 ? exponent : -exponent ;
193 bits = 0;
194 if ( shift )
195 {
196 if ( accum.exp )
197 {
198 accum.exp++;
199 poly_div2((long long *)&accum.sigl);
200 }
201 while ( shift )
202 {
203 poly_div2((long long *)&accum.sigl);
204 if ( shift & 1)
205 accum.sigh |= 0x80000000;
206 shift >>= 1;
207 bits++;
208 }
209 }
210
211
212 accum.exp += bits + EXP_BIAS - 1;
213
214 reg_move(&accum, result);
215 normalize(result);
216
217 return;
218 }
219
220
221
222
223
224
225 int poly_l2p1(REG *arg, REG *result)
226 {
227 char sign = 0;
228 long long Xsq;
229 REG arg_pl1, denom, accum, local_arg, poly_arg;
230
231
232 sign = arg->sign;
233
234 reg_add(arg, &CONST_1, &arg_pl1);
235
236 if ( (arg_pl1.sign) | (arg_pl1.tag) )
237 {
238 return 1;
239 }
240
241 reg_add(&CONST_1, &arg_pl1, &denom);
242 reg_div(arg, &denom, &local_arg);
243 local_arg.sign = 0;
244
245
246
247
248 if ( local_arg.exp >= EXP_BIAS - 3 )
249 {
250 if ( (local_arg.exp > EXP_BIAS - 3) ||
251 (local_arg.sigh > (unsigned)0xafb0ccc0) )
252 {
253
254 poly_l2(&arg_pl1, result); return 0;
255 }
256 }
257
258
259 reg_move(&local_arg, &poly_arg);
260
261
262 shrx((unsigned *)&(poly_arg.sigl), -(poly_arg.exp - EXP_BIAS + 3));
263
264 mul64((long long *)&(poly_arg.sigl), (long long *)&(poly_arg.sigl), &Xsq);
265 poly_div16(&Xsq);
266
267
268 polynomial(&(accum.sigl), (unsigned *)&Xsq, lterms, HIPOWER-1);
269
270 accum.tag = TW_Valid;
271 accum.sign = SIGN_POS;
272
273
274 accum.exp = EXP_BIAS - 1;
275 normalize(&accum);
276
277 reg_u_mul(&local_arg, &accum, &accum);
278 accum.exp -= EXP_BIAS - 1;
279
280 reg_u_add(&local_arg, &accum, result);
281
282
283 result->exp++;
284
285 result->sign = sign;
286
287 return 0;
288 }