Commit 760a49ee authored by Martin Reinecke's avatar Martin Reinecke

sync with pocketfft

parent 1783e797
...@@ -406,17 +406,10 @@ struct util // hack to avoid duplicate symbols ...@@ -406,17 +406,10 @@ struct util // hack to avoid duplicate symbols
size_t res=1; size_t res=1;
while ((n&1)==0) while ((n&1)==0)
{ res=2; n>>=1; } { res=2; n>>=1; }
for (size_t x=3; x*x<=n; x+=2)
size_t limit=size_t(sqrt(double(n)+0.01)); while ((n%x)==0)
for (size_t x=3; x<=limit; x+=2) { res=x; n/=x; }
while ((n/x)*x==n)
{
res=x;
n/=x;
limit=size_t(sqrt(double(n)+0.01));
}
if (n>1) res=n; if (n>1) res=n;
return res; return res;
} }
...@@ -427,17 +420,13 @@ struct util // hack to avoid duplicate symbols ...@@ -427,17 +420,13 @@ struct util // hack to avoid duplicate symbols
double result=0.; double result=0.;
while ((n&1)==0) while ((n&1)==0)
{ result+=2; n>>=1; } { result+=2; n>>=1; }
for (size_t x=3; x*x<=n; x+=2)
size_t limit=size_t(sqrt(double(n)+0.01)); while ((n%x)==0)
for (size_t x=3; x<=limit; x+=2) {
while ((n/x)*x==n) result+= (x<=5) ? double(x) : lfp*double(x); // penalize larger prime factors
{ n/=x;
result+= (x<=5) ? double(x) : lfp*double(x); // penalize larger prime factors }
n/=x;
limit=size_t(sqrt(double(n)+0.01));
}
if (n>1) result+=(n<=5) ? double(n) : lfp*double(n); if (n>1) result+=(n<=5) ? double(n) : lfp*double(n);
return result*double(ni); return result*double(ni);
} }
...@@ -1151,16 +1140,11 @@ template<bool bwd, typename T> void pass_all(T c[], T0 fct) ...@@ -1151,16 +1140,11 @@ template<bool bwd, typename T> void pass_all(T c[], T0 fct)
add_factor(2); add_factor(2);
swap(fact[0].fct, fact.back().fct); swap(fact[0].fct, fact.back().fct);
} }
size_t maxl = size_t(sqrt(double(len)))+1; for (size_t divisor=3; divisor*divisor<=len; divisor+=2)
for (size_t divisor=3; (len>1)&&(divisor<maxl); divisor+=2) while ((len%divisor)==0)
if ((len%divisor)==0)
{ {
while ((len%divisor)==0) add_factor(divisor);
{ len/=divisor;
add_factor(divisor);
len/=divisor;
}
maxl=size_t(sqrt(double(len)))+1;
} }
if (len>1) add_factor(len); if (len>1) add_factor(len);
} }
...@@ -1942,16 +1926,11 @@ template<typename T> void radbg(size_t ido, size_t ip, size_t l1, ...@@ -1942,16 +1926,11 @@ template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
add_factor(2); add_factor(2);
swap(fact[0].fct, fact.back().fct); swap(fact[0].fct, fact.back().fct);
} }
size_t maxl = size_t(sqrt(double(len)))+1; for (size_t divisor=3; divisor*divisor<=len; divisor+=2)
for (size_t divisor=3; (len>1)&&(divisor<maxl); divisor+=2) while ((len%divisor)==0)
if ((len%divisor)==0)
{ {
while ((len%divisor)==0) add_factor(divisor);
{ len/=divisor;
add_factor(divisor);
len/=divisor;
}
maxl=size_t(sqrt(double(len)))+1;
} }
if (len>1) add_factor(len); if (len>1) add_factor(len);
} }
...@@ -2127,8 +2106,8 @@ template<typename T0> class pocketfft_c ...@@ -2127,8 +2106,8 @@ template<typename T0> class pocketfft_c
: len(length) : len(length)
{ {
if (length==0) throw runtime_error("zero-length FFT requested"); if (length==0) throw runtime_error("zero-length FFT requested");
if ((length<50) || size_t tmp = (length<50) ? 0 : util::largest_prime_factor(length);
(double(util::largest_prime_factor(length))<=sqrt(double(length)))) if (tmp*tmp <= length)
{ {
packplan=unique_ptr<cfftp<T0>>(new cfftp<T0>(length)); packplan=unique_ptr<cfftp<T0>>(new cfftp<T0>(length));
return; return;
...@@ -2167,8 +2146,8 @@ template<typename T0> class pocketfft_r ...@@ -2167,8 +2146,8 @@ template<typename T0> class pocketfft_r
: len(length) : len(length)
{ {
if (length==0) throw runtime_error("zero-length FFT requested"); if (length==0) throw runtime_error("zero-length FFT requested");
if ((length<50) size_t tmp = (length<50) ? 0 : util::largest_prime_factor(length);
|| (double(util::largest_prime_factor(length))<=sqrt(double(length)))) if (tmp*tmp <= length)
{ {
packplan=unique_ptr<rfftp<T0>>(new rfftp<T0>(length)); packplan=unique_ptr<rfftp<T0>>(new rfftp<T0>(length));
return; return;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment