181 #endif
182 {
183 int jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
184 double z,fw,f[20],fq[20],q[20];
185
186 /* initialize jk*/
187 jk = init_jk[prec];
188 jp = jk;
189
190 /* determine jx,jv,q0, note that 3>q0 */
191 jx = nx-1;
192 jv = (e0-3)/24; if(jv<0) jv=0;
193 q0 = e0-24*(jv+1);
194
195 /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
196 j = jv-jx; m = jx+jk;
197 for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
198
199 /* compute q[0],q[1],...q[jk] */
200 for (i=0;i<=jk;i++) {
201 for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
202 }
203
204 jz = jk;
205 recompute:
206 /* distill q[] into iq[] reversingly */
207 for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
208 fw = (double)((int)(twon24* z));
209 iq[i] = (int)(z-two24*fw);
210 z = q[j-1]+fw;
211 }
212
213 /* compute n */
214 z = scalbn(z,q0); /* actual value of z */
215 z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */
216 n = (int) z;
217 z -= (double)n;
218 ih = 0;
219 if(q0>0) { /* need iq[jz-1] to determine n */
220 i = (iq[jz-1]>>(24-q0)); n += i;
221 iq[jz-1] -= i<<(24-q0);
|
181 #endif
182 {
183 int jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
184 double z,fw,f[20],fq[20],q[20];
185
186 /* initialize jk*/
187 jk = init_jk[prec];
188 jp = jk;
189
190 /* determine jx,jv,q0, note that 3>q0 */
191 jx = nx-1;
192 jv = (e0-3)/24; if(jv<0) jv=0;
193 q0 = e0-24*(jv+1);
194
195 /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
196 j = jv-jx; m = jx+jk;
197 for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
198
199 /* compute q[0],q[1],...q[jk] */
200 for (i=0;i<=jk;i++) {
201 for (j=0,fw=0.0;j<=jx;j++) {
202 fw += x[j]*f[jx+i-j]; q[i] = fw;
203 }
204 }
205
206 jz = jk;
207 recompute:
208 /* distill q[] into iq[] reversingly */
209 for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
210 fw = (double)((int)(twon24* z));
211 iq[i] = (int)(z-two24*fw);
212 z = q[j-1]+fw;
213 }
214
215 /* compute n */
216 z = scalbn(z,q0); /* actual value of z */
217 z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */
218 n = (int) z;
219 z -= (double)n;
220 ih = 0;
221 if(q0>0) { /* need iq[jz-1] to determine n */
222 i = (iq[jz-1]>>(24-q0)); n += i;
223 iq[jz-1] -= i<<(24-q0);
|