1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
|
/* Copyright (C) 1993 Hongjiu Lu
This file is part of the Linux C Library.
The Linux C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public License as
published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.
The Linux C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Library General Public License for more details. */
#include <errno.h>
#include <math.h>
#include <fpu_control.h>
double pow (x,y)
double x,y;
{
int negative;
__volatile__ unsigned short cw, saved_cw; /* 80387 control word */
if (y == 0.0) return 1.0;
if (x == 0.0) return (y < 0.0) ? __infnan(EDOM) : 0.0;
if (y == 1.0) return x;
if (x < 0.0) {
long long tmp;
tmp = (long long) y;
/* Even or odd */
negative = tmp & 1;
/* That is only valid if |y| < 2^64. */
if (y != (double) tmp) {
return __infnan( EDOM);
}
x = -x;
} else {
negative = 0;
}
/*
* Inline assembler implements the following algorithm:
*
* 1 - x' = y *log2(x)
* 2 - find x1, x2 where x' = x1 + x2 and |x2| <= 0.5
* 3 - x = 2^x2 scaled by x1
*/
__asm__ __volatile__ ("fnstcw %0" : "=m" (cw) : );
saved_cw = cw;
cw &= 0xf3ff; /* force round to nearest */
cw |= 0x003f; /* turn off exceptions (see ANSI/ISO 9989 section 7.5.1) */
__asm__ __volatile__ ("fldcw %0" : : "m" (cw));
__asm__ __volatile__ ("fyl2x;fstl %2;frndint;fstl %%st(2);fsubrp;f2xm1;"
"fld1;faddp;fscale"
: "=t" (x) : "0" (x), "u" (y));
/* Restore the control word */
__asm__ __volatile__ ("fldcw %0" : : "m" (saved_cw));
/* Check for results that can not be represented (kludgy): */
if (__isinf(x) || __isnan(x))
errno = ERANGE;
return (negative) ? -x : x;
}
|
