2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001 Free Software Foundation
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.1 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 Lesser 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:slowpow.c */
23 /* FUNCTION:slowpow */
25 /*FILES NEEDED:mpa.h */
26 /* mpa.c mpexp.c mplog.c halfulp.c */
28 /* Given two IEEE double machine numbers y,x , routine computes the */
29 /* correctly rounded (to nearest) value of x^y. Result calculated by */
30 /* multiplication (in halfulp.c) or if result isn't accurate enough */
31 /* then routine converts x and y into multi-precision doubles and */
32 /* calls to mpexp routine */
33 /*************************************************************************/
36 #include "math_private.h"
38 void __mpexp(mp_no
*x
, mp_no
*y
, int p
);
39 void __mplog(mp_no
*x
, mp_no
*y
, int p
);
41 double __halfulp(double x
,double y
);
43 double __slowpow(double x
, double y
, double z
) {
45 mp_no mpx
, mpy
, mpz
,mpw
,mpp
,mpr
,mpr1
;
46 static const mp_no eps
= {-3,{1.0,4.0}};
49 res
= __halfulp(x
,y
); /* halfulp() returns -10 or x^y */
50 if (res
>= 0) return res
; /* if result was really computed by halfulp */
51 /* else, if result was not really computed by halfulp */
52 p
= 10; /* p=precision */
56 __mplog(&mpx
, &mpz
, p
); /* log(x) = z */
57 __mul(&mpy
,&mpz
,&mpw
,p
); /* y * z =w */
58 __mpexp(&mpw
, &mpp
, p
); /* e^w =pp */
59 __add(&mpp
,&eps
,&mpr
,p
); /* pp+eps =r */
60 __mp_dbl(&mpr
, &res
, p
);
61 __sub(&mpp
,&eps
,&mpr1
,p
); /* pp -eps =r1 */
62 __mp_dbl(&mpr1
, &res1
, p
); /* converting into double precision */
63 if (res
== res1
) return res
;
65 p
= 32; /* if we get here result wasn't calculated exactly, continue */
66 __dbl_mp(x
,&mpx
,p
); /* for more exact calculation */
69 __mplog(&mpx
, &mpz
, p
); /* log(c)=z */
70 __mul(&mpy
,&mpz
,&mpw
,p
); /* y*z =w */
71 __mpexp(&mpw
, &mpp
, p
); /* e^w=pp */
72 __mp_dbl(&mpp
, &res
, p
); /* converting into double precision */