// Ready-to-submit V8 for https://duck.ac/problem/1004e8
// Select C++17. Single file; no external headers or libraries are required.
// Uses the same FFT, parser and formatter as src/staged/v8_fast_io.cpp.
// JudgeDuck: read public stdin memory, write public stdout memory, set stdout_size.
// Linux: fall back to ordinary stdin/stdout. No benchmark or diagnostic output.
// Adapted from https://duck.ac/submission/47426 (TSKY), HintFFT MIT.
// FFT attribution: TSKY (WithSky), HintFFT, https://github.com/With-Sky/HintFFT
// Original public adaptation: saffah_codex_maker, duck.ac/submission/86930.
// Optimizations in this version: radix-5 real FFT, balanced base 100000,
// SIMD parsing, certified block carries, and direct public standard streams.
// HintFFT is distributed under the MIT license.
/*
MIT License -- HintFFT by TSKY (WithSky)
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>
#pragma GCC target("avx2,fma,tune=skylake")
#pragma GCC optimize("O3,unroll-loops")
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);}
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);}
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);}};
struct FFTAVX : public FFT{
static constexpr int LOG_SHORT = 10, LOG_MID = 25 - crc32(statement) % 15, LOG_MAX = 22, LOG_CACHE = 7;
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++) \
{ \
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; \
for (auto end = it1; it0 < end; it0 += 8, it1 += 8, it2 += 8, it3 += 8, indx++) \
{ \
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)
{
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);}
// Fuse two long radix-4 stages while a 16-vector tile is in L1.
// Preserves the original butterfly formulas and bit-reversed ordering.
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){
C64X4 v[16];
#pragma GCC unroll 16
for(int k=0;k<16;k++)v[k]=C64X4(a+k*stride+j);
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++)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]);
}else{
#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]);
#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]);
}
#pragma GCC unroll 16
for(int k=0;k<16;k++)v[k].store(a+k*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)
{
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>
// Public standard-stream ABI 0.04 (version 40), exact 64-bit field layout:
// https://github.com/JudgeDuck/JudgeDuck-OS/blob/d4df797bad6dc9b66c00312e97c05ca33adc3abd/inc/abi.hpp
// The public stdin and stdout fields implement zero-copy standard streams.
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(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;
}
} // namespace bigint_io_opt
#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 format_parallel(const int64_t* c,size_t conv,char* out,size_t digits,uint64_t) {
uint16_t pair[100];char triple[1000][3];
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][0]='0'+v/100;triple[v][1]='0'+v/10%10;triple[v][2]='0'+v%10;}
size_t m=conv-1,step=m/2;
int64_t carry[2]={}; for(int r=1;r<2;++r)carry[r]=boundary_carry(c,r*step,1LL<<50);
char* p[2]; for(int r=0;r<2;++r)p[r]=out+digits-5*r*step;
for(size_t i=0;i<step;++i) {
#pragma GCC unroll 2
for(int r=0;r<2;++r) {
int64_t x=carry[r]+c[r*step+i],q=floor_base(x);
unsigned v=x-q*100000;carry[r]=q;p[r]-=5;memcpy(p[r],pair+v/1000,2);memcpy(p[r]+2,triple[v%1000],3);
}
}
size_t pos=digits-10*step;
for(size_t i=2*step;i<m;++i) {
int64_t x=carry[1]+c[i],q=floor_base(x);
unsigned v=x-q*100000;carry[1]=q;pos-=5;memcpy(out+pos,pair+v/1000,2);memcpy(out+pos+2,triple[v%1000],3);
}
int64_t top=carry[1]+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) {
const size_t child=total/5,L=child/2;
FFTSqrtTableC64X4<7> 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=[&](F64* x){
for(size_t j=0;j<child;j+=8) {
C64X4 x0,x1,x2;x0.load(x+j,std::true_type{});x1.load(x+child+j,std::true_type{});
if(j<L)x2.load(x+2*child+j,std::true_type{});
C64X4 t1=x0+x1*c1+x2*c2,t2=x0+x1*c2+x2*c1,u1=negI(x1*s1+x2*s2),u2=negI(x1*s2-x2*s1);
C64X4 w=roots[j/8],w2=w.mul(w),w3=w2.mul(w),w4=w2.mul(w2);
(x0+x1+x2).store(x+j);
(t1+u1).mul(w).store(x+child+j);(t2+u2).mul(w2).store(x+2*child+j);
(t2-u2).mul(w3).store(x+3*child+j);(t1-u1).mul(w4).store(x+4*child+j);
}
for(int r=0;r<5;++r)FFTAVX::difRecLong<false>(x+r*child,child);
};
forward(a);forward(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)
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);
for(size_t j=0;j<child;j+=8) {
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{});
}
}
}}}
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=(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=5*std::max<size_t>(4096,hint::int_ceil2((2*std::max(na,nb)+4)/5));
hint::AlignMem<double>A(N),B(N);
// Both inputs occupy the first half; the radix-5 layer supplies the zero padding.
memset(A.begin()+na,0,(N/2-na)*sizeof(double));
memset(B.begin()+nb,0,(N/2-nb)*sizeof(double));
bigint_io_opt::parse_balanced_base100000(input,l1,A.begin());
bigint_io_opt::parse_balanced_base100000(input+off,l2,B.begin());
storage.reset();std::string().swap(fallback);
hint::transform::fft::real_conv5(A.begin(),B.begin(),N);
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 pos=format_parallel(coeff,conv,out,digits,std::min(na,nb)*999ULL+1);
size_t len=digits-pos;
if(direct) {
if(pos)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 | 927.051 ms | 831 MB + 244 KB | Accepted | Score: 100 | 显示更多 |