← Back to branch summary

~ubuntu-branches/ubuntu/trusty/mingw-w64/trusty

« back to all changes in this revision

Viewing changes to mingw-w64-crt/math/powif.c

  • Committer: Bazaar Package Importer
  • Author(s): Robert Millan
  • Date: 2010-01-25 22:58:25 UTC
  • Revision ID: james.westby@ubuntu.com-20100125225825-eomk8xw36mwwg1yl
Tags: upstream-0~20100125
Import upstream version 0~20100125

Show diffs side-by-side

added added

removed removed

Lines of Context:
mingw-w64-crt/math/powif.c
 
1
/**
 
2
 * This file has no copyright assigned and is placed in the Public Domain.
 
3
 * This file is part of the w64 mingw-runtime package.
 
4
 * No warranty is given; refer to the file DISCLAIMER.PD within this package.
 
5
 */
 
6
#include "cephes_mconf.h"
 
7
 
 
8
#ifndef _SET_ERRNO
 
9
#define _SET_ERRNO(x)
 
10
#endif
 
11
float __powif (float x, int nn);
 
12
 
 
13
float __powif (float x, int nn)
 
14
{
 
15
        int n, e, sign, asign, lx;
 
16
        float w, y, s;
 
17
 
 
18
        /* See pow.c for these tests.  */
 
19
        if (x == 0.0F)
 
20
        {
 
21
                if (nn == 0)
 
22
                        return (1.0F );
 
23
                else if (nn < 0)
 
24
                        return (INFINITYF);
 
25
                else
 
26
                {
 
27
                        if (nn & 1)
 
28
                                return (x);
 
29
                        else
 
30
                                return (0.0);
 
31
                }
 
32
        }
 
33
 
 
34
        if (nn == 0)
 
35
                return (1.0);
 
36
 
 
37
        if (nn == -1)
 
38
                return (1.0/x);
 
39
 
 
40
        if (x < 0.0)
 
41
        {
 
42
                asign = -1;
 
43
                x = -x;
 
44
        }
 
45
        else
 
46
                asign = 0;
 
47
 
 
48
        if (nn < 0)
 
49
        {
 
50
                sign = -1;
 
51
                n = -nn;
 
52
        }
 
53
        else
 
54
        {
 
55
                sign = 1;
 
56
                n = nn;
 
57
        }
 
58
 
 
59
        /* Even power will be positive. */
 
60
        if ((n & 1) == 0)
 
61
                asign = 0;
 
62
 
 
63
        /* Overflow detection */
 
64
 
 
65
        /* Calculate approximate logarithm of answer */
 
66
        s = frexpf(x, &lx);
 
67
        e = (lx - 1)*n;
 
68
        if ((e == 0) || (e > 64) || (e < -64))
 
69
        {
 
70
                s = (s - 7.0710678118654752e-1) / (s +  7.0710678118654752e-1);
 
71
                s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2F;
 
72
        }
 
73
        else
 
74
        {
 
75
                s = LOGE2F * e;
 
76
        }
 
77
 
 
78
        if (s > MAXLOGF)
 
79
        {
 
80
                mtherr("__powif", OVERFLOW);
 
81
                _SET_ERRNO(ERANGE);
 
82
                y = INFINITYF;
 
83
                goto done;
 
84
        }
 
85
 
 
86
#if DENORMAL
 
87
        if (s < MINLOGF)
 
88
        {
 
89
                y = 0.0;
 
90
                goto done;
 
91
        }
 
92
 
 
93
        /* Handle tiny denormal answer, but with less accuracy
 
94
         * since roundoff error in 1.0/x will be amplified.
 
95
         * The precise demarcation should be the gradual underflow threshold.
 
96
         */
 
97
        if ((s < (-MAXLOGF+2.0)) && (sign < 0))
 
98
        {
 
99
                x = 1.0/x;
 
100
                sign = -sign;
 
101
        }
 
102
#else
 
103
        /* do not produce denormal answer */
 
104
        if (s < -MAXLOGF)
 
105
                return (0.0);
 
106
#endif
 
107
 
 
108
        /* First bit of the power */
 
109
        if (n & 1)
 
110
                y = x;
 
111
        else
 
112
                y = 1.0;
 
113
 
 
114
        w = x;
 
115
        n >>= 1;
 
116
        while (n)
 
117
        {
 
118
                w = w * w;      /* arg to the 2-to-the-kth power */
 
119
                if (n & 1)      /* if that bit is set, then include in product */
 
120
                        y *= w;
 
121
                n >>= 1;
 
122
        }
 
123
 
 
124
        if (sign < 0)
 
125
                y = 1.0/y;
 
126
 
 
127
done:
 
128
 
 
129
        if (asign)
 
130
        {
 
131
                /* odd power of negative number */
 
132
                if (y == 0.0)
 
133
                        y = NEGZEROF;
 
134
                else
 
135
                        y = -y;
 
136
        }
 
137
        return (y);
 
138
}
 
139