1 /*
2  * ====================================================
3  * Copyright (C) 1998, 2002 by Red Hat Inc. All rights reserved.
4  *
5  * Permission to use, copy, modify, and distribute this
6  * software is freely granted, provided that this notice
7  * is preserved.
8  * ====================================================
9  */
10 
11 #if !defined(_SOFT_FLOAT)
12 
13 /*
14 Fast version of pow using Intel float instructions.
15 
16    double _f_pow (double x, double y);
17 
18 Function calculates x to power of y.
19 The function optimizes the case where x is >0.0 and y is finite.
20 In such a case, there is no error checking or setting of errno.
21 All other cases defer to normal pow() function which will
22 set errno as normal.
23 */
24 
25 #include <math.h>
26 #include <ieeefp.h>
27 #include "f_math.h"
28 
_f_pow(double x,double y)29 double _f_pow (double x, double y)
30 {
31   /* following sequence handles the majority of cases for pow() */
32   if (x > 0.0 && check_finite(y))
33     {
34       double result;
35       /* calculate x ** y as 2 ** (y log2(x)).  On Intel, can only
36          raise 2 to an integer or a small fraction, thus, we have
37          to perform two steps 2**integer portion * 2**fraction. */
38       __asm__("fyl2x; fld %%st; frndint; fsub %%st,%%st(1);"\
39            "fxch; fchs; f2xm1; fld1; faddp; fxch; fld1; fscale; fstp %%st(1);"\
40            "fmulp" : "=t" (result) : "0" (x), "u" (y) : "st(1)" );
41       return result;
42     }
43   else /* all other strange cases, defer to normal pow() */
44     return pow (x,y);
45 }
46 
47 #endif
48