@@ -94,6 +94,41 @@ static const npy_uint64 MAGIC64[] = {0x5555555555555555ull, 0x3333333333333333ul
94
94
* classification instead of calling npy_isinf/npy_isnan: we should have some
95
95
* macros for this, though, instead of doing it manually
96
96
*/
97
+ #ifndef HAVE_LOG2
98
+ NPY_INPLACE double npy_log2 (double x )
99
+ {
100
+ #ifdef HAVE_FREXP
101
+ if (!npy_isfinite (x ) || x <= 0. ) {
102
+ /* special value result */
103
+ return npy_log (x );
104
+ }
105
+ else {
106
+ /*
107
+ * fallback implementation copied from python3.4 math.log2
108
+ * provides int(log(2**i)) == i for i 1-64 in default rounding mode.
109
+ *
110
+ * We want log2(m * 2**e) == log(m) / log(2) + e. Care is needed when
111
+ * x is just greater than 1.0: in that case e is 1, log(m) is negative,
112
+ * and we get significant cancellation error from the addition of
113
+ * log(m) / log(2) to e. The slight rewrite of the expression below
114
+ * avoids this problem.
115
+ */
116
+ int e ;
117
+ double m = frexp (x , & e );
118
+ if (x >= 1.0 ) {
119
+ return log (2.0 * m ) / log (2.0 ) + (e - 1 );
120
+ }
121
+ else {
122
+ return log (m ) / log (2.0 ) + e ;
123
+ }
124
+ }
125
+ #else
126
+ /* does not provide int(log(2**i)) == i */
127
+ return NPY_LOG2E * npy_log (x );
128
+ #endif
129
+ }
130
+ #endif
131
+
97
132
/*
98
133
*
99
134
* sin, cos, tan
@@ -140,10 +175,12 @@ static const npy_uint64 MAGIC64[] = {0x5555555555555555ull, 0x3333333333333333ul
140
175
#define WORKAROUND_APPLE_TRIG_BUG 0
141
176
#endif
142
177
178
+ /* mandatory C99 functions */
179
+
143
180
/**begin repeat1
144
181
* #kind = sin,cos,tan,sinh,cosh,tanh,fabs,floor,ceil,rint,trunc,sqrt,log10,
145
- * log,exp,expm1,asin,acos,atan,asinh,acosh,atanh,log1p,exp2,log2 #
146
- * #TRIG_WORKAROUND = WORKAROUND_APPLE_TRIG_BUG*3, 0*22 #
182
+ * log,exp,expm1,asin,acos,atan,asinh,acosh,atanh,log1p,exp2#
183
+ * #TRIG_WORKAROUND = WORKAROUND_APPLE_TRIG_BUG*3, 0*21 #
147
184
*/
148
185
NPY_INPLACE @type @ npy_ @kind @@c @(@type @ x )
149
186
{
@@ -155,6 +192,19 @@ NPY_INPLACE @type@ npy_@kind@@c@(@type@ x)
155
192
return NPY__FP_SFX (@kind @)(x );
156
193
}
157
194
195
+ /**end repeat1**/
196
+ /* Optional C99 functions */
197
+ /**begin repeat1
198
+ * #kind = log2#
199
+ * #KIND = LOG2#
200
+ */
201
+ #ifdef HAVE_ @KIND @@C @
202
+ NPY_INPLACE @type @ npy_ @kind @@c @(@type @ x )
203
+ {
204
+ return NPY__FP_SFX (@kind @)(x );
205
+ }
206
+ #endif
207
+
158
208
/**end repeat1**/
159
209
160
210
#undef WORKAROUND_APPLE_TRIG_BUG
0 commit comments