]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/dbl-64/mpexp.c
3 * IBM Accurate Mathematical Library
4 * Copyright (c) International Business Machines Corp., 2001
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20 /*************************************************************************/
21 /* MODULE_NAME:mpexp.c */
23 /* FUNCTIONS: mpexp */
25 /* FILES NEEDED: mpa.h endian.h mpexp.h */
28 /* Multi-Precision exponential function subroutine */
29 /* ( for p >= 4, 2**(-55) <= abs(x) <= 1024 ). */
30 /*************************************************************************/
36 /* Multi-Precision exponential function subroutine (for p >= 4, */
37 /* 2**(-55) <= abs(x) <= 1024). */
38 void __mpexp(mp_no
*x
, mp_no
*y
, int p
) {
42 static const int np
[33] = {0,0,0,0,3,3,4,4,5,4,4,5,5,5,6,6,6,6,6,6,
43 6,6,6,6,7,7,7,7,8,8,8,8,8};
44 static const int m1p
[33]= {0,0,0,0,17,23,23,28,27,38,42,39,43,47,43,47,50,54,
45 57,60,64,67,71,74,68,71,74,77,70,73,76,78,81};
46 static const int m1np
[7][18] = {
47 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
48 { 0, 0, 0, 0,36,48,60,72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
49 { 0, 0, 0, 0,24,32,40,48,56,64,72, 0, 0, 0, 0, 0, 0, 0},
50 { 0, 0, 0, 0,17,23,29,35,41,47,53,59,65, 0, 0, 0, 0, 0},
51 { 0, 0, 0, 0, 0, 0,23,28,33,38,42,47,52,57,62,66, 0, 0},
52 { 0, 0, 0, 0, 0, 0, 0, 0,27, 0, 0,39,43,47,51,55,59,63},
53 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,43,47,50,54}};
54 mp_no mpone
= {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
55 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
56 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
57 mp_no mpk
= {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
58 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
59 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
60 mp_no mps
,mpak
,mpt1
,mpt2
;
62 /* Choose m,n and compute a=2**(-m) */
63 n
= np
[p
]; m1
= m1p
[p
]; a
= twomm1
[p
].d
;
64 for (i
=0; i
<EX
; i
++) a
*= RADIXI
;
65 for ( ; i
>EX
; i
--) a
*= RADIX
;
66 b
= X
[1]*RADIXI
; m2
= 24*EX
;
67 for (; b
<HALF
; m2
--) { a
*= TWO
; b
*= TWO
; }
69 for (i
=2; i
<=p
; i
++) { if (X
[i
]!=ZERO
) break; }
70 if (i
==p
+1) { m2
--; a
*= TWO
; }
74 for (i
=n
-1; i
>0; i
--,n
--) { if (m1np
[i
][p
]+m2
>0) break; }
77 /* Compute s=x*2**(-m). Put result in mps */
79 __mul(x
,&mpt1
,&mps
,p
);
81 /* Evaluate the polynomial. Put result in mpt2 */
82 mpone
.e
=1; mpone
.d
[0]=ONE
; mpone
.d
[1]=ONE
;
83 mpk
.e
= 1; mpk
.d
[0] = ONE
; mpk
.d
[1]=nn
[n
].d
;
84 __dvd(&mps
,&mpk
,&mpt1
,p
);
85 __add(&mpone
,&mpt1
,&mpak
,p
);
86 for (k
=n
-1; k
>1; k
--) {
87 __mul(&mps
,&mpak
,&mpt1
,p
);
89 __dvd(&mpt1
,&mpk
,&mpt2
,p
);
90 __add(&mpone
,&mpt2
,&mpak
,p
);
92 __mul(&mps
,&mpak
,&mpt1
,p
);
93 __add(&mpone
,&mpt1
,&mpt2
,p
);
95 /* Raise polynomial value to the power of 2**m. Put result in y */
96 for (k
=0,j
=0; k
<m
; ) {
97 __mul(&mpt2
,&mpt2
,&mpt1
,p
); k
++;
98 if (k
==m
) { j
=1; break; }
99 __mul(&mpt1
,&mpt1
,&mpt2
,p
); k
++;
101 if (j
) __cpy(&mpt1
,y
,p
);
102 else __cpy(&mpt2
,y
,p
);