// ===== 引用区（本账号 saffah_v41f1017_agg4）=====
// [1] duck.ac 用户 saffah_cc_v41_agg1，提交 #130085 <https://duck.ac/submission/130085>
// （1004e7 实时 T=57.424 ms 榜首）。本文件正文 = 该公开提交的逐字节副本 + 下方所述唯一改动；
// 其自带的完整上游引用链（saffah_codex_6a_agg3 #128192/#128168/#128146/#128142、
// saffah_cc_v41_agg1 #128350/#127228/#127097 等账号与 URL、内容说明）原样保留在下方。
// 该公开提交无独立许可证声明；按站点公开条款具名引用，未声称为原创。
// [2] duck.ac/problem/1004e7 题面：a、b 各 10,000,000 位，传统题 stdin→stdout。
// 合规：以上均为 duck.ac 公开提交；正文中以注释形式保留其全部上游引用。
// ===== 思路 =====
// 单变量试验：把 fwd4 顶层分片的两条软预取距离 +8192 改为 +16384（两处同步）。
// 目的：在真实判题机（i3-8100, GCC 9.3）上测量该参数；其余逐字节未动，输出应与底盘逐字节同。
// ======================
#include <complex>
#include <iostream>
#include <type_traits>
#include <cstdint>
#include <climits>
#include <cstring>
#include <immintrin.h>
#pragma GCC target("avx2,fma,tune=skylake")
#pragma GCC optimize("O3,no-unroll-loops","rename-registers","align-loops=32","align-functions=32")
#pragma GCC optimize("prefetch-loop-arrays")
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));}
class Float64X4{
public:
using F64 = double;
using F64X4 = Float64X4;
Float64X4() : data(_mm256_setzero_pd()) {}
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(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);
auto real_i64=_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(real,magic)),bits);
auto imag_i64=_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(imag,magic)),bits);
return C64X4(__m256d(real_i64),__m256d(imag_i64));}
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;
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){
const Float theta = -HINT_2PI * factor / theta_divider;
for (size_t begin = 0, index = 0; begin < OMEGA_LEN * 2; begin += stride * 2)
{
for (size_t j = 0; j < stride; j++, index++)
{
table[begin + j] = std::cos(theta * index);
table[begin + j + stride] = std::sin(theta * index);}}}
constexpr const Float &operator[](size_t index) const{
return table[index];}};
void initOmegaX4(F64 *arr, size_t fft_len, int table_len, int factor){
table_len /= 4;
const F64 theta = -HINT_2PI * factor / fft_len;
auto arrx4 = reinterpret_cast<C64X4 *>(arr);
arr[0] = 1, arr[4] = 0;
arr[1] = std::cos(theta), arr[5] = std::sin(theta);
arr[2] = std::cos(theta * 2), arr[6] = std::sin(theta * 2);
arr[3] = std::cos(theta * 3), arr[7] = std::sin(theta * 3);
for (size_t begin = 1; begin < table_len; begin *= 2){
size_t nth = begin * 4;
C64X4 unit;
unit.set1(std::cos(theta * nth), std::sin(theta * nth));
for (size_t i = 0; i < begin; i++)
{
arrx4[begin + i] = arrx4[i].mul(unit);}}}
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){
static_assert(std::is_same<Float, Float64>::value);
size_t len = size_t(1) << LOG_BEGIN, cpx_len = len / DIV;
auto it = getBeginLog(LOG_BEGIN);
initOmegaX4(it, len, cpx_len, factor);
for (int log_len = LOG_BEGIN + 1; log_len <= LOG_END; log_len++)
{
len = size_t(1) << log_len, cpx_len = len / DIV;
Float theta = -HINT_2PI * factor / len;
auto it = getBeginLog(log_len), it_last = getBeginLog(log_len - 1);
C64X4 unit(std::cos(theta), std::sin(theta));
for (auto end = it + cpx_len * 2; it < end; it += 16, it_last += 8)
{
C64X4 omega0, omega1;
omega0.load(it_last);
omega1 = omega0.mul(unit);
transpose64_2X4(omega0.real, omega1.real);
transpose64_2X4(omega0.imag, omega1.imag);
omega0.store(it), omega1.store(it + 8);}}}
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;
size_t low_len = CACHE_LEN * C4_COUNT, high_len = table_len / low_len;
if (high != nullptr)
{
delete[] high;}
high = new C64[high_len];
auto p = reinterpret_cast<F64 *>(&low[0]);
initOmegaX4(p, fft_len, low_len, factor);
const F64 theta = -HINT_2PI * factor / fft_len;
high[0] = C64(1, 0);
for (size_t begin = 1; begin < high_len; begin *= 2)
{
C64 unit = std::polar<F64>(1.0, theta * begin * low_len);
for (size_t i = 0; i < begin; i++)
{
high[i + begin] = high[i] * unit;}}}
C64X4 operator[](size_t i) const{
C64X4 hi;
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 dotjit { extern int on; extern Float64 *a, *b; extern size_t child; extern int cid; void leaf(Float64* p, size_t n); void set_child(int r);
}
struct FFTAVX : public FFT{
static constexpr int LOG_SHORT = 9, LOG_MID = 17, LOG_MAX = 0x16, 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;
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);}}}
static void difIter(Float64 in_out[], size_t float_len){
size_t fft_len = float_len / 2;
C64X4 c0, c1, c2, c3;
size_t stride = fft_len / 2;
auto it0 = in_out, it1 = it0 + stride, it2 = it1 + stride, it3 = it2 + stride;
for (size_t rank = 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){
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 = omega1.mul(omega1), c1 = c1.mul(omega2); \
c3 = c3.mul(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+192),_MM_HINT_T0); \
_mm_prefetch((const char*)(it1+256),_MM_HINT_T0); \
_mm_prefetch((const char*)(it2+320),_MM_HINT_T0); \
_mm_prefetch((const char*)(it3+384),_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 = omega1.mul(omega1), c1 = c1.mul(omega2); \
c3 = c3.mul(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+192),_MM_HINT_T0); \
_mm_prefetch((const char*)(it1+256),_MM_HINT_T0); \
_mm_prefetch((const char*)(it2+320),_MM_HINT_T0); \
_mm_prefetch((const char*)(it3+384),_MM_HINT_T0); \
C64X4 c0 = it0, c1 = it1, c2 = it2, c3 = it3, omega1, omega2; \
omega1 = table[indx], c2 = c2.mulConj(omega1); \
omega2 = omega1.mul(omega1), c1 = c1.mulConj(omega2); \
c3 = c3.mulConj(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 (dotjit::on) dotjit::leaf(in_out, float_len);
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);
}
};
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);
butterfly(v[4*q],v[4*q+1],v[4*q+2],v[4*q+3],t1[j/8]);
}
#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++){
butterfly(v[4*q],v[4*q+1],v[4*q+2],v[4*q+3],t1[j/8]);
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);
}
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;}
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::TableSqrt FFTAVX::sqrt_table_1(1);
constexpr uint32_t bitrev32(uint32_t n){
constexpr uint32_t mask55 = 0x55555555;
constexpr uint32_t mask33 = 0x33333333;
constexpr uint32_t mask0f = 0x0f0f0f0f;
constexpr uint32_t maskff = 0x00ff00ff;
n = ((n & mask55) << 1) | ((n >> 1) & mask55);
n = ((n & mask33) << 2) | ((n >> 2) & mask33);
n = ((n & mask0f) << 4) | ((n >> 4) & mask0f);
n = ((n & maskff) << 8) | ((n >> 8) & maskff);
return (n << 16) | (n >> 16);}
constexpr uint32_t bitrev(uint32_t n, int len){
return bitrev32(n) >> (32 - len);}
template <int MAX_LOG_LEN, int DIV>
class BinRevTableC64X4HP{
public:
static constexpr int LOG_BLOCK = 2, BLOCK = 1 << LOG_BLOCK;
static constexpr size_t MAX_LEN = size_t(1) << MAX_LOG_LEN;
struct Unit{
C64 units[MAX_LOG_LEN]{};
F64 block[BLOCK * 2]{};
Unit()
{
constexpr F64 factor = F64(1) / DIV;
for (int i = 0; i < MAX_LOG_LEN; i++)
{
units[i] = getOmega(size_t(1) << (i + 1), 1, factor);}
block[0] = 1, block[BLOCK] = 0;
for (int i = 1; i < BLOCK; i++)
{
C64 omega = getOmega(BLOCK, bitrev(i, LOG_BLOCK), factor);
block[i] = omega.real(), block[i + BLOCK] = omega.imag();}}};
BinRevTableC64X4HP() : index(0), pop(0){
std::memcpy(table, unit_table.block, sizeof(unit_table.block));}
void reset(size_t i = 0){
if (i == 0)
{
pop = 0, index = i;
return;}
pop = 1, index = i / BLOCK;
int zero = __builtin_ctzll(index);
auto fp = reinterpret_cast<const F64 *>(&unit_table.units[zero + 2]);
table[1].load1(fp, fp + 1);
table[1] = table[1].mul(table[0]);}
C64X4 iterate(){
C64X4 res = table[pop], unit4;
index++;
int zero = __builtin_ctzll(index);
auto fp = reinterpret_cast<const F64 *>(&unit_table.units[zero + 2]);
unit4.load1(fp, fp + 1);
pop -= zero;
table[pop + 1] = table[pop].mul(unit4);
pop++;
return res;}
static C64 getOmega(size_t n, size_t index, F64 factor = 1){
F64 theta = -HINT_2PI * index / n;
return std::polar<F64>(1, theta * factor);}
private:
alignas(64) static const Unit unit_table;
alignas(64) C64X4 table[MAX_LOG_LEN];
size_t index;
int pop;
int log_max_iter, log_fft_len;};
template <int MAX_LOG_LEN, int DIV>
alignas(64) const typename BinRevTableC64X4HP<MAX_LOG_LEN, DIV>::Unit BinRevTableC64X4HP<MAX_LOG_LEN, DIV>::unit_table;
template <size_t RI_DIFF = 1, typename FloatTy>
inline void dot_rfft(FloatTy *inout0, FloatTy *inout1, const FloatTy *in0, const FloatTy *in1,
const std::complex<FloatTy> &omega, const FloatTy inv = 1){
using Complex = std::complex<FloatTy>;
auto addConj = [](Complex c0, Complex c1){ return Complex(c0.real() + c1.real(), c0.imag() - c1.imag()); };
Complex x0(inout0[0], inout0[RI_DIFF]), x1(inout1[0], inout1[RI_DIFF]),
y0(in0[0], in0[RI_DIFF]), y1(in1[0], in1[RI_DIFF]);
auto t0 = x0 * y0, t1 = x1 * y1, xy0 = addConj(x0, x1), xy1 = addConj(y0, y1);
auto t2 = xy0 * xy1;
y1 = addConj(t0, t1);
x1 = (y1 + y1 - t2) * omega * omega;
const auto inv2 = inv + inv;
x0 = (t2 - x1) * inv, x1 = Complex(t0.real() - t1.real(), t0.imag() + t1.imag()) * inv2;
Complex out0 = x0 + x1, out1(x0.real() - x1.real(), x1.imag() - x0.imag());
inout0[0] = out0.real(), inout0[RI_DIFF] = out0.imag();
inout1[0] = out1.real(), inout1[RI_DIFF] = out1.imag();}
inline void dot_rfftX4(F64 *inout0, F64 *inout1, const F64 *in0, const F64 *in1, const C64X4 &omega, const F64X4 &inv){
auto addConj = [](const C64X4 &x0, const C64X4 &x1){
return C64X4(x0.real + x1.real, x0.imag - x1.imag);};
C64X4 x0 = inout0, x1 = inout1, y0 = in0, y1 = in1;
x1 = x1.reverse();
y1 = y1.reverse();
C64X4 t0 = x0.mul(y0), t1 = x1.mul(y1);
C64X4 xy0 = addConj(x0, x1), xy1 = addConj(y0, y1);
C64X4 t2 = xy0.mul(xy1);
y1 = addConj(t0, t1);
x1 = (y1 + y1 - t2).mul(omega);
const F64X4 inv2 = inv + inv;
x0 = (t2 - x1) * inv, x1 = C64X4(t0.real - t1.real, t0.imag + t1.imag) * inv2;
C64X4 out0 = x0 + x1, out1(x0.real - x1.real, x1.imag - x0.imag);
out0.store(inout0), out1.reverse().store(inout1);}
inline void real_dot_binrev4(Float64 in_out[], Float64 in[], size_t float_len,Float64 scale=1){
Float64 inv = 2.0 * scale / float_len;{
auto r0 = in_out[0], i0 = in_out[4], r1 = in[0], i1 = in[4];
transform2(r0, i0);
transform2(r1, i1);
r0 *= r1, i0 *= i1;
transform2(r0, i0);
in_out[0] = r0 * 0.5 * inv, in_out[4] = i0 * 0.5 * inv;}
auto temp = C64(in_out[1], in_out[5]) * C64(in[1], in[5]) * inv;
in_out[1] = temp.real(), in_out[5] = temp.imag();
inv /= 4;
dot_rfft<4>(&in_out[2], &in_out[3], &in[2], &in[3], C64(COS_PI_8, -COS_PI_8), inv);
constexpr Float64 COS_16_1 = 0.92387953251128675612818318939;
constexpr Float64 SIN_16_1 = 0.38268343236508977172845998403;
dot_rfft<4>(&in_out[8], &in_out[11], &in[8], &in[11], C64(COS_16_1, -SIN_16_1), inv);
dot_rfft<4>(&in_out[9], &in_out[10], &in[9], &in[10], C64(-SIN_16_1, -COS_16_1), inv);
const Float64X4 inv4 = F64X4(0.5 * scale / float_len);
BinRevTableC64X4HP<28, 1> table;
for (size_t begin = 16; begin < float_len; begin *= 2){
table.reset(begin / 2);
auto it0 = in_out + begin, it1 = it0 + begin - 8, it2 = in + begin, it3 = it2 + begin - 8;
for (; it0 < it1; it0 += 8, it1 -= 8, it2 += 8, it3 -= 8)
{
_mm_prefetch((const char*)(it0+256),_MM_HINT_T0);
_mm_prefetch((const char*)(it1-256),_MM_HINT_T0);
_mm_prefetch((const char*)(it2+256),_MM_HINT_T0);
_mm_prefetch((const char*)(it3-256),_MM_HINT_T0);
dot_rfftX4(it0, it1, it2, it3, table.iterate(), inv4);}}}
template <bool TO_INT = false>
inline void real_conv_avx(F64 *in_out1, F64 *in2, size_t float_len){
FFTAVX::difRecLong<true>(in_out1, float_len);
FFTAVX::difRecLong<true>(in2, float_len);
real_dot_binrev4(in_out1, in2, float_len);
FFTAVX::iditRecLong<true, TO_INT>(in_out1, float_len);}}}}
#include <cstdio>
#include <string>
#include <memory>
#include <sys/stat.h>
#include <unistd.h>
#include <sys/auxv.h>
#include <cstddef>
namespace duck_public_stdio {
constexpr unsigned long AT_DUCK = 0x6b637564;
struct DuckInfo_t {
uint64_t abi_version;
const char *stdin_ptr;
uint64_t stdin_size;
char *stdout_ptr;
uint64_t stdout_limit;
uint64_t stdout_size;
} __attribute__((packed));
static_assert(sizeof(void*)==8 && sizeof(DuckInfo_t)==48,"64-bit Duck ABI required");
static_assert(offsetof(DuckInfo_t,stdin_ptr)==8 && offsetof(DuckInfo_t,stdin_size)==16,
"Public stdin fields must match ABI 0.04");
static_assert(offsetof(DuckInfo_t,stdout_ptr)==24 && offsetof(DuckInfo_t,stdout_size)==40,
"Public stdout fields must match ABI 0.04");
}
#include <cstddef>
#include <cstdint>
#include <cstring>
namespace bigint_io_opt {
inline void parse_balanced_base100000_dual(const char* s0,std::size_t len0,double* out0,
const char* s1,std::size_t len1,double* out1) {
size_t k0=0,k1=0; int carry0=0,carry1=0;
__m128i prev0=_mm_setzero_si128(),prev1=_mm_setzero_si128();
const __m128i zeros=_mm_set1_epi8('0'),w1=_mm_set1_epi16(0x010a),w2=_mm_set1_epi32(0x00010064);
const __m128i sh=_mm_setr_epi8(11,12,13,14,6,7,8,9,-128,-128,-128,-128,-128,-128,-128,-128);
const __m128i tail=_mm_setr_epi8(15,-128,-128,-128,10,-128,-128,-128,-128,-128,-128,-128,-128,-128,-128,-128);
const __m128i threshold=_mm_set1_epi32(49999),minusbase=_mm_set1_epi32(-100000);
const __m256i Z256=_mm256_set1_epi8((char)'0'),W1v=_mm256_set1_epi16(0x010a),W2v=_mm256_set1_epi32(0x00010064);
const __m256i SHv=_mm256_broadcastsi128_si256(sh),TAv=_mm256_broadcastsi128_si256(tail);
const __m256i THv=_mm256_set1_epi32(49999),MBv=_mm256_set1_epi32(-100000);
__m256i prevA=_mm256_setzero_si256(),prevB=_mm256_setzero_si256();
bool used256=false;
while(len0>=46 && len1>=46) {
__m256i A0=_mm256_set_m128i(_mm_loadu_si128((const __m128i*)(s0+len0-36)),
_mm_loadu_si128((const __m128i*)(s0+len0-16)));
__m256i B0=_mm256_set_m128i(_mm_loadu_si128((const __m128i*)(s0+len0-46)),
_mm_loadu_si128((const __m128i*)(s0+len0-26)));
__m256i A1=_mm256_set_m128i(_mm_loadu_si128((const __m128i*)(s1+len1-36)),
_mm_loadu_si128((const __m128i*)(s1+len1-16)));
__m256i B1=_mm256_set_m128i(_mm_loadu_si128((const __m128i*)(s1+len1-46)),
_mm_loadu_si128((const __m128i*)(s1+len1-26)));
__m256i x0=_mm256_unpacklo_epi64(_mm256_shuffle_epi8(A0,SHv),_mm256_shuffle_epi8(B0,SHv));
__m256i y0=_mm256_unpacklo_epi64(_mm256_shuffle_epi8(A0,TAv),_mm256_shuffle_epi8(B0,TAv));
__m256i x1=_mm256_unpacklo_epi64(_mm256_shuffle_epi8(A1,SHv),_mm256_shuffle_epi8(B1,SHv));
__m256i y1=_mm256_unpacklo_epi64(_mm256_shuffle_epi8(A1,TAv),_mm256_shuffle_epi8(B1,TAv));
x0=_mm256_madd_epi16(_mm256_maddubs_epi16(x0,W1v),W2v);
x1=_mm256_madd_epi16(_mm256_maddubs_epi16(x1,W1v),W2v);
x0=_mm256_sub_epi32(_mm256_add_epi32(_mm256_add_epi32(_mm256_slli_epi32(x0,3),_mm256_slli_epi32(x0,1)),y0),_mm256_set1_epi32(533328));
x1=_mm256_sub_epi32(_mm256_add_epi32(_mm256_add_epi32(_mm256_slli_epi32(x1,3),_mm256_slli_epi32(x1,1)),y1),_mm256_set1_epi32(533328));
__m256i g0=_mm256_cmpgt_epi32(x0,THv);
__m256i g1=_mm256_cmpgt_epi32(x1,THv);
__m256i P0=_mm256_permute2x128_si256(prevA,g0,0x21);
__m256i P1=_mm256_permute2x128_si256(prevB,g1,0x21);
__m256i c0=_mm256_alignr_epi8(g0,P0,12);
__m256i c1=_mm256_alignr_epi8(g1,P1,12);
x0=_mm256_add_epi32(_mm256_sub_epi32(x0,c0),_mm256_and_si256(g0,MBv));
x1=_mm256_add_epi32(_mm256_sub_epi32(x1,c1),_mm256_and_si256(g1,MBv));
prevA=g0;prevB=g1; used256=true;
_mm256_stream_pd(out0+k0,_mm256_cvtepi32_pd(_mm256_castsi256_si128(x0)));
_mm256_stream_pd(out0+k0+4,_mm256_cvtepi32_pd(_mm256_extracti128_si256(x0,1)));
k0+=8;len0-=40;
_mm256_stream_pd(out1+k1,_mm256_cvtepi32_pd(_mm256_castsi256_si128(x1)));
_mm256_stream_pd(out1+k1+4,_mm256_cvtepi32_pd(_mm256_extracti128_si256(x1,1)));
k1+=8;len1-=40;
}
if(used256){prev0=_mm256_extracti128_si256(prevA,1);prev1=_mm256_extracti128_si256(prevB,1);}
while(len0>=26 && len1>=26) {
__m128i a0=_mm_sub_epi8(_mm_loadu_si128((const __m128i*)(s0+len0-16)),zeros);
__m128i b0=_mm_sub_epi8(_mm_loadu_si128((const __m128i*)(s0+len0-26)),zeros);
__m128i a1=_mm_sub_epi8(_mm_loadu_si128((const __m128i*)(s1+len1-16)),zeros);
__m128i b1=_mm_sub_epi8(_mm_loadu_si128((const __m128i*)(s1+len1-26)),zeros);
__m128i x0=_mm_unpacklo_epi64(_mm_shuffle_epi8(a0,sh),_mm_shuffle_epi8(b0,sh));
__m128i y0=_mm_unpacklo_epi64(_mm_shuffle_epi8(a0,tail),_mm_shuffle_epi8(b0,tail));
__m128i x1=_mm_unpacklo_epi64(_mm_shuffle_epi8(a1,sh),_mm_shuffle_epi8(b1,sh));
__m128i y1=_mm_unpacklo_epi64(_mm_shuffle_epi8(a1,tail),_mm_shuffle_epi8(b1,tail));
x0=_mm_madd_epi16(_mm_maddubs_epi16(x0,w1),w2);
x1=_mm_madd_epi16(_mm_maddubs_epi16(x1,w1),w2);
x0=_mm_add_epi32(_mm_add_epi32(_mm_slli_epi32(x0,3),_mm_slli_epi32(x0,1)),y0);
x1=_mm_add_epi32(_mm_add_epi32(_mm_slli_epi32(x1,3),_mm_slli_epi32(x1,1)),y1);
__m128i g0=_mm_cmpgt_epi32(x0,threshold);
__m128i g1=_mm_cmpgt_epi32(x1,threshold);
x0=_mm_add_epi32(_mm_sub_epi32(x0,_mm_alignr_epi8(g0,prev0,12)),_mm_and_si128(g0,minusbase));
x1=_mm_add_epi32(_mm_sub_epi32(x1,_mm_alignr_epi8(g1,prev1,12)),_mm_and_si128(g1,minusbase));
prev0=g0;prev1=g1;
_mm256_stream_pd(out0+k0,_mm256_cvtepi32_pd(x0));k0+=4;len0-=20;
_mm256_stream_pd(out1+k1,_mm256_cvtepi32_pd(x1));k1+=4;len1-=20;
}
_mm_sfence();
carry0=_mm_extract_epi32(prev0,3)&1;carry1=_mm_extract_epi32(prev1,3)&1;
while(len0>=5) { int v=0;for(int j=5;j>0;--j)v=v*10+s0[len0-j]-'0';v+=carry0; carry0=v>=50000;out0[k0++]=v-carry0*100000;len0-=5; }
{ int v=0;for(std::size_t j=0;j<len0;++j)v=v*10+s0[j]-'0';out0[k0]=v+carry0; }
while(len1>=5) { int v=0;for(int j=5;j>0;--j)v=v*10+s1[len1-j]-'0';v+=carry1; carry1=v>=50000;out1[k1++]=v-carry1*100000;len1-=5; }
{ int v=0;for(std::size_t j=0;j<len1;++j)v=v*10+s1[j]-'0';out1[k1]=v+carry1; }
}
inline void parse_balanced_base100000(const char* s,std::size_t len,double* out) {
size_t k=0;int carry=0;
const __m128i zeros=_mm_set1_epi8('0'),w1=_mm_set1_epi16(0x010a),w2=_mm_set1_epi32(0x00010064);
const __m128i sh=_mm_setr_epi8(11,12,13,14,6,7,8,9,-128,-128,-128,-128,-128,-128,-128,-128);
const __m128i tail=_mm_setr_epi8(15,-128,-128,-128,10,-128,-128,-128,-128,-128,-128,-128,-128,-128,-128,-128);
const __m128i threshold=_mm_set1_epi32(49999),minusbase=_mm_set1_epi32(-100000);
while(len>=26) {
__m128i a=_mm_sub_epi8(_mm_loadu_si128((const __m128i*)(s+len-16)),zeros);
__m128i b=_mm_sub_epi8(_mm_loadu_si128((const __m128i*)(s+len-26)),zeros);
__m128i x=_mm_unpacklo_epi64(_mm_shuffle_epi8(a,sh),_mm_shuffle_epi8(b,sh));
__m128i y=_mm_unpacklo_epi64(_mm_shuffle_epi8(a,tail),_mm_shuffle_epi8(b,tail));
x=_mm_madd_epi16(_mm_maddubs_epi16(x,w1),w2);
x=_mm_add_epi32(_mm_add_epi32(_mm_slli_epi32(x,3),_mm_slli_epi32(x,1)),y);
x=_mm_add_epi32(x,_mm_cvtsi32_si128(carry));
__m128i g=_mm_cmpgt_epi32(x,threshold),p=_mm_cmpeq_epi32(x,threshold);
if(!_mm_testz_si128(p,p)){
g=_mm_or_si128(g,_mm_and_si128(p,_mm_slli_si128(g,4)));
p=_mm_and_si128(p,_mm_slli_si128(p,4));
g=_mm_or_si128(g,_mm_and_si128(p,_mm_slli_si128(g,8)));
}
x=_mm_add_epi32(_mm_sub_epi32(x,_mm_slli_si128(g,4)),_mm_and_si128(g,minusbase));
carry=_mm_extract_epi32(g,3)&1;
_mm256_store_pd(out+k,_mm256_cvtepi32_pd(x));k+=4;len-=20;
}
while(len>=5) {
int v=0;for(int j=5;j>0;--j)v=v*10+s[len-j]-'0';v+=carry;
carry=v>=50000;out[k++]=v-carry*100000;len-=5;
}
int v=0;for(size_t j=0;j<len;++j)v=v*10+s[j]-'0';out[k]=v+carry;
}
}
#include <cerrno>
static bool write_all_stdout(const char* data,size_t remaining) {
while(remaining){
const ssize_t sent=write(STDOUT_FILENO,data,remaining);
if(sent>0){data+=sent;remaining-=(size_t)sent;}
else if(sent<0 && errno==EINTR)continue;
else return false;
}
return true;
}
static inline int64_t floor_base(int64_t v) {
return int64_t((uint64_t(v)+4000000000000000000ULL)/100000ULL)-40000000000000LL;
}
static inline int64_t boundary_carry(const int64_t* c,size_t end,int64_t bound) {
size_t width=6;
for(;;) {
size_t start=end>width?end-width:0;
int64_t lo=start?-bound:0,hi=start?bound:0;
for(size_t j=start;j<end;++j) {lo=floor_base(lo+c[j]);hi=floor_base(hi+c[j]);}
if(lo==hi)return lo;width*=2;
}
}
static size_t output_leading_zeros(const int64_t* c,size_t conv,size_t digits){
size_t m=conv-1;if(m<2)return 0;
if(c[m]!=0)return 0;
size_t t=m-1;if(c[t]==0)return 0;
const size_t step=m/2;
if(step*2!=m)return 0;
if(5*step*2!=digits)return 0;
if(boundary_carry(c,m,1LL<<50)+c[m]!=0)return 0;
int64_t x=c[t]+boundary_carry(c,t,1LL<<50);
int64_t v=x%100000;if(v<0)v+=100000;
if(v<=0)return 0;
size_t dd=0;for(int64_t y=v;y;y/=10)++dd;
size_t L=5*t+dd;
if(L>=digits)return 0;
size_t sh=digits-L;
if(sh>=5)return 0;
return sh;}
__attribute__((target("bmi2"))) static size_t format_parallel(const int64_t* c,size_t conv,char* out,size_t digits,uint64_t) {
uint16_t pair[100];uint32_t triple[1000];
for(unsigned v=0;v<100;++v)pair[v]=('0'+v/10)|(('0'+v%10)<<8);
for(unsigned v=0;v<1000;++v){triple[v]=((uint32_t)('0'+v/100)<<8)|((uint32_t)('0'+v/10%10)<<16)|((uint32_t)('0'+v%10)<<24);}
enum { NC = 2 };
size_t m=conv-1,step=m/NC;
uint64_t S0=400000000000000000ULL;
uint64_t S1=(uint64_t)boundary_carry(c,step,1LL<<50)+400000000000000000ULL;
char* p0=out+digits;
char* p1=out+digits-5*step;
const int64_t* cp0=c;
const int64_t* cp1=c+step;
const uint64_t MAGIC=12089258196146291748ULL, KOFF=399996000000000000ULL;
size_t nblk=step/8, lc=nblk;
if(nblk){
__asm__ __volatile__(
"movq %[koff], %%rax\n\t"
"movq %[S0], %%r8\n\t"
"movq %[S1], %%r10\n\t"
"movq %[cp0], %%r12\n\t"
"movq %[cp1], %%r13\n\t"
"movq %[p0], %%r14\n\t"
"movq %[p1], %%r15\n\t"
"leaq %[tb], %%rbx\n\t"
"leaq %[pb], %%rdi\n\t"
".p2align 4,,10\n"
"1:\n\t"
"prefetcht0 1024(%%r12)\n\t"
"prefetcht0 1024(%%r13)\n\t"
"addq 0(%%r12), %%r8\n\t"
"mulx %%r8, %%rsi, %%r9\n\t"
"shrq $16, %%r9\n\t"
"imull $100000, %%r9d, %%ecx\n\t"
"subl %%ecx, %%r8d\n\t"
"imulq $274877907, %%r8, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r8d\n\t"
"movl (%%rbx,%%r8,4), %%esi\n\t"
"movl %%esi, -4(%%r14)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r8d\n\t"
"movw %%r8w, -5(%%r14)\n\t"
"addq %%rax, %%r9\n\t"
"addq 0(%%r13), %%r10\n\t"
"mulx %%r10, %%rsi, %%r11\n\t"
"shrq $16, %%r11\n\t"
"imull $100000, %%r11d, %%ecx\n\t"
"subl %%ecx, %%r10d\n\t"
"imulq $274877907, %%r10, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r10d\n\t"
"movl (%%rbx,%%r10,4), %%esi\n\t"
"movl %%esi, -4(%%r15)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r10d\n\t"
"movw %%r10w, -5(%%r15)\n\t"
"addq %%rax, %%r11\n\t"
"addq 8(%%r12), %%r9\n\t"
"mulx %%r9, %%rsi, %%r8\n\t"
"shrq $16, %%r8\n\t"
"imull $100000, %%r8d, %%ecx\n\t"
"subl %%ecx, %%r9d\n\t"
"imulq $274877907, %%r9, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r9d\n\t"
"movl (%%rbx,%%r9,4), %%esi\n\t"
"movl %%esi, -9(%%r14)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r9d\n\t"
"movw %%r9w, -10(%%r14)\n\t"
"addq %%rax, %%r8\n\t"
"addq 8(%%r13), %%r11\n\t"
"mulx %%r11, %%rsi, %%r10\n\t"
"shrq $16, %%r10\n\t"
"imull $100000, %%r10d, %%ecx\n\t"
"subl %%ecx, %%r11d\n\t"
"imulq $274877907, %%r11, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r11d\n\t"
"movl (%%rbx,%%r11,4), %%esi\n\t"
"movl %%esi, -9(%%r15)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r11d\n\t"
"movw %%r11w, -10(%%r15)\n\t"
"addq %%rax, %%r10\n\t"
"addq 16(%%r12), %%r8\n\t"
"mulx %%r8, %%rsi, %%r9\n\t"
"shrq $16, %%r9\n\t"
"imull $100000, %%r9d, %%ecx\n\t"
"subl %%ecx, %%r8d\n\t"
"imulq $274877907, %%r8, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r8d\n\t"
"movl (%%rbx,%%r8,4), %%esi\n\t"
"movl %%esi, -14(%%r14)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r8d\n\t"
"movw %%r8w, -15(%%r14)\n\t"
"addq %%rax, %%r9\n\t"
"addq 16(%%r13), %%r10\n\t"
"mulx %%r10, %%rsi, %%r11\n\t"
"shrq $16, %%r11\n\t"
"imull $100000, %%r11d, %%ecx\n\t"
"subl %%ecx, %%r10d\n\t"
"imulq $274877907, %%r10, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r10d\n\t"
"movl (%%rbx,%%r10,4), %%esi\n\t"
"movl %%esi, -14(%%r15)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r10d\n\t"
"movw %%r10w, -15(%%r15)\n\t"
"addq %%rax, %%r11\n\t"
"addq 24(%%r12), %%r9\n\t"
"mulx %%r9, %%rsi, %%r8\n\t"
"shrq $16, %%r8\n\t"
"imull $100000, %%r8d, %%ecx\n\t"
"subl %%ecx, %%r9d\n\t"
"imulq $274877907, %%r9, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r9d\n\t"
"movl (%%rbx,%%r9,4), %%esi\n\t"
"movl %%esi, -19(%%r14)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r9d\n\t"
"movw %%r9w, -20(%%r14)\n\t"
"addq %%rax, %%r8\n\t"
"addq 24(%%r13), %%r11\n\t"
"mulx %%r11, %%rsi, %%r10\n\t"
"shrq $16, %%r10\n\t"
"imull $100000, %%r10d, %%ecx\n\t"
"subl %%ecx, %%r11d\n\t"
"imulq $274877907, %%r11, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r11d\n\t"
"movl (%%rbx,%%r11,4), %%esi\n\t"
"movl %%esi, -19(%%r15)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r11d\n\t"
"movw %%r11w, -20(%%r15)\n\t"
"addq %%rax, %%r10\n\t"
"addq 32(%%r12), %%r8\n\t"
"mulx %%r8, %%rsi, %%r9\n\t"
"shrq $16, %%r9\n\t"
"imull $100000, %%r9d, %%ecx\n\t"
"subl %%ecx, %%r8d\n\t"
"imulq $274877907, %%r8, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r8d\n\t"
"movl (%%rbx,%%r8,4), %%esi\n\t"
"movl %%esi, -24(%%r14)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r8d\n\t"
"movw %%r8w, -25(%%r14)\n\t"
"addq %%rax, %%r9\n\t"
"addq 32(%%r13), %%r10\n\t"
"mulx %%r10, %%rsi, %%r11\n\t"
"shrq $16, %%r11\n\t"
"imull $100000, %%r11d, %%ecx\n\t"
"subl %%ecx, %%r10d\n\t"
"imulq $274877907, %%r10, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r10d\n\t"
"movl (%%rbx,%%r10,4), %%esi\n\t"
"movl %%esi, -24(%%r15)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r10d\n\t"
"movw %%r10w, -25(%%r15)\n\t"
"addq %%rax, %%r11\n\t"
"addq 40(%%r12), %%r9\n\t"
"mulx %%r9, %%rsi, %%r8\n\t"
"shrq $16, %%r8\n\t"
"imull $100000, %%r8d, %%ecx\n\t"
"subl %%ecx, %%r9d\n\t"
"imulq $274877907, %%r9, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r9d\n\t"
"movl (%%rbx,%%r9,4), %%esi\n\t"
"movl %%esi, -29(%%r14)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r9d\n\t"
"movw %%r9w, -30(%%r14)\n\t"
"addq %%rax, %%r8\n\t"
"addq 40(%%r13), %%r11\n\t"
"mulx %%r11, %%rsi, %%r10\n\t"
"shrq $16, %%r10\n\t"
"imull $100000, %%r10d, %%ecx\n\t"
"subl %%ecx, %%r11d\n\t"
"imulq $274877907, %%r11, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r11d\n\t"
"movl (%%rbx,%%r11,4), %%esi\n\t"
"movl %%esi, -29(%%r15)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r11d\n\t"
"movw %%r11w, -30(%%r15)\n\t"
"addq %%rax, %%r10\n\t"
"addq 48(%%r12), %%r8\n\t"
"mulx %%r8, %%rsi, %%r9\n\t"
"shrq $16, %%r9\n\t"
"imull $100000, %%r9d, %%ecx\n\t"
"subl %%ecx, %%r8d\n\t"
"imulq $274877907, %%r8, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r8d\n\t"
"movl (%%rbx,%%r8,4), %%esi\n\t"
"movl %%esi, -34(%%r14)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r8d\n\t"
"movw %%r8w, -35(%%r14)\n\t"
"addq %%rax, %%r9\n\t"
"addq 48(%%r13), %%r10\n\t"
"mulx %%r10, %%rsi, %%r11\n\t"
"shrq $16, %%r11\n\t"
"imull $100000, %%r11d, %%ecx\n\t"
"subl %%ecx, %%r10d\n\t"
"imulq $274877907, %%r10, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r10d\n\t"
"movl (%%rbx,%%r10,4), %%esi\n\t"
"movl %%esi, -34(%%r15)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r10d\n\t"
"movw %%r10w, -35(%%r15)\n\t"
"addq %%rax, %%r11\n\t"
"addq 56(%%r12), %%r9\n\t"
"mulx %%r9, %%rsi, %%r8\n\t"
"shrq $16, %%r8\n\t"
"imull $100000, %%r8d, %%ecx\n\t"
"subl %%ecx, %%r9d\n\t"
"imulq $274877907, %%r9, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r9d\n\t"
"movl (%%rbx,%%r9,4), %%esi\n\t"
"movl %%esi, -39(%%r14)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r9d\n\t"
"movw %%r9w, -40(%%r14)\n\t"
"addq %%rax, %%r8\n\t"
"addq 56(%%r13), %%r11\n\t"
"mulx %%r11, %%rsi, %%r10\n\t"
"shrq $16, %%r10\n\t"
"imull $100000, %%r10d, %%ecx\n\t"
"subl %%ecx, %%r11d\n\t"
"imulq $274877907, %%r11, %%rcx\n\t"
"shrq $38, %%rcx\n\t"
"imull $1000, %%ecx, %%esi\n\t"
"subl %%esi, %%r11d\n\t"
"movl (%%rbx,%%r11,4), %%esi\n\t"
"movl %%esi, -39(%%r15)\n\t"
"movzwl (%%rdi,%%rcx,2), %%r11d\n\t"
"movw %%r11w, -40(%%r15)\n\t"
"addq %%rax, %%r10\n\t"
"addq $64, %%r12\n\t"
"addq $64, %%r13\n\t"
"subq $40, %%r14\n\t"
"subq $40, %%r15\n\t"
"decq %[lc]\n\t"
"jne 1b\n\t"
"movq %%r8, %[S0]\n\t"
"movq %%r10, %[S1]\n\t"
"movq %%r12, %[cp0]\n\t"
"movq %%r13, %[cp1]\n\t"
"movq %%r14, %[p0]\n\t"
"movq %%r15, %[p1]\n\t"
: [S0]"+m"(S0), [S1]"+m"(S1), [cp0]"+m"(cp0), [cp1]"+m"(cp1), [p0]"+m"(p0), [p1]"+m"(p1), [lc]"+m"(lc)
: [magic]"d"(MAGIC), [tb]"m"(triple), [pb]"m"(pair), [koff]"m"(KOFF)
: "rax","rbx","rcx","rsi","rdi","r8","r9","r10","r11","r12","r13","r14","r15","cc","memory");
}
size_t i=nblk*8;
for(;i<step;++i) {
uint64_t t=(uint64_t)c[i]+S0;
unsigned long long hh;_mulx_u64(t,MAGIC,&hh);
uint64_t q=hh>>16;S0=q+399996000000000000ULL;
uint32_t v=(uint32_t)t-(uint32_t)(q*100000ULL);
uint32_t hi=v/1000u,lo=v-hi*1000u;
char* d=p0-5;p0=d;
*(uint32_t*)(d+1)=triple[lo];
memcpy(d,pair+hi,2);
}
for(size_t i2=0;i2<step-(nblk*8);++i2) {
size_t i3=nblk*8+i2;
uint64_t t=(uint64_t)c[step+i3]+S1;
unsigned long long hh;_mulx_u64(t,MAGIC,&hh);
uint64_t q=hh>>16;S1=q+399996000000000000ULL;
uint32_t v=(uint32_t)t-(uint32_t)(q*100000ULL);
uint32_t hi=v/1000u,lo=v-hi*1000u;
char* d=p1-5;p1=d;
*(uint32_t*)(d+1)=triple[lo];
memcpy(d,pair+hi,2);
}
size_t pos=digits-5*step*NC;
int64_t carry1=(int64_t)(S1-400000000000000000ULL);
for(size_t i=(size_t)NC*step;i<m;++i) {
int64_t x=carry1+c[i],q=floor_base(x);
unsigned v=x-q*100000;carry1=q;pos-=5;
*(uint32_t*)(out+pos+1)=triple[v%1000];
memcpy(out+pos,pair+v/1000,2);
}
int64_t top=carry1+c[m];
while(top>0) {out[--pos]='0'+top%10;top/=10;}
if(pos==digits)out[--pos]='0';
while(pos+1<digits&&out[pos]=='0')++pos;
return pos;
}
namespace hint { namespace transform { namespace fft {
static void real_conv5(F64* a,F64* b,size_t total,size_t used_a,size_t used_b) {
const size_t child=total/5,L=child/2;
FFTSqrtTableC64X4<6> roots(total/2,5,1);
const F64X4 c1(0.30901699437494742410229341718282),c2(-0.80901699437494742410229341718282),
s1(0.95105651629515357211643933337938),s2(0.58778525229247312916870595463907),zero(0.0);
auto negI=[&](const C64X4& z){return C64X4(z.imag,zero-z.real);};
auto posI=[&](const C64X4& z){return C64X4(zero-z.imag,z.real);};
auto forward_stage=[&](auto has2_tag,F64* x,size_t j0,size_t jend){
constexpr bool HAS2=decltype(has2_tag)::value;
for(size_t j=j0;j<jend;j+=8) {
_mm_prefetch((const char*)(x+j+512),_MM_HINT_T0);
_mm_prefetch((const char*)(x+child+j+512),_MM_HINT_T0);
if constexpr(HAS2)_mm_prefetch((const char*)(x+2*child+j+512),_MM_HINT_T0);
C64X4 x0,x1,x2;
x0.load(x+j,std::true_type{});
x1.load(x+child+j,std::true_type{});
if constexpr(HAS2)x2.load(x+2*child+j,std::true_type{});
C64X4 t1,t2,u1,u2;
if constexpr(HAS2) {
t1=x0+x1*c1+x2*c2;
t2=x0+x1*c2+x2*c1;
u1=negI(x1*s1+x2*s2);
u2=negI(x1*s2-x2*s1);
} else {
t1=x0+x1*c1;
t2=x0+x1*c2;
u1=negI(x1*s1);
u2=negI(x1*s2);
}
C64X4 w=roots[j/8],w2=w.mul(w),w3=w2.mul(w),w4=w2.mul(w2);
if constexpr(HAS2)(x0+x1+x2).store(x+j);
else (x0+x1).store(x+j);
(t1+u1).mul(w).store(x+child+j);
(t2+u2).mul(w2).store_nt(x+2*child+j);
(t2-u2).mul(w3).store_nt(x+3*child+j);
(t1-u1).mul(w4).store_nt(x+4*child+j);
}
};
auto forward=[&](F64* x,size_t used){
size_t active2=used>2*child?used-2*child:0;
size_t split=std::min(child,(active2+7)&~size_t(7));
forward_stage(std::true_type{},x,0,split);
forward_stage(std::false_type{},x,split,child);
_mm_sfence();
for(int r=0;r<5;++r)FFTAVX::difRecLong<false>(x+r*child,child);
};
forward(a,used_a);forward(b,used_b);
real_dot_binrev4(a,b,child,0.2);
const F64X4 inv(0.5/total);
for(int r=1;r<=2;++r) {
double theta=-HINT_2PI*r/(5*L);C64X4 offset(std::cos(theta),std::sin(theta));
BinRevTableC64X4HP<28,1> table;
for(size_t j=0;j<child;j+=8) {
if(__builtin_expect(j + 192 < child,1)) {
_mm_prefetch((const char*)(a+r*child+j+192),_MM_HINT_T0);
_mm_prefetch((const char*)(a+(6-r)*child-8-j-192),_MM_HINT_T0);
_mm_prefetch((const char*)(b+r*child+j+192),_MM_HINT_T0);
_mm_prefetch((const char*)(b+(6-r)*child-8-j-192),_MM_HINT_T0);
}
dot_rfftX4(a+r*child+j,a+(6-r)*child-8-j,b+r*child+j,b+(6-r)*child-8-j,table.iterate().mul(offset),inv);
}
}
for(int r=0;r<5;++r)FFTAVX::iditRecLong<false,false>(a+r*child,child);
const size_t conv0=used_a+used_b-1;
size_t jlim=0;
if(conv0>4*child){ jlim=(conv0-4*child<child)?(conv0-4*child):child; jlim&=~size_t(7); }
for(size_t j=0;j<jlim;j+=8) {
if(__builtin_expect(j + 496 < child,1)) {
_mm_prefetch((const char*)(a+j+256),_MM_HINT_T0);
_mm_prefetch((const char*)(a+child+j+316),_MM_HINT_T0);
_mm_prefetch((const char*)(a+2*child+j+376),_MM_HINT_T0);
_mm_prefetch((const char*)(a+3*child+j+436),_MM_HINT_T0);
_mm_prefetch((const char*)(a+4*child+j+496),_MM_HINT_T0);
}
C64X4 w=roots[j/8],w2=w.mul(w),w3=w2.mul(w),w4=w2.mul(w2);
C64X4 x0=a+j,x1=C64X4(a+child+j).mulConj(w),x2=C64X4(a+2*child+j).mulConj(w2),
x3=C64X4(a+3*child+j).mulConj(w3),x4=C64X4(a+4*child+j).mulConj(w4);
C64X4 sum1=x1+x4,sum2=x2+x3,diff1=x1-x4,diff2=x2-x3;
C64X4 t1=x0+sum1*c1+sum2*c2,t2=x0+sum1*c2+sum2*c1,u1=posI(diff1*s1+diff2*s2),u2=posI(diff1*s2-diff2*s1);
(x0+sum1+sum2).transToI64(std::true_type{}).store(a+j,std::true_type{});
(t1+u1).transToI64(std::true_type{}).store(a+child+j,std::true_type{});
(t2+u2).transToI64(std::true_type{}).store(a+2*child+j,std::true_type{});
(t2-u2).transToI64(std::true_type{}).store(a+3*child+j,std::true_type{});
(t1-u1).transToI64(std::true_type{}).store(a+4*child+j,std::true_type{});
}
for(size_t j=jlim;j<child;j+=8) {
if(__builtin_expect(j + 496 < child,1)) {
_mm_prefetch((const char*)(a+j+256),_MM_HINT_T0);
_mm_prefetch((const char*)(a+child+j+316),_MM_HINT_T0);
_mm_prefetch((const char*)(a+2*child+j+376),_MM_HINT_T0);
_mm_prefetch((const char*)(a+3*child+j+436),_MM_HINT_T0);
}
C64X4 w=roots[j/8],w2=w.mul(w),w3=w2.mul(w),w4=w2.mul(w2);
C64X4 x0=a+j,x1=C64X4(a+child+j).mulConj(w),x2=C64X4(a+2*child+j).mulConj(w2),
x3=C64X4(a+3*child+j).mulConj(w3),x4=C64X4(a+4*child+j).mulConj(w4);
C64X4 sum1=x1+x4,sum2=x2+x3,diff1=x1-x4,diff2=x2-x3;
C64X4 t1=x0+sum1*c1+sum2*c2,t2=x0+sum1*c2+sum2*c1,u1=posI(diff1*s1+diff2*s2),u2=posI(diff1*s2-diff2*s1);
(x0+sum1+sum2).transToI64(std::true_type{}).store(a+j,std::true_type{});
(t1+u1).transToI64(std::true_type{}).store(a+child+j,std::true_type{});
(t2+u2).transToI64(std::true_type{}).store(a+2*child+j,std::true_type{});
(t2-u2).transToI64(std::true_type{}).store(a+3*child+j,std::true_type{});
}
}
}}}
#include <vector>
#include <complex>
#include <cmath>
namespace r4 {
using namespace hint::transform::fft;
static constexpr size_t RSYNC = 1024;
struct alignas(64) SplitTwiddles {
double real[1024], imag[1024];
SplitTwiddles(){for(size_t i=0;i<1024;i++){size_t r=0;for(size_t j=0;j<10;j++)r=(r<<1)|((i>>j)&1);double th=-3.14159265358979323846*2.0*(double)r/2048.0;real[i]=std::cos(th);imag[i]=std::sin(th);}}
};
static inline const SplitTwiddles& ts_table(){ static const SplitTwiddles t; return t; }
static inline size_t brev(size_t v, size_t w){ size_t r=0; for(size_t i=0;i<w;i++) r=(r<<1)|((v>>i)&1); return r; }
using F64=double; using C64X4=hint::Complex64X4; using F64X4=hint::Float64X4;
static inline C64X4 negI(const C64X4& z){ return C64X4(z.imag, F64X4(0.0)-z.real); }
static inline C64X4 posI(const C64X4& z){ return C64X4(F64X4(0.0)-z.imag, z.real); }
static void mid23(F64* a, size_t cs){
const size_t sa = cs>>2, sb = cs>>4;
const C64X4* TA = reinterpret_cast<const C64X4*>(FFTAVX::multi_table_1.getBegin(cs>>1));
const C64X4* TB = reinterpret_cast<const C64X4*>(FFTAVX::multi_table_1.getBegin(cs>>3));
for(size_t j0=0;j0<sb;j0+=8){
C64X4 A[4][4];
for(int k2=0;k2<4;k2++){
const size_t j = (size_t)k2*(sb>>3) + (j0>>3);
F64* pj = a + ((size_t)k2*sb + j0);
_mm_prefetch((const char*)(pj+sa+192+8*k2),_MM_HINT_T0);
_mm_prefetch((const char*)(pj+2*sa+200+8*k2),_MM_HINT_T0);
_mm_prefetch((const char*)(pj+3*sa+208+8*k2),_MM_HINT_T0);
C64X4 c0(pj), c1(pj+sa), c2(pj+2*sa), c3(pj+3*sa);
FFTAVX::dif4(c0.real,c0.imag,c1.real,c1.imag,c2.real,c2.imag,c3.real,c3.imag);
C64X4 w = TA[j], 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;
}
C64X4 w = TB[j0>>3], w2 = w.mul(w), w3 = w2.mul(w);
for(int k1=0;k1<4;k1++){
C64X4 c0=A[0][k1], c1=A[1][k1], c2=A[2][k1], c3=A[3][k1];
FFTAVX::dif4(c0.real,c0.imag,c1.real,c1.imag,c2.real,c2.imag,c3.real,c3.imag);
c2 = c2.mul(w); c1 = c1.mul(w2); c3 = c3.mul(w3);
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++)
FFTAVX::difRecMid<false>(a + (size_t)k1*sa + (size_t)k2*sb, sb);
}
static void fwd4(F64* p, size_t child, FFTSqrtTableC64X4<6>& roots, size_t used){
const size_t cs = child >> 2;
const bool zq = (used <= (child<<1));
const size_t flen = child >> 1;
if(flen <= FFTAVX::MID_LEN || flen >= (size_t(1)<<22)){
for(size_t j=0;j<child;j+=8){
C64X4 x0; x0.load(p+j, std::true_type{});
C64X4 x1; x1.load(p+child+j, std::true_type{});
C64X4 x2; x2.load(p+2*child+j, std::true_type{});
C64X4 x3; x3.load(p+3*child+j, std::true_type{});
C64X4 t0=x0+x2,t1=x0-x2,t2=x1+x3,t3=x1-x3;
C64X4 y0=t0+t2,y2=t0-t2,y1=t1+negI(t3),y3=t1+posI(t3);
C64X4 w=roots[j/8],w2=w.mul(w),w3=w2.mul(w);
y0.store(p+j); y1.mul(w).store(p+child+j); y2.mul(w2).store_nt(p+2*child+j); y3.mul(w3).store_nt(p+3*child+j);
}
_mm_sfence();
for(int r=0;r<4;++r) FFTAVX::difRecLong<false>(p+r*child,child);
return;
}
const auto &t2 = FFTAVX::sqrt_table_1[hint::hint_log2(flen)];
constexpr size_t SL = 512;
for(size_t j0=0;j0<cs;j0+=SL){
size_t Tj = cs-j0; if(Tj>SL) Tj=SL;
for(int s=0;s<4;++s){
F64 *b0=p+(size_t)s*cs+j0, *b1=b0+child, *b2=b1+child, *b3=b2+child;
for(size_t j=0;j<Tj;j+=8){
_mm_prefetch((const char*)(b0+j+16384),_MM_HINT_T0);
_mm_prefetch((const char*)(b1+j+16384),_MM_HINT_T0);
C64X4 x0; x0.load(b0+j, std::true_type{});
C64X4 x1; x1.load(b1+j, std::true_type{});
C64X4 y0,y1,y2,y3;
if(!zq){
C64X4 x2; x2.load(b2+j, std::true_type{});
C64X4 x3; x3.load(b3+j, std::true_type{});
C64X4 r0=x0+x2,r1=x0-x2,r2=x1+x3,r3=x1-x3;
y0=r0+r2;y2=r0-r2;y1=r1+negI(r3);y3=r1+posI(r3);
} else {
y0=x0+x1; y1=x0+negI(x1); y2=x0-x1; y3=x0+posI(x1);
}
C64X4 w=roots[(s*cs+j0+j)/8],w2=w.mul(w),w3=w2.mul(w);
y0.store(b0+j); y1.mul(w).store(b1+j); y2.mul(w2).store(b2+j); y3.mul(w3).store(b3+j);
}
}
for(int r=0;r<4;++r){
F64* q=p+(size_t)r*child+j0;
for(size_t j=0;j<Tj;j+=8){
_mm_prefetch((const char*)(q+j+192),_MM_HINT_T0);
_mm_prefetch((const char*)(q+cs+j+192),_MM_HINT_T0);
_mm_prefetch((const char*)(q+2*cs+j+192),_MM_HINT_T0);
_mm_prefetch((const char*)(q+3*cs+j+192),_MM_HINT_T0);
C64X4 c0(q+j), c1(q+cs+j), c2(q+2*cs+j), c3(q+3*cs+j), omega1, omega2;
FFTAVX::dif4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
omega1 = t2[(j0+j)/8], c2 = c2.mul(omega1);
omega2 = omega1.mul(omega1), c1 = c1.mul(omega2);
c3 = c3.mul(omega2.mul(omega1));
c0.store(q+j), c1.store(q+cs+j), c2.store(q+2*cs+j), c3.store(q+3*cs+j);
}
}
}
_mm_sfence();
for(int r=0;r<4;++r) for(int s=0;s<4;++s){
F64* q = p+(size_t)r*child+(size_t)s*cs;
if(cs>=(size_t(1)<<15) && cs<=(size_t(1)<<18)) mid23(q, cs);
else FFTAVX::difRecMid<false>(q, cs);
}
}
template<bool TO_INT=false>
static void inv4(F64* p, F64* pb, size_t child, FFTSqrtTableC64X4<6>& roots, double scale=1.0, bool seam=false){
if(seam){ dotjit::a=p; dotjit::b=pb; dotjit::child=child; }
for(int r=0;r<4;++r){ dotjit::cid=r; dotjit::set_child(r); dotjit::on=seam?1:0;
FFTAVX::iditRecLong<false,false>(p+r*child,child); }
dotjit::on=0;
for(size_t j=0;j<child;j+=8){
_mm_prefetch((const char*)(p+j+256),_MM_HINT_T0);
_mm_prefetch((const char*)(p+child+j+264),_MM_HINT_T0);
_mm_prefetch((const char*)(p+2*child+j+272),_MM_HINT_T0);
_mm_prefetch((const char*)(p+3*child+j+280),_MM_HINT_T0);
C64X4 w=roots[j/8],w2=w.mul(w),w3=w2.mul(w);
C64X4 x0=p+j,x1=C64X4(p+child+j).mulConj(w),x2=C64X4(p+2*child+j).mulConj(w2),x3=C64X4(p+3*child+j).mulConj(w3);
C64X4 t0=x0+x2,t1=x0-x2,t2=x1+x3,t3=x1-x3;
if constexpr(TO_INT){
const F64X4 sc(scale);
(((t0+t2))* sc).transToI64(std::true_type{}).store(p+j, std::true_type{});
(((t1+posI(t3)))* sc).transToI64(std::true_type{}).store(p+child+j, std::true_type{});
(((t0-t2))* sc).transToI64(std::true_type{}).store(p+2*child+j, std::true_type{});
(((t1-posI(t3)))* sc).transToI64(std::true_type{}).store(p+3*child+j, std::true_type{});
} else {
(t0+t2).store(p+j); (t1+posI(t3)).store(p+child+j); (t0-t2).store(p+2*child+j); (t1-posI(t3)).store(p+3*child+j);
}
}
}
static void dot_vec(F64* a, F64* b, size_t child, size_t total){
const size_t R=total, Lc=child/2;
const double PI=3.14159265358979323846;
real_dot_binrev4(a, b, child, (double)Lc);
size_t L=0; while(((size_t)1<<L) < Lc) ++L;
std::vector<std::complex<double>> unit(L);
for(size_t j=0;j<L;j++){ double th=-2*PI*(double)((size_t)1<<j)/(double)(2.0*(double)Lc);
unit[j]=std::complex<double>(std::cos(th),std::sin(th)); }
const F64X4 half(0.5);
const C64X4 negHalfI(F64X4(0.0), F64X4(-0.5));
auto conjC=[&](const C64X4& v){ return C64X4(v.real, F64X4(0.0)-v.imag); };
auto mulI =[&](const C64X4& v){ return C64X4(F64X4(0.0)-v.imag, v.real); };
auto run_child=[&](int r, bool selfpair){
std::complex<double> cur(1,0); size_t m=0;
double th=-2*PI*(double)r/(double)R; std::complex<double> off(std::cos(th),std::sin(th));
const auto& TS=ts_table();
const size_t steps = selfpair ? (child/2)/8 : child/8;
F64* it1=(selfpair? a+3*child-8 : a+(5-r)*child-8);
F64* jt1=(selfpair? b+3*child-8 : b+(5-r)*child-8);
F64* it0=a+r*child; F64* jt0=b+r*child;
for(size_t begin=0; begin<steps; begin+=RSYNC/4){
std::complex<double> sc=std::polar(1.0, -3.14159265358979323846*2.0*(double)brev(begin/(RSYNC/4), L-10)/(double)(2*Lc));
std::complex<double> segc=sc*off;
const __m256d sr=_mm256_set1_pd(segc.real()), si=_mm256_set1_pd(segc.imag());
const size_t extent=std::min(steps-begin, RSYNC/4);
for(size_t st=0; st<extent; st++, it0+=8, it1-=8, jt0+=8, jt1-=8){
_mm_prefetch((const char*)(it0+256),_MM_HINT_T0);
_mm_prefetch((const char*)(it1-256),_MM_HINT_T0);
_mm_prefetch((const char*)(jt0+256),_MM_HINT_T0);
_mm_prefetch((const char*)(jt1-256),_MM_HINT_T0);
__m256d tr=_mm256_load_pd(TS.real+4*st), ti=_mm256_load_pd(TS.imag+4*st);
C64X4 w(F64X4(_mm256_sub_pd(_mm256_mul_pd(tr,sr),_mm256_mul_pd(ti,si))),
F64X4(_mm256_add_pd(_mm256_mul_pd(tr,si),_mm256_mul_pd(ti,sr))));
C64X4 x0 = it0, x1 = C64X4(it1).reverse();
C64X4 y0 = jt0, y1 = C64X4(jt1).reverse();
C64X4 cx1=conjC(x1), cy1=conjC(y1);
C64X4 Ea=(x0+cx1)*half, Oa=(x0-cx1).mul(negHalfI);
C64X4 Eb=(y0+cy1)*half, Ob=(y0-cy1).mul(negHalfI);
C64X4 wo=w.mul(Oa), wb=w.mul(Ob);
C64X4 Xak=Ea+wo, Xakc=Ea-wo, Xbk=Eb+wb, Xbkc=Eb-wb;
C64X4 Pk=Xak.mul(Xbk), Pkc=Xakc.mul(Xbkc);
C64X4 Ep=(Pk+Pkc)*half;
C64X4 Op=(Pk-Pkc).mulConj(w)*half;
C64X4 Zp=Ep+mulI(Op);
C64X4 Zq=conjC(Ep-mulI(Op));
Zp.store(it0); Zq.reverse().store(it1); } } };
run_child(1,false);
run_child(2,true);
}
static void pack_inplace(F64* p, size_t total, size_t used){
(void)used;
for(size_t j=0;j<total;j+=8){
C64X4 v; v.load(p+j, std::true_type{});
v.store(p+j); }
}
static void unpack_int64(F64* p, size_t total, double scale){
const F64X4 sc(scale);
int64_t* o=reinterpret_cast<int64_t*>(p);
for(size_t j=0;j<total;j+=8){
C64X4 v(p+j);
v = v * sc;
v.transToI64(std::true_type{}).store(o+j, std::true_type{}); }
}
static void real_conv4(F64* a, F64* b, size_t total, size_t used_a, size_t used_b){
const size_t child=total/4;
FFTSqrtTableC64X4<6> roots(total/2,4,1);
fwd4(a,child,roots,used_a); fwd4(b,child,roots,used_b);
const bool seam=(child>=(size_t)4096)&&((child&(child-1))==0);
if(!seam) dot_vec(a,b,child,total);
inv4<true>(a,b,child,roots,2.0/(double)total,seam);
}
}
namespace hint { namespace transform { namespace fft { namespace dotjit {
int on = 0;
Float64 *a = 0, *b = 0;
size_t child = 0;
int cid = 0;
size_t sw_R = 0, sw_Lc = 0, sw_L = 0;
double sw_offr = 0, sw_offi = 0;
size_t cached_segment = (size_t)-1;
double cached_segment_real = 0, cached_segment_imag = 0;
void set_child(int r){
cached_segment = (size_t)-1;
sw_R = 4 * child; sw_Lc = child / 2;
size_t L = 0; while (((size_t)1 << L) < sw_Lc) ++L;
sw_L = L;
const double PI = 3.14159265358979323846;
double th = -2 * PI * (double)r / (double)sw_R;
sw_offr = std::cos(th); sw_offi = std::sin(th);
}
static BinRevTableC64X4HP<28, 1>& tab0ref(){ static BinRevTableC64X4HP<28, 1> t; return t; }
static size_t tab0_level = (size_t)-1;
static inline C64X4 conjv(const C64X4& v){ return C64X4(v.real, F64X4(0.0) - v.imag); }
static inline C64X4 muli(const C64X4& v){ return C64X4(F64X4(0.0) - v.imag, v.real); }
static void seam0_head(Float64* A, Float64* B, size_t ch){
const Float64 scale = (Float64)(ch / 2);
Float64 inv = 2.0 * scale / (Float64)ch;
auto r0 = A[0], i0 = A[4], r1 = B[0], i1 = B[4];
transform2(r0, i0);
transform2(r1, i1);
r0 *= r1, i0 *= i1;
transform2(r0, i0);
A[0] = r0 * 0.5 * inv, A[4] = i0 * 0.5 * inv;
auto temp = C64(A[1], A[5]) * C64(B[1], B[5]) * inv;
A[1] = temp.real(), A[5] = temp.imag();
inv /= 4;
dot_rfft<4>(&A[2], &A[3], &B[2], &B[3], C64(COS_PI_8, -COS_PI_8), inv);
constexpr Float64 COS_16_1 = 0.92387953251128675612818318939;
constexpr Float64 SIN_16_1 = 0.38268343236508977172845998403;
dot_rfft<4>(&A[8], &A[11], &B[8], &B[11], C64(COS_16_1, -SIN_16_1), inv);
dot_rfft<4>(&A[9], &A[10], &B[9], &B[10], C64(-SIN_16_1, -COS_16_1), inv);
const Float64X4 inv4 = F64X4(0.5 * scale / (Float64)ch);
BinRevTableC64X4HP<28, 1> tb;
for (size_t begin = 16; begin < 512; begin *= 2){
tb.reset(begin / 2);
auto it0 = A + begin, it1 = it0 + begin - 8, it2 = B + begin, it3 = it2 + begin - 8;
for (; it0 < it1; it0 += 8, it1 -= 8, it2 += 8, it3 -= 8)
{
_mm_prefetch((const char*)(it0+256),_MM_HINT_T0);
_mm_prefetch((const char*)(it1-256),_MM_HINT_T0);
_mm_prefetch((const char*)(it2+256),_MM_HINT_T0);
_mm_prefetch((const char*)(it3-256),_MM_HINT_T0);
dot_rfftX4(it0, it1, it2, it3, tb.iterate(), inv4);}}
}
static void seam0_level(Float64* A, Float64* B, size_t ch, size_t s, size_t k){
const Float64 scale = (Float64)(ch / 2);
const Float64X4 inv4 = F64X4(0.5 * scale / (Float64)ch);
const size_t JB = s >> 9;
const size_t t0 = 64 * (k - JB);
if (s != tab0_level){ tab0ref().reset(s / 2); tab0_level = s; }
auto it0 = A + s + 8 * t0, it1 = A + 2 * s - 8 - 8 * t0;
auto it2 = B + s + 8 * t0, it3 = B + 2 * s - 8 - 8 * t0;
const size_t np = (t0 + 64 <= (s >> 4)) ? 64 : ((s >> 4) - t0);
for (size_t i = 0; i < np; ++i){
_mm_prefetch((const char*)(it0+256),_MM_HINT_T0);
_mm_prefetch((const char*)(it1-256),_MM_HINT_T0);
_mm_prefetch((const char*)(it2+256),_MM_HINT_T0);
_mm_prefetch((const char*)(it3-256),_MM_HINT_T0);
dot_rfftX4(it0, it1, it2, it3, tab0ref().iterate(), inv4);
it0 += 8; it1 -= 8; it2 += 8; it3 -= 8;}
}
static void seam_pair(int cc, size_t k){
(void)cc;
const auto& TS = r4::ts_table();
const size_t t0 = 64 * k;
F64 *it0, *it1, *jt0, *jt1;
if (cid == 1){
it0 = a + child + 512 * k; it1 = a + 4 * child - 8 - 512 * k;
jt0 = b + child + 512 * k; jt1 = b + 4 * child - 8 - 512 * k;}
else {
it0 = a + 2 * child + 512 * k; it1 = a + 3 * child - 8 - 512 * k;
jt0 = b + 2 * child + 512 * k; jt1 = b + 3 * child - 8 - 512 * k;}
const F64X4 half(0.5);
const C64X4 negHalfI(F64X4(0.0), F64X4(-0.5));
size_t step = t0, left = 64;
while (left){
const size_t s2 = step >> 8, st = step & 255;
if (s2 != cached_segment) {
std::complex<double> sc = std::polar(1.0, -3.14159265358979323846 * 2.0 * (double)r4::brev(s2, sw_L - 10) / (double)(2 * sw_Lc));
std::complex<double> segc = sc * std::complex<double>(sw_offr, sw_offi);
cached_segment_real = segc.real(); cached_segment_imag = segc.imag();
cached_segment = s2;
}
const __m256d sr = _mm256_set1_pd(cached_segment_real), si = _mm256_set1_pd(cached_segment_imag);
size_t n = 256 - st; if (n > left) n = left;
for (size_t u = 0; u < n; ++u){
_mm_prefetch((const char*)(it0+256),_MM_HINT_T0);
_mm_prefetch((const char*)(it1-256),_MM_HINT_T0);
_mm_prefetch((const char*)(jt0+256),_MM_HINT_T0);
_mm_prefetch((const char*)(jt1-256),_MM_HINT_T0);
__m256d tr = _mm256_load_pd(TS.real + 4 * (st + u)), ti = _mm256_load_pd(TS.imag + 4 * (st + u));
C64X4 w(F64X4(_mm256_sub_pd(_mm256_mul_pd(tr, sr), _mm256_mul_pd(ti, si))),
F64X4(_mm256_add_pd(_mm256_mul_pd(tr, si), _mm256_mul_pd(ti, sr))));
C64X4 x0 = it0, x1 = C64X4(it1).reverse();
C64X4 y0 = jt0, y1 = C64X4(jt1).reverse();
C64X4 cx1 = conjv(x1), cy1 = conjv(y1);
C64X4 Ea = (x0 + cx1) * half, Oa = (x0 - cx1).mul(negHalfI);
C64X4 Eb = (y0 + cy1) * half, Ob = (y0 - cy1).mul(negHalfI);
C64X4 wo = w.mul(Oa), wb = w.mul(Ob);
C64X4 Xak = Ea + wo, Xakc = Ea - wo, Xbk = Eb + wb, Xbkc = Eb - wb;
C64X4 Pk = Xak.mul(Xbk), Pkc = Xakc.mul(Xbkc);
C64X4 Ep = (Pk + Pkc) * half;
C64X4 Op = (Pk - Pkc).mulConj(w) * half;
C64X4 Zp = Ep + muli(Op);
C64X4 Zq = conjv(Ep - muli(Op));
Zp.store(it0); Zq.reverse().store(it1);
it0 += 8; it1 -= 8; jt0 += 8; jt1 -= 8;}
step += n; left -= n;}
}
static void seam_cell(int cc, size_t k){
if (cc == 0){
if (k == 0){ seam0_head(a, b, child); return; }
size_t JB = 1; while ((JB << 1) <= k) JB <<= 1;
if (k < JB + ((JB + 1) >> 1)) seam0_level(a, b, child, JB << 9, k);
return;}
if (cc == 3) return;
if (cc == 2 && k >= (child >> 9) / 2) return;
seam_pair(cc == 2 ? 2 : 1, k);
}
void leaf(Float64* p, size_t n){
for (size_t u = 0; u < n; u += 512){
const size_t k = (size_t)((p + u) - (a + (size_t)cid * child)) >> 9;
seam_cell(cid, k);}
}
}}}}
int main(){
struct stat st{};
std::unique_ptr<char[]> storage;
std::string fallback;
const char* input=nullptr;
size_t input_size=0;
auto* duck=reinterpret_cast<duck_public_stdio::DuckInfo_t*>(getauxval(duck_public_stdio::AT_DUCK));
if(duck && duck->abi_version==40 && duck->stdin_ptr){
input=duck->stdin_ptr;input_size=duck->stdin_size;
}else if(fstat(STDIN_FILENO,&st)==0 && st.st_size>0){
size_t capacity=(size_t)st.st_size;
storage.reset(new char[capacity]);
while(input_size<capacity){
size_t got=fread(storage.get()+input_size,1,capacity-input_size,stdin);
if(!got)break;
input_size+=got;
}
input=storage.get();
}else{
char block[1<<16];size_t got;
while((got=fread(block,1,sizeof block,stdin)))fallback.append(block,got);
input=fallback.data();input_size=fallback.size();
}
const char* newline=(input_size>10000000 && input[10000000]=='\n') ? input+10000000 : (const char*)memchr(input,'\n',input_size);
if(!newline)return 2;
size_t cut=(size_t)(newline-input);
size_t l1=cut;while(l1&&input[l1-1]<'0')--l1;
size_t off=cut+1;while(off<input_size&&input[off]<'0')++off;
size_t l2=input_size-off;while(l2&&input[off+l2-1]<'0')--l2;
size_t na=l1/5+1,nb=l2/5+1,conv=na+nb-1,N=std::max<size_t>(16384,hint::int_ceil2(na+nb+4));
hint::AlignMem<double>A(N),B(N+256);
bigint_io_opt::parse_balanced_base100000_dual(input,l1,A.begin(),input+off,l2,B.begin()+256);
storage.reset();std::string().swap(fallback);
{
double* pa=A.begin(); double* pb=B.begin()+256;
const size_t PG=512;
for(size_t i=na;i<N;i+=PG) pa[i]=0.0;
for(size_t i=nb;i<N;i+=PG) pb[i]=0.0;
}
r4::real_conv4(A.begin(),B.begin()+256,N,na,nb);
int64_t* coeff=(int64_t*)A.begin();
const size_t digits=l1+l2;
bool direct=duck && duck->abi_version==40 && duck->stdout_ptr && duck->stdout_limit>=digits+1;
char* out=direct?duck->stdout_ptr:(char*)B.begin();
size_t lead=direct?output_leading_zeros(coeff,conv,digits):0;
char* dst=lead?out-lead:out;
size_t pos=format_parallel(coeff,conv,dst,digits,std::min(na,nb)*999ULL+1);
size_t len=digits-pos;
if(direct) {
if(pos&&!lead)memmove(out,out+pos,len);
out[len]='\n';duck->stdout_size=len+1;
} else {
if(!write_all_stdout(out+pos,len) || !write_all_stdout("\n",1))return 3;
}
return 0;
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 57.775 ms | 84 MB + 356 KB | Accepted | Score: 100 | 显示更多 |