// DUCK_PUBLIC_SURROGATE_BENCHMARK
// Independent public arrays only; no production generator, input or seeds.
#define BENCH_CASE 0
#define BENCH_H 8192
#include <bits/stdc++.h>
#include <x86intrin.h>
// FFT-only extraction from batch convpluse8 reference / Duck #141068 HintFFT.
// Original TSKY (WithSky) HintFFT MIT license is retained below.
// Uses the certified direct-root constructor from the batch reference.
#include <bits/stdc++.h>
#include <immintrin.h>
/*
MIT License
Copyright (c) 2023 With-Sky
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
constexpr auto statement = "Do NOT modify this statement! "
"https://github.com/With-Sky/HintFFT "
"TSKY (WithSky)";
#include <complex>
#include <iostream>
#include <type_traits>
#include <cstdint>
#include <climits>
#include <cstring>
#include <immintrin.h>
#define NI64 C64X4(hint::NoInitTag{})
#pragma GCC target("avx2,fma,tune=skylake")
#pragma GCC optimize("O3,tree-vectorize")
#pragma GCC optimize("align-functions=32")
#pragma GCC optimize("ira-loop-pressure")
#pragma GCC optimize("sched-pressure")
namespace hint
{
template <typename T, size_t ALIGN = 64>
class AlignMem{
public:
using Ptr = T *;
using ConstPtr = const T *;
~AlignMem(){
if (ptr){
_mm_free(ptr);}};
AlignMem() : ptr(nullptr), len(0) {}
AlignMem(size_t n) : ptr(reinterpret_cast<Ptr>(_mm_malloc(n * sizeof(T), ALIGN))), len(n) {}
AlignMem(const AlignMem &) = delete;
AlignMem &operator=(const AlignMem &) = delete;
T &operator[](size_t i){
return ptr[i];}
const T &operator[](size_t i) const{
return ptr[i];}
Ptr begin(){
return ptr;}
Ptr end(){
return ptr + len;}
ConstPtr begin() const{
return ptr;}
ConstPtr end() const{
return ptr + len;}
size_t size() const{
return len;}
private:
T *ptr;
size_t len;};
template <typename YMM>
inline void transpose64_2X4(YMM &row0, YMM &row1){
auto t0 = _mm256_unpacklo_pd(__m256d(row0), __m256d(row1));
auto t1 = _mm256_unpackhi_pd(__m256d(row0), __m256d(row1));
row0 = YMM(_mm256_permute2f128_pd(t0, t1, 0x20));
row1 = YMM(_mm256_permute2f128_pd(t0, t1, 0x31));}
template <typename YMM>
inline void transpose64_4X2(YMM &row0, YMM &row1){
auto t0 = _mm256_permute2f128_pd(__m256d(row0), __m256d(row1), 0x20);
auto t1 = _mm256_permute2f128_pd(__m256d(row0), __m256d(row1), 0x31);
row0 = YMM(_mm256_unpacklo_pd(t0, t1));
row1 = YMM(_mm256_unpackhi_pd(t0, t1));}
template <typename YMM>
inline void transpose64_4X4(YMM &row0, YMM &row1, YMM &row2, YMM &row3){
auto t0 = _mm256_unpacklo_pd(__m256d(row0), __m256d(row1));
auto t1 = _mm256_unpackhi_pd(__m256d(row0), __m256d(row1));
auto t2 = _mm256_unpacklo_pd(__m256d(row2), __m256d(row3));
auto t3 = _mm256_unpackhi_pd(__m256d(row2), __m256d(row3));
row0 = YMM(_mm256_permute2f128_pd(t0, t2, 0x20));
row1 = YMM(_mm256_permute2f128_pd(t1, t3, 0x20));
row2 = YMM(_mm256_permute2f128_pd(t0, t2, 0x31));
row3 = YMM(_mm256_permute2f128_pd(t1, t3, 0x31));}
struct NoInitTag {};
class Float64X4{
public:
using F64 = double;
using F64X4 = Float64X4;
Float64X4() : data(_mm256_setzero_pd()) {}
Float64X4(NoInitTag) {}
Float64X4(__m256d in_data) : data(in_data) {}
Float64X4(F64 in_data) : data(_mm256_set1_pd(in_data)) {}
Float64X4(const F64 *in_data) : data(_mm256_load_pd(in_data)) {}
F64X4 operator+(const F64X4 &other) const{
return _mm256_add_pd(data, other.data);}
F64X4 operator-(const F64X4 &other) const{
return _mm256_sub_pd(data, other.data);}
F64X4 operator*(const F64X4 &other) const{
return _mm256_mul_pd(data, other.data);}
F64X4 operator/(const F64X4 &other) const{
return _mm256_div_pd(data, other.data);}
static F64X4 fmadd(const F64X4 &a, const F64X4 &b, const F64X4 &c){
return _mm256_fmadd_pd(a.data, b.data, c.data);}
static F64X4 fmsub(const F64X4 &a, const F64X4 &b, const F64X4 &c){
return _mm256_fmsub_pd(a.data, b.data, c.data);}
template <int N>
F64X4 permute4x64() const{
return _mm256_permute4x64_pd(data, N);}
F64X4 reverse() const{
return permute4x64<0b00011011>();}
void load(const F64 *p){
data = _mm256_load_pd(p);}
void load1(const F64 *p){
data = _mm256_broadcast_sd(p);}
void store(F64 *p) const{
_mm256_store_pd(p, data);}
void store_nt(F64 *p) const{
_mm256_stream_pd(p, data);}
operator __m256d() const{
return data;}
__m256i toI64X4() const{
constexpr uint64_t mask = (uint64_t(1) << 52) - 1;
constexpr uint64_t offset = (uint64_t(1) << 10) - 1;
const __m256i f64bits = _mm256_castpd_si256(data);
__m256i tail = _mm256_and_si256(f64bits, _mm256_set1_epi64x(mask));
tail = _mm256_or_si256(tail, _mm256_set1_epi64x(mask + 1));
__m256i exp = _mm256_srli_epi64(f64bits, 52);
exp = _mm256_sub_epi64(_mm256_set1_epi64x(offset + 52), exp);
return _mm256_srlv_epi64(tail, exp);}
private:
__m256d data;};
struct Complex64X4{
using C64X4 = Complex64X4;
using F64X4 = Float64X4;
using F64 = double;
Complex64X4() {}
Complex64X4(NoInitTag) : real(NoInitTag{}), imag(NoInitTag{}) {}
Complex64X4(F64X4 real, F64X4 imag) : real(real), imag(imag) {}
Complex64X4(const F64 *p) : real(p), imag(p + 4) {}
Complex64X4(const F64 *p_real, const F64 *p_imag) : real(p_real), imag(p_imag) {}
C64X4 operator+(const C64X4 &other) const{
return C64X4(real + other.real, imag + other.imag);}
C64X4 operator-(const C64X4 &other) const{
return C64X4(real - other.real, imag - other.imag);}
C64X4 operator*(const F64X4 &other) const{
return C64X4(real * other, imag * other);}
C64X4 mul(const C64X4 &other) const{
const F64X4 ii = imag * other.imag;
const F64X4 ri = real * other.imag;
const F64X4 r = F64X4::fmsub(real, other.real, ii);
const F64X4 i = F64X4::fmadd(imag, other.real, ri);
return C64X4(r, i);}
C64X4 mulConj(const C64X4 &other) const{
const F64X4 ii = imag * other.imag;
const F64X4 ri = real * other.imag;
const F64X4 r = F64X4::fmadd(real, other.real, ii);
const F64X4 i = F64X4::fmsub(imag, other.real, ri);
return C64X4(r, i);}
C64X4 reverse() const{
return C64X4(real.reverse(), imag.reverse());}
void set1(F64 real_in, F64 imag_in){
real = F64X4(real_in);
imag = F64X4(imag_in);}
template <typename T>
void load(const T *p, std::false_type){
this->load(p);}
template <typename T>
void load(const T *p, std::true_type){
this->load(p);
*this = this->toRRIIPermu();}
template <typename T>
void load(const T *p){
real.load(reinterpret_cast<const F64 *>(p));
imag.load(reinterpret_cast<const F64 *>(p) + 4);}
void load1(const F64 *real_p, const F64 *imag_p){
real.load1(real_p);
imag.load1(imag_p);}
template <typename T>
void store(T *p, std::false_type) const{
this->store(p);}
template <typename T>
void store(T *p, std::true_type) const{
this->toRIRIPermu().store(p);}
template <typename T>
void store(T *p) const{
real.store(reinterpret_cast<F64 *>(p));
imag.store(reinterpret_cast<F64 *>(p) + 4);}
template <typename T>
void store_nt(T *p) const{
real.store_nt(reinterpret_cast<F64 *>(p));
imag.store_nt(reinterpret_cast<F64 *>(p) + 4);}
C64X4 toRIRIPermu() const{
C64X4 res = *this;
transpose64_2X4(res.real, res.imag);
return res;}
C64X4 toRRIIPermu() const{
C64X4 res = *this;
transpose64_4X2(res.real, res.imag);
return res;}
C64X4 transToI64(std::false_type) const{
return *this;}
C64X4 transToI64(std::true_type) const{
const __m256d magic=_mm256_set1_pd(6755399441055744.0);
const __m256i bits=_mm256_castpd_si256(magic);return C64X4(_mm256_castsi256_pd(_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(real,magic)),bits)),_mm256_castsi256_pd(_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(imag,magic)),bits)));}
F64X4 real, imag;};
using Float32 = float;
using Float64 = double;
constexpr Float64 HINT_PI = 3.141592653589793238462643;
constexpr Float64 HINT_2PI = HINT_PI * 2;
constexpr Float64 COS_PI_8 = 0.707106781186547524400844;
template <typename T>
constexpr T int_floor2(T n){
constexpr int bits = sizeof(n) * 8;
for (int i = 1; i < bits; i *= 2){
n |= (n >> i);}
return (n >> 1) + 1;}
template <typename T>
constexpr T int_ceil2(T n){
constexpr int bits = sizeof(n) * 8;
n--;
for (int i = 1; i < bits; i *= 2){
n |= (n >> i);}
return n + 1;}
template <typename IntTy>
constexpr bool is_2pow(IntTy n){
return n != 0 && (n & (n - 1)) == 0;}
template <typename T>
constexpr int hint_log2(T n){
constexpr int bits = sizeof(n) * 8;
int l = -1, r = bits;
while ((l + 1) != r){
int mid = (l + r) / 2;
if ((T(1) << mid) > n){
r = mid;}
else{
l = mid;}}
return l;}
constexpr uint32_t crc32(const char *str){
uint32_t crc = 0xFFFFFFFF;
while (*str != '\0'){
crc ^= *str;
for (int i = 0; i < 8; ++i){
crc = (crc >> 1) ^ (0 - (crc & 1)) & 0xEDB88320;}
str++;}
return ~crc;}
namespace transform{
template <typename T>
inline void transform2(T &sum, T &diff){
T temp0 = sum, temp1 = diff;
sum = temp0 + temp1;
diff = temp0 - temp1;}
namespace fft{
using F64 = Float64;
using C64 = std::complex<F64>;
using F64X4 = Float64X4;
using C64X4 = Complex64X4;
// Direct long-double Taylor root generation. No long trigonometric recurrence
// and no dependency on the platform accuracy of sinl/cosl. Runtime FFT remains double.
inline C64 certifiedRoot(size_t n,size_t index){
index%=n;
const size_t q=(4*index+n/2)/n;
const int64_t d=int64_t(4*index)-int64_t(q*n);
constexpr long double pi=3.141592653589793238462643383279502884L;
const long double x=(static_cast<long double>(d)/static_cast<long double>(n))*(pi/2);
const long double xx=x*x;
long double sn=x,cs=1,st=x,ct=1;
for(unsigned k=1;k<=18;++k){
st*=(-xx)/static_cast<long double>((2*k)*(2*k+1));
ct*=(-xx)/static_cast<long double>((2*k-1)*(2*k));
sn+=st;cs+=ct;
}
long double re,im;
switch(q&3){case 0:re=cs;im=sn;break;case 1:re=-sn;im=cs;break;case 2:re=-cs;im=-sn;break;default:re=sn;im=-cs;break;}
return C64(static_cast<double>(re),-static_cast<double>(im));
}
template <typename Float, size_t OMEGA_LEN>
class TableFix{
alignas(64) Float table[OMEGA_LEN * 2];
public:
TableFix(size_t theta_divider, size_t factor, size_t stride){
for(size_t begin=0,index=0;begin<OMEGA_LEN*2;begin+=stride*2)
for(size_t j=0;j<stride;++j,++index){C64 z=certifiedRoot(theta_divider,index*factor);table[begin+j]=Float(z.real());table[begin+j+stride]=Float(z.imag());}
}
constexpr const Float &operator[](size_t index) const{
return table[index];}};
void initOmegaX4(F64 *arr, size_t fft_len, int table_len, int factor){
for(int j=0;j<table_len;++j){C64 z=certifiedRoot(fft_len,size_t(j)*factor);size_t p=size_t(j/4)*8+j%4;arr[p]=z.real();arr[p+4]=z.imag();}
}
template <typename Float, int LOG_BEGIN, int LOG_END, int DIV>
class TableFixMulti{
static_assert(LOG_END >= LOG_BEGIN);
static_assert(is_2pow(DIV));
static constexpr size_t TABLE_CPX_LEN = (size_t(1) << (LOG_END + 1)) / DIV;
alignas(64) Float table[TABLE_CPX_LEN * 2];
public:
TableFixMulti(size_t factor, size_t stride = 4){
initBottomUp(factor, stride);}
void initBottomUp(size_t factor, size_t stride){
(void)stride;static_assert(std::is_same<Float,Float64>::value);
for(int log_len=LOG_BEGIN;log_len<=LOG_END;++log_len){size_t len=size_t(1)<<log_len;initOmegaX4(getBeginLog(log_len),len,len/DIV,factor);}
}
constexpr const Float *getBeginLog(int log_rank) const{
return getBegin(size_t(1) << log_rank);}
constexpr Float *getBeginLog(int log_rank){
return getBegin(size_t(1) << log_rank);}
constexpr const Float *getBegin(size_t rank) const{
return &table[rank * 2 / DIV];}
constexpr Float *getBegin(size_t rank){
return &table[rank * 2 / DIV];}};
template <int CACHE_LOG_LEN>
class FFTSqrtTableC64X4{
public:
using F64 = double;
using C64 = std::complex<double>;
using C64X4 = hint::Complex64X4;
static constexpr size_t CACHE_LEN = size_t(1) << CACHE_LOG_LEN;
static constexpr size_t MASK = CACHE_LEN - 1;
static constexpr size_t C4_COUNT = sizeof(C64X4) / sizeof(C64);
~FFTSqrtTableC64X4(){
if (high)
{
delete[] high;}}
FFTSqrtTableC64X4() {}
FFTSqrtTableC64X4(size_t fft_len, int len_div, int factor){
init(fft_len, len_div, factor);}
void init(size_t fft_len, int len_div, int factor){
size_t table_len=fft_len/len_div,low_len=CACHE_LEN*C4_COUNT,high_len=table_len/low_len;
if(high!=nullptr)delete[]high;high=new C64[high_len];
initOmegaX4(reinterpret_cast<F64*>(&low[0]),fft_len,low_len,factor);
for(size_t i=0;i<high_len;++i)high[i]=certifiedRoot(fft_len,i*low_len*factor);
}
static constexpr size_t LOGV=CACHE_LOG_LEN,MASKV=MASK;
const C64X4* lowPtr() const{return low;}
const C64* highPtr() const{return high;}
C64X4 operator[](size_t i) const{
C64X4 hi(NoInitTag{});
auto p = reinterpret_cast<const F64 *>(&high[i >> CACHE_LOG_LEN]);
hi.load1(p, p + 1);
return low[i & MASK].mul(hi);}
private:
alignas(64) C64X4 low[CACHE_LEN];
C64 *high = nullptr;};
template <int DIV, int LOG_BEGIN, int LOG_MAX, int CACHE_LOG_LEN>
class FFTTableSqrt{
using TableLong = FFTSqrtTableC64X4<CACHE_LOG_LEN>;
static constexpr size_t SHORT_LEN = size_t(1) << LOG_BEGIN;
static constexpr size_t TABLE_LEN = LOG_MAX - LOG_BEGIN + 1;
public:
FFTTableSqrt(int factor){
for (int i = 0; i < TABLE_LEN; i++)
{
size_t fft_len = SHORT_LEN << i;
table[i].init(fft_len, DIV, factor);}}
const TableLong &operator[](int log_len) const{
log_len -= LOG_BEGIN;
return table[log_len];}
TableLong &operator[](int log_len){
log_len -= LOG_BEGIN;
return table[log_len];}
private:
TableLong table[TABLE_LEN];};
struct FFT{
template <typename Float>
static void trans2MulI(Float &r0, Float &i0, Float &r1, Float &i1){
auto temp = r1;
r1 = r0 + i1;
r0 = r0 - i1;
i1 = i0 - temp;
i0 = i0 + temp;}
template <typename Float>
static void trans2MulNegI(Float &r0, Float &i0, Float &r1, Float &i1){
auto temp = r1;
r1 = r0 - i1;
r0 = r0 + i1;
i1 = i0 + temp;
i0 = i0 - temp;}
template <typename Float>
static void dif4(Float &r0, Float &i0, Float &r1, Float &i1, Float &r2, Float &i2, Float &r3, Float &i3){
difSplit(r0, i0, r1, i1, r2, i2, r3, i3);
transform2(r0, r1);
transform2(i0, i1);}
template <typename Float>
static void idit4(Float &r0, Float &i0, Float &r1, Float &i1, Float &r2, Float &i2, Float &r3, Float &i3){
transform2(r0, r1);
transform2(i0, i1);
iditSplit(r0, i0, r1, i1, r2, i2, r3, i3);}
template <typename Float>
static void difSplit(Float &r0, Float &i0, Float &r1, Float &i1, Float &r2, Float &i2, Float &r3, Float &i3){
transform2(r0, r2);
transform2(i0, i2);
transform2(r1, r3);
transform2(i1, i3);
trans2MulNegI(r2, i2, r3, i3);}
template <typename Float>
static void iditSplit(Float &r0, Float &i0, Float &r1, Float &i1, Float &r2, Float &i2, Float &r3, Float &i3){
transform2(r2, r3);
transform2(i2, i3);
transform2(r0, r2);
transform2(i0, i2);
trans2MulI(r1, i1, r3, i3);}};
namespace paired_leaf { extern bool enabled; void run(F64*,size_t); }
static bool g_skip_leaf=false;
struct FFTAVX : public FFT{
static constexpr int LOG_SHORT = 11, LOG_MID = 17, LOG_MAX = 23, LOG_CACHE = 8;
static constexpr size_t SHORT_LEN = size_t(1) << LOG_SHORT, MID_LEN = size_t(1) << LOG_MID, MAX_LEN = size_t(1) << LOG_MAX;
using TableFix4 = const TableFix<Float64, 4>;
using TableFix8 = const TableFix<Float64, 8>;
using TableMulti1 = const TableFixMulti<Float64, LOG_SHORT + 1, LOG_MID, 4>;
using TableMulti2 = const TableFixMulti<Float64, 6, LOG_SHORT + 1, 4>;
using TableMulti3 = const TableFixMulti<Float64, 6, LOG_SHORT, 4>;
using TableSqrt = const FFTTableSqrt<4, LOG_MID + 1, LOG_MAX, LOG_CACHE>;
static TableFix4 table_8, table_16_1, table_16_3;
static TableFix8 table_32_1, table_32_3;
static TableMulti2 multi_table_2;
static TableMulti3 multi_table_3;
static TableMulti1 multi_table_1, multi_table_mid2, multi_table_mid3;
static TableSqrt sqrt_table_1;
static constexpr const Float64 *it8 = &table_8[0], *it16_1 = &table_16_1[0], *it16_3 = &table_16_3[0], *it32_1 = &table_32_1[0], *it32_3 = &table_32_3[0];
static void dif4x4(F64X4 &r0, F64X4 &i0, F64X4 &r1, F64X4 &i1, F64X4 &r2, F64X4 &i2, F64X4 &r3, F64X4 &i3){
transpose64_4X4(r0, r1, r2, r3);
transpose64_4X4(i0, i1, i2, i3);
dif4(r0, i0, r1, i1, r2, i2, r3, i3);
transpose64_4X4(r0, r1, r2, r3);
transpose64_4X4(i0, i1, i2, i3);}
static void idit4x4(F64X4 &r0, F64X4 &i0, F64X4 &r1, F64X4 &i1, F64X4 &r2, F64X4 &i2, F64X4 &r3, F64X4 &i3){
transpose64_4X4(r0, r1, r2, r3);
transpose64_4X4(i0, i1, i2, i3);
idit4(r0, i0, r1, i1, r2, i2, r3, i3);
transpose64_4X4(r0, r1, r2, r3);
transpose64_4X4(i0, i1, i2, i3);}
static void dif8x2(C64X4 &c0, C64X4 &c1, C64X4 &c2, C64X4 &c3){
C64X4 omega(it8);
transform2(c0, c1);
transform2(c2, c3);
c1 = c1.mul(omega), c3 = c3.mul(omega);
dif4x4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);}
static void idit8x2(C64X4 &c0, C64X4 &c1, C64X4 &c2, C64X4 &c3){
C64X4 omega(it8);
idit4x4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c1 = c1.mulConj(omega), c3 = c3.mulConj(omega);
transform2(c0, c1);
transform2(c2, c3);}
static void dif16(Float64 in_out[]){
auto p = reinterpret_cast<C64X4 *>(in_out);
C64X4 c0 = p[0], c1 = p[1], c2 = p[2], c3 = p[3];
dif4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c1 = c1.mul(C64X4(it8)), c2 = c2.mul(C64X4(it16_1)), c3 = c3.mul(C64X4(it16_3));
dif4x4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
p[0] = c0, p[1] = c1, p[2] = c2, p[3] = c3;}
static void idit16(Float64 in_out[]){
auto p = reinterpret_cast<C64X4 *>(in_out);
C64X4 c0 = p[0], c1 = p[1], c2 = p[2], c3 = p[3], omega;
idit4x4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c1 = c1.mulConj(C64X4(it8)), c2 = c2.mulConj(C64X4(it16_1)), c3 = c3.mulConj(C64X4(it16_3));
idit4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
p[0] = c0, p[1] = c1, p[2] = c2, p[3] = c3;}
static void dif32(Float64 in_out[]){
auto p = reinterpret_cast<C64X4 *>(in_out);
C64X4 c0 = p[0], c1 = p[2], c2 = p[4], c3 = p[6];
difSplit(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c2 = c2.mul(C64X4(it32_1)), c3 = c3.mul(C64X4(it32_3));
p[0] = c0, p[2] = c1, p[4] = c2, p[6] = c3;
c0 = p[1], c1 = p[3], c2 = p[5], c3 = p[7];
difSplit(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c2 = c2.mul(C64X4(it32_1 + 8)), c3 = c3.mul(C64X4(it32_3 + 8));
p[1] = c0, p[3] = c1, c0 = p[4], c1 = p[6];
dif8x2(c0, c2, c1, c3);
p[4] = c0, p[5] = c2, p[6] = c1, p[7] = c3;
dif16(in_out);}
static void idit32(Float64 in_out[]){
idit16(in_out);
auto p = reinterpret_cast<C64X4 *>(in_out);
C64X4 c0 = p[4], c1 = p[5], c2 = p[6], c3 = p[7];
idit8x2(c0, c1, c2, c3);
p[5] = c1, p[7] = c3, c1 = p[0], c3 = p[2];
c0 = c0.mulConj(C64X4(it32_1)), c2 = c2.mulConj(C64X4(it32_3));
iditSplit(c1.real, c1.imag, c3.real, c3.imag, c0.real, c0.imag, c2.real, c2.imag);
p[0] = c1, p[2] = c3, p[4] = c0, p[6] = c2;
c0 = p[1], c1 = p[3], c2 = p[5], c3 = p[7];
c2 = c2.mulConj(C64X4(it32_1 + 8)), c3 = c3.mulConj(C64X4(it32_3 + 8));
iditSplit(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
p[1] = c0, p[3] = c1, p[5] = c2, p[7] = c3;}
template <typename F32, typename F16>
static void fftTiny(Float64 in_out[], size_t float_len, F32 &&func32, F16 &&func16){
if (hint_log2(float_len / 2) % 2 == 0)
{
for (auto end = in_out + float_len; in_out < end; in_out += 32)
{
func16(in_out);}}
else
{
for (auto end = in_out + float_len; in_out < end; in_out += 64)
{
func32(in_out);}}}
#define E8Z3_FUSE_MINRANK 256u
static void difTop2(Float64 *in_out, size_t float_len, size_t R){
const size_t S1 = R / 2, S2 = R / 8;
for (Float64 *begin = in_out, *bend = in_out + float_len; begin < bend; begin += R * 2) {
const Float64 *t1a = multi_table_2.getBegin(R * 2), *t2a = multi_table_2.getBegin(R), *t3a = multi_table_3.getBegin(R);
const Float64 *t1b = multi_table_2.getBegin(R / 2), *t2b = multi_table_2.getBegin(R / 4), *t3b = multi_table_3.getBegin(R / 4);
for (size_t j2 = 0; j2 < S2; j2 += 8) {
C64X4 v[4][4]={NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64};
for (int q = 0; q < 4; q++) {
const size_t off = (size_t)q * S2 + j2;
C64X4 c0(begin + off), c1(begin + S1 + off), c2(begin + 2 * S1 + off), c3(begin + 3 * S1 + off);
dif4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c1 = c1.mul(C64X4(t2a + off)); c2 = c2.mul(C64X4(t1a + off)); c3 = c3.mul(C64X4(t3a + off));
v[0][q] = c0; v[1][q] = c1; v[2][q] = c2; v[3][q] = c3;
}
for (int s = 0; s < 4; s++) {
C64X4 c0 = v[s][0], c1 = v[s][1], c2 = v[s][2], c3 = v[s][3];
dif4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c1 = c1.mul(C64X4(t2b + j2)); c2 = c2.mul(C64X4(t1b + j2)); c3 = c3.mul(C64X4(t3b + j2));
c0.store(begin + (size_t)s * S1 + j2);
c1.store(begin + (size_t)s * S1 + S2 + j2);
c2.store(begin + (size_t)s * S1 + 2 * S2 + j2);
c3.store(begin + (size_t)s * S1 + 3 * S2 + j2);
}
}
}
}
static void difIter(Float64 in_out[], size_t float_len){
size_t fft_len = float_len / 2;
C64X4 c0(NoInitTag{}), c1(NoInitTag{}), c2(NoInitTag{}), c3(NoInitTag{});
size_t stride = fft_len / 2;
auto it0 = in_out, it1 = it0 + stride, it2 = it1 + stride, it3 = it2 + stride;
if (fft_len >= E8Z3_FUSE_MINRANK) { difTop2(in_out, float_len, fft_len); }
for (size_t rank = fft_len >= E8Z3_FUSE_MINRANK ? fft_len / 16 : fft_len; rank >= 64; rank /= 4)
{
stride = rank / 2;
for (auto begin = in_out, end = in_out + float_len; begin < end; begin += rank * 2)
{
auto table1 = multi_table_2.getBegin(rank * 2), table2 = multi_table_2.getBegin(rank), table3 = multi_table_3.getBegin(rank);
it0 = begin, it1 = it0 + stride, it2 = it1 + stride, it3 = it2 + stride;
for (; it0 < begin + stride; it0 += 8, it1 += 8, it2 += 8, it3 += 8, table1 += 8, table2 += 8, table3 += 8)
{
c0 = it0, c1 = it1, c2 = it2, c3 = it3;
dif4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c1 = c1.mul(C64X4(table2)), c2 = c2.mul(C64X4(table1)), c3 = c3.mul(C64X4(table3));
c0.store(it0), c1.store(it1), c2.store(it2), c3.store(it3);
}}}
fftTiny(in_out, float_len, dif32, dif16);}
static void iditIter(Float64 in_out[], size_t float_len){
if(paired_leaf::enabled)paired_leaf::run(in_out,float_len);
size_t fft_len = float_len / 2;
size_t rank = hint_log2(fft_len) % 2 == 0 ? 64 : 128;
fftTiny(in_out, float_len, idit32, idit16);
for (; rank <= fft_len; rank *= 4)
{
const size_t stride = rank / 2;
for (auto begin = in_out, end = in_out + float_len; begin < end; begin += rank * 2)
{
auto table1 = multi_table_2.getBegin(rank * 2), table2 = multi_table_2.getBegin(rank), table3 = multi_table_3.getBegin(rank);
auto it0 = begin, it1 = it0 + stride, it2 = it1 + stride, it3 = it2 + stride;
for (; it0 < begin + stride; it0 += 8, it1 += 8, it2 += 8, it3 += 8, table1 += 8, table2 += 8, table3 += 8)
{
C64X4 c0 = it0, c1 = it1, c2 = it2, c3 = it3;
c1 = c1.mulConj(C64X4(table2)), c2 = c2.mulConj(C64X4(table1)), c3 = c3.mulConj(C64X4(table3));
idit4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c0.store(it0), c1.store(it1), c2.store(it2), c3.store(it3);
}}}}
#define difLayer(dif_func, in_out, stride, table) \
do \
{ \
auto it0 = in_out, it1 = in_out + stride, it2 = it1 + stride, it3 = it2 + stride; \
size_t indx = 0; \
for (auto end = it1; FROM_RIRI_PERM && it0 < end; it0 += 8, it1 += 8, it2 += 8, it3 += 8, indx++) \
{ \
C64X4 c0, c1, c2, c3, omega1, omega2; \
c0.load(it0, FromRIRI{}), c1.load(it1, FromRIRI{}), c2 = c0,c3 = c1; \
transform2(c0,c1); \
trans2MulNegI(c2.real,c2.imag,c3.real,c3.imag); \
omega1 = table[indx], c2 = c2.mul(omega1); \
omega2 = fft_len<=MID_LEN?C64X4(multi_table_mid2.getBegin(fft_len)+8*indx):omega1.mul(omega1), c1 = c1.mul(omega2); \
c3 = c3.mul(fft_len<=MID_LEN?C64X4(multi_table_mid3.getBegin(fft_len)+8*indx):omega2.mul(omega1)); \
c0.store(it0), c1.store(it1), c2.store(it2), c3.store(it3); \
} \
for (auto end = it1; (!FROM_RIRI_PERM) && it0 < end; it0 += 8, it1 += 8, it2 += 8, it3 += 8, indx++) \
{ \
_mm_prefetch((const char*)(it0+64),_MM_HINT_T0); \
_mm_prefetch((const char*)(it1+64),_MM_HINT_T0); \
_mm_prefetch((const char*)(it2+64),_MM_HINT_T0); \
_mm_prefetch((const char*)(it3+64),_MM_HINT_T0); \
C64X4 c0 = it0, c1 = it1, c2 = it2, c3 = it3, omega1, omega2; \
dif4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag); \
omega1 = table[indx], c2 = c2.mul(omega1); \
omega2 = fft_len<=MID_LEN?C64X4(multi_table_mid2.getBegin(fft_len)+8*indx):omega1.mul(omega1), c1 = c1.mul(omega2); \
c3 = c3.mul(fft_len<=MID_LEN?C64X4(multi_table_mid3.getBegin(fft_len)+8*indx):omega2.mul(omega1)); \
c0.store(it0), c1.store(it1), c2.store(it2), c3.store(it3); \
} \
dif_func(in_out, stride); \
dif_func(in_out + stride, stride); \
dif_func(in_out + stride * 2, stride); \
dif_func(in_out + stride * 3, stride); \
} while (0)
#define iditLayer(idit_func, in_out, stride, table) \
do \
{ \
idit_func(in_out, stride); \
idit_func(in_out + stride, stride); \
idit_func(in_out + stride * 2, stride); \
idit_func(in_out + stride * 3, stride); \
auto it0 = in_out, it1 = in_out + stride, it2 = it1 + stride, it3 = it2 + stride; \
size_t indx = 0; \
_Pragma("GCC unroll 4") \
for (auto end = it1; it0 < end; it0 += 8, it1 += 8, it2 += 8, it3 += 8, indx++) \
{ \
_mm_prefetch((const char*)(it0+64),_MM_HINT_T0); \
_mm_prefetch((const char*)(it1+64),_MM_HINT_T0); \
_mm_prefetch((const char*)(it2+64),_MM_HINT_T0); \
_mm_prefetch((const char*)(it3+64),_MM_HINT_T0); \
C64X4 c0 = it0, c1 = it1, c2 = it2, c3 = it3, omega1, omega2; \
omega1 = table[indx], c2 = c2.mulConj(omega1); \
omega2 = fft_len<=MID_LEN?C64X4(multi_table_mid2.getBegin(fft_len)+8*indx):omega1.mul(omega1), c1 = c1.mulConj(omega2); \
c3 = c3.mulConj(fft_len<=MID_LEN?C64X4(multi_table_mid3.getBegin(fft_len)+8*indx):omega2.mul(omega1)); \
idit4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag); \
c0 = c0.transToI64(ToI64{}), c1 = c1.transToI64(ToI64{}), c2 = c2.transToI64(ToI64{}), c3 = c3.transToI64(ToI64{}); \
c0.store(it0, ToRIRI{}), c1.store(it1, ToRIRI{}), c2.store(it2, ToRIRI{}), c3.store(it3, ToRIRI{}); \
} \
} while (0)
template <bool FROM_RIRI_PERM = false>
static void difRecMid(Float64 in_out[], size_t float_len){
const size_t fft_len = float_len / 2;
if (fft_len <= SHORT_LEN)
{
difIter(in_out, float_len);
return;}
using FromRIRI = std::integral_constant<bool, FROM_RIRI_PERM>;
auto table1 = reinterpret_cast<const C64X4 *>(multi_table_1.getBegin(fft_len));
const size_t stride = float_len / 4;
difLayer(difRecMid, in_out, stride, table1);}
template <bool TO_RIRI_PERM = false, bool TO_INT64 = false>
static void iditRecMid(Float64 in_out[], size_t float_len){
const size_t fft_len = float_len / 2;
if (fft_len <= SHORT_LEN)
{
if(!g_skip_leaf)iditIter(in_out, float_len);
return;}
using ToRIRI = std::integral_constant<bool, TO_RIRI_PERM>;
using ToI64 = std::integral_constant<bool, TO_INT64>;
const size_t stride = float_len / 4;
auto table1 = reinterpret_cast<const C64X4 *>(multi_table_1.getBegin(fft_len));
iditLayer(iditRecMid, in_out, stride, table1);}
static void fusedLong(Float64* a,size_t len,bool inverse){
size_t stride=len/16;
const auto &t0=sqrt_table_1[hint_log2(len/2)];
const auto &t1=sqrt_table_1[hint_log2(len/8)];
if(inverse)for(int k=0;k<16;k++)iditRecLong<false,false>(a+k*stride,stride);
for(size_t j=0;j<stride;j+=8){
if(j>=64 && j+416<stride){
for(int k=0;k<16;k++) _mm_prefetch((const char*)(a+(size_t)k*stride+j-64+32*k),_MM_HINT_T0);
}
C64X4 v[16];
auto butterfly=[&](C64X4 &c0,C64X4 &c1,C64X4 &c2,C64X4 &c3,const C64X4 &w){
C64X4 w2=w.mul(w),w3=w2.mul(w);
if(inverse){
c1=c1.mulConj(w2);c2=c2.mulConj(w);c3=c3.mulConj(w3);
idit4(c0.real,c0.imag,c1.real,c1.imag,c2.real,c2.imag,c3.real,c3.imag);
}else{
dif4(c0.real,c0.imag,c1.real,c1.imag,c2.real,c2.imag,c3.real,c3.imag);
c1=c1.mul(w2);c2=c2.mul(w);c3=c3.mul(w3);
}
};
const C64X4 w=t1[j/8],w2=w.mul(w),w3=w2.mul(w);
auto butterflyCommon=[&](C64X4 &c0,C64X4 &c1,C64X4 &c2,C64X4 &c3){
if(inverse){
c1=c1.mulConj(w2);c2=c2.mulConj(w);c3=c3.mulConj(w3);
idit4(c0.real,c0.imag,c1.real,c1.imag,c2.real,c2.imag,c3.real,c3.imag);
}else{
dif4(c0.real,c0.imag,c1.real,c1.imag,c2.real,c2.imag,c3.real,c3.imag);
c1=c1.mul(w2);c2=c2.mul(w);c3=c3.mul(w3);
}
};
if(inverse){
#pragma GCC unroll 4
for(int q=0;q<4;q++){
v[4*q]=C64X4(a+(4*q)*stride+j);
v[4*q+1]=C64X4(a+(4*q+1)*stride+j);
v[4*q+2]=C64X4(a+(4*q+2)*stride+j);
v[4*q+3]=C64X4(a+(4*q+3)*stride+j);
butterflyCommon(v[4*q],v[4*q+1],v[4*q+2],v[4*q+3]);
}
#pragma GCC unroll 4
for(int k=0;k<4;k++){
butterfly(v[k],v[k+4],v[k+8],v[k+12],t0[(k*stride+j)/8]);
v[k].store(a+k*stride+j);
v[k+4].store(a+(k+4)*stride+j);
v[k+8].store(a+(k+8)*stride+j);
v[k+12].store(a+(k+12)*stride+j);
}
}else{
#pragma GCC unroll 4
for(int k=0;k<4;k++){
v[k]=C64X4(a+k*stride+j);
v[k+4]=C64X4(a+(k+4)*stride+j);
v[k+8]=C64X4(a+(k+8)*stride+j);
v[k+12]=C64X4(a+(k+12)*stride+j);
butterfly(v[k],v[k+4],v[k+8],v[k+12],t0[(k*stride+j)/8]);
}
#pragma GCC unroll 4
for(int q=0;q<4;q++){
butterflyCommon(v[4*q],v[4*q+1],v[4*q+2],v[4*q+3]);
v[4*q].store(a+(4*q)*stride+j);
v[4*q+1].store(a+(4*q+1)*stride+j);
v[4*q+2].store(a+(4*q+2)*stride+j);
v[4*q+3].store(a+(4*q+3)*stride+j);
}
}
}
if(!inverse)for(int k=0;k<16;k++)difRecLong<false>(a+k*stride,stride);
}
static void difFuseAB(Float64* a, size_t float_len){
const size_t fft_len = float_len/2;
const size_t sa = float_len/4;
const size_t sb = sa/4;
const auto &TA = sqrt_table_1[hint_log2(fft_len)];
const C64X4 *TB = reinterpret_cast<const C64X4 *>(multi_table_1.getBegin(sa/2));
for(size_t j0=0;j0<sb;j0+=8){
C64X4 A[4][4]={NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64,NI64};
const size_t jb=j0>>3;
const C64X4* talp=TA.lowPtr()+(jb&TA.MASKV);
const C64X4 talv=*talp;
const size_t tahbase=(jb>>TA.LOGV),tahstep=(sb>>3)>>TA.LOGV;
for(int k2=0;k2<4;k2++){
F64* pj = a + (size_t)k2*sb + j0;
_mm_prefetch((const char*)(pj+sa+128+8*k2),_MM_HINT_T0);
_mm_prefetch((const char*)(pj+2*sa+136+8*k2),_MM_HINT_T0);
_mm_prefetch((const char*)(pj+3*sa+144+8*k2),_MM_HINT_T0);
C64X4 c0(pj), c1(pj+sa), c2(pj+2*sa), c3(pj+3*sa);
dif4(c0.real,c0.imag,c1.real,c1.imag,c2.real,c2.imag,c3.real,c3.imag);
const C64* tahp=TA.highPtr()+tahbase+(size_t)k2*tahstep;
C64X4 tah(NoInitTag{});{const F64* hq=reinterpret_cast<const F64*>(tahp);tah.load1(hq,hq+1);}
C64X4 w = talv.mul(tah), w2 = w.mul(w);
c2 = c2.mul(w); c1 = c1.mul(w2); c3 = c3.mul(w2.mul(w));
A[k2][0]=c0; A[k2][1]=c1; A[k2][2]=c2; A[k2][3]=c3;
}
const C64X4 wB=TB[j0>>3],wB2=wB.mul(wB),wB3=wB2.mul(wB);
for(int k1=0;k1<4;k1++){
C64X4 c0=A[0][k1], c1=A[1][k1], c2=A[2][k1], c3=A[3][k1];
dif4(c0.real,c0.imag,c1.real,c1.imag,c2.real,c2.imag,c3.real,c3.imag);
c2=c2.mul(wB);c1=c1.mul(wB2);c3=c3.mul(wB3);
F64* q = a + (size_t)k1*sa + j0;
c0.store(q); c1.store(q+sb); c2.store(q+2*sb); c3.store(q+3*sb);
}
}
for(int k1=0;k1<4;k1++) for(int k2=0;k2<4;k2++)
difRecMid<false>(a + (size_t)k1*sa + (size_t)k2*sb, sb);
}
template <bool FROM_RIRI_PERM = false>
static void difRecLong(Float64 in_out[], size_t float_len){
if constexpr(!FROM_RIRI_PERM) if(float_len/2>=(1ULL<<22)){fusedLong(in_out,float_len,false);return;}
const size_t fft_len = float_len / 2;
if (fft_len <= MID_LEN)
{
difRecMid<FROM_RIRI_PERM>(in_out, float_len);
return;}
if constexpr(!FROM_RIRI_PERM)
if(float_len>(1ULL<<18) && float_len<=(1ULL<<20)){difFuseAB(in_out,float_len);return;}
using FromRIRI = std::integral_constant<bool, FROM_RIRI_PERM>;
const auto &table1 = sqrt_table_1[hint_log2(fft_len)];
const size_t stride = float_len / 4;
difLayer(difRecLong, in_out, stride, table1);}
template <bool TO_RIRI_PERM = false, bool TO_INT64 = false>
static void iditRecLong(Float64 in_out[], size_t float_len){
if constexpr(!TO_RIRI_PERM && !TO_INT64) if(float_len/2>=(1ULL<<22)){fusedLong(in_out,float_len,true);return;}
const size_t fft_len = float_len / 2;
if (fft_len <= MID_LEN)
{
iditRecMid<TO_RIRI_PERM, TO_INT64>(in_out, float_len);
return;}
using ToRIRI = std::integral_constant<bool, TO_RIRI_PERM>;
using ToI64 = std::integral_constant<bool, TO_INT64>;
const size_t stride = float_len / 4;
const auto &table1 = sqrt_table_1[hint_log2(fft_len)];
iditLayer(iditRecLong, in_out, stride, table1);}};
#undef difLayer
#undef iditLayer
constexpr int FFTAVX::LOG_SHORT, FFTAVX::LOG_MID, FFTAVX::LOG_MAX, FFTAVX::LOG_CACHE;
constexpr size_t FFTAVX::SHORT_LEN, FFTAVX::MID_LEN, FFTAVX::MAX_LEN;
FFTAVX::TableFix4 FFTAVX::table_8(8, 1, 4), FFTAVX::table_16_1(16, 1, 4), FFTAVX::table_16_3(16, 3, 4);
FFTAVX::TableFix8 FFTAVX::table_32_1(32, 1, 4), FFTAVX::table_32_3(32, 3, 4);
FFTAVX::TableMulti2 FFTAVX::multi_table_2(2);
FFTAVX::TableMulti3 FFTAVX::multi_table_3(3);
FFTAVX::TableMulti1 FFTAVX::multi_table_1(1), FFTAVX::multi_table_mid2(2), FFTAVX::multi_table_mid3(3);
FFTAVX::TableSqrt FFTAVX::sqrt_table_1(1);
namespace paired_leaf {bool enabled=false;inline void run(F64*,size_t){std::abort();}}
}}}
#pragma GCC push_options
#pragma GCC target("avx2,fma")
#pragma GCC optimize("O3")
namespace packed_fft {
using u32=uint32_t;using C=std::complex<double>;using FFT=hint::transform::fft::FFTAVX;
inline unsigned at(unsigned i){return 8*(i/4)+i%4;}
inline C get(const double*a,unsigned i){unsigned k=at(i);return {a[k],a[k+4]};}
inline void put(double*a,unsigned i,C v){unsigned k=at(i);a[k]=v.real();a[k+4]=v.imag();}
inline unsigned negslot(unsigned i){return i?i^((1u<<(31-__builtin_clz(i)))-1):0;}
struct Result{bool certified=false;long double bound=0;double observed_max_fraction=0;};
// out[t] = sum_i a[i+t]*b[i], 0<=t<h, modulo 2^32.
// A has 2h entries. Refusal leaves out unchanged. Each accepted integer is
// recovered by a posteriori norm error certificate, not a random-input guess.
class MiddlePlan {
unsigned h,n;hint::AlignMem<double>x,y,z;uint64_t normB=0;
using CV=hint::Complex64X4;using FV=hint::Float64X4;
static CV conjugate(CV v){v.imag=FV(_mm256_xor_pd(__m256d(v.imag),_mm256_set1_pd(-0.)));return v;}
static CV product(CV X,CV Xc,CV Y,CV Yc){return (X+Xc).mul(Y+Yc)*FV(.25)+(X.mul(Y)-Xc.mul(Yc))*FV(.5);}
public:
explicit MiddlePlan(const u32*b,unsigned size):h(size),n(2*size),x(4*size),y(4*size),z(4*size){
assert(h>=8&&h<=(1u<<20)&&(h&(h-1))==0);std::fill(y.begin(),y.begin()+2*n,0.);
for(unsigned i=0;i<h;++i){int lo=int16_t(b[i]),hi=int16_t((b[i]+32768u)>>16);put(y.begin(),h-1-i,{double(lo),double(hi)});normB+=int64_t(lo)*lo+int64_t(hi)*hi;}
FFT::difRecLong<false>(y.begin(),2*n);
}
Result apply(const u32*a,u32*out){
uint64_t normA=0;for(unsigned i=0;i<n;++i){int lo=int16_t(a[i]),hi=int16_t((a[i]+32768u)>>16);put(x.begin(),i,{double(lo),double(hi)});normA+=int64_t(lo)*lo+int64_t(hi)*hi;}
FFT::difRecLong<false>(x.begin(),2*n);long double normZ2=0;
for(unsigned i=0;i<4;++i){unsigned j=negslot(i);C X=get(x.begin(),i),Xc=std::conj(get(x.begin(),j)),Y=get(y.begin(),i),Yc=std::conj(get(y.begin(),j));C Z=.25*(X+Xc)*(Y+Yc)+.5*(X*Y-Xc*Yc);put(z.begin(),i,Z);normZ2+=static_cast<long double>(Z.real())*Z.real()+static_cast<long double>(Z.imag())*Z.imag();}
__m256d norms=_mm256_setzero_pd();for(unsigned block=4;block<n;block*=2)for(unsigned i=block;i<block*2;i+=4){unsigned j=3*block-4-i;CV X(x.begin()+2*i),Xc=conjugate(CV(x.begin()+2*j).reverse()),Y(y.begin()+2*i),Yc=conjugate(CV(y.begin()+2*j).reverse());CV Z=product(X,Xc,Y,Yc);Z.store(z.begin()+2*i);norms=_mm256_fmadd_pd(__m256d(Z.real),__m256d(Z.real),norms);norms=_mm256_fmadd_pd(__m256d(Z.imag),__m256d(Z.imag),norms);}
alignas(32)double sums[4];_mm256_store_pd(sums,norms);for(double value:sums)normZ2+=value;
constexpr long double u=0x1p-53L;long double op=(n<=(1u<<17)?5:64)*(unsigned(__builtin_ctz(n))+1),gamma=op*u/(1-op*u),normprod=sqrtl(static_cast<long double>(normA)*normB);long double normZ=sqrtl(normZ2)*(1+32*n*u);
long double error=(4*gamma+2*gamma*gamma+128*u*(1+gamma)*(1+gamma))*normprod+gamma*normZ/sqrtl(n);Result result;result.bound=error;if(!(error<.49L))return result;
if(_MM_GET_ROUNDING_MODE()!=_MM_ROUND_NEAREST)return result;
FFT::iditRecLong<false,false>(z.begin(),2*n);const __m256d scale=_mm256_set1_pd(1./n),magic=_mm256_set1_pd(0x1.8p52);const __m256i mb=_mm256_castpd_si256(magic),perm=_mm256_setr_epi32(0,2,4,6,1,3,5,7);unsigned t=0;
auto scalar=[&](unsigned atout){C v=get(z.begin(),h-1+atout)/double(n);out[atout]=u32(llround(v.real()))+(u32(llround(v.imag()))<<16);};scalar(t++);
for(;t+4<=h;t+=4){unsigned pos=2*(h-1+t);__m256d re=_mm256_mul_pd(_mm256_load_pd(z.begin()+pos),scale),im=_mm256_mul_pd(_mm256_load_pd(z.begin()+pos+4),scale);__m256i ri=_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(re,magic)),mb),ii=_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(im,magic)),mb);__m128i low=_mm256_castsi256_si128(_mm256_permutevar8x32_epi32(ri,perm)),cross=_mm256_castsi256_si128(_mm256_permutevar8x32_epi32(ii,perm));_mm_storeu_si128((__m128i*)(out+t),_mm_add_epi32(low,_mm_slli_epi32(cross,16)));}
for(;t<h;++t)scalar(t);result.certified=true;return result;
}
};
inline Result middle(const u32*a,const u32*b,u32*out,unsigned h){MiddlePlan plan(b,h);return plan.apply(a,out);}
}
#pragma GCC pop_options
#pragma GCC push_options
#pragma GCC target("avx2,fma")
#pragma GCC optimize("O3")
namespace packed_fft {
struct PreparedA{
unsigned h,n;uint64_t normA=0;hint::AlignMem<double>spectrum;
PreparedA(const u32*a,unsigned size):h(size),n(2*size),spectrum(4*size){assert(h>=8&&h<=(1u<<20)&&(h&(h-1))==0);for(unsigned i=0;i<n;++i){int lo=int16_t(a[i]),hi=int16_t((a[i]+32768u)>>16);put(spectrum.begin(),i,{double(lo),double(hi)});normA+=int64_t(lo)*lo+int64_t(hi)*hi;}FFT::difRecLong<false>(spectrum.begin(),2*n);}
};
class AReuseWorkspace{
unsigned h,n;hint::AlignMem<double>y,z;
using CV=hint::Complex64X4;using FV=hint::Float64X4;
static CV conjugate(CV v){v.imag=FV(_mm256_xor_pd(__m256d(v.imag),_mm256_set1_pd(-0.)));return v;}
static CV product(CV X,CV Xc,CV Y,CV Yc){return (X+Xc).mul(Y+Yc)*FV(.25)+(X.mul(Y)-Xc.mul(Yc))*FV(.5);}
public:
explicit AReuseWorkspace(unsigned size):h(size),n(2*size),y(4*size),z(4*size){}
Result apply(const PreparedA&pa,const u32*b,u32*out){assert(pa.h==h);uint64_t normA=pa.normA,normB=0;std::fill(y.begin(),y.begin()+2*n,0.);for(unsigned i=0;i<h;++i){int lo=int16_t(b[i]),hi=int16_t((b[i]+32768u)>>16);put(y.begin(),h-1-i,{double(lo),double(hi)});normB+=int64_t(lo)*lo+int64_t(hi)*hi;}FFT::difRecLong<false>(y.begin(),2*n);
long double normZ2=0;
for(unsigned i=0;i<4;++i){unsigned j=negslot(i);C X=get(pa.spectrum.begin(),i),Xc=std::conj(get(pa.spectrum.begin(),j)),Y=get(y.begin(),i),Yc=std::conj(get(y.begin(),j));C Z=.25*(X+Xc)*(Y+Yc)+.5*(X*Y-Xc*Yc);put(z.begin(),i,Z);normZ2+=static_cast<long double>(Z.real())*Z.real()+static_cast<long double>(Z.imag())*Z.imag();}
__m256d norms=_mm256_setzero_pd();for(unsigned block=4;block<n;block*=2)for(unsigned i=block;i<block*2;i+=4){unsigned j=3*block-4-i;CV X(pa.spectrum.begin()+2*i),Xc=conjugate(CV(pa.spectrum.begin()+2*j).reverse()),Y(y.begin()+2*i),Yc=conjugate(CV(y.begin()+2*j).reverse());CV Z=product(X,Xc,Y,Yc);Z.store(z.begin()+2*i);norms=_mm256_fmadd_pd(__m256d(Z.real),__m256d(Z.real),norms);norms=_mm256_fmadd_pd(__m256d(Z.imag),__m256d(Z.imag),norms);}
alignas(32)double sums[4];_mm256_store_pd(sums,norms);for(double value:sums)normZ2+=value;
constexpr long double u=0x1p-53L;long double op=(n<=(1u<<17)?5:64)*(unsigned(__builtin_ctz(n))+1),gamma=op*u/(1-op*u),normprod=sqrtl(static_cast<long double>(normA)*normB);long double normZ=sqrtl(normZ2)*(1+32*n*u);
long double error=(4*gamma+2*gamma*gamma+128*u*(1+gamma)*(1+gamma))*normprod+gamma*normZ/sqrtl(n);Result result;result.bound=error;if(!(error<.49L))return result;
if(_MM_GET_ROUNDING_MODE()!=_MM_ROUND_NEAREST)return result;
FFT::iditRecLong<false,false>(z.begin(),2*n);const __m256d scale=_mm256_set1_pd(1./n),magic=_mm256_set1_pd(0x1.8p52);const __m256i mb=_mm256_castpd_si256(magic),perm=_mm256_setr_epi32(0,2,4,6,1,3,5,7);unsigned t=0;
auto scalar=[&](unsigned atout){C v=get(z.begin(),h-1+atout)/double(n);out[atout]=u32(llround(v.real()))+(u32(llround(v.imag()))<<16);};scalar(t++);
for(;t+4<=h;t+=4){unsigned pos=2*(h-1+t);__m256d re=_mm256_mul_pd(_mm256_load_pd(z.begin()+pos),scale),im=_mm256_mul_pd(_mm256_load_pd(z.begin()+pos+4),scale);__m256i ri=_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(re,magic)),mb),ii=_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(im,magic)),mb);__m128i low=_mm256_castsi256_si128(_mm256_permutevar8x32_epi32(ri,perm)),cross=_mm256_castsi256_si128(_mm256_permutevar8x32_epi32(ii,perm));_mm_storeu_si128((__m128i*)(out+t),_mm_add_epi32(low,_mm_slli_epi32(cross,16)));}
for(;t<h;++t)scalar(t);result.certified=true;return result;
}
};
}
#pragma GCC pop_options
#include <memory>
#include <vector>
#include <cstring>
#include <algorithm>
namespace tiled_mod {
using u32=uint32_t;using u64=uint64_t;
struct State{
const u32*a;unsigned n,H;packed_fft::AReuseWorkspace work;
std::vector<std::unique_ptr<packed_fft::PreparedA>> cache;
std::vector<u32>input,D,tmp;
u64 accepted=0,refused=0,prepared=0,tiles=0;long double max_bound=0;
State(unsigned size,const u32*data,unsigned h):a(data),n(size),H(h),work(h),cache((size+h-1)/h+2),input(2*h),D(h),tmp(h){}
const packed_fft::PreparedA&get(unsigned block){
auto&slot=cache[block];if(!slot){unsigned start=block*H,count=start<n?std::min(2*H,n-start):0;
std::fill(input.begin(),input.end(),0);if(count)std::copy(a+start,a+start+count,input.begin());
slot=std::make_unique<packed_fft::PreparedA>(input.data(),H);++prepared;
}return *slot;
}
// Exact all-q cross: k in [l,m), j in [m,r), q>=1, i=q*j+k<n.
void cross(const u32*b,u32*c,unsigned l,unsigned m,unsigned r){
if(m>=n-l||l>=m||m>=r)return;r=std::min(r,n);unsigned wk=std::min(m-l,n-l-m),last=n-1-l,qmax=last/m;
std::vector<std::pair<unsigned,unsigned>> ranges;
for(unsigned q=1;q<=qmax;++q){unsigned lo=u64(q)*m/H,hi=std::min<u64>(u64(q)*(r-1),last)/H;
if(!ranges.empty()&&lo<=ranges.back().second+1)ranges.back().second=hi;else ranges.push_back({lo,hi});
}
for(auto range:ranges)for(unsigned block=range.first;block<=range.second;++block){
unsigned dlo=block*H,dhi=std::min(last,dlo+H-1),qlo=std::max(1u,(dlo+r-2)/(r-1)),qhi=dhi/m;bool any=false;
std::fill(D.begin(),D.end(),0);
for(unsigned q=qlo;q<=qhi;++q){unsigned first=std::max(m,(dlo+q-1)/q),end=std::min(r-1,dhi/q);
if(first<=end)any=true;for(unsigned j=first;j<=end;++j)D[q*j-dlo]+=b[j];
}
if(!any)continue;++tiles;
unsigned kb0=l/H,kb1=(l+wk-1)/H;
for(unsigned kb=kb0;kb<=kb1;++kb){unsigned ast=(block+kb)*H;if(ast>=n)continue;
unsigned klo=std::max(l,kb*H),khi=std::min(l+wk,(kb+1)*H),offset=klo-kb*H;
auto result=work.apply(get(block+kb),D.data(),tmp.data());max_bound=std::max(max_bound,result.bound);
if(result.certified){++accepted;for(unsigned k=klo;k<khi;++k)c[k]+=tmp[k-kb*H];}
else {++refused;for(unsigned i=0;i<H&&dlo+i+klo<n;++i){u32 value=D[i];if(!value)continue;unsigned end=std::min(khi,n-dlo-i);for(unsigned k=klo;k<end;++k)c[k]+=a[dlo+i+k]*value;}}
}
}
}
};
}
namespace tiled_mod {
struct Cross{unsigned l,m,r;};
// Schedule small output intervals by their A window. PreparedA is constructed
// once per nonempty bucket and released immediately, so its storage is O(H),
// independent of n. Each task's D tile is generated from its code-only ranges.
struct StreamState{
const u32*a;unsigned n,H;u64 accepted=0,refused=0,prepared=0,tasks=0;long double max_bound=0;
StreamState(unsigned size,const u32*data,unsigned h):a(data),n(size),H(h){}
void run(const u32*b,u32*c,const std::vector<Cross>&nodes){
std::vector<std::vector<unsigned>> buckets((n+H-1)/H+2);
for(unsigned id=0;id<nodes.size();++id){auto node=nodes[id];unsigned l=node.l,m=node.m,r=std::min(node.r,n);
if(l>=m||m>=r||m>=n-l)continue;unsigned wk=std::min(m-l,n-l-m),kb=l/H,last=n-1-l;
if((l+wk-1)/H!=kb)throw std::invalid_argument("stream output crosses an H boundary");
unsigned qmax=last/m;std::vector<std::pair<unsigned,unsigned>> ranges;
for(unsigned q=1;q<=qmax;++q){unsigned lo=u64(q)*m/H,hi=std::min<u64>(u64(q)*(r-1),last)/H;
if(!ranges.empty()&&lo<=ranges.back().second+1)ranges.back().second=hi;else ranges.push_back({lo,hi});
}
for(auto range:ranges)for(unsigned block=range.first;block<=range.second;++block)if(u64(block+kb)*H<n){buckets[block+kb].push_back(id);++tasks;}
}
packed_fft::AReuseWorkspace work(H);std::vector<u32>input(2*H),D(H),tmp(H);
for(unsigned ast=0;ast<buckets.size();++ast){if(buckets[ast].empty())continue;unsigned start=ast*H,count=std::min(2*H,n-start);
std::fill(input.begin(),input.end(),0);std::copy(a+start,a+start+count,input.begin());packed_fft::PreparedA pa(input.data(),H);++prepared;
for(unsigned id:buckets[ast]){auto node=nodes[id];unsigned l=node.l,m=node.m,r=std::min(node.r,n),wk=std::min(m-l,n-l-m),kb=l/H,block=ast-kb,dlo=block*H,dhi=std::min(n-1-l,dlo+H-1);bool any=false;
std::fill(D.begin(),D.end(),0);unsigned qlo=std::max(1u,(dlo+r-2)/(r-1)),qhi=dhi/m;
for(unsigned q=qlo;q<=qhi;++q){unsigned first=std::max(m,(dlo+q-1)/q),end=std::min(r-1,dhi/q);if(first<=end)any=true;for(unsigned j=first;j<=end;++j)D[q*j-dlo]+=b[j];}
if(!any)continue;auto result=work.apply(pa,D.data(),tmp.data());max_bound=std::max(max_bound,result.bound);
if(result.certified){++accepted;for(unsigned k=l;k<l+wk;++k)c[k]+=tmp[k-kb*H];}
else {++refused;for(unsigned i=0;i<H&&dlo+i+l<n;++i){u32 value=D[i];if(!value)continue;unsigned end=std::min(l+wk,n-dlo-i);for(unsigned k=l;k<end;++k)c[k]+=a[dlo+i+k]*value;}}
}
}
}
};
}
namespace partial_fft_inplace {
using packed_fft::C;using packed_fft::get;using packed_fft::put;using packed_fft::FFT;
class PartialInverse {
using CV=hint::Complex64X4;using FV=hint::Float64X4;
unsigned n,start,len,base;
std::vector<std::unique_ptr<hint::AlignMem<double>>> roots;
public:
PartialInverse(unsigned size,unsigned first,unsigned count):n(size),start(first),len(count),base(16){
assert(len&&start+len<=n);while(base<len)base*=2;assert(base<=n);
for(unsigned step=base;step<n;step*=2){auto table=std::make_unique<hint::AlignMem<double>>(2*base);for(unsigned t=0;t<len;++t)put(table->begin(),(start+t)&(base-1),std::conj(hint::transform::fft::certifiedRoot(2*step,(start+t)%(2*step))));roots.push_back(std::move(table));}
}
unsigned size()const{return len;}
unsigned block_size()const{return base;}
// Input is the Hint binary-bit-reversed spectrum in four-lane AoSoA layout.
// Return contains precisely outputs start..start+len, unnormalized, in AoSoA.
const double*apply(double*spectrum){
for(unsigned block=0;block<n;block+=base)FFT::iditRecLong<false,false>(spectrum+2*block,2*base);
unsigned rank=0;for(unsigned step=base;step<n;step*=2,++rank){const double*omega=roots[rank]->begin();for(unsigned block=0;block<n;block+=2*step){double*left=spectrum+2*block;const double*right=spectrum+2*(block+step);unsigned t=0;while(t<len){unsigned p=(start+t)&(base-1);if(!(p&3)&&t+4<=len){CV x(left+2*p),y(right+2*p),w(omega+2*p);(x+y.mul(w)).store(left+2*p);t+=4;}else{put(left,p,get(left,p)+get(right,p)*get(omega,p));++t;}}}}
return spectrum;
}
};
}
#pragma GCC push_options
#pragma GCC target("avx2,fma")
#pragma GCC optimize("O3")
namespace packed_ab_partial_inplace {
using packed_fft::PreparedA;using packed_fft::Result;using packed_fft::u32;using packed_fft::C;using packed_fft::FFT;using packed_fft::get;using packed_fft::put;using packed_fft::negslot;
class ABWorkspace{
unsigned h,n;uint64_t normB=0;hint::AlignMem<double>y,z;
std::map<std::pair<unsigned,unsigned>,std::unique_ptr<partial_fft_inplace::PartialInverse>> plans;
using CV=hint::Complex64X4;using FV=hint::Float64X4;
static CV conjugate(CV v){v.imag=FV(_mm256_xor_pd(__m256d(v.imag),_mm256_set1_pd(-0.)));return v;}
static CV product(CV X,CV Xc,CV Y,CV Yc){return (X+Xc).mul(Y+Yc)*FV(.25)+(X.mul(Y)-Xc.mul(Yc))*FV(.5);}
public:
uint64_t pruned=0,full=0;
explicit ABWorkspace(unsigned size):h(size),n(2*size),y(4*size),z(4*size){}
void set_b(const u32*b){normB=0;std::fill(y.begin(),y.begin()+2*n,0.);for(unsigned i=0;i<h;++i){int lo=int16_t(b[i]),hi=int16_t((b[i]+32768u)>>16);put(y.begin(),h-1-i,{double(lo),double(hi)});normB+=int64_t(lo)*lo+int64_t(hi)*hi;}FFT::difRecLong<false>(y.begin(),2*n);}
Result apply(const PreparedA&pa,u32*out,unsigned offset=0,unsigned count=0){if(!count)count=h;assert(pa.h==h&&offset+count<=h);uint64_t normA=pa.normA;long double normZ2=0;
for(unsigned i=0;i<4;++i){unsigned j=negslot(i);C X=get(pa.spectrum.begin(),i),Xc=std::conj(get(pa.spectrum.begin(),j)),Y=get(y.begin(),i),Yc=std::conj(get(y.begin(),j));C Z=.25*(X+Xc)*(Y+Yc)+.5*(X*Y-Xc*Yc);put(z.begin(),i,Z);normZ2+=static_cast<long double>(Z.real())*Z.real()+static_cast<long double>(Z.imag())*Z.imag();}
__m256d norms=_mm256_setzero_pd();
{unsigned i=4;CV X(pa.spectrum.begin()+2*i),Xc=conjugate(X.reverse()),Y(y.begin()+2*i),Yc=conjugate(Y.reverse());CV Z=product(X,Xc,Y,Yc);Z.store(z.begin()+2*i);norms=_mm256_fmadd_pd(__m256d(Z.real),__m256d(Z.real),norms);norms=_mm256_fmadd_pd(__m256d(Z.imag),__m256d(Z.imag),norms);}
for(unsigned block=8;block<n;block*=2)for(unsigned i=block;i<block+block/2;i+=4){unsigned j=3*block-4-i;CV X(pa.spectrum.begin()+2*i),Xc=conjugate(CV(pa.spectrum.begin()+2*j).reverse()),Y(y.begin()+2*i),Yc=conjugate(CV(y.begin()+2*j).reverse());CV T=(X+Xc).mul(Y+Yc)*FV(.25),S=(X.mul(Y)-Xc.mul(Yc))*FV(.5),Z=T+S,W=conjugate(T-S).reverse();Z.store(z.begin()+2*i);W.store(z.begin()+2*j);norms=_mm256_fmadd_pd(__m256d(Z.real),__m256d(Z.real),norms);norms=_mm256_fmadd_pd(__m256d(Z.imag),__m256d(Z.imag),norms);norms=_mm256_fmadd_pd(__m256d(W.real),__m256d(W.real),norms);norms=_mm256_fmadd_pd(__m256d(W.imag),__m256d(W.imag),norms);}
alignas(32)double sums[4];_mm256_store_pd(sums,norms);for(double value:sums)normZ2+=value;
constexpr long double u=0x1p-53L;long double op=(n<=(1u<<17)?5:64)*(unsigned(__builtin_ctz(n))+1),gamma=op*u/(1-op*u),normprod=sqrtl(static_cast<long double>(normA)*normB);long double normZ=sqrtl(normZ2)*(1+32*n*u);
long double error=(4*gamma+2*gamma*gamma+128*u*(1+gamma)*(1+gamma))*normprod+gamma*normZ/sqrtl(n);Result result;result.bound=error;if(!(error<.49L))return result;
if(_MM_GET_ROUNDING_MODE()!=_MM_ROUND_NEAREST)return result;
const double*inverse=z.begin();unsigned indexbase=h-1+offset,indexmask=~0u;bool usepartial=false;
if(count<=h/4){long double opI=(n<=(1u<<17)?7:64)*(unsigned(__builtin_ctz(n))+1),gammaI=opI*u/(1-opI*u);long double ep=(4*gamma+2*gamma*gamma+128*u*(1+gamma)*(1+gamma))*normprod+gammaI*normZ/sqrtl(n);if(ep<.49L){auto key=std::make_pair(offset,count);auto&plan=plans[key];if(!plan)plan=std::make_unique<partial_fft_inplace::PartialInverse>(n,h-1+offset,count);inverse=plan->apply(z.begin());indexmask=plan->block_size()-1;indexbase=(h-1+offset)&indexmask;result.bound=ep;usepartial=true;++pruned;}}
if(!usepartial){FFT::iditRecLong<false,false>(z.begin(),2*n);++full;}
const __m256d scale=_mm256_set1_pd(1./n),magic=_mm256_set1_pd(0x1.8p52);const __m256i mb=_mm256_castpd_si256(magic),perm=_mm256_setr_epi32(0,2,4,6,1,3,5,7);unsigned t=0;
auto scalar=[&](unsigned atout){C v=get(inverse,(indexbase+atout)&indexmask)/double(n);out[offset+atout]=u32(llround(v.real()))+(u32(llround(v.imag()))<<16);};
while(t<count&&((indexbase+t)&3))scalar(t++);
for(;t+4<=count;t+=4){unsigned pos=2*((indexbase+t)&indexmask);__m256d re=_mm256_mul_pd(_mm256_load_pd(inverse+pos),scale),im=_mm256_mul_pd(_mm256_load_pd(inverse+pos+4),scale);__m256i ri=_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(re,magic)),mb),ii=_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(im,magic)),mb);__m128i low=_mm256_castsi256_si128(_mm256_permutevar8x32_epi32(ri,perm)),cross=_mm256_castsi256_si128(_mm256_permutevar8x32_epi32(ii,perm));_mm_storeu_si128((__m128i*)(out+offset+t),_mm_add_epi32(low,_mm_slli_epi32(cross,16)));}
for(;t<count;++t)scalar(t);result.certified=true;return result;
}
};
}
#pragma GCC pop_options
#include <bits/stdc++.h>
#include <immintrin.h>
#pragma GCC target("avx2,fma")
#pragma GCC optimize("O3")
namespace sparse_native {
using u32=uint32_t;struct Group{u32 u,b;};
inline void add(unsigned n,const u32*a,unsigned dlo,const u32*D,unsigned H,u32*c,unsigned klo,unsigned khi){
std::vector<Group>groups;for(unsigned u=0;u<H&&dlo+u+klo<n;++u)if(D[u])groups.push_back({u,D[u]});if(groups.empty())return;
if(n<128){for(const auto&g:groups)for(unsigned k=klo;k<std::min(khi,n-dlo-g.u);++k)c[k]+=a[dlo+g.u+k]*g.b;return;}
for(unsigned k0=klo;k0<khi&&dlo+groups.front().u+k0<n;k0+=64){
__m256i z0=_mm256_loadu_si256((const __m256i*)(c+k0+0)), z1=_mm256_loadu_si256((const __m256i*)(c+k0+8)), z2=_mm256_loadu_si256((const __m256i*)(c+k0+16)), z3=_mm256_loadu_si256((const __m256i*)(c+k0+24)), z4=_mm256_loadu_si256((const __m256i*)(c+k0+32)), z5=_mm256_loadu_si256((const __m256i*)(c+k0+40)), z6=_mm256_loadu_si256((const __m256i*)(c+k0+48)), z7=_mm256_loadu_si256((const __m256i*)(c+k0+56));size_t stop=0;
if(k0+64<=khi&&dlo+k0+64<=n){unsigned maxu=n-dlo-k0-64;if(groups.back().u<=maxu)stop=groups.size();else stop=std::upper_bound(groups.begin(),groups.end(),maxu,[](u32 u,const Group&g){return u<g.u;})-groups.begin();}
for(size_t at=0;at<stop;++at){const auto&g=groups[at];const u32*x=a+dlo+g.u+k0;__m256i bj=_mm256_set1_epi32(g.b);
z0=_mm256_add_epi32(z0,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+0)),bj));
z1=_mm256_add_epi32(z1,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+8)),bj));
z2=_mm256_add_epi32(z2,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+16)),bj));
z3=_mm256_add_epi32(z3,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+24)),bj));
z4=_mm256_add_epi32(z4,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+32)),bj));
z5=_mm256_add_epi32(z5,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+40)),bj));
z6=_mm256_add_epi32(z6,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+48)),bj));
z7=_mm256_add_epi32(z7,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+56)),bj));
}
for(size_t at=stop;at<groups.size();++at){const auto&g=groups[at];if(dlo+g.u+k0>=n)break;unsigned cnt=std::min(khi-k0,n-dlo-g.u-k0);const u32*x=a+dlo+g.u+k0;__m256i bj=_mm256_set1_epi32(g.b);
#define MASKADD(V,Z) if(cnt>V*8){__m256i mask=_mm256_cmpgt_epi32(_mm256_set1_epi32(cnt),_mm256_setr_epi32(V*8,V*8+1,V*8+2,V*8+3,V*8+4,V*8+5,V*8+6,V*8+7));Z=_mm256_add_epi32(Z,_mm256_mullo_epi32(_mm256_maskload_epi32((const int*)(x+V*8),mask),bj));}
MASKADD(0,z0)
MASKADD(1,z1)
MASKADD(2,z2)
MASKADD(3,z3)
MASKADD(4,z4)
MASKADD(5,z5)
MASKADD(6,z6)
MASKADD(7,z7)
#undef MASKADD
}
_mm256_storeu_si256((__m256i*)(c+k0+0),z0);
_mm256_storeu_si256((__m256i*)(c+k0+8),z1);
_mm256_storeu_si256((__m256i*)(c+k0+16),z2);
_mm256_storeu_si256((__m256i*)(c+k0+24),z3);
_mm256_storeu_si256((__m256i*)(c+k0+32),z4);
_mm256_storeu_si256((__m256i*)(c+k0+40),z5);
_mm256_storeu_si256((__m256i*)(c+k0+48),z6);
_mm256_storeu_si256((__m256i*)(c+k0+56),z7);
}
}
inline uint64_t count(unsigned n,unsigned dlo,const u32*D,unsigned H,unsigned klo,unsigned khi){if(dlo+klo>=n)return 0;if(uint64_t(dlo)+H+khi<=n){unsigned nonzero=0,u=0;for(;u+8<=H;u+=8){unsigned zmask=_mm256_movemask_ps(_mm256_castsi256_ps(_mm256_cmpeq_epi32(_mm256_loadu_si256((const __m256i*)(D+u)),_mm256_setzero_si256())));nonzero+=8-__builtin_popcount(zmask);}for(;u<H;++u)nonzero+=bool(D[u]);return uint64_t(nonzero)*(khi-klo);}uint64_t value=0;unsigned use=std::min(H,n-dlo-klo);for(unsigned u=0;u<use;++u)value+=uint64_t(D[u]!=0)*(std::min(khi,n-dlo-u)-klo);return value;}
}
namespace tiled_mod {
// An output interval may span several H blocks. Tasks are still scheduled by
// their A window; rejected transforms share all smaller A windows in that same
// 2H region. No input-dependent spectrum survives the current A bucket.
class RecursiveStreamState {
struct Task{unsigned node,kb;};
const u32*a;unsigned n,H,root_start=0;
std::map<unsigned,std::unique_ptr<packed_ab_partial_inplace::ABWorkspace>> work;
std::map<unsigned,std::vector<u32>> output;
std::map<u64,std::unique_ptr<packed_fft::PreparedA>> local_a;
std::map<u64,u64> local_norm;
static u64 key(unsigned h,unsigned start){return (u64(h)<<32)|start;}
u64 norm(unsigned h,unsigned start){
auto id=key(h,start);auto found=local_norm.find(id);if(found!=local_norm.end())return found->second;
u64 value=0;for(unsigned i=start;i<std::min(n,start+2*h);++i){int lo=int16_t(a[i]),hi=int16_t((a[i]+32768u)>>16);value+=int64_t(lo)*lo+int64_t(hi)*hi;}
local_norm.emplace(id,value);return value;
}
const packed_fft::PreparedA&prepared_a(unsigned h,unsigned start){
assert(start>=root_start&&u64(start)+2*h<=u64(root_start)+2*H);
auto&slot=local_a[key(h,start)];if(!slot){
if(start+2*h<=n)slot=std::make_unique<packed_fft::PreparedA>(a+start,h);
else{std::vector<u32>pad(2*h);if(start<n)std::copy(a+start,a+n,pad.begin());slot=std::make_unique<packed_fft::PreparedA>(pad.data(),h);}
++prepared;
}return *slot;
}
void add_segment(unsigned dlo,const u32*db,unsigned h,u32*c,unsigned klo,unsigned khi){
if(klo>=khi||dlo+klo>=n)return;
// Positive-q MOD cross outputs are below n/2, so C64's preserved lanes
// remain inside the allocated C array. This helper is not a general range op.
assert(khi<=n/2);
u64 products=sparse_native::count(n,dlo,db,h,klo,khi);
if(h<=64||products<=u64(h)*250){++native_tasks;native_products+=products;sparse_native::add(n,a,dlo,db,h,c,klo,khi);return;}
u64 normB=0;for(unsigned i=0;i<h;++i){int lo=int16_t(db[i]),hi=int16_t((db[i]+32768u)>>16);normB+=int64_t(lo)*lo+int64_t(hi)*hi;}
constexpr long double u=0x1p-53L;unsigned N=2*h;
long double op=(N<=(1u<<17)?5:64)*(unsigned(__builtin_ctz(N))+1),gamma=op*u/(1-op*u);
long double coefficient=4*gamma+2*gamma*gamma+128*u*(1+gamma)*(1+gamma);
bool b_ready=false;
for(unsigned kb=klo/h;kb<=(khi-1)/h;++kb){
unsigned start=dlo+kb*h,kl=std::max(klo,kb*h),kh=std::min(khi,(kb+1)*h);if(start>=n)continue;
bool retry=coefficient*sqrtl(static_cast<long double>(norm(h,start))*normB)>=.49L;
if(retry)++pre_refused;
if(!retry){
auto&ws=work[h];if(!ws)ws=std::make_unique<packed_ab_partial_inplace::ABWorkspace>(h);
if(!b_ready){ws->set_b(db);b_ready=true;}
auto&tmp=output[h];tmp.resize(h);
auto result=ws->apply(prepared_a(h,start),tmp.data(),kl-kb*h,kh-kl);
if(result.certified){++accepted;for(unsigned k=kl;k<kh;++k)c[k]+=tmp[k-kb*h];continue;}
}
++refused;
add_segment(dlo,db,h/2,c,kl,kh);
add_segment(dlo+h/2,db+h/2,h/2,c,kl,kh);
}
}
public:
u64 accepted=0,refused=0,pre_refused=0,prepared=0,tasks=0,native_tasks=0,native_products=0;
RecursiveStreamState(unsigned size,const u32*data,unsigned h):a(data),n(size),H(h){assert(h>=8&&(h&(h-1))==0);}
void run(const u32*b,u32*c,const std::vector<Cross>&nodes){
std::vector<std::vector<Task>> buckets((n+H-1)/H+2);
for(unsigned id=0;id<nodes.size();++id){auto node=nodes[id];unsigned l=node.l,m=node.m,r=std::min(node.r,n);
if(l>=m||m>=r||m>=n-l)continue;unsigned wk=std::min(m-l,n-l-m),last=n-1-l;
std::vector<std::pair<unsigned,unsigned>> ranges;
for(unsigned q=1;q<=last/m;++q){unsigned lo=u64(q)*m/H,hi=std::min<u64>(u64(q)*(r-1),last)/H;
if(!ranges.empty()&&lo<=ranges.back().second+1)ranges.back().second=hi;else ranges.push_back({lo,hi});
}
for(auto range:ranges)for(unsigned block=range.first;block<=range.second;++block)
for(unsigned kb=l/H;kb<=(l+wk-1)/H;++kb)if(u64(block+kb)*H<n){buckets[block+kb].push_back({id,kb});++tasks;}
}
std::vector<u32>D(H);
for(unsigned ast=0;ast<buckets.size();++ast){if(buckets[ast].empty())continue;
root_start=ast*H;local_a.clear();local_norm.clear();
for(auto task:buckets[ast]){auto node=nodes[task.node];unsigned l=node.l,m=node.m,r=std::min(node.r,n),wk=std::min(m-l,n-l-m),dlo=(ast-task.kb)*H,dhi=std::min(n-1-l,dlo+H-1);
std::fill(D.begin(),D.end(),0);unsigned qlo=std::max(1u,(dlo+r-2)/(r-1)),qhi=dhi/m;bool any=false;if(qlo>qhi)continue;
if(qhi-qlo+1<r-m){
for(unsigned q=qlo;q<=qhi;++q){unsigned first=std::max(m,(dlo+q-1)/q),end=std::min(r-1,dhi/q);if(first<=end)any=true;for(unsigned j=first;j<=end;++j)D[q*j-dlo]+=b[j];}
}else{
for(unsigned j=m;j<r;++j){unsigned first=std::max(1u,(dlo+j-1)/j),end=dhi/j;if(first<=end)any=true;for(unsigned q=first;q<=end;++q)D[q*j-dlo]+=b[j];}
}
if(!any)continue;
unsigned kl=std::max(l,task.kb*H),kh=std::min(l+wk,(task.kb+1)*H);add_segment(dlo,D.data(),H,c,kl,kh);
}
}
local_a.clear();local_norm.clear();
}
};
}
#include <bits/stdc++.h>
#include <immintrin.h>
#pragma GCC target("avx2")
#pragma GCC optimize("O3")
namespace mod_filtered {
using u32=uint32_t;using u64=uint64_t;
struct Group{u32 d,b;};
template<unsigned DBlock=2048>
void cross_filtered(unsigned n,const u32*a,const u32*b,u32*c,unsigned l,unsigned r,unsigned qmin=1,unsigned qmax=~0u){
unsigned m=l+(r-l)/2,jend=std::min(r,n-l);qmin=std::max(1u,qmin);
if(m<=l||m>=jend||qmin>qmax)return;
if(n<128){for(unsigned j=m;j<jend;++j)for(unsigned q=qmin;q<=std::min(qmax,(n-1-l)/j);++q)for(unsigned k=l;k<std::min(m,n-q*j);++k)c[k]+=a[q*j+k]*b[j];return;}
const unsigned limit=std::min(m,n-m);std::vector<u32>weights(DBlock);std::vector<Group>groups;groups.reserve(DBlock);
for(unsigned lo=0;lo<n-l;lo+=DBlock){unsigned hi=std::min(n-l,lo+DBlock),qb=std::max(qmin,(lo+jend-2)/(jend-1)),qe=std::min(qmax,(hi-1)/m);if(qb>qe)continue;
std::fill(weights.begin(),weights.end(),0u);
if(qe-qb+1<jend-m){for(unsigned q=qb;q<=qe;++q){unsigned jb=std::max(m,(lo+q-1)/q),je=std::min(jend-1,(hi-1)/q);for(unsigned j=jb;j<=je;++j)weights[q*j-lo]+=b[j];}}
else {unsigned je=std::min(jend-1,(hi-1)/qmin);for(unsigned j=m;j<=je;++j){unsigned bq=std::max(qmin,(lo+j-1)/j),eq=std::min(qmax,(hi-1)/j);if(bq<=eq)for(unsigned d=bq*j;d<=eq*j;d+=j)weights[d-lo]+=b[j];}}
groups.clear();for(unsigned u=0;u<hi-lo;++u)if(weights[u])groups.push_back({lo+u,weights[u]});if(groups.empty())continue;
unsigned kend=std::min(limit,n-groups.front().d);
for(unsigned k0=l;k0<kend;k0+=64){
__m256i z0=_mm256_loadu_si256((const __m256i*)(c+k0+0)), z1=_mm256_loadu_si256((const __m256i*)(c+k0+8)), z2=_mm256_loadu_si256((const __m256i*)(c+k0+16)), z3=_mm256_loadu_si256((const __m256i*)(c+k0+24)), z4=_mm256_loadu_si256((const __m256i*)(c+k0+32)), z5=_mm256_loadu_si256((const __m256i*)(c+k0+40)), z6=_mm256_loadu_si256((const __m256i*)(c+k0+48)), z7=_mm256_loadu_si256((const __m256i*)(c+k0+56));
size_t stop=0;if(k0+64<=m){u32 maxd=n-k0-64;if(groups.back().d<=maxd)stop=groups.size();else{stop=std::upper_bound(groups.begin(),groups.end(),maxd,[](u32 d,const Group&g){return d<g.d;})-groups.begin();}}
for(size_t at=0;at<stop;++at){const Group&g=groups[at];const u32*x=a+g.d+k0;__m256i bj=_mm256_set1_epi32(g.b);
z0=_mm256_add_epi32(z0,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+0)),bj));
z1=_mm256_add_epi32(z1,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+8)),bj));
z2=_mm256_add_epi32(z2,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+16)),bj));
z3=_mm256_add_epi32(z3,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+24)),bj));
z4=_mm256_add_epi32(z4,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+32)),bj));
z5=_mm256_add_epi32(z5,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+40)),bj));
z6=_mm256_add_epi32(z6,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+48)),bj));
z7=_mm256_add_epi32(z7,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x+56)),bj));
}
for(size_t at=stop;at<groups.size();++at){const Group&g=groups[at];if(g.d+k0>=n)break;unsigned cnt=std::min(m-k0,n-g.d-k0);const u32*x=a+g.d+k0;__m256i bj=_mm256_set1_epi32(g.b);
#define MASKADD(V,Z) if(cnt>V*8){__m256i mask=_mm256_cmpgt_epi32(_mm256_set1_epi32(cnt),_mm256_setr_epi32(V*8,V*8+1,V*8+2,V*8+3,V*8+4,V*8+5,V*8+6,V*8+7));Z=_mm256_add_epi32(Z,_mm256_mullo_epi32(_mm256_maskload_epi32((const int*)(x+V*8),mask),bj));}
MASKADD(0,z0)
MASKADD(1,z1)
MASKADD(2,z2)
MASKADD(3,z3)
MASKADD(4,z4)
MASKADD(5,z5)
MASKADD(6,z6)
MASKADD(7,z7)
#undef MASKADD
}
_mm256_storeu_si256((__m256i*)(c+k0+0),z0);
_mm256_storeu_si256((__m256i*)(c+k0+8),z1);
_mm256_storeu_si256((__m256i*)(c+k0+16),z2);
_mm256_storeu_si256((__m256i*)(c+k0+24),z3);
_mm256_storeu_si256((__m256i*)(c+k0+32),z4);
_mm256_storeu_si256((__m256i*)(c+k0+40),z5);
_mm256_storeu_si256((__m256i*)(c+k0+48),z6);
_mm256_storeu_si256((__m256i*)(c+k0+56),z7);
}
}
}
}
#ifndef DUCK_SHA256_HPP
#define DUCK_SHA256_HPP
#include <stddef.h>
#include <stdint.h>
#include <string.h>
// Portable SHA-256; input is the exact byte sequence, with no native-word hashing.
namespace duck_sha256 {
static inline uint32_t rotate(uint32_t x,unsigned n) {return (x>>n)|(x<<(32-n));}
static void compress(uint32_t h[8],const unsigned char block[64]) {
static const uint32_t k[64]={
0x428a2f98u,0x71374491u,0xb5c0fbcfu,0xe9b5dba5u,0x3956c25bu,0x59f111f1u,0x923f82a4u,0xab1c5ed5u,
0xd807aa98u,0x12835b01u,0x243185beu,0x550c7dc3u,0x72be5d74u,0x80deb1feu,0x9bdc06a7u,0xc19bf174u,
0xe49b69c1u,0xefbe4786u,0x0fc19dc6u,0x240ca1ccu,0x2de92c6fu,0x4a7484aau,0x5cb0a9dcu,0x76f988dau,
0x983e5152u,0xa831c66du,0xb00327c8u,0xbf597fc7u,0xc6e00bf3u,0xd5a79147u,0x06ca6351u,0x14292967u,
0x27b70a85u,0x2e1b2138u,0x4d2c6dfcu,0x53380d13u,0x650a7354u,0x766a0abbu,0x81c2c92eu,0x92722c85u,
0xa2bfe8a1u,0xa81a664bu,0xc24b8b70u,0xc76c51a3u,0xd192e819u,0xd6990624u,0xf40e3585u,0x106aa070u,
0x19a4c116u,0x1e376c08u,0x2748774cu,0x34b0bcb5u,0x391c0cb3u,0x4ed8aa4au,0x5b9cca4fu,0x682e6ff3u,
0x748f82eeu,0x78a5636fu,0x84c87814u,0x8cc70208u,0x90befffau,0xa4506cebu,0xbef9a3f7u,0xc67178f2u};
uint32_t w[64];
for(unsigned i=0;i<16;++i)w[i]=(uint32_t(block[4*i])<<24)|(uint32_t(block[4*i+1])<<16)|(uint32_t(block[4*i+2])<<8)|block[4*i+3];
for(unsigned i=16;i<64;++i){uint32_t a=w[i-15],b=w[i-2];w[i]=w[i-16]+(rotate(a,7)^rotate(a,18)^(a>>3))+w[i-7]+(rotate(b,17)^rotate(b,19)^(b>>10));}
uint32_t a=h[0],b=h[1],c=h[2],d=h[3],e=h[4],f=h[5],g=h[6],v=h[7];
for(unsigned i=0;i<64;++i){
uint32_t t1=v+(rotate(e,6)^rotate(e,11)^rotate(e,25))+((e&f)^(~e&g))+k[i]+w[i];
uint32_t t2=(rotate(a,2)^rotate(a,13)^rotate(a,22))+((a&b)^(a&c)^(b&c));
v=g;g=f;f=e;e=d+t1;d=c;c=b;b=a;a=t1+t2;
}
h[0]+=a;h[1]+=b;h[2]+=c;h[3]+=d;h[4]+=e;h[5]+=f;h[6]+=g;h[7]+=v;
}
inline void digest(const void*data,size_t length,unsigned char out[32]) {
uint32_t h[8]={0x6a09e667u,0xbb67ae85u,0x3c6ef372u,0xa54ff53au,0x510e527fu,0x9b05688cu,0x1f83d9abu,0x5be0cd19u};
const unsigned char*p=(const unsigned char*)data;
size_t whole=length/64;
for(size_t i=0;i<whole;++i)compress(h,p+64*i);
unsigned char tail[128]={};size_t remain=length%64;
if(remain)memcpy(tail,p+64*whole,remain);
tail[remain]=0x80;size_t padded=remain<56?64:128;
uint64_t bits=uint64_t(length)*8;
for(unsigned i=0;i<8;++i)tail[padded-1-i]=(unsigned char)(bits>>(8*i));
compress(h,tail);if(padded==128)compress(h,tail+64);
for(unsigned i=0;i<8;++i)for(unsigned j=0;j<4;++j)out[4*i+j]=(unsigned char)(h[i]>>(24-8*j));
}
}
#endif
struct PublicRng{uint64_t state;unsigned next(){uint64_t z=(state+=0x9e3779b97f4a7c15ULL);z=(z^(z>>30))*0xbf58476d1ce4e5b9ULL;z=(z^(z>>27))*0x94d049bb133111ebULL;return unsigned(z^(z>>31));}};
static uint64_t ticks(){_mm_lfence();return __rdtsc();}
int main(){constexpr unsigned n=10000000,H=BENCH_H;std::vector<unsigned>a(n),b(n),c(n),gold(n);PublicRng ar{0x1357214865430012ULL},br{0x2143658730024589ULL};uint64_t t0=ticks();for(unsigned i=0;i<n;++i){a[i]=ar.next();b[i]=br.next();}
std::vector<tiled_mod::Cross>nodes;
if(BENCH_CASE==0)nodes={{0,8192,16384},{196608,212992,229376}};
else nodes={{98304,131072,163840},{1000000,1001024,1002048}};
uint64_t t1=ticks();tiled_mod::RecursiveStreamState state(n,a.data(),H);state.run(b.data(),c.data(),nodes);uint64_t t2=ticks();
for(auto node:nodes)mod_filtered::cross_filtered(n,a.data(),b.data(),gold.data(),node.l,node.r);
uint64_t t3=ticks();if(c!=gold)__builtin_trap();unsigned char hash[32];duck_sha256::digest(c.data(),size_t(n)*4,hash);uint64_t t4=ticks();
for(unsigned char x:hash)printf("%02x",x);puts("");
printf("G%llu F%llu N%llu V%llu ok1 a%llu r%llu p%llu t%llu nt%llu\n",(unsigned long long)(t1-t0),(unsigned long long)(t2-t1),(unsigned long long)(t3-t2),(unsigned long long)(t4-t3),(unsigned long long)state.accepted,(unsigned long long)state.refused,(unsigned long long)state.prepared,(unsigned long long)state.tasks,(unsigned long long)state.native_tasks);
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 5.415 s | 156 MB + 768 KB | Accepted | Score: 100 | 显示更多 |