On 24.12.2014 13:07, Sturla Molden wrote:
> On 24/12/14 04:33, Robert McGibbon wrote:
>> Alex Griffing pointed out on github that this feature was recently added
>> to scipy in https://github.com/scipy/scipy/pull/3144. Sweet!
>
>
> I use different padsize search than the one in SciPy. It would be
> interesting to see which is faster.
>
hm this is a brute force search, probably fast enough but slower than
scipy's code (if it also were cython)
I also ported it to C a while back, so that could be used for numpy if
speed is an issue. Code attached.
static inline unsigned long get_msb(unsigned long a)
{
/* dumb integer msb, (63 - __builtin_clzl(a)) would be faster */
unsigned long msb = 0;
while (a >>= 1) {
msb++;
}
return msb;
}
/* ---------------------------------------------------------------------------*/
/**
* @brief get next 5-smooth number next to a certain value
* @param a lower bound of result
*
* 5-smooth number is a number with only prime factors 2, 3 or 5.
* This can be used to get minimal padding sizes that still allow efficient fft
* transforms.
*/
/* ---------------------------------------------------------------------------*/
size_t visir_get_next_regular(size_t a)
{
/* fftw can also deal efficiently with factors of 7 and 13 but it needs
* testing that the speed of these algorithms will not cancel out the gain
* of smaller input sizes */
if (a <= 6)
return a;
/* is already power of 2 */
if ((a & (a - 1)) == 0)
return a;
/* overflow */
if (5 > SIZE_MAX / a)
return a;
size_t match = SIZE_MAX;
size_t p5 = 1;
while(p5 < a) {
size_t p35 = p5;
while (p35 < a) {
/* ceil division */
size_t quotient = a % p35 != 0 ? a / p35 + 1 : a / p35;
/* next power of two of quotient */
size_t p2 = 2 << get_msb(quotient - 1);
size_t n = p2 * p35;
if (n == a)
return n;
else if (n < match)
match = n;
p35 *= 3;
if (p35 == a)
return p35;
}
if (p35 < match)
match = p35;
p5 *= 5;
if (p5 == a)
return p5;
}
if (p5 < match)
match = p5;
return match;
}
_______________________________________________
NumPy-Discussion mailing list
[email protected]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
