OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libquadmath/] [math/] [rem_pio2q.c] - Blame information for rev 820

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 740 jeremybenn
#include "quadmath-imp.h"
2
#include <math.h>
3
 
4
 
5
/* @(#)k_rem_pio2.c 5.1 93/09/24 */
6
/*
7
 * ====================================================
8
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
9
 *
10
 * Developed at SunPro, a Sun Microsystems, Inc. business.
11
 * Permission to use, copy, modify, and distribute this
12
 * software is freely granted, provided that this notice
13
 * is preserved.
14
 * ====================================================
15
 */
16
 
17
/*
18
 * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
19
 * double x[],y[]; int e0,nx,prec; int ipio2[];
20
 *
21
 * __kernel_rem_pio2 return the last three digits of N with
22
 *              y = x - N*pi/2
23
 * so that |y| < pi/2.
24
 *
25
 * The method is to compute the integer (mod 8) and fraction parts of
26
 * (2/pi)*x without doing the full multiplication. In general we
27
 * skip the part of the product that are known to be a huge integer (
28
 * more accurately, = 0 mod 8 ). Thus the number of operations are
29
 * independent of the exponent of the input.
30
 *
31
 * (2/pi) is represented by an array of 24-bit integers in ipio2[].
32
 *
33
 * Input parameters:
34
 *      x[]     The input value (must be positive) is broken into nx
35
 *              pieces of 24-bit integers in double precision format.
36
 *              x[i] will be the i-th 24 bit of x. The scaled exponent
37
 *              of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
38
 *              match x's up to 24 bits.
39
 *
40
 *              Example of breaking a double positive z into x[0]+x[1]+x[2]:
41
 *                      e0 = ilogb(z)-23
42
 *                      z  = scalbn(z,-e0)
43
 *              for i = 0,1,2
44
 *                      x[i] = floor(z)
45
 *                      z    = (z-x[i])*2**24
46
 *
47
 *
48
 *      y[]     ouput result in an array of double precision numbers.
49
 *              The dimension of y[] is:
50
 *                      24-bit  precision       1
51
 *                      53-bit  precision       2
52
 *                      64-bit  precision       2
53
 *                      113-bit precision       3
54
 *              The actual value is the sum of them. Thus for 113-bit
55
 *              precision, one may have to do something like:
56
 *
57
 *              long double t,w,r_head, r_tail;
58
 *              t = (long double)y[2] + (long double)y[1];
59
 *              w = (long double)y[0];
60
 *              r_head = t+w;
61
 *              r_tail = w - (r_head - t);
62
 *
63
 *      e0      The exponent of x[0]
64
 *
65
 *      nx      dimension of x[]
66
 *
67
 *      prec    an integer indicating the precision:
68
 *                      0        24  bits (single)
69
 *                      1       53  bits (double)
70
 *                      2       64  bits (extended)
71
 *                      3       113 bits (quad)
72
 *
73
 *      ipio2[]
74
 *              integer array, contains the (24*i)-th to (24*i+23)-th
75
 *              bit of 2/pi after binary point. The corresponding
76
 *              floating value is
77
 *
78
 *                      ipio2[i] * 2^(-24(i+1)).
79
 *
80
 * External function:
81
 *      double scalbn(), floor();
82
 *
83
 *
84
 * Here is the description of some local variables:
85
 *
86
 *      jk      jk+1 is the initial number of terms of ipio2[] needed
87
 *              in the computation. The recommended value is 2,3,4,
88
 *              6 for single, double, extended,and quad.
89
 *
90
 *      jz      local integer variable indicating the number of
91
 *              terms of ipio2[] used.
92
 *
93
 *      jx      nx - 1
94
 *
95
 *      jv      index for pointing to the suitable ipio2[] for the
96
 *              computation. In general, we want
97
 *                      ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
98
 *              is an integer. Thus
99
 *                      e0-3-24*jv >= 0 or (e0-3)/24 >= jv
100
 *              Hence jv = max(0,(e0-3)/24).
101
 *
102
 *      jp      jp+1 is the number of terms in PIo2[] needed, jp = jk.
103
 *
104
 *      q[]     double array with integral value, representing the
105
 *              24-bits chunk of the product of x and 2/pi.
106
 *
107
 *      q0      the corresponding exponent of q[0]. Note that the
108
 *              exponent for q[i] would be q0-24*i.
109
 *
110
 *      PIo2[]  double precision array, obtained by cutting pi/2
111
 *              into 24 bits chunks.
112
 *
113
 *      f[]     ipio2[] in floating point
114
 *
115
 *      iq[]    integer array by breaking up q[] in 24-bits chunk.
116
 *
117
 *      fq[]    final product of x*(2/pi) in fq[0],..,fq[jk]
118
 *
119
 *      ih      integer. If >0 it indicates q[] is >= 0.5, hence
120
 *              it also indicates the *sign* of the result.
121
 *
122
 */
123
 
124
/*
125
 * Constants:
126
 * The hexadecimal values are the intended ones for the following
127
 * constants. The decimal values may be used, provided that the
128
 * compiler will convert from decimal to binary accurately enough
129
 * to produce the hexadecimal values shown.
130
 */
131
 
132
 
133
static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
134
 
135
static const double PIo2[] = {
136
  1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
137
  7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
138
  5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
139
  3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
140
  1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
141
  1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
142
  2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
143
  2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
144
};
145
 
146
static const double
147
  zero   = 0.0,
148
  one    = 1.0,
149
  two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
150
  twon24  =  5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
151
 
152
 
153
static int
154
__quadmath_kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec, const int32_t *ipio2)
155
{
156
        int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
157
        double z,fw,f[20],fq[20],q[20];
158
 
159
    /* initialize jk*/
160
        jk = init_jk[prec];
161
        jp = jk;
162
 
163
    /* determine jx,jv,q0, note that 3>q0 */
164
        jx =  nx-1;
165
        jv = (e0-3)/24; if(jv<0) jv=0;
166
        q0 =  e0-24*(jv+1);
167
 
168
    /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
169
        j = jv-jx; m = jx+jk;
170
        for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
171
 
172
    /* compute q[0],q[1],...q[jk] */
173
        for (i=0;i<=jk;i++) {
174
            for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
175
        }
176
 
177
        jz = jk;
178
recompute:
179
    /* distill q[] into iq[] reversingly */
180
        for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
181
            fw    =  (double)((int32_t)(twon24* z));
182
            iq[i] =  (int32_t)(z-two24*fw);
183
            z     =  q[j-1]+fw;
184
        }
185
 
186
    /* compute n */
187
        z  = scalbn(z,q0);              /* actual value of z */
188
        z -= 8.0*floor(z*0.125);                /* trim off integer >= 8 */
189
        n  = (int32_t) z;
190
        z -= (double)n;
191
        ih = 0;
192
        if(q0>0) {       /* need iq[jz-1] to determine n */
193
            i  = (iq[jz-1]>>(24-q0)); n += i;
194
            iq[jz-1] -= i<<(24-q0);
195
            ih = iq[jz-1]>>(23-q0);
196
        }
197
        else if(q0==0) ih = iq[jz-1]>>23;
198
        else if(z>=0.5) ih=2;
199
 
200
        if(ih>0) {       /* q > 0.5 */
201
            n += 1; carry = 0;
202
            for(i=0;i<jz ;i++) { /* compute 1-q */
203
                j = iq[i];
204
                if(carry==0) {
205
                    if(j!=0) {
206
                        carry = 1; iq[i] = 0x1000000- j;
207
                    }
208
                } else  iq[i] = 0xffffff - j;
209
            }
210
            if(q0>0) {           /* rare case: chance is 1 in 12 */
211
                switch(q0) {
212
                case 1:
213
                   iq[jz-1] &= 0x7fffff; break;
214
                case 2:
215
                   iq[jz-1] &= 0x3fffff; break;
216
                }
217
            }
218
            if(ih==2) {
219
                z = one - z;
220
                if(carry!=0) z -= scalbn(one,q0);
221
            }
222
        }
223
 
224
    /* check if recomputation is needed */
225
        if(z==zero) {
226
            j = 0;
227
            for (i=jz-1;i>=jk;i--) j |= iq[i];
228
            if(j==0) { /* need recomputation */
229
                for(k=1;iq[jk-k]==0;k++);   /* k = no. of terms needed */
230
 
231
                for(i=jz+1;i<=jz+k;i++) {   /* add q[jz+1] to q[jz+k] */
232
                    f[jx+i] = (double) ipio2[jv+i];
233
                    for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
234
                    q[i] = fw;
235
                }
236
                jz += k;
237
                goto recompute;
238
            }
239
        }
240
 
241
    /* chop off zero terms */
242
        if(z==0.0) {
243
            jz -= 1; q0 -= 24;
244
            while(iq[jz]==0) { jz--; q0-=24;}
245
        } else { /* break z into 24-bit if necessary */
246
            z = scalbn(z,-q0);
247
            if(z>=two24) {
248
                fw = (double)((int32_t)(twon24*z));
249
                iq[jz] = (int32_t)(z-two24*fw);
250
                jz += 1; q0 += 24;
251
                iq[jz] = (int32_t) fw;
252
            } else iq[jz] = (int32_t) z ;
253
        }
254
 
255
    /* convert integer "bit" chunk to floating-point value */
256
        fw = scalbn(one,q0);
257
        for(i=jz;i>=0;i--) {
258
            q[i] = fw*(double)iq[i]; fw*=twon24;
259
        }
260
 
261
    /* compute PIo2[0,...,jp]*q[jz,...,0] */
262
        for(i=jz;i>=0;i--) {
263
            for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
264
            fq[jz-i] = fw;
265
        }
266
 
267
    /* compress fq[] into y[] */
268
        switch(prec) {
269
            case 0:
270
                fw = 0.0;
271
                for (i=jz;i>=0;i--) fw += fq[i];
272
                y[0] = (ih==0)? fw: -fw;
273
                break;
274
            case 1:
275
            case 2:
276
                fw = 0.0;
277
                for (i=jz;i>=0;i--) fw += fq[i];
278
                y[0] = (ih==0)? fw: -fw;
279
                fw = fq[0]-fw;
280
                for (i=1;i<=jz;i++) fw += fq[i];
281
                y[1] = (ih==0)? fw: -fw;
282
                break;
283
            case 3:     /* painful */
284
                for (i=jz;i>0;i--) {
285
#if __FLT_EVAL_METHOD__ != 0
286
                    volatile
287
#endif
288
                    double fv = (double)(fq[i-1]+fq[i]);
289
                    fq[i]  += fq[i-1]-fv;
290
                    fq[i-1] = fv;
291
                }
292
                for (i=jz;i>1;i--) {
293
#if __FLT_EVAL_METHOD__ != 0
294
                    volatile
295
#endif
296
                    double fv = (double)(fq[i-1]+fq[i]);
297
                    fq[i]  += fq[i-1]-fv;
298
                    fq[i-1] = fv;
299
                }
300
                for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
301
                if(ih==0) {
302
                    y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
303
                } else {
304
                    y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
305
                }
306
        }
307
        return n&7;
308
}
309
 
310
 
311
 
312
 
313
 
314
/* Quad-precision floating point argument reduction.
315
   Copyright (C) 1999 Free Software Foundation, Inc.
316
   This file is part of the GNU C Library.
317
   Contributed by Jakub Jelinek <jj@ultra.linux.cz>
318
 
319
   The GNU C Library is free software; you can redistribute it and/or
320
   modify it under the terms of the GNU Lesser General Public
321
   License as published by the Free Software Foundation; either
322
   version 2.1 of the License, or (at your option) any later version.
323
 
324
   The GNU C Library is distributed in the hope that it will be useful,
325
   but WITHOUT ANY WARRANTY; without even the implied warranty of
326
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
327
   Lesser General Public License for more details.
328
 
329
   You should have received a copy of the GNU Lesser General Public
330
   License along with the GNU C Library; if not, write to the Free
331
   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
332
   02111-1307 USA.  */
333
 
334
/*
335
 * Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi
336
 */
337
static const int32_t two_over_pi[] = {
338
0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62,
339
0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a,
340
0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129,
341
0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41,
342
0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8,
343
0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf,
344
0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5,
345
0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08,
346
0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3,
347
0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880,
348
0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b,
349
0x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6,
350
0xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2,
351
0xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35,
352
0xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30,
353
0x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c,
354
0x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4,
355
0xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770,
356
0xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7,
357
0xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19,
358
0xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522,
359
0x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16,
360
0xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6,
361
0x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e,
362
0x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48,
363
0xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3,
364
0xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf,
365
0xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55,
366
0x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612,
367
0xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929,
368
0x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec,
369
0xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b,
370
0x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c,
371
0x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4,
372
0x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb,
373
0xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc,
374
0x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c,
375
0x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f,
376
0xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5,
377
0x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437,
378
0x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b,
379
0x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea,
380
0x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad,
381
0x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3,
382
0x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3,
383
0xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717,
384
0x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f,
385
0xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61,
386
0xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db,
387
0xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51,
388
0xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0,
389
0x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c,
390
0x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6,
391
0x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc,
392
0xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed,
393
0x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328,
394
0x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d,
395
0x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0,
396
0xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b,
397
0x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4,
398
0xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3,
399
0xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f,
400
0x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad,
401
0x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b,
402
0x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4,
403
0xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761,
404
0x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31,
405
0x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30,
406
0xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262,
407
0x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e,
408
0xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1,
409
0xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c,
410
0x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4,
411
0xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08,
412
0x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196,
413
0xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9,
414
0x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4,
415
0x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc,
416
0x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c,
417
0xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0,
418
0x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c,
419
0xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0,
420
0x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac,
421
0xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22,
422
0x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893,
423
0x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7,
424
0xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5,
425
0x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f,
426
0xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4,
427
0x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf,
428
0x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b,
429
0x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2,
430
0x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138,
431
0xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e,
432
0xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569,
433
0x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34,
434
0xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9,
435
0x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d,
436
0xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f,
437
0xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855,
438
0x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569,
439
0xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b,
440
0x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe,
441
0x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41,
442
0x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49,
443
0xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f,
444
0xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110,
445
0x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8,
446
0x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365,
447
0xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a,
448
0x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270,
449
0x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5,
450
0x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616,
451
0x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b,
452
0x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0,
453
0xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb,
454
0x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a,
455
0x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e,
456
0xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa,
457
0x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5,
458
0xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0,
459
0xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2,
460
0xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886,
461
0x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142,
462
0xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba,
463
0xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4,
464
0xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708,
465
0x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555,
466
0xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3,
467
0x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55,
468
0xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58,
469
0xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5,
470
0x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c,
471
0x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe,
472
0x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b,
473
0x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8,
474
0x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005,
475
0x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7,
476
0xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50,
477
0x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604,
478
0x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643,
479
0x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485,
480
0xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d,
481
0xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6,
482
0xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2,
483
0x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02,
484
0xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3,
485
0x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412,
486
0xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274,
487
0x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755,
488
0xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849,
489
0xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce,
490
0xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5,
491
0x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba,
492
0xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6,
493
0xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d,
494
0x7b7b89, 0x483d38,
495
};
496
 
497
static const __float128 c[] = {
498
/* 93 bits of pi/2 */
499
#define PI_2_1 c[0]
500
 1.57079632679489661923132169155131424e+00Q, /* 3fff921fb54442d18469898cc5100000 */
501
 
502
/* pi/2 - PI_2_1 */
503
#define PI_2_1t c[1]
504
 8.84372056613570112025531863263659260e-29Q, /* 3fa1c06e0e68948127044533e63a0106 */
505
};
506
 
507
 
508
int32_t
509
__quadmath_rem_pio2q (__float128 x, __float128 *y)
510
{
511
  __float128 z, w, t;
512
  double tx[8];
513
  int64_t exp, n, ix, hx;
514
  uint64_t lx;
515
 
516
  GET_FLT128_WORDS64 (hx, lx, x);
517
  ix = hx & 0x7fffffffffffffffLL;
518
  if (ix <= 0x3ffe921fb54442d1LL)       /* x in <-pi/4, pi/4> */
519
    {
520
      y[0] = x;
521
      y[1] = 0;
522
      return 0;
523
    }
524
 
525
  if (ix < 0x40002d97c7f3321dLL)        /* |x| in <pi/4, 3pi/4) */
526
    {
527
      if (hx > 0)
528
        {
529
          /* 113 + 93 bit PI is ok */
530
          z = x - PI_2_1;
531
          y[0] = z - PI_2_1t;
532
          y[1] = (z - y[0]) - PI_2_1t;
533
          return 1;
534
        }
535
      else
536
        {
537
          /* 113 + 93 bit PI is ok */
538
          z = x + PI_2_1;
539
          y[0] = z + PI_2_1t;
540
          y[1] = (z - y[0]) + PI_2_1t;
541
          return -1;
542
        }
543
    }
544
 
545
  if (ix >= 0x7fff000000000000LL)       /* x is +=oo or NaN */
546
    {
547
      y[0] = x - x;
548
      y[1] = y[0];
549
      return 0;
550
    }
551
 
552
  /* Handle large arguments.
553
     We split the 113 bits of the mantissa into 5 24bit integers
554
     stored in a double array.  */
555
  exp = (ix >> 48) - 16383 - 23;
556
 
557
  /* This is faster than doing this in floating point, because we
558
     have to convert it to integers anyway and like this we can keep
559
     both integer and floating point units busy.  */
560
  tx [0] = (double)(((ix >> 25) & 0x7fffff) | 0x800000);
561
  tx [1] = (double)((ix >> 1) & 0xffffff);
562
  tx [2] = (double)(((ix << 23) | (lx >> 41)) & 0xffffff);
563
  tx [3] = (double)((lx >> 17) & 0xffffff);
564
  tx [4] = (double)((lx << 7) & 0xffffff);
565
 
566
  n = __quadmath_kernel_rem_pio2 (tx, tx + 5, exp,
567
                                  ((lx << 7) & 0xffffff) ? 5 : 4,
568
                                  3, two_over_pi);
569
 
570
  /* The result is now stored in 3 double values, we need to convert it into
571
     two __float128 values.  */
572
  t = (__float128) tx [6] + (__float128) tx [7];
573
  w = (__float128) tx [5];
574
 
575
  if (hx >= 0)
576
    {
577
      y[0] = w + t;
578
      y[1] = t - (y[0] - w);
579
      return n;
580
    }
581
  else
582
    {
583
      y[0] = -(w + t);
584
      y[1] = -t - (y[0] + w);
585
      return -n;
586
    }
587
}

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.