From 760a49ee73a94161cbb8a772b807ae77dacb4e28 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Fri, 21 Jun 2019 15:41:46 +0200 Subject: [PATCH] sync with pocketfft --- pocketfft_hdronly.h | 63 +++++++++++++++------------------------------ 1 file changed, 21 insertions(+), 42 deletions(-) diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h index 2662a01..359375c 100644 --- a/pocketfft_hdronly.h +++ b/pocketfft_hdronly.h @@ -406,17 +406,10 @@ struct util // hack to avoid duplicate symbols size_t res=1; while ((n&1)==0) { res=2; n>>=1; } - - size_t limit=size_t(sqrt(double(n)+0.01)); - for (size_t x=3; x<=limit; x+=2) - while ((n/x)*x==n) - { - res=x; - n/=x; - limit=size_t(sqrt(double(n)+0.01)); - } + for (size_t x=3; x*x<=n; x+=2) + while ((n%x)==0) + { res=x; n/=x; } if (n>1) res=n; - return res; } @@ -427,17 +420,13 @@ struct util // hack to avoid duplicate symbols double result=0.; while ((n&1)==0) { result+=2; n>>=1; } - - size_t limit=size_t(sqrt(double(n)+0.01)); - 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; - limit=size_t(sqrt(double(n)+0.01)); - } + for (size_t x=3; x*x<=n; x+=2) + while ((n%x)==0) + { + result+= (x<=5) ? double(x) : lfp*double(x); // penalize larger prime factors + n/=x; + } if (n>1) result+=(n<=5) ? double(n) : lfp*double(n); - return result*double(ni); } @@ -1151,16 +1140,11 @@ template void pass_all(T c[], T0 fct) add_factor(2); swap(fact[0].fct, fact.back().fct); } - size_t maxl = size_t(sqrt(double(len)))+1; - for (size_t divisor=3; (len>1)&&(divisor1) add_factor(len); } @@ -1942,16 +1926,11 @@ template void radbg(size_t ido, size_t ip, size_t l1, add_factor(2); swap(fact[0].fct, fact.back().fct); } - size_t maxl = size_t(sqrt(double(len)))+1; - for (size_t divisor=3; (len>1)&&(divisor1) add_factor(len); } @@ -2127,8 +2106,8 @@ template class pocketfft_c : len(length) { if (length==0) throw runtime_error("zero-length FFT requested"); - if ((length<50) || - (double(util::largest_prime_factor(length))<=sqrt(double(length)))) + size_t tmp = (length<50) ? 0 : util::largest_prime_factor(length); + if (tmp*tmp <= length) { packplan=unique_ptr>(new cfftp(length)); return; @@ -2167,8 +2146,8 @@ template class pocketfft_r : len(length) { if (length==0) throw runtime_error("zero-length FFT requested"); - if ((length<50) - || (double(util::largest_prime_factor(length))<=sqrt(double(length)))) + size_t tmp = (length<50) ? 0 : util::largest_prime_factor(length); + if (tmp*tmp <= length) { packplan=unique_ptr>(new rfftp(length)); return; -- GitLab