#define SIMD_PARSE
#define SIMD_ROOT_BUILD
#define SPARSE_FORWARD_TOP
#define REALFFT
#define SMALL_ROOTS
#define VECTOR_REAL_PRODUCT
#define ALGEBRA_PRODUCT
#pragma GCC optimize("O3,unroll-loops,omit-frame-pointer")
#pragma GCC target("arch=skylake")
#ifndef DUCK_FASTIO_H
#define DUCK_FASTIO_H
typedef unsigned long duck_u64;
typedef long duck_i64;
typedef struct {
duck_u64 abi_version;
const char *stdin_ptr;
duck_u64 stdin_size;
char *stdout_ptr;
duck_u64 stdout_limit;
duck_u64 stdout_size;
char *stderr_ptr;
duck_u64 stderr_limit;
duck_u64 stderr_size;
const char *ib_ptr;
duck_u64 ib_limit;
char *ob_ptr;
duck_u64 ob_limit;
duck_u64 tsc_frequency;
} __attribute__((packed)) DuckInfo;
static __attribute__((always_inline)) inline DuckInfo *duck_info(long argc, char **argv) {
char **p = argv + argc + 1;
while (*p) ++p;
duck_u64 *aux = (duck_u64 *)(p + 1);
while (aux[0]) {
if (aux[0] == 0x6b637564UL) return (DuckInfo *)aux[1];
aux += 2;
}
return (DuckInfo *)0;
}
static __attribute__((always_inline)) inline duck_u64 duck_read_u64(const char **cursor) {
const char *p = *cursor;
while ((unsigned char)(*p - '0') > 9) ++p;
duck_u64 value = 0;
do {
value = value * 10 + (unsigned char)(*p - '0');
++p;
} while ((unsigned char)(*p - '0') <= 9);
*cursor = p;
return value;
}
static __attribute__((always_inline)) inline duck_i64 duck_read_i64(const char **cursor) {
const char *p = *cursor;
while (*p != '-' && (unsigned char)(*p - '0') > 9) ++p;
int negative = *p == '-';
p += negative;
duck_u64 value = 0;
do {
value = value * 10 + (unsigned char)(*p - '0');
++p;
} while ((unsigned char)(*p - '0') <= 9);
*cursor = p;
return negative ? -(duck_i64)value : (duck_i64)value;
}
static __attribute__((always_inline)) inline char *duck_write_u64(char *out, duck_u64 value) {
char tmp[24];
unsigned n = 0;
do {
tmp[n++] = (char)('0' + value % 10);
value /= 10;
} while (value);
do *out++ = tmp[--n]; while (n);
return out;
}
static __attribute__((always_inline)) inline char *duck_write_i64(char *out, duck_i64 value) {
if (value < 0) {
*out++ = '-';
return duck_write_u64(out, (duck_u64)(-value));
}
return duck_write_u64(out, (duck_u64)value);
}
static __attribute__((always_inline, noreturn)) inline void duck_exit(void) {
__asm__ volatile("mov $60,%%eax;xor %%edi,%%edi;syscall" ::: "rax", "rdi", "rcx", "r11", "memory");
__builtin_unreachable();
}
#endif
#include <immintrin.h>
#include <math.h>
#ifdef PARALLEL_FORWARD
#include <pthread.h>
#endif
#ifdef LOCAL_TEST
#include <stdio.h>
#endif
struct C { double r, i; };
#ifdef BASE5
enum { N=1<<19, H=N>>1, DIGITS=1000000, LIMBS=DIGITS/5, BASE=100000 };
#else
enum { N=1<<19, H=N>>1, DIGITS=1000000, LIMBS=DIGITS/4, BASE=10000 };
#endif
alignas(32) static C z[N];
alignas(32) static unsigned digit4[10000];
#ifdef SMALL_ROOTS
alignas(32) static C roots16[4],roots64[16],roots256[64],roots1024[256],roots4096[1024],roots16384[4096];
alignas(32) static C roots16_2[4],roots64_2[16],roots256_2[64],roots1024_2[256],roots4096_2[1024],roots16384_2[4096];
alignas(32) static C roots16_3[4],roots64_3[16],roots256_3[64],roots1024_3[256],roots4096_3[1024],roots16384_3[4096];
#ifdef PRODUCT_ROOT_TABLE
alignas(32) static C product_roots[4096];
#endif
#ifdef SIMD_ROOT_BUILD
static inline __m256d root_cmul(__m256d a,__m256d b){
return _mm256_fmaddsub_pd(_mm256_movedup_pd(a),b,_mm256_mul_pd(_mm256_permute_pd(a,15),_mm256_permute_pd(b,5)));
}
#endif
static void build_small_roots(){
C *tabs[6]={roots16,roots64,roots256,roots1024,roots4096,roots16384},*tabs2[6]={roots16_2,roots64_2,roots256_2,roots1024_2,roots4096_2,roots16384_2},*tabs3[6]={roots16_3,roots64_3,roots256_3,roots1024_3,roots4096_3,roots16384_3};unsigned lens[6]={16,64,256,1024,4096,16384};
#ifdef SIMD_ROOT_BUILD
for(unsigned t=0;t<6;++t){
double a=-2.0*3.141592653589793238462643383279502884/lens[t];C step={cos(a),sin(a)},step2={step.r*step.r-step.i*step.i,2*step.r*step.i};
__m256d rw=_mm256_setr_pd(1,0,step.r,step.i),advance=_mm256_setr_pd(step2.r,step2.i,step2.r,step2.i);
for(unsigned j=0;j<lens[t]/4;j+=2){
if(j&&(j&63u)==0){__m256d sq=_mm256_mul_pd(rw,rw),mag=_mm256_add_pd(_mm256_movedup_pd(sq),_mm256_permute_pd(sq,15));rw=_mm256_mul_pd(rw,_mm256_sub_pd(_mm256_set1_pd(1.5),_mm256_mul_pd(_mm256_set1_pd(.5),mag)));}
__m256d rw2=root_cmul(rw,rw),rw3=root_cmul(rw2,rw);
_mm256_store_pd((double*)(tabs[t]+j),rw);_mm256_store_pd((double*)(tabs2[t]+j),rw2);_mm256_store_pd((double*)(tabs3[t]+j),rw3);
rw=root_cmul(rw,advance);
}
}
#else
for(unsigned t=0;t<6;++t){double a=-2.0*3.141592653589793238462643383279502884/lens[t];C step={cos(a),sin(a)},w={1,0};for(unsigned j=0;j<lens[t]/4;++j){if(j&&(j&63u)==0){double s=1.5-.5*(w.r*w.r+w.i*w.i);w.r*=s;w.i*=s;}C w2={w.r*w.r-w.i*w.i,2*w.r*w.i};tabs[t][j]=w;tabs2[t][j]=w2;tabs3[t][j]={w2.r*w.r-w2.i*w.i,w2.r*w.i+w2.i*w.r};w={w.r*step.r-w.i*step.i,w.r*step.i+w.i*step.r};}}
#endif
#ifdef PRODUCT_ROOT_TABLE
{double a=-2.0*3.141592653589793238462643383279502884/N;C step={cos(a),sin(a)},w={1,0};for(unsigned j=0;j<4096;++j){if(j&&(j&63u)==0){double s=1.5-.5*(w.r*w.r+w.i*w.i);w.r*=s;w.i*=s;}product_roots[j]=w;w={w.r*step.r-w.i*step.i,w.r*step.i+w.i*step.r};}}
#endif
}
#endif
#ifdef ALL_ROOTS
alignas(32) static C roots_all[87380];
static void build_all_roots(){
for(unsigned len=16;len<=H;len<<=2){unsigned q=len/4,off=(q-4)/3;double a=-2.0*3.141592653589793238462643383279502884/len;C step={cos(a),sin(a)},w={1,0};
for(unsigned j=0;j<q;++j){if(j&&(j&63u)==0){double s=1.5-.5*(w.r*w.r+w.i*w.i);w.r*=s;w.i*=s;}roots_all[off+j]=w;w={w.r*step.r-w.i*step.i,w.r*step.i+w.i*step.r};}}
}
#endif
static inline __m256d cmul(__m256d a,__m256d b){
__m256d ar=_mm256_movedup_pd(a),ai=_mm256_permute_pd(a,15),bs=_mm256_permute_pd(b,5);
return _mm256_fmaddsub_pd(ar,b,_mm256_mul_pd(ai,bs));
}
static inline void square_cube(__m256d w,__m256d *w2,__m256d *w3){
__m256d ar=_mm256_movedup_pd(w),ai=_mm256_permute_pd(w,15),sq=_mm256_mul_pd(w,w);
__m256d rr=_mm256_movedup_pd(sq),ii=_mm256_permute_pd(sq,15),ri=_mm256_mul_pd(ar,ai);
*w2=_mm256_blend_pd(_mm256_sub_pd(rr,ii),_mm256_add_pd(ri,ri),10);
__m256d threeii=_mm256_add_pd(ii,_mm256_add_pd(ii,ii)),threerr=_mm256_add_pd(rr,_mm256_add_pd(rr,rr));
*w3=_mm256_blend_pd(_mm256_mul_pd(ar,_mm256_sub_pd(rr,threeii)),_mm256_mul_pd(ai,_mm256_sub_pd(threerr,ii)),10);
}
static inline __m256d conjv(__m256d a){
return _mm256_xor_pd(a,_mm256_castsi256_pd(_mm256_setr_epi64x(0,0x8000000000000000ULL,0,0x8000000000000000ULL)));
}
static inline __m256d minus_i(__m256d a){
return _mm256_xor_pd(_mm256_permute_pd(a,5),_mm256_castsi256_pd(_mm256_setr_epi64x(0,0x8000000000000000ULL,0,0x8000000000000000ULL)));
}
static inline __m256d plus_i(__m256d a){
return _mm256_xor_pd(_mm256_permute_pd(a,5),_mm256_castsi256_pd(_mm256_setr_epi64x(0x8000000000000000ULL,0,0x8000000000000000ULL,0)));
}
static inline C mul(C a,C b){return {a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r};}
static inline C norm(C a){double s=1.5-.5*(a.r*a.r+a.i*a.i);return {a.r*s,a.i*s};}
static void radix2_forward(){
const double ang=-2.0*3.141592653589793238462643383279502884/N;
C step={cos(ang),sin(ang)},w={1,0};
for(unsigned j=0;j<H;j+=2){
if(j&&(j&63u)==0)w=norm(w);
C w0=w,w1=mul(w0,step);w=mul(w1,step);
__m256d a=_mm256_load_pd((double*)(z+j)),b=_mm256_load_pd((double*)(z+H+j));
_mm256_store_pd((double*)(z+j),_mm256_add_pd(a,b));
_mm256_store_pd((double*)(z+H+j),cmul(_mm256_sub_pd(a,b),_mm256_setr_pd(w0.r,w0.i,w1.r,w1.i)));
}
}
static void radix2_inverse(){
const double ang=-2.0*3.141592653589793238462643383279502884/N;
C step={cos(ang),sin(ang)},w={1,0};
const __m256d scale=_mm256_set1_pd(1.0/N);
for(unsigned j=0;j<H;j+=2){
if(j&&(j&63u)==0)w=norm(w);
C w0=w,w1=mul(w,step);w=mul(w1,step);
__m256d a=_mm256_load_pd((double*)(z+j));
__m256d b=cmul(_mm256_load_pd((double*)(z+H+j)),_mm256_setr_pd(w0.r,-w0.i,w1.r,-w1.i));
_mm256_store_pd((double*)(z+j),_mm256_mul_pd(_mm256_add_pd(a,b),scale));
_mm256_store_pd((double*)(z+H+j),_mm256_mul_pd(_mm256_sub_pd(a,b),scale));
}
}
static __attribute__((unused)) void forward_old(){
radix2_forward();
for(unsigned len=N/2;len>=8;len>>=2){
unsigned q=len>>2;
double ang=-2.0*3.141592653589793238462643383279502884/len;
C step={cos(ang),sin(ang)};
for(unsigned p=0;p<N;p+=len){
C w={1,0};
for(unsigned j=0;j<q;j+=2){
if(j&&(j&63u)==0)w=norm(w);
C w0=w,w1=mul(w,step);w=mul(w1,step);
__m256d rw=_mm256_setr_pd(w0.r,w0.i,w1.r,w1.i);
__m256d rw2=cmul(rw,rw),rw3=cmul(rw2,rw);
__m256d a=_mm256_load_pd((double*)(z+p+j));
__m256d b=_mm256_load_pd((double*)(z+p+q+j));
__m256d c=_mm256_load_pd((double*)(z+p+2*q+j));
__m256d d=_mm256_load_pd((double*)(z+p+3*q+j));
__m256d t0=_mm256_add_pd(a,c),t1=_mm256_sub_pd(a,c);
__m256d t2=_mm256_add_pd(b,d),t3=minus_i(_mm256_sub_pd(b,d));
_mm256_store_pd((double*)(z+p+j),_mm256_add_pd(t0,t2));
_mm256_store_pd((double*)(z+p+q+j),cmul(_mm256_add_pd(t1,t3),rw));
_mm256_store_pd((double*)(z+p+2*q+j),cmul(_mm256_sub_pd(t0,t2),rw2));
_mm256_store_pd((double*)(z+p+3*q+j),cmul(_mm256_sub_pd(t1,t3),rw3));
}
}
}
for(unsigned p=0;p<N;p+=4){
C a=z[p],b=z[p+1],c=z[p+2],d=z[p+3];
C t0={a.r+c.r,a.i+c.i},t1={a.r-c.r,a.i-c.i};
C t2={b.r+d.r,b.i+d.i},x={b.r-d.r,b.i-d.i},t3={x.i,-x.r};
z[p]={t0.r+t2.r,t0.i+t2.i};z[p+1]={t1.r+t3.r,t1.i+t3.i};
z[p+2]={t0.r-t2.r,t0.i-t2.i};z[p+3]={t1.r-t3.r,t1.i-t3.i};
}
}
static __attribute__((unused)) void inverse_old(){
for(unsigned p=0;p<N;p+=4){
C a=z[p],b=z[p+1],c=z[p+2],d=z[p+3];
C t0={a.r+c.r,a.i+c.i},t1={a.r-c.r,a.i-c.i};
C t2={b.r+d.r,b.i+d.i},x={b.r-d.r,b.i-d.i},t3={-x.i,x.r};
z[p]={t0.r+t2.r,t0.i+t2.i};z[p+1]={t1.r+t3.r,t1.i+t3.i};
z[p+2]={t0.r-t2.r,t0.i-t2.i};z[p+3]={t1.r-t3.r,t1.i-t3.i};
}
for(unsigned len=16;len<=N/2;len<<=2){
unsigned q=len>>2;
double ang=-2.0*3.141592653589793238462643383279502884/len;
C step={cos(ang),sin(ang)};
for(unsigned p=0;p<N;p+=len){
C w={1,0};
for(unsigned j=0;j<q;j+=2){
if(j&&(j&63u)==0)w=norm(w);
C w0=w,w1=mul(w,step);w=mul(w1,step);
__m256d rw=conjv(_mm256_setr_pd(w0.r,w0.i,w1.r,w1.i));
__m256d rw2=cmul(rw,rw),rw3=cmul(rw2,rw);
__m256d a=_mm256_load_pd((double*)(z+p+j));
__m256d b=cmul(_mm256_load_pd((double*)(z+p+q+j)),rw);
__m256d c=cmul(_mm256_load_pd((double*)(z+p+2*q+j)),rw2);
__m256d d=cmul(_mm256_load_pd((double*)(z+p+3*q+j)),rw3);
__m256d t0=_mm256_add_pd(a,c),t1=_mm256_sub_pd(a,c);
__m256d t2=_mm256_add_pd(b,d),t3=plus_i(_mm256_sub_pd(b,d));
_mm256_store_pd((double*)(z+p+j),_mm256_add_pd(t0,t2));
_mm256_store_pd((double*)(z+p+q+j),_mm256_add_pd(t1,t3));
_mm256_store_pd((double*)(z+p+2*q+j),_mm256_sub_pd(t0,t2));
_mm256_store_pd((double*)(z+p+3*q+j),_mm256_sub_pd(t1,t3));
}
}
}
radix2_inverse();
}
template<unsigned M> static inline void forward_block(C *x){
constexpr unsigned q=M/4;
constexpr double ang=-2.0*3.141592653589793238462643383279502884/M;
C step={cos(ang),sin(ang)},w={1,0};
for(unsigned j=0;j<q;j+=2){
if(j&&(j&63u)==0)w=norm(w);
C w0=w,w1=mul(w,step);w=mul(w1,step);
__m256d rw=_mm256_setr_pd(w0.r,w0.i,w1.r,w1.i),rw2=cmul(rw,rw),rw3=cmul(rw2,rw);
__m256d a=_mm256_load_pd((double*)(x+j)),b=_mm256_load_pd((double*)(x+q+j));
__m256d c=_mm256_load_pd((double*)(x+2*q+j)),d=_mm256_load_pd((double*)(x+3*q+j));
__m256d t0=_mm256_add_pd(a,c),t1=_mm256_sub_pd(a,c),t2=_mm256_add_pd(b,d),t3=minus_i(_mm256_sub_pd(b,d));
_mm256_store_pd((double*)(x+j),_mm256_add_pd(t0,t2));
_mm256_store_pd((double*)(x+q+j),cmul(_mm256_add_pd(t1,t3),rw));
_mm256_store_pd((double*)(x+2*q+j),cmul(_mm256_sub_pd(t0,t2),rw2));
_mm256_store_pd((double*)(x+3*q+j),cmul(_mm256_sub_pd(t1,t3),rw3));
}
}
template<unsigned M> static inline void inverse_block(C *x){
constexpr unsigned q=M/4;
constexpr double ang=-2.0*3.141592653589793238462643383279502884/M;
C step={cos(ang),sin(ang)},w={1,0};
for(unsigned j=0;j<q;j+=2){
if(j&&(j&63u)==0)w=norm(w);
C w0=w,w1=mul(w,step);w=mul(w1,step);
__m256d rw=conjv(_mm256_setr_pd(w0.r,w0.i,w1.r,w1.i)),rw2=cmul(rw,rw),rw3=cmul(rw2,rw);
__m256d a=_mm256_load_pd((double*)(x+j)),b=cmul(_mm256_load_pd((double*)(x+q+j)),rw);
__m256d c=cmul(_mm256_load_pd((double*)(x+2*q+j)),rw2),d=cmul(_mm256_load_pd((double*)(x+3*q+j)),rw3);
__m256d t0=_mm256_add_pd(a,c),t1=_mm256_sub_pd(a,c),t2=_mm256_add_pd(b,d),t3=plus_i(_mm256_sub_pd(b,d));
_mm256_store_pd((double*)(x+j),_mm256_add_pd(t0,t2));
_mm256_store_pd((double*)(x+q+j),_mm256_add_pd(t1,t3));
_mm256_store_pd((double*)(x+2*q+j),_mm256_sub_pd(t0,t2));
_mm256_store_pd((double*)(x+3*q+j),_mm256_sub_pd(t1,t3));
}
}
template<unsigned M> static inline void inverse_block_round(C *x){
constexpr unsigned q=M/4;
constexpr double ang=-2.0*3.141592653589793238462643383279502884/M;
C step={cos(ang),-sin(ang)},step2=mul(step,step);
__m256d rw=_mm256_setr_pd(1,0,step.r,step.i),advance=_mm256_setr_pd(step2.r,step2.i,step2.r,step2.i);
const __m256d magic=_mm256_set1_pd(0x1p52);
const __m256i bias=_mm256_set1_epi64x(0x4330000000000000ULL);
for(unsigned j=0;j<q;j+=2){
if(j&&(j&63u)==0){__m256d sq=_mm256_mul_pd(rw,rw),mag=_mm256_add_pd(_mm256_movedup_pd(sq),_mm256_permute_pd(sq,15));rw=_mm256_mul_pd(rw,_mm256_sub_pd(_mm256_set1_pd(1.5),_mm256_mul_pd(_mm256_set1_pd(.5),mag)));}
__m256d rw2=cmul(rw,rw),rw3=cmul(rw2,rw);
__m256d a=_mm256_load_pd((double*)(x+j)),b=cmul(_mm256_load_pd((double*)(x+q+j)),rw);
__m256d c=cmul(_mm256_load_pd((double*)(x+2*q+j)),rw2),d=cmul(_mm256_load_pd((double*)(x+3*q+j)),rw3);
__m256d t0=_mm256_add_pd(a,c),t1=_mm256_sub_pd(a,c),t2=_mm256_add_pd(b,d),t3=plus_i(_mm256_sub_pd(b,d));
__m256d y0=_mm256_add_pd(t0,t2),y1=_mm256_add_pd(t1,t3),y2=_mm256_sub_pd(t0,t2),y3=_mm256_sub_pd(t1,t3);
_mm256_store_si256((__m256i*)(x+j),_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(y0,magic)),bias));
_mm256_store_si256((__m256i*)(x+q+j),_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(y1,magic)),bias));
_mm256_store_si256((__m256i*)(x+2*q+j),_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(y2,magic)),bias));
_mm256_store_si256((__m256i*)(x+3*q+j),_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(y3,magic)),bias));
rw=cmul(rw,advance);
}
}
#ifdef ALL_ROOTS
template<unsigned M> static inline C *all_roots(){return roots_all+(M/4-4)/3;}
template<unsigned M> static inline void forward_block_alltab(C*x){
constexpr unsigned q=M/4;C *rt=all_roots<M>();
for(unsigned j=0;j<q;j+=2){__m256d rw=_mm256_load_pd((double*)(rt+j)),rw2=cmul(rw,rw),rw3=cmul(rw2,rw);
__m256d a=_mm256_load_pd((double*)(x+j)),b=_mm256_load_pd((double*)(x+q+j)),c=_mm256_load_pd((double*)(x+2*q+j)),d=_mm256_load_pd((double*)(x+3*q+j));
__m256d t0=_mm256_add_pd(a,c),t1=_mm256_sub_pd(a,c),t2=_mm256_add_pd(b,d),t3=minus_i(_mm256_sub_pd(b,d));
_mm256_store_pd((double*)(x+j),_mm256_add_pd(t0,t2));_mm256_store_pd((double*)(x+q+j),cmul(_mm256_add_pd(t1,t3),rw));
_mm256_store_pd((double*)(x+2*q+j),cmul(_mm256_sub_pd(t0,t2),rw2));_mm256_store_pd((double*)(x+3*q+j),cmul(_mm256_sub_pd(t1,t3),rw3));}
}
template<unsigned M> static inline void inverse_block_alltab(C*x){
constexpr unsigned q=M/4;C *rt=all_roots<M>();
for(unsigned j=0;j<q;j+=2){__m256d rw=conjv(_mm256_load_pd((double*)(rt+j))),rw2=cmul(rw,rw),rw3=cmul(rw2,rw);
__m256d a=_mm256_load_pd((double*)(x+j)),b=cmul(_mm256_load_pd((double*)(x+q+j)),rw),c=cmul(_mm256_load_pd((double*)(x+2*q+j)),rw2),d=cmul(_mm256_load_pd((double*)(x+3*q+j)),rw3);
__m256d t0=_mm256_add_pd(a,c),t1=_mm256_sub_pd(a,c),t2=_mm256_add_pd(b,d),t3=plus_i(_mm256_sub_pd(b,d));
_mm256_store_pd((double*)(x+j),_mm256_add_pd(t0,t2));_mm256_store_pd((double*)(x+q+j),_mm256_add_pd(t1,t3));
_mm256_store_pd((double*)(x+2*q+j),_mm256_sub_pd(t0,t2));_mm256_store_pd((double*)(x+3*q+j),_mm256_sub_pd(t1,t3));}
}
#endif
#ifdef SMALL_ROOTS
template<unsigned M> static inline C *small_roots();
template<unsigned M> static inline C *small_roots2();
template<unsigned M> static inline C *small_roots3();
template<> inline C *small_roots<16>(){return roots16;}
template<> inline C *small_roots<64>(){return roots64;}
template<> inline C *small_roots<256>(){return roots256;}
template<> inline C *small_roots<1024>(){return roots1024;}
template<> inline C *small_roots<4096>(){return roots4096;}
template<> inline C *small_roots<16384>(){return roots16384;}
template<> inline C *small_roots2<16>(){return roots16_2;}
template<> inline C *small_roots2<64>(){return roots64_2;}
template<> inline C *small_roots2<256>(){return roots256_2;}
template<> inline C *small_roots2<1024>(){return roots1024_2;}
template<> inline C *small_roots2<4096>(){return roots4096_2;}
template<> inline C *small_roots2<16384>(){return roots16384_2;}
template<> inline C *small_roots3<16>(){return roots16_3;}
template<> inline C *small_roots3<64>(){return roots64_3;}
template<> inline C *small_roots3<256>(){return roots256_3;}
template<> inline C *small_roots3<1024>(){return roots1024_3;}
template<> inline C *small_roots3<4096>(){return roots4096_3;}
template<> inline C *small_roots3<16384>(){return roots16384_3;}
template<unsigned M> static inline void forward_block_tab(C*x){
constexpr unsigned q=M/4;C *rt=small_roots<M>(),*rt2=small_roots2<M>(),*rt3=small_roots3<M>();
for(unsigned j=0;j<q;j+=2){__m256d rw=_mm256_load_pd((double*)(rt+j)),rw2=_mm256_load_pd((double*)(rt2+j)),rw3=_mm256_load_pd((double*)(rt3+j));
__m256d a=_mm256_load_pd((double*)(x+j)),b=_mm256_load_pd((double*)(x+q+j)),c=_mm256_load_pd((double*)(x+2*q+j)),d=_mm256_load_pd((double*)(x+3*q+j));
__m256d t0=_mm256_add_pd(a,c),t1=_mm256_sub_pd(a,c),t2=_mm256_add_pd(b,d),t3=minus_i(_mm256_sub_pd(b,d));
_mm256_store_pd((double*)(x+j),_mm256_add_pd(t0,t2));_mm256_store_pd((double*)(x+q+j),cmul(_mm256_add_pd(t1,t3),rw));
_mm256_store_pd((double*)(x+2*q+j),cmul(_mm256_sub_pd(t0,t2),rw2));_mm256_store_pd((double*)(x+3*q+j),cmul(_mm256_sub_pd(t1,t3),rw3));}
}
template<unsigned M> static inline void inverse_block_tab(C*x){
constexpr unsigned q=M/4;C *rt=small_roots<M>(),*rt2=small_roots2<M>(),*rt3=small_roots3<M>();
for(unsigned j=0;j<q;j+=2){__m256d rw=conjv(_mm256_load_pd((double*)(rt+j))),rw2=conjv(_mm256_load_pd((double*)(rt2+j))),rw3=conjv(_mm256_load_pd((double*)(rt3+j)));
__m256d a=_mm256_load_pd((double*)(x+j)),b=cmul(_mm256_load_pd((double*)(x+q+j)),rw),c=cmul(_mm256_load_pd((double*)(x+2*q+j)),rw2),d=cmul(_mm256_load_pd((double*)(x+3*q+j)),rw3);
__m256d t0=_mm256_add_pd(a,c),t1=_mm256_sub_pd(a,c),t2=_mm256_add_pd(b,d),t3=plus_i(_mm256_sub_pd(b,d));
_mm256_store_pd((double*)(x+j),_mm256_add_pd(t0,t2));_mm256_store_pd((double*)(x+q+j),_mm256_add_pd(t1,t3));
_mm256_store_pd((double*)(x+2*q+j),_mm256_sub_pd(t0,t2));_mm256_store_pd((double*)(x+3*q+j),_mm256_sub_pd(t1,t3));}
}
#endif
static inline void forward4(C*x){C a=x[0],b=x[1],c=x[2],d=x[3],t0={a.r+c.r,a.i+c.i},t1={a.r-c.r,a.i-c.i},t2={b.r+d.r,b.i+d.i},v={b.r-d.r,b.i-d.i},t3={v.i,-v.r};x[0]={t0.r+t2.r,t0.i+t2.i};x[1]={t1.r+t3.r,t1.i+t3.i};x[2]={t0.r-t2.r,t0.i-t2.i};x[3]={t1.r-t3.r,t1.i-t3.i};}
static inline void inverse4(C*x){C a=x[0],b=x[1],c=x[2],d=x[3],t0={a.r+c.r,a.i+c.i},t1={a.r-c.r,a.i-c.i},t2={b.r+d.r,b.i+d.i},v={b.r-d.r,b.i-d.i},t3={-v.i,v.r};x[0]={t0.r+t2.r,t0.i+t2.i};x[1]={t1.r+t3.r,t1.i+t3.i};x[2]={t0.r-t2.r,t0.i-t2.i};x[3]={t1.r-t3.r,t1.i-t3.i};}
static inline void forward4x2(C*x){
__m256d r0=_mm256_load_pd((double*)x),r1=_mm256_load_pd((double*)(x+2)),r2=_mm256_load_pd((double*)(x+4)),r3=_mm256_load_pd((double*)(x+6));
__m256d a=_mm256_permute2f128_pd(r0,r2,0x20),b=_mm256_permute2f128_pd(r0,r2,0x31),c=_mm256_permute2f128_pd(r1,r3,0x20),d=_mm256_permute2f128_pd(r1,r3,0x31);
__m256d t0=_mm256_add_pd(a,c),t1=_mm256_sub_pd(a,c),t2=_mm256_add_pd(b,d),t3=minus_i(_mm256_sub_pd(b,d));
a=_mm256_add_pd(t0,t2);b=_mm256_add_pd(t1,t3);c=_mm256_sub_pd(t0,t2);d=_mm256_sub_pd(t1,t3);
_mm256_store_pd((double*)x,_mm256_permute2f128_pd(a,b,0x20));_mm256_store_pd((double*)(x+2),_mm256_permute2f128_pd(c,d,0x20));
_mm256_store_pd((double*)(x+4),_mm256_permute2f128_pd(a,b,0x31));_mm256_store_pd((double*)(x+6),_mm256_permute2f128_pd(c,d,0x31));
}
static inline void inverse4x2(C*x){
__m256d r0=_mm256_load_pd((double*)x),r1=_mm256_load_pd((double*)(x+2)),r2=_mm256_load_pd((double*)(x+4)),r3=_mm256_load_pd((double*)(x+6));
__m256d a=_mm256_permute2f128_pd(r0,r2,0x20),b=_mm256_permute2f128_pd(r0,r2,0x31),c=_mm256_permute2f128_pd(r1,r3,0x20),d=_mm256_permute2f128_pd(r1,r3,0x31);
__m256d t0=_mm256_add_pd(a,c),t1=_mm256_sub_pd(a,c),t2=_mm256_add_pd(b,d),t3=plus_i(_mm256_sub_pd(b,d));
a=_mm256_add_pd(t0,t2);b=_mm256_add_pd(t1,t3);c=_mm256_sub_pd(t0,t2);d=_mm256_sub_pd(t1,t3);
_mm256_store_pd((double*)x,_mm256_permute2f128_pd(a,b,0x20));_mm256_store_pd((double*)(x+2),_mm256_permute2f128_pd(c,d,0x20));
_mm256_store_pd((double*)(x+4),_mm256_permute2f128_pd(a,b,0x31));_mm256_store_pd((double*)(x+6),_mm256_permute2f128_pd(c,d,0x31));
}
template<unsigned M> struct RecFFT{
#ifdef ALL_ROOTS
static void forward(C*x){forward_block_alltab<M>(x);constexpr unsigned q=M/4;RecFFT<q>::forward(x);RecFFT<q>::forward(x+q);RecFFT<q>::forward(x+2*q);RecFFT<q>::forward(x+3*q);}
static void inverse(C*x){constexpr unsigned q=M/4;RecFFT<q>::inverse(x);RecFFT<q>::inverse(x+q);RecFFT<q>::inverse(x+2*q);RecFFT<q>::inverse(x+3*q);inverse_block_alltab<M>(x);}
#else
static void forward(C*x){forward_block<M>(x);constexpr unsigned q=M/4;RecFFT<q>::forward(x);RecFFT<q>::forward(x+q);RecFFT<q>::forward(x+2*q);RecFFT<q>::forward(x+3*q);}
static void inverse(C*x){constexpr unsigned q=M/4;RecFFT<q>::inverse(x);RecFFT<q>::inverse(x+q);RecFFT<q>::inverse(x+2*q);RecFFT<q>::inverse(x+3*q);inverse_block<M>(x);}
#endif
};
template<> struct RecFFT<256>{
#ifdef ALL_ROOTS
static void forward(C*x){forward_block_alltab<256>(x);for(unsigned p=0;p<256;p+=64)forward_block_alltab<64>(x+p);for(unsigned p=0;p<256;p+=16)forward_block_alltab<16>(x+p);for(unsigned p=0;p<256;p+=8)forward4x2(x+p);}
static void inverse(C*x){for(unsigned p=0;p<256;p+=8)inverse4x2(x+p);for(unsigned p=0;p<256;p+=16)inverse_block_alltab<16>(x+p);for(unsigned p=0;p<256;p+=64)inverse_block_alltab<64>(x+p);inverse_block_alltab<256>(x);}
#elif defined(SMALL_ROOTS)
static void forward(C*x){forward_block_tab<256>(x);for(unsigned p=0;p<256;p+=64)forward_block_tab<64>(x+p);for(unsigned p=0;p<256;p+=16)forward_block_tab<16>(x+p);for(unsigned p=0;p<256;p+=8)forward4x2(x+p);}
static void inverse(C*x){for(unsigned p=0;p<256;p+=8)inverse4x2(x+p);for(unsigned p=0;p<256;p+=16)inverse_block_tab<16>(x+p);for(unsigned p=0;p<256;p+=64)inverse_block_tab<64>(x+p);inverse_block_tab<256>(x);}
#else
static void forward(C*x){forward_block<256>(x);for(unsigned p=0;p<256;p+=64)forward_block<64>(x+p);for(unsigned p=0;p<256;p+=16)forward_block<16>(x+p);for(unsigned p=0;p<256;p+=8)forward4x2(x+p);}
static void inverse(C*x){for(unsigned p=0;p<256;p+=8)inverse4x2(x+p);for(unsigned p=0;p<256;p+=16)inverse_block<16>(x+p);for(unsigned p=0;p<256;p+=64)inverse_block<64>(x+p);inverse_block<256>(x);}
#endif
};
#if defined(SMALL_ROOTS) && !defined(ALL_ROOTS)
template<> struct RecFFT<1024>{
static void forward(C*x){forward_block_tab<1024>(x);for(unsigned k=0;k<4;++k)RecFFT<256>::forward(x+256*k);}
static void inverse(C*x){for(unsigned k=0;k<4;++k)RecFFT<256>::inverse(x+256*k);inverse_block_tab<1024>(x);}
};
template<> struct RecFFT<4096>{
static void forward(C*x){forward_block_tab<4096>(x);for(unsigned k=0;k<4;++k)RecFFT<1024>::forward(x+1024*k);}
static void inverse(C*x){for(unsigned k=0;k<4;++k)RecFFT<1024>::inverse(x+1024*k);inverse_block_tab<4096>(x);}
};
template<> struct RecFFT<16384>{
static void forward(C*x){forward_block_tab<16384>(x);for(unsigned k=0;k<4;++k)RecFFT<4096>::forward(x+4096*k);}
static void inverse(C*x){for(unsigned k=0;k<4;++k)RecFFT<4096>::inverse(x+4096*k);inverse_block_tab<16384>(x);}
};
#endif
#ifdef SPARSE_FORWARD_TOP
static void forward_sparse_top(C *x){
constexpr unsigned q=H/4,active=LIMBS/2-q;
constexpr double ang=-2.0*3.141592653589793238462643383279502884/H;
C step={cos(ang),sin(ang)},step2=mul(step,step);
__m256d rw=_mm256_setr_pd(1,0,step.r,step.i),advance=_mm256_setr_pd(step2.r,step2.i,step2.r,step2.i);
unsigned j=0;
for(;j<active;j+=2){
if(j&&(j&63u)==0){__m256d sq=_mm256_mul_pd(rw,rw),mag=_mm256_add_pd(_mm256_movedup_pd(sq),_mm256_permute_pd(sq,15));rw=_mm256_mul_pd(rw,_mm256_sub_pd(_mm256_set1_pd(1.5),_mm256_mul_pd(_mm256_set1_pd(.5),mag)));}
__m256d rw2=cmul(rw,rw),rw3=cmul(rw2,rw),av=_mm256_load_pd((double*)(x+j)),bv=_mm256_load_pd((double*)(x+q+j)),ib=minus_i(bv);
_mm256_store_pd((double*)(x+j),_mm256_add_pd(av,bv));
_mm256_store_pd((double*)(x+q+j),cmul(_mm256_add_pd(av,ib),rw));
_mm256_store_pd((double*)(x+2*q+j),cmul(_mm256_sub_pd(av,bv),rw2));
_mm256_store_pd((double*)(x+3*q+j),cmul(_mm256_sub_pd(av,ib),rw3));
rw=cmul(rw,advance);
}
for(;j<q;j+=2){
if(j&&(j&63u)==0){__m256d sq=_mm256_mul_pd(rw,rw),mag=_mm256_add_pd(_mm256_movedup_pd(sq),_mm256_permute_pd(sq,15));rw=_mm256_mul_pd(rw,_mm256_sub_pd(_mm256_set1_pd(1.5),_mm256_mul_pd(_mm256_set1_pd(.5),mag)));}
__m256d rw2=cmul(rw,rw),rw3=cmul(rw2,rw),av=_mm256_load_pd((double*)(x+j));
_mm256_store_pd((double*)(x+j),av);
_mm256_store_pd((double*)(x+q+j),cmul(av,rw));
_mm256_store_pd((double*)(x+2*q+j),cmul(av,rw2));
_mm256_store_pd((double*)(x+3*q+j),cmul(av,rw3));
rw=cmul(rw,advance);
}
for(unsigned k=0;k<4;++k)RecFFT<q>::forward(x+k*q);
}
#endif
#ifdef REALFFT
static inline void forward_twiddled4(C*x,unsigned j,unsigned q,__m256d rw,__m256d rw2,__m256d rw3){
__m256d a=_mm256_load_pd((double*)(x+j)),b=_mm256_load_pd((double*)(x+q+j));
__m256d c=_mm256_load_pd((double*)(x+2*q+j)),d=_mm256_load_pd((double*)(x+3*q+j));
__m256d t0=_mm256_add_pd(a,c),t1=_mm256_sub_pd(a,c),t2=_mm256_add_pd(b,d),t3=minus_i(_mm256_sub_pd(b,d));
_mm256_store_pd((double*)(x+j),_mm256_add_pd(t0,t2));
_mm256_store_pd((double*)(x+q+j),cmul(_mm256_add_pd(t1,t3),rw));
_mm256_store_pd((double*)(x+2*q+j),cmul(_mm256_sub_pd(t0,t2),rw2));
_mm256_store_pd((double*)(x+3*q+j),cmul(_mm256_sub_pd(t1,t3),rw3));
}
template<unsigned M> static inline void forward_block_pair(C*x,C*y){
constexpr unsigned q=M/4;constexpr double ang=-2.0*3.141592653589793238462643383279502884/M;
C step={cos(ang),sin(ang)},w={1,0};
for(unsigned j=0;j<q;j+=2){
if(j&&(j&63u)==0)w=norm(w);C w0=w,w1=mul(w,step);w=mul(w1,step);
__m256d rw=_mm256_setr_pd(w0.r,w0.i,w1.r,w1.i),rw2=cmul(rw,rw),rw3=cmul(rw2,rw);
forward_twiddled4(x,j,q,rw,rw2,rw3);forward_twiddled4(y,j,q,rw,rw2,rw3);
}
}
template<unsigned M> struct RecFFTPair{
static void forward(C*x,C*y){forward_block_pair<M>(x,y);constexpr unsigned q=M/4;RecFFTPair<q>::forward(x,y);RecFFTPair<q>::forward(x+q,y+q);RecFFTPair<q>::forward(x+2*q,y+2*q);RecFFTPair<q>::forward(x+3*q,y+3*q);}
};
template<> struct RecFFTPair<256>{
static void forward(C*x,C*y){
forward_block_pair<256>(x,y);for(unsigned p=0;p<256;p+=64)forward_block_pair<64>(x+p,y+p);for(unsigned p=0;p<256;p+=16)forward_block_pair<16>(x+p,y+p);
for(unsigned p=0;p<256;p+=8){forward4x2(x+p);forward4x2(y+p);}
}
};
static void forward_pair_top(C*x,C*y){
forward_block_pair<H>(x,y);constexpr unsigned q=H/4,r=q/4;
for(unsigned k=0;k<4;++k){C *a=x+k*q,*b=y+k*q;forward_block_pair<q>(a,b);for(unsigned j=0;j<4;++j)RecFFT<r>::forward(a+j*r);for(unsigned j=0;j<4;++j)RecFFT<r>::forward(b+j*r);}
}
#endif
#ifdef RADIX8
static inline void dft8_forward(__m256d&x0,__m256d&x1,__m256d&x2,__m256d&x3,__m256d&x4,__m256d&x5,__m256d&x6,__m256d&x7){
const __m256d cw1=_mm256_setr_pd(.7071067811865475244,-.7071067811865475244,.7071067811865475244,-.7071067811865475244);
const __m256d cw3=_mm256_setr_pd(-.7071067811865475244,-.7071067811865475244,-.7071067811865475244,-.7071067811865475244);
__m256d a0=_mm256_add_pd(x0,x4),b0=_mm256_sub_pd(x0,x4),a1=_mm256_add_pd(x1,x5),b1=cmul(_mm256_sub_pd(x1,x5),cw1);
__m256d a2=_mm256_add_pd(x2,x6),b2=minus_i(_mm256_sub_pd(x2,x6)),a3=_mm256_add_pd(x3,x7),b3=cmul(_mm256_sub_pd(x3,x7),cw3);
__m256d e0=_mm256_add_pd(a0,a2),e1=_mm256_sub_pd(a0,a2),e2=_mm256_add_pd(a1,a3),e3=minus_i(_mm256_sub_pd(a1,a3));
__m256d o0=_mm256_add_pd(b0,b2),o1=_mm256_sub_pd(b0,b2),o2=_mm256_add_pd(b1,b3),o3=minus_i(_mm256_sub_pd(b1,b3));
x0=_mm256_add_pd(e0,e2);x2=_mm256_add_pd(e1,e3);x4=_mm256_sub_pd(e0,e2);x6=_mm256_sub_pd(e1,e3);
x1=_mm256_add_pd(o0,o2);x3=_mm256_add_pd(o1,o3);x5=_mm256_sub_pd(o0,o2);x7=_mm256_sub_pd(o1,o3);
}
static inline void dft8_inverse(__m256d&x0,__m256d&x1,__m256d&x2,__m256d&x3,__m256d&x4,__m256d&x5,__m256d&x6,__m256d&x7){
const __m256d cw1=_mm256_setr_pd(.7071067811865475244,.7071067811865475244,.7071067811865475244,.7071067811865475244);
const __m256d cw3=_mm256_setr_pd(-.7071067811865475244,.7071067811865475244,-.7071067811865475244,.7071067811865475244);
__m256d a0=_mm256_add_pd(x0,x4),b0=_mm256_sub_pd(x0,x4),a1=_mm256_add_pd(x1,x5),b1=cmul(_mm256_sub_pd(x1,x5),cw1);
__m256d a2=_mm256_add_pd(x2,x6),b2=plus_i(_mm256_sub_pd(x2,x6)),a3=_mm256_add_pd(x3,x7),b3=cmul(_mm256_sub_pd(x3,x7),cw3);
__m256d e0=_mm256_add_pd(a0,a2),e1=_mm256_sub_pd(a0,a2),e2=_mm256_add_pd(a1,a3),e3=plus_i(_mm256_sub_pd(a1,a3));
__m256d o0=_mm256_add_pd(b0,b2),o1=_mm256_sub_pd(b0,b2),o2=_mm256_add_pd(b1,b3),o3=plus_i(_mm256_sub_pd(b1,b3));
x0=_mm256_add_pd(e0,e2);x2=_mm256_add_pd(e1,e3);x4=_mm256_sub_pd(e0,e2);x6=_mm256_sub_pd(e1,e3);
x1=_mm256_add_pd(o0,o2);x3=_mm256_add_pd(o1,o3);x5=_mm256_sub_pd(o0,o2);x7=_mm256_sub_pd(o1,o3);
}
template<unsigned M> static inline void forward_block8(C*x){
constexpr unsigned q=M/8;constexpr double ang=-2.0*3.141592653589793238462643383279502884/M;
C step={cos(ang),sin(ang)},w={1,0};
for(unsigned j=0;j<q;j+=2){
if(j&&(j&63u)==0)w=norm(w);C w0=w,w1=mul(w,step);w=mul(w1,step);
__m256d rw=_mm256_setr_pd(w0.r,w0.i,w1.r,w1.i),wp[8];wp[1]=rw;
for(unsigned k=2;k<8;++k)wp[k]=cmul(wp[k-1],rw);
__m256d x0=_mm256_load_pd((double*)(x+j)),x1=_mm256_load_pd((double*)(x+q+j)),x2=_mm256_load_pd((double*)(x+2*q+j)),x3=_mm256_load_pd((double*)(x+3*q+j));
__m256d x4=_mm256_load_pd((double*)(x+4*q+j)),x5=_mm256_load_pd((double*)(x+5*q+j)),x6=_mm256_load_pd((double*)(x+6*q+j)),x7=_mm256_load_pd((double*)(x+7*q+j));
dft8_forward(x0,x1,x2,x3,x4,x5,x6,x7);
_mm256_store_pd((double*)(x+j),x0);_mm256_store_pd((double*)(x+q+j),cmul(x1,wp[1]));_mm256_store_pd((double*)(x+2*q+j),cmul(x2,wp[2]));_mm256_store_pd((double*)(x+3*q+j),cmul(x3,wp[3]));
_mm256_store_pd((double*)(x+4*q+j),cmul(x4,wp[4]));_mm256_store_pd((double*)(x+5*q+j),cmul(x5,wp[5]));_mm256_store_pd((double*)(x+6*q+j),cmul(x6,wp[6]));_mm256_store_pd((double*)(x+7*q+j),cmul(x7,wp[7]));
}
}
template<unsigned M> static inline void inverse_block8(C*x){
constexpr unsigned q=M/8;constexpr double ang=-2.0*3.141592653589793238462643383279502884/M;
C step={cos(ang),sin(ang)},w={1,0};
for(unsigned j=0;j<q;j+=2){
if(j&&(j&63u)==0)w=norm(w);C w0=w,w1=mul(w,step);w=mul(w1,step);
__m256d rw=conjv(_mm256_setr_pd(w0.r,w0.i,w1.r,w1.i)),wp[8];wp[1]=rw;
for(unsigned k=2;k<8;++k)wp[k]=cmul(wp[k-1],rw);
__m256d x0=_mm256_load_pd((double*)(x+j)),x1=cmul(_mm256_load_pd((double*)(x+q+j)),wp[1]),x2=cmul(_mm256_load_pd((double*)(x+2*q+j)),wp[2]),x3=cmul(_mm256_load_pd((double*)(x+3*q+j)),wp[3]);
__m256d x4=cmul(_mm256_load_pd((double*)(x+4*q+j)),wp[4]),x5=cmul(_mm256_load_pd((double*)(x+5*q+j)),wp[5]),x6=cmul(_mm256_load_pd((double*)(x+6*q+j)),wp[6]),x7=cmul(_mm256_load_pd((double*)(x+7*q+j)),wp[7]);
dft8_inverse(x0,x1,x2,x3,x4,x5,x6,x7);
_mm256_store_pd((double*)(x+j),x0);_mm256_store_pd((double*)(x+q+j),x1);_mm256_store_pd((double*)(x+2*q+j),x2);_mm256_store_pd((double*)(x+3*q+j),x3);
_mm256_store_pd((double*)(x+4*q+j),x4);_mm256_store_pd((double*)(x+5*q+j),x5);_mm256_store_pd((double*)(x+6*q+j),x6);_mm256_store_pd((double*)(x+7*q+j),x7);
}
}
static inline C add(C a,C b){return {a.r+b.r,a.i+b.i};}
static inline C sub(C a,C b){return {a.r-b.r,a.i-b.i};}
static inline C mi(C a){return {a.i,-a.r};}
static inline C pi(C a){return {-a.i,a.r};}
static inline void forward8(C*x){
constexpr double s=.7071067811865475244;
C a0=add(x[0],x[4]),b0=sub(x[0],x[4]),a1=add(x[1],x[5]),v1=sub(x[1],x[5]),b1={(v1.r+v1.i)*s,(v1.i-v1.r)*s};
C a2=add(x[2],x[6]),b2=mi(sub(x[2],x[6])),a3=add(x[3],x[7]),v3=sub(x[3],x[7]),b3={(v3.i-v3.r)*s,-(v3.r+v3.i)*s};
C e0=add(a0,a2),e1=sub(a0,a2),e2=add(a1,a3),e3=mi(sub(a1,a3));
C o0=add(b0,b2),o1=sub(b0,b2),o2=add(b1,b3),o3=mi(sub(b1,b3));
x[0]=add(e0,e2);x[2]=add(e1,e3);x[4]=sub(e0,e2);x[6]=sub(e1,e3);
x[1]=add(o0,o2);x[3]=add(o1,o3);x[5]=sub(o0,o2);x[7]=sub(o1,o3);
}
static inline void inverse8(C*x){
constexpr double s=.7071067811865475244;
C a0=add(x[0],x[4]),b0=sub(x[0],x[4]),a1=add(x[1],x[5]),v1=sub(x[1],x[5]),b1={(v1.r-v1.i)*s,(v1.r+v1.i)*s};
C a2=add(x[2],x[6]),b2=pi(sub(x[2],x[6])),a3=add(x[3],x[7]),v3=sub(x[3],x[7]),b3={-(v3.r+v3.i)*s,(v3.r-v3.i)*s};
C e0=add(a0,a2),e1=sub(a0,a2),e2=add(a1,a3),e3=pi(sub(a1,a3));
C o0=add(b0,b2),o1=sub(b0,b2),o2=add(b1,b3),o3=pi(sub(b1,b3));
x[0]=add(e0,e2);x[2]=add(e1,e3);x[4]=sub(e0,e2);x[6]=sub(e1,e3);
x[1]=add(o0,o2);x[3]=add(o1,o3);x[5]=sub(o0,o2);x[7]=sub(o1,o3);
}
template<unsigned M> struct RecFFT8{
static void forward(C*x){forward_block8<M>(x);constexpr unsigned q=M/8;for(unsigned k=0;k<8;++k)RecFFT8<q>::forward(x+k*q);}
static void inverse(C*x){constexpr unsigned q=M/8;for(unsigned k=0;k<8;++k)RecFFT8<q>::inverse(x+k*q);inverse_block8<M>(x);}
};
template<> struct RecFFT8<512>{
static void forward(C*x){forward_block8<512>(x);for(unsigned p=0;p<512;p+=64)forward_block8<64>(x+p);for(unsigned p=0;p<512;p+=8)forward8(x+p);}
static void inverse(C*x){for(unsigned p=0;p<512;p+=8)inverse8(x+p);for(unsigned p=0;p<512;p+=64)inverse_block8<64>(x+p);inverse_block8<512>(x);}
};
static void forward_radix8(){radix2_forward();RecFFT8<H>::forward(z);RecFFT8<H>::forward(z+H);}
static void inverse_radix8(){RecFFT8<H>::inverse(z);RecFFT8<H>::inverse(z+H);radix2_inverse();}
#endif
static void forward(){radix2_forward();RecFFT<H>::forward(z);RecFFT<H>::forward(z+H);}
static void inverse(){RecFFT<H>::inverse(z);RecFFT<H>::inverse(z+H);radix2_inverse();}
static inline unsigned parse4(const char *p){return (unsigned)(p[0]-'0')*1000u+(unsigned)(p[1]-'0')*100u+(unsigned)(p[2]-'0')*10u+(unsigned)(p[3]-'0');}
#ifdef BASE5
static inline unsigned parse5(const char *p){return (unsigned)(p[0]-'0')*10000u+parse4(p+1);}
#endif
static inline void spectral_pair(unsigned p,unsigned q){
if(p==q){z[p]={z[p].r*z[p].i,0};return;}
C a=z[p],b={z[q].r,-z[q].i};C av={(a.r+b.r)*.5,(a.i+b.i)*.5};C bv={(a.i-b.i)*.5,(b.r-a.r)*.5};C v=mul(av,bv);z[p]=v;z[q]={v.r,-v.i};
}
static void spectral_range(unsigned p,unsigned q){
const __m256d half=_mm256_set1_pd(.5);
while(p+1<q-1){
__m256d a=_mm256_load_pd((double*)(z+p));
__m256d b=_mm256_permute4x64_pd(_mm256_loadu_pd((double*)(z+q-1)),0x4e);b=conjv(b);
__m256d av=_mm256_mul_pd(_mm256_add_pd(a,b),half),bv=_mm256_mul_pd(minus_i(_mm256_sub_pd(a,b)),half);
__m256d v=cmul(av,bv);_mm256_store_pd((double*)(z+p),v);
v=_mm256_permute4x64_pd(conjv(v),0x4e);_mm256_storeu_pd((double*)(z+q-1),v);
p+=2;q-=2;
}
while(p<=q){spectral_pair(p,q);++p;--q;}
}
#ifdef RADIX8
static inline unsigned rev_chunks18(unsigned x){
unsigned lo=((x&7u)<<6)|(x&56u)|((x>>6)&7u),y=x>>9;
unsigned hi=((y&7u)<<6)|(y&56u)|((y>>6)&7u);
return (lo<<9)|hi;
}
static inline unsigned radix8_conjugate_index(unsigned p){
unsigned f=(p>=H) | (rev_chunks18(p&(H-1))<<1);
f=(-f)&(N-1);
return ((f&1)*H)|rev_chunks18(f>>1);
}
static void spectral_radix8(){
for(unsigned p=0;p<N;++p){unsigned q=radix8_conjugate_index(p);if(p<=q)spectral_pair(p,q);}
}
#endif
#ifdef REALFFT
static inline unsigned rev4_18(unsigned x){
return ((x&3u)<<16)|((x&12u)<<12)|((x&48u)<<8)|((x&192u)<<4)|(x&768u)|
((x&3072u)>>4)|((x&12288u)>>8)|((x&49152u)>>12)|((x&196608u)>>16);
}
static inline C real_spectrum(C u,C v,C w){
v.i=-v.i;
C s={(u.r+v.r)*.5,(u.i+v.i)*.5},d={(u.r-v.r)*.5,(u.i-v.i)*.5},t=mul(w,d);
return {s.r+t.i,s.i-t.r};
}
static inline C pack_spectrum(C a,C bh,C w,double scale){
bh.i=-bh.i;
C s={(a.r+bh.r)*.5,(a.i+bh.i)*.5},d={(a.r-bh.r)*.5,(a.i-bh.i)*.5};
C o=mul(d,{w.r,-w.i});
return {(s.r-o.i)*scale,(s.i+o.r)*scale};
}
#ifdef VECTOR_REAL_PRODUCT
static inline __m256d real_spectrum_v(__m256d u,__m256d v,__m256d w){
v=conjv(v);__m256d s=_mm256_mul_pd(_mm256_add_pd(u,v),_mm256_set1_pd(.5));
__m256d d=_mm256_mul_pd(_mm256_sub_pd(u,v),_mm256_set1_pd(.5));return _mm256_add_pd(s,minus_i(cmul(w,d)));
}
static inline __m256d pack_spectrum_v(__m256d a,__m256d bh,__m256d w){
bh=conjv(bh);__m256d hs=_mm256_set1_pd(.5/H);
__m256d s=_mm256_mul_pd(_mm256_add_pd(a,bh),hs),d=_mm256_mul_pd(_mm256_sub_pd(a,bh),hs);
return _mm256_add_pd(s,plus_i(cmul(d,conjv(w))));
}
#endif
#ifdef PARALLEL_FORWARD
static void *real_forward_worker(void *p){RecFFT<H>::forward((C*)p);return 0;}
#endif
#ifdef PARALLEL_INVERSE
static void *real_inverse_worker(void *p){
C *x=(C*)p;constexpr unsigned q=H/4;
RecFFT<q>::inverse(x);RecFFT<q>::inverse(x+q);return 0;
}
static void real_inverse_parallel(C *x){
constexpr unsigned q=H/4;pthread_t worker;
int threaded=pthread_create(&worker,0,real_inverse_worker,x+2*q)==0;
RecFFT<q>::inverse(x);RecFFT<q>::inverse(x+q);
if(threaded)pthread_join(worker,0);
else {RecFFT<q>::inverse(x+2*q);RecFFT<q>::inverse(x+3*q);}
inverse_block<H>(x);
}
#endif
#ifdef LOCAL_TEST
static void solve(DuckInfo *di){
#else
static __attribute__((noreturn)) void solve(DuckInfo *di){
#endif
const char *in=di->stdin_ptr;C *a=z,*b=z+H;
#ifdef ALL_ROOTS
build_all_roots();
#elif defined(SMALL_ROOTS)
build_small_roots();
#endif
for(unsigned j=0;j<LIMBS/2;++j){
#if defined(BASE5) && !defined(BALANCED5)
a[j].r=parse5(in+DIGITS-5-10*j);a[j].i=parse5(in+DIGITS-10-10*j);
b[j].r=parse5(in+2*DIGITS+1-5-10*j);b[j].i=parse5(in+2*DIGITS+1-10-10*j);
#elif !defined(BASE5)
#ifdef SIMD_PARSE
const char *pa=in+DIGITS-8-8*j,*pb=in+2*DIGITS+1-8-8*j;
__m128i raw=_mm_unpacklo_epi64(_mm_loadl_epi64((const __m128i*)pa),_mm_loadl_epi64((const __m128i*)pb));
raw=_mm_sub_epi8(raw,_mm_set1_epi8('0'));
__m128i p2=_mm_maddubs_epi16(raw,_mm_set1_epi16(0x010a));
__m128i p4=_mm_madd_epi16(p2,_mm_set1_epi32(0x00010064));
__m256d v=_mm256_cvtepi32_pd(_mm_shuffle_epi32(p4,0xb1));
_mm_store_pd((double*)(a+j),_mm256_castpd256_pd128(v));
_mm_store_pd((double*)(b+j),_mm256_extractf128_pd(v,1));
#else
a[j].r=parse4(in+DIGITS-4-8*j);a[j].i=parse4(in+DIGITS-8-8*j);
b[j].r=parse4(in+2*DIGITS+1-4-8*j);b[j].i=parse4(in+2*DIGITS+1-8-8*j);
#endif
#endif
}
#ifdef BALANCED5
int ca=0,cb=0;for(unsigned j=0;j<LIMBS;++j){int x=(int)parse5(in+DIGITS-5-5*j)+ca,y=(int)parse5(in+2*DIGITS+1-5-5*j)+cb;ca=x>=50000;cb=y>=50000;x-=ca*BASE;y-=cb*BASE;if(j&1){a[j/2].i=x;b[j/2].i=y;}else{a[j/2].r=x;b[j/2].r=y;}}
a[LIMBS/2].r=ca;b[LIMBS/2].r=cb;
#endif
#ifdef PROFILE_PARSE
di->stdout_size=0;duck_exit();
#endif
#ifdef PARALLEL_FORWARD
pthread_t worker;int threaded=pthread_create(&worker,0,real_forward_worker,b)==0;
RecFFT<H>::forward(a);if(threaded)pthread_join(worker,0);else RecFFT<H>::forward(b);
#else
#ifdef SPARSE_FORWARD_TOP
forward_sparse_top(a);forward_sparse_top(b);
#else
RecFFT<H>::forward(a);RecFFT<H>::forward(b);
#endif
#endif
#ifdef PROFILE_FORWARD
di->stdout_size=0;duck_exit();
#endif
C a0=a[0],b0=b[0];double c0=(a0.r+a0.i)*(b0.r+b0.i),ch=(a0.r-a0.i)*(b0.r-b0.i);
const double scale=1.0/H;a[0]={(c0+ch)*.5*scale,(c0-ch)*.5*scale};
const double ang=-2.0*3.141592653589793238462643383279502884/N;C jump[9];
for(unsigned t=0;t<9;++t){
int delta=1<<(16-2*t);for(unsigned j=0;j<t;++j)delta-=3<<(16-2*j);
jump[t]={cos(ang*delta),sin(ang*delta)};
}
for(unsigned h=1;h<H;h<<=2){
unsigned p=h,q=4*h-1,k=rev4_18(p);C w;
#ifdef PRODUCT_ROOT_TABLE
if(h>=64)w=product_roots[k];else w={cos(ang*k),sin(ang*k)};
#else
w={cos(ang*k),sin(ang*k)};
#endif
#ifdef VECTOR_REAL_PRODUCT
if(h==1){
for(;p<=q;++p,--q){C wh={-w.r,w.i};C ak=real_spectrum(a[p],a[q],w),ah=real_spectrum(a[q],a[p],wh),bk=real_spectrum(b[p],b[q],w),bh=real_spectrum(b[q],b[p],wh);C ck=mul(ak,bk),c_h=mul(ah,bh);a[p]=pack_spectrum(ck,c_h,w,scale);if(p!=q)a[q]=pack_spectrum(c_h,ck,wh,scale);unsigned x=p,t=0;while((x&3u)==3u){x>>=2;++t;}w=mul(w,jump[t]);}
continue;
}
for(;p<q;p+=2,q-=2){
if((p&63u)==0){
#ifdef PRODUCT_NORMALIZE_ONLY
w=norm(w);
#elif defined(PRODUCT_ROOT_TABLE)
w=product_roots[rev4_18(p)];
#else
unsigned rk=rev4_18(p);w={cos(ang*rk),sin(ang*rk)};
#endif
}
unsigned x,t;C w1=mul(w,jump[0]);
__m256d vw=_mm256_setr_pd(w.r,w.i,w1.r,w1.i),vwh=_mm256_xor_pd(vw,_mm256_castsi256_pd(_mm256_setr_epi64x(0x8000000000000000ULL,0,0x8000000000000000ULL,0)));
__m256d ap=_mm256_load_pd((double*)(a+p)),aq=_mm256_permute4x64_pd(_mm256_load_pd((double*)(a+q-1)),0x4e);
__m256d bp=_mm256_load_pd((double*)(b+p)),bq=_mm256_permute4x64_pd(_mm256_load_pd((double*)(b+q-1)),0x4e);
#ifdef ALGEBRA_PRODUCT
__m256d one=_mm256_setr_pd(1,0,1,0),vs=_mm256_set1_pd(1.0/H),w2=cmul(vw,vw),f=_mm256_mul_pd(_mm256_add_pd(one,w2),_mm256_set1_pd(.25));
__m256d da=_mm256_sub_pd(ap,conjv(aq)),db=_mm256_sub_pd(bp,conjv(bq));
__m256d corr=cmul(f,cmul(da,db));
__m256d op=_mm256_mul_pd(_mm256_sub_pd(cmul(ap,bp),corr),vs);
__m256d oq=_mm256_mul_pd(_mm256_sub_pd(cmul(aq,bq),conjv(corr)),vs);
#else
__m256d ak=real_spectrum_v(ap,aq,vw),ah=real_spectrum_v(aq,ap,vwh),bk=real_spectrum_v(bp,bq,vw),bh=real_spectrum_v(bq,bp,vwh);
__m256d ck=cmul(ak,bk),ch=cmul(ah,bh),op=pack_spectrum_v(ck,ch,vw),oq=pack_spectrum_v(ch,ck,vwh);
#endif
_mm256_store_pd((double*)(a+p),op);_mm256_store_pd((double*)(a+q-1),_mm256_permute4x64_pd(oq,0x4e));
x=p+1;t=0;while((x&3u)==3u){x>>=2;++t;}w=mul(w1,jump[t]);
}
#else
for(;p<=q;++p,--q){
if((p&63u)==0){
#ifdef PRODUCT_NORMALIZE_ONLY
w=norm(w);
#elif defined(PRODUCT_ROOT_TABLE)
w=product_roots[rev4_18(p)];
#else
unsigned rk=rev4_18(p);w={cos(ang*rk),sin(ang*rk)};
#endif
}C wh={-w.r,w.i};
C ak=real_spectrum(a[p],a[q],w),ah=real_spectrum(a[q],a[p],wh);
C bk=real_spectrum(b[p],b[q],w),bh=real_spectrum(b[q],b[p],wh);
C ck=mul(ak,bk),c_h=mul(ah,bh);
a[p]=pack_spectrum(ck,c_h,w,scale);if(p!=q)a[q]=pack_spectrum(c_h,ck,wh,scale);
unsigned x=p,t=0;while((x&3u)==3u){x>>=2;++t;}w=mul(w,jump[t]);
}
#endif
}
#ifdef PROFILE_PRODUCT
di->stdout_size=0;duck_exit();
#endif
#ifdef PARALLEL_INVERSE
real_inverse_parallel(a);
#else
for(unsigned k=0;k<4;++k)RecFFT<H/4>::inverse(a+k*(H/4));
inverse_block_round<H>(a);
#endif
#ifdef PROFILE_INVERSE
di->stdout_size=0;duck_exit();
#endif
for(unsigned v=0;v<10000;++v)digit4[v]=(unsigned)('0'+v/1000)|((unsigned)('0'+v/100%10)<<8)|((unsigned)('0'+v/10%10)<<16)|((unsigned)('0'+v%10)<<24);
char *out=di->stdout_ptr;double *coef=(double*)a;
#ifdef BALANCED5
long long carry=0;for(unsigned j=0;j<2*LIMBS;++j){long long v=(long long)(coef[j]+(coef[j]>=0?.5:-.5))+carry,q=v/BASE,r=v-q*BASE;if(r<0){r+=BASE;--q;}carry=q;unsigned pos=2*DIGITS-5-5*j;out[pos]=(char)('0'+r/10000);*(unsigned*)(out+pos+1)=digit4[r%10000];}
unsigned skip=0;while(skip<4&&out[skip]=='0')++skip;
#elif defined(BASE5)
unsigned long long carry=0;
for(unsigned j=0;j<2*LIMBS;++j){unsigned long long v=(unsigned long long)(coef[j]+.5)+carry;carry=v/BASE;unsigned r=v-carry*BASE,pos=2*DIGITS-5-5*j;out[pos]=(char)('0'+r/10000);*(unsigned*)(out+pos+1)=digit4[r%10000];}
unsigned skip=0;while(skip<4&&out[skip]=='0')++skip;
#else
unsigned long long *icoef=(unsigned long long*)a;
unsigned long long carry=0;
for(unsigned j=0;j<2*LIMBS;j+=2){
unsigned long long v0=icoef[j];
unsigned long long v1=icoef[j+1];
unsigned long long pair=v0+carry+v1*BASE;
carry=pair/100000000ULL;pair-=carry*100000000ULL;
unsigned hi=pair/BASE,lo=pair-hi*BASE;
*(unsigned long long*)(out+4*(2*LIMBS-2-j))=(unsigned long long)digit4[hi]|((unsigned long long)digit4[lo]<<32);
}
unsigned skip=0;while(skip<3&&out[skip]=='0')++skip;
#endif
if(skip)for(unsigned j=0;j<2*DIGITS-skip;++j)out[j]=out[j+skip];
out[2*DIGITS-skip]='\n';di->stdout_size=2*DIGITS+1-skip;
#ifdef LOCAL_TEST
return;
#else
duck_exit();
#endif
}
#else
#ifdef LOCAL_TEST
static void solve(DuckInfo *di){
#else
static __attribute__((noreturn)) void solve(DuckInfo *di){
#endif
const char *in=di->stdin_ptr;
#ifdef LOCAL_IMPULSE
z[1].r=1;
#elif defined(LOCAL_PACKIMPULSE)
z[1]={1,1};
#elif defined(LOCAL_ONES)
for(unsigned i=0;i<LIMBS;++i)z[i]={1,1};
#elif defined(LOCAL_REALONES)
for(unsigned i=0;i<LIMBS;++i)z[i].r=1;
#else
for(unsigned i=0;i<LIMBS;++i){
z[i].r=parse4(in+DIGITS-4-4*i);
z[i].i=parse4(in+2*DIGITS+1-4-4*i);
}
#endif
#ifdef PROFILE_PARSE
di->stdout_size=0;duck_exit();
#endif
#ifdef RADIX8
forward_radix8();
#else
forward();
#endif
#ifdef PROFILE_FORWARD
di->stdout_size=0;duck_exit();
#endif
#ifdef INSPECT_FORWARD
for(unsigned i=0;i<32;++i)fprintf(stderr,"%u %.9f %.9f\n",i,z[i].r,z[i].i);
for(unsigned i=1;i<32;++i){unsigned p=i*(N/32);fprintf(stderr,"P %u %.9f %.9f\n",p,z[p].r,z[p].i);}
return;
#endif
#ifdef LOCAL_REALONES
for(unsigned p=0;p<N;++p)z[p]=mul(z[p],z[p]);
#else
#ifndef ROUNDTRIP
#ifdef RADIX8
spectral_radix8();
#else
spectral_pair(0,0);
for(unsigned h=1;h<H;h<<=2)spectral_range(h,4*h-1);
spectral_range(H,N-1);
#endif
#endif
#endif
#ifdef PROFILE_PRODUCT
di->stdout_size=0;duck_exit();
#endif
#ifdef RADIX8
inverse_radix8();
#else
inverse();
#endif
#ifdef PROFILE_INVERSE
di->stdout_size=0;duck_exit();
#endif
#ifdef LOCAL_DEBUG
for(unsigned i=0;i<16;++i)fprintf(stderr,"%u %.6f %.6f\n",i,z[i].r,z[i].i);
#endif
for(unsigned v=0;v<10000;++v)digit4[v]=(unsigned)('0'+v/1000)|((unsigned)('0'+v/100%10)<<8)|((unsigned)('0'+v/10%10)<<16)|((unsigned)('0'+v%10)<<24);
char *out=di->stdout_ptr;unsigned long long carry=0;
for(unsigned i=0;i<2*LIMBS;++i){unsigned long long v=(unsigned long long)(z[i].r+.5)+carry;carry=v/BASE;((unsigned*)out)[2*LIMBS-1-i]=digit4[v-carry*BASE];}
unsigned skip=0;while(skip<3&&out[skip]=='0')++skip;
if(skip)for(unsigned i=0;i<2*DIGITS-skip;++i)out[i]=out[i+skip];
out[2*DIGITS-skip]='\n';di->stdout_size=2*DIGITS+1-skip;
#ifdef LOCAL_TEST
return;
#else
duck_exit();
#endif
}
#endif
#ifndef LOCAL_TEST
extern "C" __attribute__((noreturn)) void __libc_start_main(void*,long argc,char **argv){solve(duck_info(argc,argv));}
int main(){}
#else
static char local_input[2000016],local_output[2000016];
int main(){
DuckInfo di={};di.stdin_ptr=local_input;di.stdout_ptr=local_output;
fread(local_input,1,sizeof(local_input),stdin);
solve(&di);
fwrite(local_output,1,di.stdout_size,stdout);
}
#endif
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 6.69 ms | 10 MB + 208 KB | Accepted | Score: 100 | 显示更多 |