// DUCK_PUBLIC_SURROGATE_BENCHMARK
// Independent public arrays only; no production seeds or data.
#include <bits/stdc++.h>
#include <x86intrin.h>
#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
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;}
}
#include <memory>
#include <vector>
#include <cstring>
#include <algorithm>
namespace tiled_ab_sparse {
// Charge a full 64-output SIMD tile for short tails when choosing a backend.
inline uint64_t native_score(uint64_t logical,unsigned width){return width&&width<64?(logical*64+width-1)/width:logical;}
using u32=uint32_t;using u64=uint64_t;
struct State{
const u32*a;unsigned n,H;packed_ab_partial_inplace::ABWorkspace work;
std::vector<std::unique_ptr<packed_fft::PreparedA>> cache;
std::vector<u32>input,D,tmp;std::vector<uint64_t>norm_blocks;
u64 accepted=0,refused=0,prepared=0,tiles=0,fallback_accepted=0,fallback_refused=0,fallback_native=0,pre_refused=0,native_tiles=0,native_products=0;uint64_t native_limit;long double max_bound=0;
State(unsigned size,const u32*data,unsigned h,uint64_t limit=500000):a(data),n(size),H(h),work(h),cache((size+h-1)/h+2),input(2*h),D(h),tmp(h),norm_blocks(cache.size()+1,UINT64_MAX),native_limit(limit){}
uint64_t normblock(unsigned block){auto&value=norm_blocks[block];if(value==UINT64_MAX){value=0;unsigned start=block*H;if(start<n)for(unsigned i=start;i<std::min(n,start+H);++i){int lo=int16_t(a[i]),hi=int16_t((a[i]+32768u)>>16);value+=int64_t(lo)*lo+int64_t(hi)*hi;}}return value;}
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;
}
void fallback(unsigned dlo,const u32*db,unsigned hs,u32*c,unsigned klo,unsigned khi){
if(dlo+klo>=n)return;
uint64_t products=sparse_native::count(n,dlo,db,hs,klo,khi);if(native_score(products,khi-klo)<=native_limit){fallback_native+=products;sparse_native::add(n,a,dlo,db,hs,c,klo,khi);return;}
if(hs<=64){for(unsigned u=0;u<hs&&dlo+u+klo<n;++u){u32 value=db[u];if(!value)continue;unsigned end=std::min(khi,n-dlo-u);fallback_native+=end-klo;for(unsigned k=klo;k<end;++k)c[k]+=a[dlo+u+k]*value;}return;}
packed_fft::MiddlePlan plan(db,hs);std::vector<u32>subinput(2*hs),subout(hs);
for(unsigned kb=klo/hs;kb<=(khi-1)/hs;++kb){unsigned ast=dlo+kb*hs;if(ast>=n)continue;unsigned kl=std::max(klo,kb*hs),kh=std::min(khi,(kb+1)*hs),count=std::min(2*hs,n-ast);std::fill(subinput.begin(),subinput.end(),0u);std::copy(a+ast,a+ast+count,subinput.begin());auto result=plan.apply(subinput.data(),subout.data());
if(result.certified){++fallback_accepted;for(unsigned k=kl;k<kh;++k)c[k]+=subout[k-kb*hs];}
else{++fallback_refused;fallback(dlo,db,hs/2,c,kl,kh);fallback(dlo+hs/2,db+hs/2,hs/2,c,kl,kh);}
}
}
// 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;uint64_t products=sparse_native::count(n,dlo,D.data(),H,l,l+wk);if(native_score(products,wk)<=native_limit){++native_tiles;native_products+=products;sparse_native::add(n,a,dlo,D.data(),H,c,l,l+wk);continue;}bool b_ready=false;uint64_t normB=0;for(u32 value:D){int lo=int16_t(value),hi=int16_t((value+32768u)>>16);normB+=int64_t(lo)*lo+int64_t(hi)*hi;}
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;
constexpr long double unit=0x1p-53L;unsigned nn=2*H;long double op=(nn<=(1u<<17)?5:64)*(unsigned(__builtin_ctz(nn))+1),gamma=op*unit/(1-op*unit);uint64_t normA=normblock(block+kb)+normblock(block+kb+1);long double first_bound=(4*gamma+2*gamma*gamma+128*unit*(1+gamma)*(1+gamma))*sqrtl(static_cast<long double>(normA)*normB);
if(first_bound>=.49L){++refused;++pre_refused;max_bound=std::max(max_bound,first_bound);fallback(dlo,D.data(),H/2,c,klo,khi);fallback(dlo+H/2,D.data()+H/2,H/2,c,klo,khi);continue;}
if(!b_ready){work.set_b(D.data());b_ready=true;}auto result=work.apply(get(block+kb),tmp.data(),offset,khi-klo);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;fallback(dlo,D.data(),H/2,c,klo,khi);fallback(dlo+H/2,D.data()+H/2,H/2,c,klo,khi);}
}
}
}
};
}
#pragma GCC push_options
#pragma GCC optimize("fp-contract=off")
#include <bits/stdc++.h>
namespace mod_plan_count {using u64=uint64_t;
inline u64 direct_count(unsigned n,unsigned l,unsigned r){u64 total=0,A=n-1-l,end=std::min<u64>(r-1,A);for(u64 j=l+1;j<=end;){u64 q=A/j,hi=std::min(end,A/q),cut=std::min(hi,u64(n)/(q+1));if(cut>=j){u64 count=cut-j+1;total+=q*((j+cut)*count/2-u64(l)*count);}u64 lo=std::max(j,cut+1);if(lo<=hi){u64 count=hi-lo+1;total+=count*(n-q*l)-(lo+hi)*count/2;}j=hi+1;}return total;}}
// Cost-model research only. No unmeasured model duration is acceptance evidence.
namespace calibrated_plan {
using u64=uint64_t;
enum Method{GROUPED=0,MERGED=1,TILED=2};
struct Cross{unsigned l,r;Method method;unsigned H=0;};struct Range{unsigned l,r;};
inline u64 term_count(unsigned n,unsigned k,unsigned first,unsigned last){u64 total=0,A=n-1-k,lastj=std::min<u64>(last,A);for(u64 j=first;j<=lastj;){u64 q=A/j,hi=std::min(lastj,A/q);total+=q*(hi-j+1);j=hi+1;}return total;}
// Upper bound on nonempty D tiles, exact whenever consecutive q*j differ <=H.
// Intervals from successive q are monotonically ordered, so union is linear in q.
inline u64 tile_cover(unsigned n,unsigned l,unsigned r,unsigned H){unsigned w=(r-l)/2,m=l+w,jlast=std::min(r-1,n-l-1);if(!w||m>jlast)return 0;u64 count=0,end=0;bool any=false;unsigned qmax=(n-1-l)/m;
for(unsigned q=1;q<=qmax;++q){u64 lo=u64(q)*m/H,hi=std::min<u64>(n-1-l,u64(q)*jlast)/H;if(!any||lo>end+1){count+=hi-lo+1;end=hi;any=true;}else if(hi>end){count+=hi-end;end=hi;}}
return count;
}
struct Costs{
double native_pair=1.157129381204629/12e9,native_record=0,group_unit=2.77e-9*1.0469382184428717,merge_unit=1.76e-9*1.063637380615047;
// Must be supplied from the actual tiled backend before interpreting seconds.
double tile_forward=30e-6,tile_inverse=49e-6,tile_record=0;
unsigned H=4096;
};
inline double tile_seconds(unsigned n,unsigned l,unsigned r,unsigned H,const Costs&cost){
unsigned w=(r-l)/2,m=l+w,wk=std::min(w,n-l-m),outputs=(l%H+wk+H-1)/H;double ratio=double(H)*__builtin_ctz(H)/(4096.*12),partial=1.;if(wk<=H/4){unsigned power=1;while(power<wk)power*=2;partial=.25+.75*std::max(4u,unsigned(__builtin_ctz(power)))/double(__builtin_ctz(2*H));}
double penalty=1.;if(H==16384){double density=double(w)/r;double refusal=density>=.38?1.:std::min(1.,3.*w/(n-l));double fallback=wk>H/2?2.2:1.28;penalty=1.+refusal*(fallback-1.);}
const double measured_scale[3]={1.3771538032523578,1.2382085846342463,1.4349022790994899};return tile_cover(n,l,r,H)*(30e-6*ratio+outputs*49e-6*ratio*partial)*penalty*measured_scale[__builtin_ctz(H)-12];
}
struct Plan{
unsigned n;Costs cost;std::vector<Range>direct;std::vector<Cross>cross;double estimate=0;
double build(unsigned l,unsigned r){if(r-l<2||2ull*l+1>=n)return 0;u64 native=mod_plan_count::direct_count(n,l,r);unsigned width=r-l,qrow=(n-1-l)/(l+1);double factor=(l==0&&width>=1024)||(l>=65536&&width>=4096&&uint64_t(width)>=128ull*qrow)?.68:1.;double nc=factor*native*cost.native_pair+cost.native_record*term_count(n,l,l+1,r-1);if(r-l<=64){direct.push_back({l,r});return nc;}
size_t ds=direct.size(),cs=cross.size();unsigned w=(r-l)/2,m=l+w;double child=build(l,m)+build(m,r);unsigned qmax=(n-1-l)/m;if(!qmax)return child;
double gu=0;for(unsigned lo=1;lo<=qmax;){unsigned len=8;while(len<u64(lo)*(w-1)+w)len*=2;unsigned hi=std::min(qmax,(len-w)/(w-1));gu+=double(hi-lo+3)/3.*len*__builtin_ctz(len);lo=hi+1;}
unsigned span=n-l-m,wk=std::min(w,span),len=8;while(len<u64(span)+wk-1)len*=2;double mg=double(len)*__builtin_ctz(len)*cost.merge_unit,gp=gu*cost.group_unit;
Method kind=gp<mg?GROUPED:MERGED;double xc=std::min(mg,gp);unsigned chosenH=0;
for(unsigned H:{4096u,8192u,16384u}){double tc=tile_seconds(n,l,r,H,cost);if(tc<xc){xc=tc;kind=TILED;chosenH=H;}}
if(nc<child+xc){direct.resize(ds);cross.resize(cs);direct.push_back({l,r});return nc;}
cross.push_back({l,r,kind,chosenH});return child+xc;
}
Plan(unsigned size,Costs model):n(size),cost(model){unsigned r=1;while(r<n)r*=2;estimate=build(0,r);}
};
// Uniform, seed-independent positive-integer coefficient bound for the subset
// of pairs assigned to the prime fields. TILED/native pairs are already mod32.
inline u64 field_pair_bound(unsigned n,const std::vector<Cross>&nodes){
std::vector<std::pair<unsigned,int64_t>>events;for(auto node:nodes){if(node.method==TILED)continue;unsigned m=node.l+(node.r-node.l)/2,end=std::min(m,n-m);if(end<=node.l)continue;u64 count=term_count(n,node.l,m,std::min(node.r-1,n-node.l-1));if(!count)continue;events.push_back({node.l,int64_t(count)});events.push_back({end,-int64_t(count)});}
std::sort(events.begin(),events.end());u64 best=0;int64_t value=0;for(size_t i=0;i<events.size();){unsigned at=events[i].first;while(i<events.size()&&events[i].first==at)value+=events[i++].second;assert(value>=0);best=std::max(best,u64(value));}return best;
}
// Exact maximum of the selected-field pair counts. Between any two start/end
// events the active node set is fixed and every floor((n-1-k)/j) decreases with
// k, hence an interval maximum occurs at its first event coordinate.
inline u64 field_pair_maximum(unsigned n,const std::vector<Cross>&nodes){
struct E{unsigned at,id;bool start;};std::vector<E>events;for(unsigned id=0;id<nodes.size();++id){auto node=nodes[id];if(node.method==TILED)continue;unsigned m=node.l+(node.r-node.l)/2,end=std::min(m,n-m);if(end<=node.l)continue;events.push_back({node.l,id,true});events.push_back({end,id,false});}
std::sort(events.begin(),events.end(),[](E a,E b){return a.at<b.at;});std::set<unsigned>active;u64 best=0;for(size_t i=0;i<events.size();){unsigned at=events[i].at;while(i<events.size()&&events[i].at==at){auto e=events[i++];if(e.start)active.insert(e.id);else active.erase(e.id);}u64 total=0;for(unsigned id:active){auto node=nodes[id];unsigned m=node.l+(node.r-node.l)/2;total+=term_count(n,at,m,node.r-1);}best=std::max(best,total);}return best;
}
}
#pragma GCC pop_options
#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
static uint64_t ticks(){_mm_lfence();return __rdtsc();}
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 void emit_sha(const unsigned*c,unsigned n){unsigned char h[32];duck_sha256::digest(c,size_t(n)*4,h);for(unsigned char x:h)printf("%02x",x);puts("");}
int main(){constexpr unsigned n=10000000,H=4096;std::vector<unsigned>a(n+64),b(n+64),c(n+64);auto t0=ticks();PublicRng ar{0x1357214865430012ULL},br{0x2143658730024589ULL};for(unsigned i=0;i<n;++i){a[i]=ar.next();b[i]=br.next();}auto t1=ticks();calibrated_plan::Plan plan(n,{});std::vector<calibrated_plan::Cross>nodes;for(auto r:plan.cross)if(r.method==calibrated_plan::TILED&&r.H==H)nodes.push_back(r);std::sort(nodes.begin(),nodes.end(),[](auto a,auto b){if(a.l!=b.l)return a.l<b.l;return a.r>b.r;});auto t2=ticks();tiled_ab_sparse::State state(n,a.data(),H,uint64_t(250)*H);for(unsigned at=0;at<(n+H-1)/H;++at)state.get(at);auto t3=ticks();for(auto r:nodes)state.cross(b.data(),c.data(),r.l,r.l+(r.r-r.l)/2,r.r);auto t4=ticks();emit_sha(c.data(),n);auto t5=ticks();printf("G%llu P%llu A%llu R%llu H%llu nodes%zu rf%llu sa%llu sr%llu nn%llu gt%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)(t5-t4),nodes.size(),(unsigned long long)state.refused,(unsigned long long)state.fallback_accepted,(unsigned long long)state.fallback_refused,(unsigned long long)state.fallback_native,(unsigned long long)state.native_tiles);}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 1.638 s | 423 MB + 644 KB | Accepted | Score: 100 | 显示更多 |