]>
Commit | Line | Data |
---|---|---|
e4d82761 UD |
1 | |
2 | /* | |
3 | * IBM Accurate Mathematical Library | |
aeb25823 AJ |
4 | * written by International Business Machines Corp. |
5 | * Copyright (C) 2001 Free Software Foundation | |
e4d82761 UD |
6 | * |
7 | * This program is free software; you can redistribute it and/or modify | |
8 | * it under the terms of the GNU Lesser General Public License as published by | |
cc7375ce | 9 | * the Free Software Foundation; either version 2.1 of the License, or |
e4d82761 | 10 | * (at your option) any later version. |
50944bca | 11 | * |
e4d82761 UD |
12 | * This program is distributed in the hope that it will be useful, |
13 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
c6c6dd48 | 15 | * GNU Lesser General Public License for more details. |
e4d82761 UD |
16 | * |
17 | * You should have received a copy of the GNU Lesser General Public License | |
59ba27a6 | 18 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
e4d82761 UD |
19 | */ |
20 | /************************************************************************/ | |
21 | /* */ | |
22 | /* MODULE_NAME:mplog.c */ | |
23 | /* */ | |
24 | /* FUNCTIONS: mplog */ | |
25 | /* */ | |
26 | /* FILES NEEDED: endian.h mpa.h mplog.h */ | |
27 | /* mpexp.c */ | |
28 | /* */ | |
29 | /* Multi-Precision logarithm function subroutine (for precision p >= 4, */ | |
30 | /* 2**(-1024) < x < 2**1024) and x is outside of the interval */ | |
31 | /* [1-2**(-54),1+2**(-54)]. Upon entry, x should be set to the */ | |
32 | /* multi-precision value of the input and y should be set into a multi- */ | |
33 | /* precision value of an approximation of log(x) with relative error */ | |
34 | /* bound of at most 2**(-52). The routine improves the accuracy of y. */ | |
35 | /* */ | |
36 | /************************************************************************/ | |
37 | #include "endian.h" | |
38 | #include "mpa.h" | |
39 | ||
ca58f1db | 40 | void __mpexp(mp_no *, mp_no *, int); |
e4d82761 | 41 | |
ca58f1db | 42 | void __mplog(mp_no *x, mp_no *y, int p) { |
e4d82761 | 43 | #include "mplog.h" |
50944bca UD |
44 | int i,m; |
45 | #if 0 | |
46 | int j,k,m1,m2,n; | |
e4d82761 | 47 | double a,b; |
50944bca | 48 | #endif |
e4d82761 UD |
49 | static const int mp[33] = {0,0,0,0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3, |
50 | 4,4,4,4,4,4,4,4,4,4,4,4,4,4}; | |
50944bca | 51 | 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, |
e4d82761 | 52 | 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, |
50944bca | 53 | 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; |
e4d82761 UD |
54 | mp_no mpt1,mpt2; |
55 | ||
56 | /* Choose m and initiate mpone */ | |
57 | m = mp[p]; mpone.e = 1; mpone.d[0]=mpone.d[1]=ONE; | |
58 | ||
59 | /* Perform m newton iterations to solve for y: exp(y)-x=0. */ | |
60 | /* The iterations formula is: y(n+1)=y(n)+(x*exp(-y(n))-1). */ | |
ca58f1db | 61 | __cpy(y,&mpt1,p); |
e4d82761 UD |
62 | for (i=0; i<m; i++) { |
63 | mpt1.d[0]=-mpt1.d[0]; | |
ca58f1db UD |
64 | __mpexp(&mpt1,&mpt2,p); |
65 | __mul(x,&mpt2,&mpt1,p); | |
66 | __sub(&mpt1,&mpone,&mpt2,p); | |
67 | __add(y,&mpt2,&mpt1,p); | |
68 | __cpy(&mpt1,y,p); | |
e4d82761 UD |
69 | } |
70 | return; | |
71 | } |