// ===== 引用区（本账号 saffah_v41f1017_agg4）=====
// [1] duck.ac 用户 saffah_cc_v41_agg1，提交 #130197 <https://duck.ac/submission/130197>
// （1005 实时 T=4.209 ms 榜首）。本文件正文 = 该公开提交的逐字节副本 + 下方所述编译指示；
// 其自带的完整上游引用链（saffah_codex_6a_agg3 #129909/#129699/#129675/#129689/#128944/
// #128715/#125838/#125697、HintFFT #19388、saffah_cc_v41_agg1 #129885 等账号与 URL、
// 复用内容说明）原样保留在下方。该公开提交无独立许可证声明；按站点公开条款具名引用。
// [2] duck.ac/problem/1005 题面：a 100,000 位、b 50,000 位，传统题 stdin→stdout。
// 合规：以上均为 duck.ac 公开提交；正文中以注释形式保留其全部上游引用。
// ===== 思路 =====
// 单变量试验：在底盘既有的 `#pragma GCC optimize("O3")` 与 `#pragma GCC target("fma")` 之后，
// 追加一条 `#pragma GCC optimize("align-functions=32","align-loops=32","rename-registers")`
// （纯编译参数，代码/算法/常数零改动）。依据：本队对同一 HintFFT 除法引擎在 1005e7 全尺度上的
// 判题机 A/B（align-functions=32 −3.3ms、align-loops=32 −5.8ms、rename-registers −10.9ms，
// 三者组合 −12.6ms ≈ −1.3%）；本题缺口为 ~1.0%。目的：以最小改动兑现该组合并取真实读数。
// 输出应与底盘逐字节相同。
// ======================
// ===== REFERENCES =====
// [1] duck.ac 账号 **saffah_codex_6a_agg3**，提交 **#129909** <https://duck.ac/submission/129909>
// （1005 实时 T = 4.287975 ms）—— **本发正文 = 该件的逐字节副本 + 下述三处改动**。
// 该件自带的全链引用块（含 #129699 / #129675 / #129689 / #128944 / #128715 / #125838 /
// #125697 / #19388(HintFFT) 等账号、URL 与复用内容说明）在下方**原样保留**（仅压缩散文，
// 账号 / URL / 复用内容行**一字未删**）。
// [2] 我方在役 **#129885** <https://duck.ac/submission/129885>（4.331376 ms）：r101a_ 的四枚
// 2-lane 融合核 + IO/打包刀 + r62c5_ 的 16 位 SWAR 与阈值 32 —— 均已含于 [1] 的血缘内。
// [3] 本发**三处新增**为本席原创（见思路段），未参考任何站外代码、论文或资料。
// 合规：以上均为 duck.ac 站内公开提交（题目页明示「你的代码将会被公开」）；按目标原文写明
// 「原作者账号 / 原始提交地址 / 所用内容」三项 ✓
// ================
// ===== 思路 =====（seat r105b_ · 1005 · 三处逐位等价的结构改动）
// 【结论】底盘上三处最大可动点，各改一处；三处**均逐位等价**（本地值闸门 stdout 逐字节相同 +
// 判题机探针的商哈希逐位相同）。
// ① [r105b-SWAR8]：abs_sub 的 SWAR 主环一轮吃的 limb 数 4 → 8（循环控制与分支减半）。
// ② [r105b-PF512]：fft_dif 蝶形环（4 条 rank 间隔流）加显式预取；距离已扫 32/128/256/512。
// ③ [r105b-EXTRACT]：试商乘法主环的元素提取由「4 次独立常量除法」改为「商链」（3 除 + 3 乘减），
// 输出四位逐位同值。
// 【判题机定价（本席静默探针，tryrun 通道，10 万位 / 5 万位真实规模，臂内极差 0.02%）】
// 底盘 #129909 = 1488.0 万拍 → 本发 = 1460.0 万拍 ⇒ **−1.89%**；同批单变量对照：① −0.99%、③ −0.77%。
// 【目的】本题为宽支（一发即翻严），靶 = 0.99T+0.001 ms = 4.246095 ms；本发为一次合并发射。
// ================
/** saffah_codex_6a_agg3 https://duck.ac/submission/129844 supplies the complete * saffah_cc_v41_agg1 https://duck.ac/submission/129885 supplies verbatim the */
/** saffah_codex_6a_agg3 https://duck.ac/submission/129832 supplies the accepted */
/** saffah_codex_6a_agg3 https://duck.ac/submission/129809 supplies the complete */
/** saffah_codex_6a_agg3 https://duck.ac/submission/129799 supplies the accepted * saffah_codex_6a_agg3 https://duck.ac/submission/129695 supplies the previously */
/** saffah_cc_v41_agg1 https://duck.ac/submission/129760 supplies the complete * saffah_codex_6a_agg3 https://duck.ac/submission/129643 and * https://duck.ac/submission/129678 supply proven normalized comparisons, */
// [1] 本发底盘 = 我方（saffah_cc_v41_agg1）在册最快件 **#129632** <https://duck.ac/submission/129632>
// #128715/#128944 既有移植注记、更上游作者链与许可说明）原位置原文字保留于下方，未删改一行。
// [2] 本发三把刀的出处 = 对手 **saffah_codex_6a_agg3** 的公开提交（本席逐行回件核过，非按 diff 猜测）：
// (a) #129699 <https://duck.ac/submission/129699>（1005 · 4.635989 ms = 本题现役 T）
// (b) #129675 <https://duck.ac/submission/129675>（1005 · 4.734 ms）
// (c) #129689 <https://duck.ac/submission/129689>（1005e7 · 1180.457756 ms = 该题现役 T）
// 许可合规：该站 FAQ 明示提交代码公开、无独立许可证 ⇒ 具名引用使用。
// [1] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#128968**
// <https://duck.ac/submission/128968>（1005 · 4.892945 ms = 我方本题在册最快件）
// [2] ★duck.ac 用户 **saffah_codex_6a_agg3**，提交 **#128944**
// <https://duck.ac/submission/128944>（1005 · 4.872047 ms = 本题榜面 T）。
// ② `Integer::string_in`：加 `__attribute__((target("avx2")))` + 32 字节 → 8 limb 的
// LICENSE：上述公开提交均**无独立许可证声明**、非来自任何外部仓库 ⇒ 按**具名引用**
// 【引用（逐条，全部为 duck.ac 公开提交，站点 FAQ 明示代码会被公开、无独立许可证 ⇒ 具名引用使用）】
// [A] 底盘 = 本账号 **saffah_cc_v41_agg1** 的在册最快件
// #128669 <https://duck.ac/submission/128669>（1005）
// #128672 <https://duck.ac/submission/128672>（1005e6）
// #128693 <https://duck.ac/submission/128693>（1005e7）
// （三件同源题族、正文逐字节同）。**沿用内容**：其正文全部（含文件头原有的全链引用与许可说明，
// [B] ★对手 **saffah_codex_6a_agg3** 提交 #128715 / #128716 / #128717
// <https://duck.ac/submission/128715> <https://duck.ac/submission/128716>
// <https://duck.ac/submission/128717>（三件逐字节同，本轮以 2.25% 反超我方 [A]）。
// [A] 底盘 = duck.ac 用户 **saffah_codex_6a_agg3**，提交 **#128206 / #128209 / #128210**
// <https://duck.ac/submission/128206>（三件逐字节同，md5 be9af3ae；1005 实时 T = 9.434486 ms）。
// 四刀 + 其 AVX2 平方环），其自带的全链引用与许可说明原样保留于下方。
// [B] 本账号 **saffah_cc_v41_agg1**：#128181 <https://duck.ac/submission/128181>（9.762335 ms，
// 许可：均为 duck.ac 公开提交（站点 FAQ：代码会被公开），无独立许可证声明 ⇒ 具名引用使用。
/** saffah_codex_6a_agg3, https://duck.ac/submission/128088, * https://duck.ac/submission/128089, https://duck.ac/submission/128090: * saffah_cc_v41_agg1, https://duck.ac/submission/128181, * https://duck.ac/submission/128184, https://duck.ac/submission/128189: * engine, with its documented inherited long-division corner limitation. */
/*- saffah_cc_v41_agg1 https://duck.ac/submission/127028, https://duck.ac/submission/127032 and https://duck.ac/submission/127038: twelve-limb long-division leaf threshold. All inherited author credits */
// [1] duck.ac 用户 saffah_codex_6a_agg3, 提交 #125838 <https://duck.ac/submission/125838>
// 用途：**本提交以此为逐字节底盘**（完整 HintFFT 分裂基 FFT + base-1000 递归除法
// [2] （保留自 [1] 的文件头）saffah_codex_maker, 提交 #125697 / #125698 / #125699
// <https://duck.ac/submission/125697>：base-1000 递归除法与 HintFFT 引擎原作者。
// [3] （保留自 [1] 的文件头）With-Sky / HintFFT, 提交 #19388
// <https://duck.ac/submission/19388>：上游 FFT 与整数运算；随件附带公开参考中已注明。
// [4] duck.ac 用户 saffah_cc_v41_agg1（本账号）, 提交 #126801 (1008) / #126784 (1011):
/*- saffah_codex_maker, https://duck.ac/submission/125697, https://duck.ac/submission/125698 and https://duck.ac/submission/125699: - With-Sky / HintFFT, https://duck.ac/submission/19388: upstream FFT and copied the complete base-1000 recursive division and HintFFT engine. integer arithmetic as attributed in the supplied public reference. */
// Adapted from With-Sky / HintFFT, Duck public submission 19388.
// Base 1000 limits floating convolution error at the new input lengths.
/*Known inherited limitation: the twelve-limb reference and this version both */
#include <algorithm>
#include <atomic>
#include <complex>
#include <future>
#include <iostream>
#include <random>
#include <string>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <immintrin.h>
#ifndef HINT_SIMD_HPP
#define HINT_SIMD_HPP
#pragma GCC optimize("O3")
#pragma GCC optimize("align-functions=32","align-loops=32","rename-registers")
#pragma GCC target("fma")
using Complex = std::complex<double>;
struct Complex2
{
__m256d data;
Complex2()
{
data = _mm256_setzero_pd();
}
Complex2(double input)
{
data = _mm256_set1_pd(input);
}
Complex2(__m256d input)
{
data = input;
}
Complex2(const Complex2 &input)
{
data = input.data;
}
Complex2(double const *ptr)
{
data = _mm256_loadu_pd(ptr);
}
Complex2(Complex a)
{
data = _mm256_broadcast_pd((__m128d *)&a);
}
Complex2(const Complex *ptr)
{
data = _mm256_loadu_pd((const double *)ptr);
}
void clr()
{
data = _mm256_setzero_pd();
}
void store(Complex *a) const
{
_mm256_storeu_pd((double *)a, data);
}
template <int M>
Complex2 element_permute() const
{
return _mm256_permute_pd(data, M);
}
Complex2 all_real() const
{
return _mm256_movedup_pd(data);
}
Complex2 all_imag() const
{
return element_permute<0XF>();
}
Complex2 swap() const
{
return element_permute<0X5>();
}
Complex2 mul_neg_i() const
{
static const __m256d subber{};
return Complex2(_mm256_addsub_pd(subber, data)).swap();
}
Complex2 operator+(Complex2 input) const
{
return _mm256_add_pd(data, input.data);
}
Complex2 operator-(Complex2 input) const
{
return _mm256_sub_pd(data, input.data);
}
Complex2 operator*(Complex2 input) const
{
auto imag = _mm256_mul_pd(all_imag().data, input.swap().data);
return _mm256_fmaddsub_pd(all_real().data, input.data, imag);
}
Complex2 operator/(Complex2 input) const
{
return _mm256_div_pd(data, input.data);
}
};
#endif
#include <cstdlib>
static std::complex<double> *g_fft_buf = nullptr;
static size_t g_fft_cap = 0;
namespace hint
{
using UINT_8 = uint8_t;
using UINT_16 = uint16_t;
using UINT_32 = uint32_t;
using UINT_64 = uint64_t;
using INT_32 = int32_t;
using INT_64 = int64_t;
using LONG = long;
using Complex = std::complex<double>;
constexpr double HINT_PI = 3.1415926535897932384626433832795;
constexpr double HINT_2PI = HINT_PI * 2;
template <typename T>
constexpr T min_2pow(T n)
{
if(n<=1)return T(1);
if constexpr(sizeof(T)<=4)
return T(1u << (32-__builtin_clz((unsigned)(n-1))));
else
return T(1ull << (64-__builtin_clzll((unsigned long long)(n-1))));
}
template <typename T>
inline void ary_clr(T *ptr, size_t len)
{
memset(ptr, 0, len * sizeof(T));
}
template <typename T>
constexpr bool is_neg(T x)
{
return x < 0;
}
template <typename T>
constexpr size_t hint_log2(T n)
{
T res = 0;
while (n > 1)
{
n /= 2;
res++;
}
return res;
}
template <typename T>
constexpr bool is_odd(T x)
{
return static_cast<bool>(x & 1);
}
template <typename T>
constexpr std::pair<T, T> div_mod(T a, T b)
{
return std::make_pair(a / b, a % b);
}
template <typename T>
void ary_copy(T *target, const T *source, size_t len)
{
if (len == 0 || target == source)
{
return;
}
if (len >= INT64_MAX)
{
throw("Ary too long\n");
}
std::memcpy(target, source, len * sizeof(T));
}
template <typename T>
inline T *ary_realloc(T *ptr, size_t len)
{
if (len * sizeof(T) < INT64_MAX)
{
ptr = static_cast<T *>(realloc(ptr, len * sizeof(T)));
}
if (ptr == nullptr)
{
throw("realloc error");
}
return ptr;
}
template <typename T>
constexpr T max_2pow(T n)
{
T res = 1;
res <<= (sizeof(T) * 8 - 1);
while (res > n)
{
res /= 2;
}
return res;
}
constexpr UINT_64 qpow(UINT_64 m, UINT_64 n, UINT_64 mod)
{
if (m <= 1)
{
return m;
}
UINT_64 result = 1;
while (n > 0)
{
if ((n & 1) != 0)
{
result = result * m % mod;
}
m = m * m % mod;
n >>= 1;
}
return result;
}
constexpr UINT_64 mod_inv(UINT_64 n, UINT_64 mod)
{
return qpow(n, mod - 2, mod);
}
template <typename T>
inline void com_ary_combine_copy(Complex *target, const T &source1, size_t len1, const T &source2, size_t len2)
{
size_t min_len = std::min(len1, len2);
size_t i = 0;
/* */
for (; i + 4 <= min_len; i += 4)
{
const __m128i va = _mm_loadu_si128((const __m128i *)(const void *)(&source1[i]));
const __m128i vb = _mm_loadu_si128((const __m128i *)(const void *)(&source2[i]));
const __m256d re = _mm256_cvtepi32_pd(va);
const __m256d im = _mm256_cvtepi32_pd(vb);
/* */
const __m256d u0 = _mm256_unpacklo_pd(re, im);
const __m256d u1 = _mm256_unpackhi_pd(re, im);
_mm256_storeu_pd((double *)(void *)(target + i), _mm256_permute2f128_pd(u0, u1, 0x20));
_mm256_storeu_pd((double *)(void *)(target + i + 2), _mm256_permute2f128_pd(u0, u1, 0x31));
}
while (i < min_len)
{
target[i] = Complex(source1[i], source2[i]);
i++;
}
const __m256d r40z = _mm256_setzero_pd();
for (; i + 4 <= len1; i += 4)
{
const __m128i va = _mm_loadu_si128((const __m128i *)(const void *)(&source1[i]));
const __m256d re = _mm256_cvtepi32_pd(va);
const __m256d u0 = _mm256_unpacklo_pd(re, r40z);
const __m256d u1 = _mm256_unpackhi_pd(re, r40z);
_mm256_storeu_pd((double *)(void *)(target + i), _mm256_permute2f128_pd(u0, u1, 0x20));
_mm256_storeu_pd((double *)(void *)(target + i + 2), _mm256_permute2f128_pd(u0, u1, 0x31));
}
while (i < len1)
{
target[i] = Complex(source1[i], 0.0);
i++;
}
for (; i + 4 <= len2; i += 4)
{
const __m128i vb = _mm_loadu_si128((const __m128i *)(const void *)(&source2[i]));
const __m256d im = _mm256_cvtepi32_pd(vb);
const __m256d u0 = _mm256_unpacklo_pd(r40z, im);
const __m256d u1 = _mm256_unpackhi_pd(r40z, im);
_mm256_storeu_pd((double *)(void *)(target + i), _mm256_permute2f128_pd(u0, u1, 0x20));
_mm256_storeu_pd((double *)(void *)(target + i + 2), _mm256_permute2f128_pd(u0, u1, 0x31));
}
while (i < len2)
{
target[i] = Complex(0.0, source2[i]);
i++;
}
}
namespace hint_transform
{
class ComplexTableY
{
private:
std::vector<std::vector<Complex>> table1;
std::vector<std::vector<Complex>> table3;
INT_32 max_log_size = 2;
INT_32 cur_log_size = 2;
static constexpr size_t FAC = 1;
ComplexTableY(const ComplexTableY &) = delete;
ComplexTableY &operator=(const ComplexTableY &) = delete;
public:
~ComplexTableY() {}
ComplexTableY(UINT_32 max_shift)
{
max_shift = std::max<size_t>(max_shift, 1);
max_log_size = max_shift;
table1.resize(max_shift + 1);
table3.resize(max_shift + 1);
table1[0] = table1[1] = table3[0] = table3[1] = std::vector<Complex>{1};
table1[2] = table3[2] = std::vector<Complex>{1};
#if TABLE_PRELOAD == 1
expand(max_shift);
#endif
}
void expand(INT_32 shift)
{
shift = std::max<INT_32>(shift, 2);
if (shift > max_log_size)
{
throw("FFT length too long for lut\n");
}
for (INT_32 i = cur_log_size + 1; i <= shift; i++)
{
size_t len = 1ull << i, vec_size = len * FAC / 4;
const double angle_scale=1.0/double(len);
table1[i].resize(vec_size);
table3[i].resize(vec_size);
table1[i][0] = table3[i][0] = Complex(1, 0);
for (size_t pos = 0; pos < vec_size / 2; pos++)
{
table1[i][pos * 2] = table1[i - 1][pos];
table3[i][pos * 2] = table3[i - 1][pos];
}
for (size_t pos = 1; pos < vec_size / 2; pos += 2)
{
double cos_theta = std::cos((HINT_2PI * pos) * angle_scale);
double sin_theta = std::sin((HINT_2PI * pos) * angle_scale);
table1[i][pos] = Complex(cos_theta, -sin_theta);
table1[i][vec_size - pos] = Complex(sin_theta, -cos_theta);
}
table1[i][vec_size / 2] = std::conj(unit_root(8, 1));
size_t odd=1;
for(;odd<vec_size/2 && odd*3<vec_size;odd+=2){
Complex tmp=table1[i][odd*3];
table3[i][odd]=tmp;
table3[i][vec_size-odd]=Complex(tmp.imag(),tmp.real());
}
for(;odd<vec_size/2;odd+=2){
Complex tmp=table1[i][len/2-odd*3];
tmp.real(-tmp.real());
table3[i][odd]=tmp;
table3[i][vec_size-odd]=Complex(tmp.imag(),tmp.real());
}
table3[i][vec_size / 2] = std::conj(unit_root(8, 3));
}
cur_log_size = std::max(cur_log_size, shift);
}
static Complex unit_root(double theta)
{
return std::polar<double>(1.0, theta);
}
static Complex unit_root(size_t m, size_t n)
{
return unit_root((HINT_2PI * n) / m);
}
Complex get_omega(UINT_32 shift, size_t n) const
{
size_t rank = 1ull << shift;
const size_t rank_ff = rank - 1, quad_n = n << 2;
size_t zone = quad_n >> shift; // 第几象限
if ((quad_n & rank_ff) == 0)
{
static constexpr Complex ONES[4] = {Complex(1, 0), Complex(0, -1), Complex(-1, 0), Complex(0, 1)};
return ONES[zone];
}
Complex tmp;
if ((zone & 2) == 0)
{
if ((zone & 1) == 0)
{
tmp = table1[shift][n];
}
else
{
tmp = table1[shift][rank / 2 - n];
tmp.real(-tmp.real());
}
}
else
{
if ((zone & 1) == 0)
{
tmp = -table1[shift][n - rank / 2];
}
else
{
tmp = table1[shift][rank - n];
tmp.imag(-tmp.imag());
}
}
return tmp;
}
Complex get_omega3(UINT_32 shift, size_t n) const
{
return table3[shift][n];
}
Complex2 get_omegaX2(UINT_32 shift, size_t n) const
{
return Complex2(table1[shift].data() + n);
}
Complex2 get_omega3X2(UINT_32 shift, size_t n) const
{
return Complex2(table3[shift].data() + n);
}
};
constexpr size_t lut_max_rank = 23;
static ComplexTableY TABLE(lut_max_rank);
template <typename T>
void binary_reverse_swap(T &ary, size_t len)
{
size_t i = 0;
for (size_t j = 1; j < len - 1; j++)
{
size_t k = len >> 1;
i ^= k;
while (k > i)
{
k >>= 1;
i ^= k;
};
if (j < i)
{
std::swap(ary[i], ary[j]);
}
}
}
template <typename SizeType = UINT_32, typename T>
void quaternary_reverse_swap(T &ary, size_t len)
{
SizeType log_n = hint_log2(len);
SizeType *rev = new SizeType[len / 4];
rev[0] = 0;
for (SizeType i = 1; i < len; i++)
{
SizeType index = (rev[i >> 2] >> 2) | ((i & 3) << (log_n - 2)); // 求rev交换数组
if (i < len / 4)
{
rev[i] = index;
}
if (i < index)
{
std::swap(ary[i], ary[index]);
}
}
delete[] rev;
}
template <typename T>
inline void fft_2point(T &sum, T &diff)
{
T tmp0 = sum;
T tmp1 = diff;
sum = tmp0 + tmp1;
diff = tmp0 - tmp1;
}
inline void fft_4point(Complex *input, size_t rank = 1)
{
Complex tmp0 = input[0];
Complex tmp1 = input[rank];
Complex tmp2 = input[rank * 2];
Complex tmp3 = input[rank * 3];
fft_2point(tmp0, tmp2);
fft_2point(tmp1, tmp3);
tmp3 = Complex(tmp3.imag(), -tmp3.real());
input[0] = tmp0 + tmp1;
input[rank] = tmp2 + tmp3;
input[rank * 2] = tmp0 - tmp1;
input[rank * 3] = tmp2 - tmp3;
}
inline void fft_dit_4point(Complex *input, size_t rank = 1)
{
Complex tmp0 = input[0];
Complex tmp1 = input[rank];
Complex tmp2 = input[rank * 2];
Complex tmp3 = input[rank * 3];
fft_2point(tmp0, tmp1);
fft_2point(tmp2, tmp3);
tmp3 = Complex(tmp3.imag(), -tmp3.real());
input[0] = tmp0 + tmp2;
input[rank] = tmp1 + tmp3;
input[rank * 2] = tmp0 - tmp2;
input[rank * 3] = tmp1 - tmp3;
}
inline void fft_dit_8point(Complex *input, size_t rank = 1)
{
Complex tmp0 = input[0];
Complex tmp1 = input[rank];
Complex tmp2 = input[rank * 2];
Complex tmp3 = input[rank * 3];
Complex tmp4 = input[rank * 4];
Complex tmp5 = input[rank * 5];
Complex tmp6 = input[rank * 6];
Complex tmp7 = input[rank * 7];
fft_2point(tmp0, tmp1);
fft_2point(tmp2, tmp3);
fft_2point(tmp4, tmp5);
fft_2point(tmp6, tmp7);
tmp3 = Complex(tmp3.imag(), -tmp3.real());
tmp7 = Complex(tmp7.imag(), -tmp7.real());
fft_2point(tmp0, tmp2);
fft_2point(tmp1, tmp3);
fft_2point(tmp4, tmp6);
fft_2point(tmp5, tmp7);
static constexpr double cos_1_8 = 0.70710678118654752440084436210485;
tmp5 = cos_1_8 * Complex(tmp5.imag() + tmp5.real(), tmp5.imag() - tmp5.real());
tmp6 = Complex(tmp6.imag(), -tmp6.real());
tmp7 = -cos_1_8 * Complex(tmp7.real() - tmp7.imag(), tmp7.real() + tmp7.imag());
input[0] = tmp0 + tmp4;
input[rank] = tmp1 + tmp5;
input[rank * 2] = tmp2 + tmp6;
input[rank * 3] = tmp3 + tmp7;
input[rank * 4] = tmp0 - tmp4;
input[rank * 5] = tmp1 - tmp5;
input[rank * 6] = tmp2 - tmp6;
input[rank * 7] = tmp3 - tmp7;
}
inline void fft_dif_4point(Complex *input, size_t rank = 1)
{
Complex tmp0 = input[0];
Complex tmp1 = input[rank];
Complex tmp2 = input[rank * 2];
Complex tmp3 = input[rank * 3];
fft_2point(tmp0, tmp2);
fft_2point(tmp1, tmp3);
tmp3 = Complex(tmp3.imag(), -tmp3.real());
input[0] = tmp0 + tmp1;
input[rank] = tmp0 - tmp1;
input[rank * 2] = tmp2 + tmp3;
input[rank * 3] = tmp2 - tmp3;
}
inline void fft_dif_8point(Complex *input, size_t rank = 1)
{
Complex tmp0 = input[0];
Complex tmp1 = input[rank];
Complex tmp2 = input[rank * 2];
Complex tmp3 = input[rank * 3];
Complex tmp4 = input[rank * 4];
Complex tmp5 = input[rank * 5];
Complex tmp6 = input[rank * 6];
Complex tmp7 = input[rank * 7];
fft_2point(tmp0, tmp4);
fft_2point(tmp1, tmp5);
fft_2point(tmp2, tmp6);
fft_2point(tmp3, tmp7);
static constexpr double cos_1_8 = 0.70710678118654752440084436210485;
tmp5 = cos_1_8 * Complex(tmp5.imag() + tmp5.real(), tmp5.imag() - tmp5.real());
tmp6 = Complex(tmp6.imag(), -tmp6.real());
tmp7 = -cos_1_8 * Complex(tmp7.real() - tmp7.imag(), tmp7.real() + tmp7.imag());
fft_2point(tmp0, tmp2);
fft_2point(tmp1, tmp3);
fft_2point(tmp4, tmp6);
fft_2point(tmp5, tmp7);
tmp3 = Complex(tmp3.imag(), -tmp3.real());
tmp7 = Complex(tmp7.imag(), -tmp7.real());
input[0] = tmp0 + tmp1;
input[rank] = tmp0 - tmp1;
input[rank * 2] = tmp2 + tmp3;
input[rank * 3] = tmp2 - tmp3;
input[rank * 4] = tmp4 + tmp5;
input[rank * 5] = tmp4 - tmp5;
input[rank * 6] = tmp6 + tmp7;
input[rank * 7] = tmp6 - tmp7;
}
inline void fft_radix2_dit_butterfly(Complex omega, Complex *input, size_t rank)
{
Complex tmp0 = input[0];
Complex tmp1 = input[rank] * omega;
input[0] = tmp0 + tmp1;
input[rank] = tmp0 - tmp1;
}
inline void fft_radix2_dif_butterfly(Complex omega, Complex *input, size_t rank)
{
Complex tmp0 = input[0];
Complex tmp1 = input[rank];
input[0] = tmp0 + tmp1;
input[rank] = (tmp0 - tmp1) * omega;
}
inline void fft_split_radix_dit_butterfly(Complex omega, Complex omega_cube,
Complex *input, size_t rank)
{
Complex tmp0 = input[0];
Complex tmp1 = input[rank];
Complex tmp2 = input[rank * 2] * omega;
Complex tmp3 = input[rank * 3] * omega_cube;
fft_2point(tmp2, tmp3);
tmp3 = Complex(tmp3.imag(), -tmp3.real());
input[0] = tmp0 + tmp2;
input[rank] = tmp1 + tmp3;
input[rank * 2] = tmp0 - tmp2;
input[rank * 3] = tmp1 - tmp3;
}
inline void fft_split_radix_dif_butterfly(Complex omega, Complex omega_cube,
Complex *input, size_t rank)
{
Complex tmp0 = input[0];
Complex tmp1 = input[rank];
Complex tmp2 = input[rank * 2];
Complex tmp3 = input[rank * 3];
fft_2point(tmp0, tmp2);
fft_2point(tmp1, tmp3);
tmp3 = Complex(tmp3.imag(), -tmp3.real());
input[0] = tmp0;
input[rank] = tmp1;
input[rank * 2] = (tmp2 + tmp3) * omega;
input[rank * 3] = (tmp2 - tmp3) * omega_cube;
}
inline void fft_split_radix_dit_butterfly(Complex2 omega, Complex2 omega_cube,
Complex *input, size_t rank)
{
Complex2 tmp0 = input;
Complex2 tmp1 = input + rank;
Complex2 tmp2 = Complex2(input + rank * 2) * omega;
Complex2 tmp3 = Complex2(input + rank * 3) * omega_cube;
fft_2point(tmp2, tmp3);
tmp3 = tmp3.mul_neg_i();
(tmp0 + tmp2).store(input);
(tmp1 + tmp3).store(input + rank);
(tmp0 - tmp2).store(input + rank * 2);
(tmp1 - tmp3).store(input + rank * 3);
}
inline void fft_split_radix_dif_butterfly(Complex2 omega, Complex2 omega_cube,
Complex *input, size_t rank)
{
Complex2 tmp0 = (input);
Complex2 tmp1 = (input + rank);
Complex2 tmp2 = (input + rank * 2);
Complex2 tmp3 = (input + rank * 3);
fft_2point(tmp0, tmp2);
fft_2point(tmp1, tmp3);
tmp3 = tmp3.mul_neg_i();
tmp0.store(input);
tmp1.store(input + rank);
((tmp2 + tmp3) * omega).store(input + rank * 2);
((tmp2 - tmp3) * omega_cube).store(input + rank * 3);
}
inline void fft_radix4_dit_butterfly(Complex omega, Complex omega_sqr, Complex omega_cube,
Complex *input, size_t rank)
{
Complex tmp0 = input[0];
Complex tmp1 = input[rank] * omega;
Complex tmp2 = input[rank * 2] * omega_sqr;
Complex tmp3 = input[rank * 3] * omega_cube;
fft_2point(tmp0, tmp2);
fft_2point(tmp1, tmp3);
tmp3 = Complex(tmp3.imag(), -tmp3.real());
input[0] = tmp0 + tmp1;
input[rank] = tmp2 + tmp3;
input[rank * 2] = tmp0 - tmp1;
input[rank * 3] = tmp2 - tmp3;
}
inline void fft_radix4_dif_butterfly(Complex omega, Complex omega_sqr, Complex omega_cube,
Complex *input, size_t rank)
{
Complex tmp0 = input[0];
Complex tmp1 = input[rank];
Complex tmp2 = input[rank * 2];
Complex tmp3 = input[rank * 3];
fft_2point(tmp0, tmp2);
fft_2point(tmp1, tmp3);
tmp3 = Complex(tmp3.imag(), -tmp3.real());
input[0] = tmp0 + tmp1;
input[rank] = (tmp2 + tmp3) * omega;
input[rank * 2] = (tmp0 - tmp1) * omega_sqr;
input[rank * 3] = (tmp2 - tmp3) * omega_cube;
}
inline void fft_conj(Complex *input, size_t fft_len, double div = 1)
{
for (size_t i = 0; i < fft_len; i++)
{
input[i] = std::conj(input[i]) / div;
}
}
inline void fft_normalize(Complex *input, size_t fft_len)
{
double len = static_cast<double>(fft_len);
for (size_t i = 0; i < fft_len; i++)
{
input[i] /= len;
}
}
/* */
#define R101A_LD(K, N) \
__m256d v##K = _mm256_insertf128_pd(_mm256_castpd128_pd256(_mm_loadu_pd((const double *)(in + (K)))), \
_mm_loadu_pd((const double *)(in + (N) + (K))), 1);
#define R101A_ST(R, K, N) \
_mm_storeu_pd((double *)(in + (K)), _mm256_castpd256_pd128(R)); \
_mm_storeu_pd((double *)(in + (N) + (K)), _mm256_extractf128_pd(R, 1));
#define R101A_SIGNIM _mm256_castsi256_pd(_mm256_setr_epi64x(0, (long long)0x8000000000000000LL, 0, (long long)0x8000000000000000LL))
#define R101A_C18 _mm256_set1_pd(0.70710678118654752440084436210485)
#define R101A_MC18 _mm256_set1_pd(-0.70710678118654752440084436210485)
__attribute__((target("avx2"))) inline void fft_dif_8point_x2(Complex *in)
{
R101A_LD(0, 8) R101A_LD(1, 8) R101A_LD(2, 8) R101A_LD(3, 8)
R101A_LD(4, 8) R101A_LD(5, 8) R101A_LD(6, 8) R101A_LD(7, 8)
__m256d s04 = _mm256_add_pd(v0, v4), d04 = _mm256_sub_pd(v0, v4);
__m256d s15 = _mm256_add_pd(v1, v5), d15 = _mm256_sub_pd(v1, v5);
__m256d s26 = _mm256_add_pd(v2, v6), d26 = _mm256_sub_pd(v2, v6);
__m256d s37 = _mm256_add_pd(v3, v7), d37 = _mm256_sub_pd(v3, v7);
const __m256d sg = R101A_SIGNIM;
{ __m256d v = d15; /* */
d15 = _mm256_mul_pd(R101A_C18, _mm256_permute_pd(_mm256_addsub_pd(_mm256_permute_pd(v, 0x5), v), 0x5)); }
d26 = _mm256_xor_pd(_mm256_permute_pd(d26, 0x5), sg); /* */
{ __m256d v = d37; /* */
d37 = _mm256_mul_pd(R101A_MC18, _mm256_addsub_pd(v, _mm256_permute_pd(v, 0x5))); }
__m256d u0 = _mm256_add_pd(s04, s26), u2 = _mm256_sub_pd(s04, s26);
__m256d u1 = _mm256_add_pd(s15, s37), u3 = _mm256_sub_pd(s15, s37);
__m256d u4 = _mm256_add_pd(d04, d26), u6 = _mm256_sub_pd(d04, d26);
__m256d u5 = _mm256_add_pd(d15, d37), u7 = _mm256_sub_pd(d15, d37);
u3 = _mm256_xor_pd(_mm256_permute_pd(u3, 0x5), sg);
u7 = _mm256_xor_pd(_mm256_permute_pd(u7, 0x5), sg);
R101A_ST(_mm256_add_pd(u0, u1), 0, 8) R101A_ST(_mm256_sub_pd(u0, u1), 1, 8)
R101A_ST(_mm256_add_pd(u2, u3), 2, 8) R101A_ST(_mm256_sub_pd(u2, u3), 3, 8)
R101A_ST(_mm256_add_pd(u4, u5), 4, 8) R101A_ST(_mm256_sub_pd(u4, u5), 5, 8)
R101A_ST(_mm256_add_pd(u6, u7), 6, 8) R101A_ST(_mm256_sub_pd(u6, u7), 7, 8)
}
__attribute__((target("avx2"))) inline void fft_dif_4point_x2(Complex *in)
{
R101A_LD(0, 4) R101A_LD(1, 4) R101A_LD(2, 4) R101A_LD(3, 4)
__m256d s02 = _mm256_add_pd(v0, v2), d02 = _mm256_sub_pd(v0, v2);
__m256d s13 = _mm256_add_pd(v1, v3), d13 = _mm256_sub_pd(v1, v3);
d13 = _mm256_xor_pd(_mm256_permute_pd(d13, 0x5), R101A_SIGNIM);
R101A_ST(_mm256_add_pd(s02, s13), 0, 4) R101A_ST(_mm256_sub_pd(s02, s13), 1, 4)
R101A_ST(_mm256_add_pd(d02, d13), 2, 4) R101A_ST(_mm256_sub_pd(d02, d13), 3, 4)
}
__attribute__((target("avx2"))) inline void fft_dit_8point_x2(Complex *in)
{
R101A_LD(0, 8) R101A_LD(1, 8) R101A_LD(2, 8) R101A_LD(3, 8)
R101A_LD(4, 8) R101A_LD(5, 8) R101A_LD(6, 8) R101A_LD(7, 8)
__m256d s01 = _mm256_add_pd(v0, v1), d01 = _mm256_sub_pd(v0, v1);
__m256d s23 = _mm256_add_pd(v2, v3), d23 = _mm256_sub_pd(v2, v3);
__m256d s45 = _mm256_add_pd(v4, v5), d45 = _mm256_sub_pd(v4, v5);
__m256d s67 = _mm256_add_pd(v6, v7), d67 = _mm256_sub_pd(v6, v7);
const __m256d sg = R101A_SIGNIM;
d23 = _mm256_xor_pd(_mm256_permute_pd(d23, 0x5), sg);
d67 = _mm256_xor_pd(_mm256_permute_pd(d67, 0x5), sg);
__m256d w0 = _mm256_add_pd(s01, s23), w2 = _mm256_sub_pd(s01, s23);
__m256d w1 = _mm256_add_pd(d01, d23), w3 = _mm256_sub_pd(d01, d23);
__m256d w4 = _mm256_add_pd(s45, s67), w6 = _mm256_sub_pd(s45, s67);
__m256d w5 = _mm256_add_pd(d45, d67), w7 = _mm256_sub_pd(d45, d67);
{ __m256d v = w5;
w5 = _mm256_mul_pd(R101A_C18, _mm256_permute_pd(_mm256_addsub_pd(_mm256_permute_pd(v, 0x5), v), 0x5)); }
w6 = _mm256_xor_pd(_mm256_permute_pd(w6, 0x5), sg);
{ __m256d v = w7;
w7 = _mm256_mul_pd(R101A_MC18, _mm256_addsub_pd(v, _mm256_permute_pd(v, 0x5))); }
R101A_ST(_mm256_add_pd(w0, w4), 0, 8) R101A_ST(_mm256_add_pd(w1, w5), 1, 8)
R101A_ST(_mm256_add_pd(w2, w6), 2, 8) R101A_ST(_mm256_add_pd(w3, w7), 3, 8)
R101A_ST(_mm256_sub_pd(w0, w4), 4, 8) R101A_ST(_mm256_sub_pd(w1, w5), 5, 8)
R101A_ST(_mm256_sub_pd(w2, w6), 6, 8) R101A_ST(_mm256_sub_pd(w3, w7), 7, 8)
}
__attribute__((target("avx2"))) inline void fft_dit_4point_x2(Complex *in)
{
R101A_LD(0, 4) R101A_LD(1, 4) R101A_LD(2, 4) R101A_LD(3, 4)
__m256d s01 = _mm256_add_pd(v0, v1), d01 = _mm256_sub_pd(v0, v1);
__m256d s23 = _mm256_add_pd(v2, v3), d23 = _mm256_sub_pd(v2, v3);
d23 = _mm256_xor_pd(_mm256_permute_pd(d23, 0x5), R101A_SIGNIM);
R101A_ST(_mm256_add_pd(s01, s23), 0, 4) R101A_ST(_mm256_add_pd(d01, d23), 1, 4)
R101A_ST(_mm256_sub_pd(s01, s23), 2, 4) R101A_ST(_mm256_sub_pd(d01, d23), 3, 4)
}
template <size_t LEN>
void fft_split_radix_dit_template(Complex *input)
{
constexpr size_t log_len = hint_log2(LEN);
constexpr size_t half_len = LEN / 2, quarter_len = LEN / 4;
fft_split_radix_dit_template<half_len>(input);
if constexpr (quarter_len == 8) fft_dit_8point_x2(input + half_len);
else if constexpr (quarter_len == 4) fft_dit_4point_x2(input + half_len);
else {
fft_split_radix_dit_template<quarter_len>(input + half_len);
fft_split_radix_dit_template<quarter_len>(input + half_len + quarter_len);
}
for (size_t i = 0; i < quarter_len; i += 2)
{
auto omega = TABLE.get_omegaX2(log_len, i);
auto omega_cube = TABLE.get_omega3X2(log_len, i);
fft_split_radix_dit_butterfly(omega, omega_cube, input + i, quarter_len);
}
}
template <>
void fft_split_radix_dit_template<0>(Complex *input) {}
template <>
void fft_split_radix_dit_template<1>(Complex *input) {}
template <>
void fft_split_radix_dit_template<2>(Complex *input)
{
fft_2point(input[0], input[1]);
}
template <>
void fft_split_radix_dit_template<4>(Complex *input)
{
fft_dit_4point(input, 1);
}
template <>
void fft_split_radix_dit_template<8>(Complex *input)
{
fft_dit_8point(input, 1);
}
template <size_t LEN>
void fft_split_radix_dif_template(Complex *input)
{
constexpr size_t log_len = hint_log2(LEN);
constexpr size_t half_len = LEN / 2, quarter_len = LEN / 4;
for (size_t i = 0; i < quarter_len; i += 2)
{
if (i + 512 < quarter_len)
{
const Complex *pf = input + i + 512;
__builtin_prefetch(pf, 1, 3);
__builtin_prefetch(pf + quarter_len, 1, 3);
__builtin_prefetch(pf + quarter_len * 2, 1, 3);
__builtin_prefetch(pf + quarter_len * 3, 1, 3);
}
auto omega = TABLE.get_omegaX2(log_len, i);
auto omega_cube = TABLE.get_omega3X2(log_len, i);
fft_split_radix_dif_butterfly(omega, omega_cube, input + i, quarter_len);
}
fft_split_radix_dif_template<half_len>(input);
if constexpr (quarter_len == 8) fft_dif_8point_x2(input + half_len);
else if constexpr (quarter_len == 4) fft_dif_4point_x2(input + half_len);
else {
fft_split_radix_dif_template<quarter_len>(input + half_len);
fft_split_radix_dif_template<quarter_len>(input + half_len + quarter_len);
}
}
template <>
void fft_split_radix_dif_template<0>(Complex *input) {}
template <>
void fft_split_radix_dif_template<1>(Complex *input) {}
template <>
void fft_split_radix_dif_template<2>(Complex *input)
{
fft_2point(input[0], input[1]);
}
template <>
void fft_split_radix_dif_template<4>(Complex *input)
{
fft_dif_4point(input, 1);
}
template <>
void fft_split_radix_dif_template<8>(Complex *input)
{
fft_dif_8point(input, 1);
}
template <size_t LEN = 1>
void fft_dit_template(Complex *input, size_t fft_len)
{
if (fft_len > LEN)
{
fft_dit_template<LEN * 2>(input, fft_len);
return;
}
TABLE.expand(hint_log2(LEN));
fft_split_radix_dit_template<LEN>(input);
}
template <>
void fft_dit_template<1 << 24>(Complex *input, size_t fft_len) {}
template <size_t LEN = 1>
void fft_dif_template(Complex *input, size_t fft_len)
{
if (fft_len > LEN)
{
fft_dif_template<LEN * 2>(input, fft_len);
return;
}
TABLE.expand(hint_log2(LEN));
fft_split_radix_dif_template<LEN>(input);
}
template <>
void fft_dif_template<1 << 24>(Complex *input, size_t fft_len) {}
/*<https://duck.ac/submission/129689>（逐行回件核，非按 diff 猜测）： [r101a-ZH1] 反向移植 saffah_codex_6a_agg3 #129689 */
template <size_t LEN = 16>
void fft_dif_zero_half(Complex *input, size_t fft_len)
{
if constexpr (LEN < (1u << 24))
{
if (fft_len > LEN)
{
fft_dif_zero_half<LEN * 2>(input, fft_len);
return;
}
constexpr size_t half = LEN / 2, quarter = LEN / 4, lg = hint_log2(LEN);
TABLE.expand(lg);
for (size_t i = 0; i < quarter; i += 2)
{
const Complex2 a(input + i), b(input + quarter + i);
const Complex2 bi = b.mul_neg_i();
((a + bi) * TABLE.get_omegaX2(lg, i)).store(input + half + i);
((a - bi) * TABLE.get_omega3X2(lg, i)).store(input + half + quarter + i);
}
fft_split_radix_dif_template<half>(input);
if constexpr (quarter == 8) fft_dif_8point_x2(input + half);
else if constexpr (quarter == 4) fft_dif_4point_x2(input + half);
else {
fft_split_radix_dif_template<quarter>(input + half);
fft_split_radix_dif_template<quarter>(input + half + quarter);
}
}
}
inline void fft_dit(Complex *input, size_t fft_len, bool bit_rev = true)
{
fft_len = max_2pow(fft_len);
if (bit_rev)
{
binary_reverse_swap(input, fft_len);
}
fft_dit_template<1>(input, fft_len);
}
inline void fft_dif(Complex *input, size_t fft_len, bool bit_rev = true)
{
fft_len = max_2pow(fft_len);
fft_dif_template<1>(input, fft_len);
if (bit_rev)
{
binary_reverse_swap(input, fft_len);
}
}
}
template <UINT_64 BASE, typename T>
constexpr void save_num_to_ary(T ary[], UINT_64 num, size_t digit)
{
size_t i = digit;
while (i > 0)
{
i--;
ary[i] = num % BASE;
num /= BASE;
}
}
UINT_64 stoui64(const std::string::const_iterator &begin,
const std::string::const_iterator &end, const UINT_32 base = 10)
{
UINT_64 result = 0;
for (auto it = begin; it < end && it - begin < 19; ++it)
{
result *= base;
char c = tolower(*it);
UINT_64 n = 0;
if (isalnum(c))
{
if (isalpha(c))
{
n = c - 'a' + 10;
}
else
{
n = c - '0';
}
}
if (n < base)
{
result += n;
}
}
return result;
}
UINT_64 stoui64(const std::string &str, const UINT_32 base = 10)
{
return stoui64(str.begin(), str.end(), base);
}
std::string ui64to_string(UINT_64 input, UINT_8 digits)
{
std::string result(digits, '0');
for (UINT_8 i = 0; i < digits; i++)
{
result[digits - i - 1] = static_cast<char>(input % 10 + '0');
input /= 10;
}
return result;
}
template <typename T = UINT_32, typename SIZE_TYPE = UINT_32>
class HintVector
{
private:
T *ary_ptr = nullptr;
SIZE_TYPE sign_n_len = 0;
SIZE_TYPE size = 0;
public:
~HintVector()
{
if (ary_ptr != nullptr)
{
delete[] ary_ptr;
ary_ptr = nullptr;
}
}
HintVector()
{
ary_ptr = nullptr;
}
HintVector(SIZE_TYPE ele_size)
{
if (ele_size > 0)
{
resize(ele_size);
}
}
HintVector(SIZE_TYPE ele_size, T ele)
{
if (ele_size > 0)
{
resize(ele_size);
change_length(ele_size);
fill_element(ele, ele_size);
}
}
HintVector(const T *ary, SIZE_TYPE len)
{
if (ary == nullptr)
{
throw("Can't copy from nullptr\n");
}
if (len > 0)
{
resize(len);
change_length(len);
ary_copy(ary_ptr, ary, len);
}
}
HintVector(const HintVector &input)
{
if (input.ary_ptr != nullptr)
{
if (ary_ptr != nullptr)
{
delete[] ary_ptr;
}
ary_ptr = nullptr;
size_t len = input.length();
resize(len);
change_length(len);
change_sign(input.sign());
ary_copy(ary_ptr, input.ary_ptr, len);
}
}
HintVector(HintVector &&input)
{
if (input.ary_ptr != nullptr)
{
size = input.size;
change_length(input.length());
change_sign(input.sign());
if (ary_ptr != nullptr)
{
delete[] ary_ptr;
}
ary_ptr = input.ary_ptr;
input.ary_ptr = nullptr;
}
}
HintVector &operator=(const HintVector &input)
{
if (this != &input)
{
if (ary_ptr != nullptr)
{
delete[] ary_ptr;
}
ary_ptr = nullptr;
size_t len = input.length();
resize(len);
change_length(len);
change_sign(input.sign());
ary_copy(ary_ptr, input.ary_ptr, len);
}
return *this;
}
HintVector &operator=(HintVector &&input)
{
if (this != &input)
{
size = input.size;
change_length(input.length());
change_sign(input.sign());
if (ary_ptr != nullptr)
{
delete[] ary_ptr;
}
ary_ptr = input.ary_ptr;
input.ary_ptr = nullptr;
}
return *this;
}
T &operator[](SIZE_TYPE index) const
{
return ary_ptr[index];
}
T *data() const
{
return ary_ptr;
}
static SIZE_TYPE size_generator(SIZE_TYPE new_size)
{
constexpr SIZE_TYPE SIZE_TYPE_BITS = sizeof(SIZE_TYPE) * 8;
constexpr SIZE_TYPE SIZE_80 = (1ull << (SIZE_TYPE_BITS - 1));
constexpr SIZE_TYPE LEN_MAX = SIZE_80 - 1;
new_size = std::min<SIZE_TYPE>(new_size, LEN_MAX);
if (new_size <= 2)
{
return 2;
}
else if (new_size <= 4)
{
return 4;
}
SIZE_TYPE size1 = min_2pow(new_size);
SIZE_TYPE size2 = size1 / 2;
size2 = size2 + size2 / 2;
if (new_size <= size2)
{
return size2;
}
else
{
return size1;
}
}
SIZE_TYPE length() const
{
constexpr SIZE_TYPE SIZE_TYPE_BITS = sizeof(SIZE_TYPE) * 8;
constexpr SIZE_TYPE SIZE_80 = (1ull << (SIZE_TYPE_BITS - 1));
constexpr SIZE_TYPE LEN_MAX = SIZE_80 - 1;
return sign_n_len & LEN_MAX;
}
bool sign() const
{
constexpr SIZE_TYPE SIZE_TYPE_BITS = sizeof(SIZE_TYPE) * 8;
constexpr SIZE_TYPE SIZE_80 = (1ull << (SIZE_TYPE_BITS - 1));
return (SIZE_80 & sign_n_len) != 0;
}
/*[r101a-PACK1] 反向移植 saffah_codex_6a_agg3 #129699：容量查询（其 #129699 新增的 3 行）。 */
SIZE_TYPE capacity() const
{
return size;
}
constexpr void resize(SIZE_TYPE new_size, SIZE_TYPE(size_func)(SIZE_TYPE) = size_generator)
{
new_size = size_func(new_size);
if (ary_ptr == nullptr)
{
size = new_size;
ary_ptr = new T[size];
change_length(0);
}
else if (size != new_size)
{
size = new_size;
ary_ptr = ary_realloc(ary_ptr, size);
change_length(length());
}
}
void change_length(SIZE_TYPE new_length)
{
constexpr SIZE_TYPE SIZE_TYPE_BITS = sizeof(SIZE_TYPE) * 8;
constexpr SIZE_TYPE SIZE_80 = (1ull << (SIZE_TYPE_BITS - 1));
constexpr SIZE_TYPE LEN_MAX = SIZE_80 - 1;
if (new_length > LEN_MAX)
{
throw("Length too long\n");
}
new_length = std::min(new_length, size);
bool sign_tmp = sign();
sign_n_len = new_length;
change_sign(sign_tmp);
}
void change_sign(bool is_sign)
{
constexpr SIZE_TYPE SIZE_TYPE_BITS = sizeof(SIZE_TYPE) * 8;
constexpr SIZE_TYPE SIZE_80 = (1ull << (SIZE_TYPE_BITS - 1));
constexpr SIZE_TYPE LEN_MAX = SIZE_80 - 1;
if ((!is_sign) || length() == 0)
{
sign_n_len = sign_n_len & LEN_MAX;
}
else
{
sign_n_len = sign_n_len | SIZE_80;
}
}
SIZE_TYPE set_true_len()
{
if (ary_ptr == nullptr)
{
return 0;
}
SIZE_TYPE t_len = length();
while (t_len > 0 && ary_ptr[t_len - 1] == 0)
{
t_len--;
}
change_length(t_len);
return length();
}
void fill_element(T ele, SIZE_TYPE count, SIZE_TYPE begin = 0)
{
if (begin >= size || ary_ptr == nullptr)
{
return;
}
if (begin + count >= size)
{
count = size - begin;
}
std::fill(ary_ptr + begin, ary_ptr + begin + count, ele);
}
};
}
namespace hint_arithm
{
using namespace hint;
using SIZE_TYPE = UINT_32;
template <typename T>
using hintvector = HintVector<T, SIZE_TYPE>;
template <typename T>
constexpr size_t ary_true_len(const T &ary, size_t len)
{
while (len > 0 && ary[len - 1] == 0)
{
len--;
}
return len;
}
template <typename T>
constexpr INT_32 abs_compare(const T ary1[], const T ary2[], size_t len1, size_t len2)
{
len1 = ary_true_len(ary1, len1);
len2 = ary_true_len(ary2, len2);
if (len1 != len2)
{
return len1 > len2 ? 1 : -1;
}
if (ary1 == ary2)
{
return 0;
}
while (len1 > 0)
{
len1--;
T num1 = ary1[len1], num2 = ary2[len1];
if (num1 != num2)
{
return num1 > num2 ? 1 : -1;
}
}
return 0;
}
template <typename T>
constexpr INT_32 abs_compare_normalized(const T ary1[], const T ary2[], size_t len1, size_t len2)
{
if (len1 != len2)
{
return len1 > len2 ? 1 : -1;
}
if (ary1 == ary2)
{
return 0;
}
while (len1 > 0)
{
len1--;
T num1 = ary1[len1], num2 = ary2[len1];
if (num1 != num2)
{
return num1 > num2 ? 1 : -1;
}
}
return 0;
}
template <UINT_64 BASE, bool is_carry = true, typename T>
constexpr void abs_mul_num(const T in[], T num, T out[], size_t len)
{
len = ary_true_len(in, len);
num %= BASE;
UINT_64 prod = 0;
if (num == 1)
{
ary_copy(out, in, len);
}
else
{
size_t i = 0;
if constexpr (BASE <= 10000)
{
/* */
constexpr UINT_64 B2 = BASE * BASE, B3 = B2 * BASE;
for (; i + 3 <= len; i += 3)
{
const UINT_64 packed = UINT_64(in[i]) + BASE * in[i + 1] + B2 * in[i + 2];
const UINT_64 value = packed * num + prod;
out[i] = T(value % BASE);
out[i + 1] = T((value / BASE) % BASE);
out[i + 2] = T((value / B2) % BASE);
prod = value / B3;
}
}
for (; i < len; i++)
{
prod += static_cast<UINT_64>(in[i]) * num;
std::tie(prod, out[i]) = div_mod<UINT_64>(prod, BASE);
}
}
if (is_carry)
{
out[len] = prod;
}
}
template <UINT_64 BASE, typename T>
constexpr INT_64 abs_div_num(const T in[], T num, T out[], size_t len)
{
size_t pos = ary_true_len(in, len);
num %= BASE;
if (num == 1)
{
ary_copy(out, in, len);
return 0;
}
UINT_64 last_rem = 0;
if constexpr (BASE <= 10000)
{
/* */
constexpr UINT_64 B2 = BASE * BASE, B3 = B2 * BASE;
while (pos >= 3)
{
pos -= 3;
const UINT_64 value = last_rem * B3 + UINT_64(in[pos + 2]) * B2
+ UINT_64(in[pos + 1]) * BASE + in[pos];
const UINT_64 q = value / num;
last_rem = value - q * num;
out[pos] = T(q % BASE);
out[pos + 1] = T((q / BASE) % BASE);
out[pos + 2] = T(q / B2);
}
}
while (pos > 0)
{
pos--;
last_rem = last_rem * BASE + in[pos];
std::tie(out[pos], last_rem) = div_mod<UINT_64>(last_rem, num);
}
return last_rem;
}
template <INT_64 BASE, bool is_carry = true, typename T>
constexpr void abs_add(const T in1[], const T in2[], T out[],
size_t len1, size_t len2)
{
if (len1 < len2)
{
std::swap(in1, in2);
std::swap(len1, len2);
}
size_t pos = 0;
UINT_64 carry = 0;
if constexpr (sizeof(T) == 4 && BASE <= 32768)
{
/* */
constexpr unsigned long long bias = 0x10000ull - (unsigned long long)BASE;
const unsigned long long biasq = bias * 0x0001000100010001ull;
const __m128i ones = _mm_set1_epi32(-1);
while (pos + 4 <= len2)
{
const __m128i a4 = _mm_loadu_si128((const __m128i *)(const void *)(in1 + pos));
const __m128i b4 = _mm_loadu_si128((const __m128i *)(const void *)(in2 + pos));
const unsigned long long va = (unsigned long long)_mm_cvtsi128_si64(_mm_packus_epi32(a4, ones));
const unsigned long long vb = (unsigned long long)_mm_cvtsi128_si64(_mm_packus_epi32(b4, ones));
const unsigned long long v = va + vb + (unsigned long long)carry + biasq;
const unsigned long long m = (v >> 15) & 0x0001000100010001ull;
const unsigned long long r = v - m * bias;
_mm_storeu_si128((__m128i *)(void *)(out + pos),
_mm_cvtepu16_epi32(_mm_cvtsi64_si128((long long)r)));
carry = 1ull - (v >> 63);
pos += 4;
}
}
while (pos < len2)
{
carry += (in1[pos] + in2[pos]);
out[pos] = carry < BASE ? carry : carry - BASE;
carry = carry < BASE ? 0 : 1;
pos++;
}
while (pos < len1 && carry > 0)
{
carry += in1[pos];
out[pos] = carry < BASE ? carry : carry - BASE;
carry = carry < BASE ? 0 : 1;
pos++;
}
ary_copy(out + pos, in1 + pos, len1 - pos);
if (is_carry)
{
out[len1] = carry % BASE;
}
}
template <INT_64 BASE, typename T>
constexpr void abs_sub(const T in1[], const T in2[], T out[],
size_t len1, size_t len2)
{
if (len1 < len2)
{
return;
}
size_t pos = 0;
INT_64 borrow = 0;
if (sizeof(T) == 4 && BASE <= 32768)
{
/* */
constexpr unsigned long long bias = 0x10000ull - (unsigned long long)BASE;
const __m128i ones = _mm_set1_epi32(-1);
/* */
while (pos + 8 <= len2)
{
const __m128i a4 = _mm_loadu_si128((const __m128i *)(const void *)(in1 + pos));
const __m128i b4 = _mm_loadu_si128((const __m128i *)(const void *)(in2 + pos));
const __m128i a4b = _mm_loadu_si128((const __m128i *)(const void *)(in1 + pos + 4));
const __m128i b4b = _mm_loadu_si128((const __m128i *)(const void *)(in2 + pos + 4));
const unsigned long long va = (unsigned long long)_mm_cvtsi128_si64(_mm_packus_epi32(a4, ones));
const unsigned long long vb = (unsigned long long)_mm_cvtsi128_si64(_mm_packus_epi32(b4, ones));
const unsigned long long va2 = (unsigned long long)_mm_cvtsi128_si64(_mm_packus_epi32(a4b, ones));
const unsigned long long vb2 = (unsigned long long)_mm_cvtsi128_si64(_mm_packus_epi32(b4b, ones));
const unsigned long long d0 = va - vb - (unsigned long long)borrow;
const unsigned long long m = (d0 >> 15) & 0x0001000100010001ull;
const unsigned long long d = d0 - m * bias;
_mm_storeu_si128((__m128i *)(void *)(out + pos),
_mm_cvtepu16_epi32(_mm_cvtsi64_si128((long long)d)));
const unsigned long long carry1 = (d0 >> 63) & 1ull;
const unsigned long long d0b = va2 - vb2 - carry1;
const unsigned long long mb = (d0b >> 15) & 0x0001000100010001ull;
const unsigned long long db = d0b - mb * bias;
_mm_storeu_si128((__m128i *)(void *)(out + pos + 4),
_mm_cvtepu16_epi32(_mm_cvtsi64_si128((long long)db)));
borrow = (INT_64)((d0b >> 63) & 1ull);
pos += 8;
}
while (pos + 4 <= len2)
{
const __m128i a4 = _mm_loadu_si128((const __m128i *)(const void *)(in1 + pos));
const __m128i b4 = _mm_loadu_si128((const __m128i *)(const void *)(in2 + pos));
unsigned long long va = (unsigned long long)_mm_cvtsi128_si64(_mm_packus_epi32(a4, ones));
unsigned long long vb = (unsigned long long)_mm_cvtsi128_si64(_mm_packus_epi32(b4, ones));
const unsigned long long d0 = va - vb - (unsigned long long)borrow;
const unsigned long long m = (d0 >> 15) & 0x0001000100010001ull;
const unsigned long long d = d0 - m * bias;
_mm_storeu_si128((__m128i *)(void *)(out + pos),
_mm_cvtepu16_epi32(_mm_cvtsi64_si128((long long)d)));
borrow = (INT_64)((d0 >> 63) & 1ull);
pos += 4;
}
borrow = borrow ? -1 : 0;
}
while (pos < len2)
{
borrow += (static_cast<INT_64>(in1[pos]) - in2[pos]);
out[pos] = (T)(borrow < 0 ? borrow + BASE : borrow);
borrow = borrow < 0 ? -1 : 0;
pos++;
}
while (pos < len1 && borrow < 0)
{
borrow += in1[pos];
out[pos] = borrow < 0 ? borrow + BASE : borrow;
borrow = borrow < 0 ? -1 : 0;
pos++;
}
ary_copy(out + pos, in1 + pos, len1 - pos);
}
template <UINT_64 BASE, typename T>
void fft_mul(const T in1[], const T in2[], T out[],
size_t len1, size_t len2)
{
using namespace hint_transform;
if (len1 == 0 || len2 == 0 || in1 == nullptr || in2 == nullptr)
{
return;
}
const size_t conv_res_len = len1 + len2 - 1; // 卷积结果长度
const size_t fft_len = min_2pow(conv_res_len); // fft长度
if (fft_len > g_fft_cap)
{
std::free(g_fft_buf);
g_fft_buf = (Complex *)std::malloc(sizeof(Complex) * fft_len);
g_fft_cap = fft_len;
}
Complex *fft_ary = g_fft_buf;
com_ary_combine_copy(fft_ary, in1, len1, in2, len2);
/* */
const size_t mx = len1 > len2 ? len1 : len2;
const bool high_zero = fft_len >= 16 && mx <= fft_len / 2;
{ // 原 new Complex[] 会逐一构造零值；复用缓冲必须显式补零
/* */
std::memset((void *)(fft_ary + mx), 0,
((high_zero ? fft_len / 2 : fft_len) - mx) * sizeof(Complex));
}
if (high_zero) fft_dif_zero_half(fft_ary, fft_len);
else fft_dif_template<1>(fft_ary, fft_len);
double inv = -1 / (2.0 * fft_len);
const __m256d scale_r=_mm256_set1_pd(inv);
const __m256d scale_i=_mm256_set1_pd(-2.0*inv);
size_t i=0;
for(;i+2<=fft_len;i+=2){
const __m256d v=_mm256_loadu_pd((const double *)(fft_ary+i));
const __m256d re=_mm256_movedup_pd(v);
const __m256d im=_mm256_permute_pd(v,15);
const __m256d real=_mm256_mul_pd(_mm256_sub_pd(_mm256_mul_pd(re,re),_mm256_mul_pd(im,im)),scale_r);
const __m256d imag=_mm256_mul_pd(_mm256_mul_pd(re,im),scale_i);
_mm256_storeu_pd((double *)(fft_ary+i),_mm256_blend_pd(real,imag,10));
}
for(;i<fft_len;++i){
Complex tmp=fft_ary[i];
fft_ary[i]=std::conj(tmp*tmp*inv);
}
fft_dit_template<1>(fft_ary, fft_len);
hint::UINT_64 carry = 0;
size_t ic = 0;
/* */
const hint::UINT_64 r62_per = (hint::UINT_64)(BASE - 1) * (BASE - 1)
* (hint::UINT_64)(BASE * BASE + BASE + 2);
if ((hint::UINT_64)((conv_res_len + 1) / 2) <= (~0ull) / r62_per)
{
const hint::UINT_64 B1 = (hint::UINT_64)BASE;
const hint::UINT_64 B2 = B1 * B1;
const hint::UINT_64 B3 = B2 * B1;
for (; ic + 3 <= conv_res_len; ic += 3)
{
const hint::UINT_64 c0 = static_cast<hint::INT_64>(fft_ary[ic ].imag() + 0.5);
const hint::UINT_64 c1 = static_cast<hint::INT_64>(fft_ary[ic + 1].imag() + 0.5);
const hint::UINT_64 c2 = static_cast<hint::INT_64>(fft_ary[ic + 2].imag() + 0.5);
const hint::UINT_64 s0 = carry + c0;
out[ic ] = (T)(s0 % B1);
out[ic + 1] = (T)((s0 / B1 + c1) % B1);
out[ic + 2] = (T)(((s0 + B1 * c1) / B2 + c2) % B1);
carry = (s0 + B1 * c1 + B2 * c2) / B3;
}
}
/* */
{
const hint::UINT_64 pair_per = (hint::UINT_64)(BASE - 1) * (BASE - 1)
* (hint::UINT_64)(BASE + 2);
if ((hint::UINT_64)((conv_res_len + 1) / 2) <= (~0ull) / pair_per)
{
for (; ic + 2 <= conv_res_len; ic += 2)
{
const hint::UINT_64 c0 = static_cast<hint::INT_64>(fft_ary[ic ].imag() + 0.5);
const hint::UINT_64 c1 = static_cast<hint::INT_64>(fft_ary[ic + 1].imag() + 0.5);
const hint::UINT_64 s0 = carry + c0;
out[ic] = (T)(s0 % BASE);
const hint::UINT_64 pair = s0 + (hint::UINT_64)BASE * c1;
out[ic + 1] = (T)((pair / BASE) % BASE);
carry = pair / ((hint::UINT_64)BASE * BASE);
}
}
}
for (; ic < conv_res_len; ic++)
{
carry += static_cast<hint::INT_64>(fft_ary[ic].imag() + 0.5);
std::tie(carry, out[ic]) = div_mod<UINT_64>(carry, BASE);
}
out[conv_res_len] = carry % BASE;
}
template <UINT_64 BASE, typename T>
constexpr T divisor_normalize(const T in[], T out[], size_t len)
{
if (in == out)
{
throw("In can't be same as out\n");
}
T multiplier = 1;
if (len == 1)
{
multiplier = (BASE - 1) / in[0];
out[0] = in[0] * multiplier;
}
else if (in[len - 1] >= (BASE / 2))
{
ary_copy(out, in, len);
}
else
{
multiplier = (BASE - 1) * BASE / (BASE * in[len - 1] + in[len - 2]);
abs_mul_num<BASE, false>(in, multiplier, out, len);
if (out[len - 1] < (BASE / 2))
{
multiplier++;
abs_add<BASE, false>(out, in, out, len, len);
}
}
return multiplier;
}
template <UINT_64 BASE, bool Prechecked=false, typename T>
void abs_long_div(T dividend[], const T divisor[], T quot[],
size_t len1, size_t len2)
{
if constexpr(!Prechecked){
len1 = ary_true_len(dividend, len1);
len2 = ary_true_len(divisor, len2);
if (divisor == nullptr || len2 == 0)
{
throw("Can't divide by zero\n");
}
if (dividend == nullptr || len1 == 0)
{
return;
}
if (abs_compare_normalized(dividend, divisor, len1, len2) < 0)
{
return;
}
}
if (len2 == 1)
{
T rem = abs_div_num<BASE>(dividend, divisor[0], quot, len1);
dividend[0] = rem;
return;
}
if (!Prechecked && divisor[len2 - 1] < (BASE / 2))
{
throw("Can't call this proc before normalize the divisor\n");
}
quot[len1 - len2] = 0;
const UINT_64 divisor_2digits = BASE * divisor[len2 - 1] + divisor[len2 - 2];
/*<https://duck.ac/submission/128944>: reciprocal of the top divisor limb, [r87e8-BP1] back-ported verbatim from saffah_codex_6a_agg3 #128944 */
const UINT_64 trial_top = divisor[len2 - 1];
const UINT_64 trial_recip = (BASE <= 65536) ? UINT_32(0xffffffffu) / UINT_32(trial_top) : 0;
(void)trial_top; (void)trial_recip;
hintvector<T> sub(len2 + 1);
/*[r101a-PACK1] 反向移植 saffah_codex_6a_agg3 #129675 <https://duck.ac/submission/129675>（逐行回件核）： */
static hintvector<UINT_64> triple_cache;
const size_t groups = (BASE <= 10000) ? len2 / 3 : 0;
if (triple_cache.capacity() < (SIZE_TYPE)groups)
{
triple_cache.resize((SIZE_TYPE)groups);
}
UINT_64 *packed = triple_cache.data();
if constexpr (BASE <= 10000)
{
for (size_t g = 0; g < groups; g++)
{
packed[g] = UINT_64(divisor[3 * g]) + BASE * divisor[3 * g + 1]
+ BASE * BASE * divisor[3 * g + 2];
}
}
auto multiply_trial = [&](T num)
{
if constexpr (BASE <= 10000)
{
constexpr UINT_64 B2 = BASE * BASE, B3 = B2 * BASE;
T *out = sub.data();
UINT_64 carry = 0;
if (num == 1)
{
ary_copy(out, divisor, len2);
}
else
{
size_t i = 0;
for (size_t g = 0; g < groups; g++, i += 3)
{
/* */
const UINT_64 value = packed[g] * num + carry;
const UINT_64 q2 = value / B3;
const UINT_64 r2 = value - q2 * B3;
const UINT_64 q1 = r2 / B2;
const UINT_64 r1 = r2 - q1 * B2;
const UINT_64 q0 = r1 / BASE;
out[i] = T(r1 - q0 * BASE);
out[i + 1] = T(q0);
out[i + 2] = T(q1);
carry = q2;
}
for (; i < len2; i++)
{
const UINT_64 value = UINT_64(divisor[i]) * num + carry;
carry = value / BASE;
out[i] = T(value - carry * BASE);
}
}
out[len2] = T(carry);
}
else
{
abs_mul_num<BASE, true, true>(divisor, num, sub.data(), len2);
}
};
while (abs_compare_normalized(dividend, divisor, len1, len2) >= 0)
{
sub.change_length(len2 + 1);
UINT_64 dividend_2digits = dividend[len1 - 1] * BASE + dividend[len1 - 2];
T quo_digit = 0;
size_t shift = len1 - len2;
if (dividend_2digits > divisor_2digits)
{
/*<https://duck.ac/submission/129699>：本支在偶 BASE 下**恒有 quo_digit == 1** [r101a-FAST1] 反向移植 saffah_codex_6a_agg3 #129699 */
if constexpr (BASE % 2 == 0)
{
sub.data()[len2] = 0; // 保持原 num==1 路径留下的进位槽状态
abs_sub<BASE>(dividend + shift, divisor, dividend + shift, len1 - shift, len2);
len1 = ary_true_len(dividend, len1);
quot[shift] = 1;
continue;
}
quo_digit = dividend_2digits / divisor_2digits;
multiply_trial(quo_digit);
sub.change_length(len2+(sub.data()[len2]!=0));
size_t sub_len = sub.length();
if (abs_compare_normalized(dividend + shift, sub.data(), len1 - shift, sub_len) < 0)
{
quo_digit--;
abs_sub<BASE>(sub.data(), divisor, sub.data(), sub_len, len2);
}
}
else if (dividend_2digits == divisor_2digits)
{
if (abs_compare_normalized(dividend + shift, divisor, len1 - shift, len2) < 0)
{
quo_digit = BASE - 1;
shift--;
multiply_trial(quo_digit);
}
else
{
quo_digit = 1;
ary_copy(sub.data(), divisor, len2);
sub.set_true_len();
}
}
else
{
if constexpr (BASE <= 65536)
{
const UINT_64 approx = (dividend_2digits * trial_recip) >> 32;
const UINT_64 rem = dividend_2digits - approx * trial_top;
quo_digit = T(approx + (rem >= trial_top));
}
else quo_digit = dividend_2digits / (divisor_2digits / BASE);
if (quo_digit >= BASE)
{
quo_digit = BASE - 1;
}
shift--;
multiply_trial(quo_digit);
sub.change_length(len2+(sub.data()[len2]!=0));
size_t sub_len = sub.length();
if (abs_compare_normalized(dividend + shift, sub.data(), len1 - shift, sub_len) < 0)
{
quo_digit--;
abs_sub<BASE>(sub.data(), divisor, sub.data(), sub_len, len2);
if (abs_compare(dividend + shift, sub.data(), len1 - shift, sub_len) < 0)
{
quo_digit--;
abs_sub<BASE>(sub.data(), divisor, sub.data(), sub_len, len2);
}
}
}
abs_sub<BASE>(dividend + shift, sub.data(), dividend + shift, len1 - shift, sub.length());
len1 = ary_true_len(dividend, len1);
quot[shift] = quo_digit;
}
}
template <UINT_64 BASE, typename T>
void abs_rec_div(T dividend[], T divisor[], hintvector<T> ",
size_t len1, size_t len2)
{
len1 = ary_true_len(dividend, len1);
len2 = ary_true_len(divisor, len2);
if (divisor == nullptr || len2 == 0)
{
throw("Can't divide by zero\n");
}
if (dividend == nullptr || len1 == 0)
{
return;
}
if (abs_compare_normalized(dividend, divisor, len1, len2) < 0)
{
return;
}
if (divisor[len2 - 1] < (BASE / 2))
{
throw("Can't call this proc before normalize the divisor\n");
}
size_t quot_len = len1 - len2 + 1;
quot.change_length(quot_len);
constexpr size_t LONG_DIV_THRESHOLD = 32;
if (len2 <= LONG_DIV_THRESHOLD) // 小于等于阈值调用长除法
{
abs_long_div<BASE,true>(dividend, divisor, quot.data(), len1, len2);
}
else if (len1 >= len2 * 2 || len1 > ((len2 + 1) / 2) * 3) // 2n/n的除法，进行两次递归
{
size_t base_len = (len1 + 3) / 4;
size_t quot_tmp_len = base_len * 3 - len2 + 2;
hintvector<T> quot_tmp(quot_tmp_len, 0);
abs_rec_div<BASE>(dividend + base_len, divisor, quot_tmp, len1 - base_len, len2);
quot_tmp_len = quot_tmp.set_true_len();
size_t dividend_len = ary_true_len(dividend, len1);
abs_rec_div<BASE>(dividend, divisor, quot, dividend_len, len2);
quot.change_length(quot_len);
quot_len = quot.set_true_len();
abs_add<BASE>(quot.data() + base_len, quot_tmp.data(), quot.data() + base_len, quot_len - base_len, quot_tmp_len);
quot.change_length(len1 - len2 + 1);
}
else
{
size_t base_len = len2 / 2;
abs_rec_div<BASE>(dividend + base_len, divisor + base_len, quot, len1 - base_len, len2 - base_len);
constexpr T ONE[1] = {1};
quot_len = quot.set_true_len();
hintvector<T> prod(base_len + quot_len, 0);
fft_mul<BASE>(divisor, quot.data(), prod.data(), base_len, quot_len);
size_t prod_len = prod.set_true_len();
len1 = ary_true_len(dividend, len1);
while (abs_compare_normalized(prod.data(), dividend, prod_len, len1) > 0)
{
abs_sub<BASE>(quot.data(), ONE, quot.data(), quot_len, 1);
abs_sub<BASE>(prod.data(), divisor, prod.data(), prod_len, base_len);
abs_add<BASE>(dividend + base_len, divisor + base_len, dividend + base_len, len1 - base_len, len2 - base_len);
quot_len = quot.set_true_len();
prod_len = prod.set_true_len();
len1 = ary_true_len(dividend, std::max(len1, len2) + 1);
}
abs_sub<BASE>(dividend, prod.data(), dividend, len1, prod_len);
}
}
template <UINT_64 BASE, typename T>
hintvector<T> abs_div(const T dividend[], T divisor[], hintvector<T> ",
size_t len1, size_t len2, bool ret_rem = true)
{
hintvector<T> normalized_divisor(len2); // 定义规则化的除数
normalized_divisor.change_length(len2);
hintvector<T> normalized_dividend(len1 + 1); // 定义规则化的被除数
normalized_dividend.change_length(len1 + 1);
T *divisor_ptr = normalized_divisor.data();
T *dividend_ptr = normalized_dividend.data();
T multiplier = divisor_normalize<BASE>(divisor, divisor_ptr, len2); // 除数规则化,获得乘数
abs_mul_num<BASE>(dividend, multiplier, dividend_ptr, len1); // 被除数规则化
len1 = normalized_dividend.set_true_len();
quot = hintvector<T>(len1 - len2 + 2, 0);
if ((!ret_rem) && (len1 + 2 < len2 * 2))
{
size_t shift = len2 * 2 - len1 - 2;
abs_rec_div<BASE>(dividend_ptr + shift, divisor_ptr + shift, quot, len1 - shift, len2 - shift);
quot.set_true_len();
return normalized_dividend;
}
abs_rec_div<BASE>(dividend_ptr, divisor_ptr, quot, len1, len2);
quot.set_true_len();
if (ret_rem)
{
len1 = normalized_dividend.set_true_len();
abs_div_num<BASE>(dividend_ptr, multiplier, dividend_ptr, len1); // 余数除以乘数得到正确的结果
normalized_dividend.set_true_len();
}
return normalized_dividend;
}}
struct DecimalLookup {
char digits[10000][4]{};
constexpr DecimalLookup(){
for(unsigned i=0;i<10000;++i){
digits[i][0]='0'+i/1000%10;
digits[i][1]='0'+i/100%10;
digits[i][2]='0'+i/10%10;
digits[i][3]='0'+i/1%10;
}
}
};
static constexpr DecimalLookup DECIMAL_LU{};
class Integer
{
public:
using DataType = hint::UINT_32;
using SizeType = hint::UINT_32;
using DataVec = hint::HintVector<DataType, SizeType>;
private:
DataVec data;
public:
static constexpr hint::UINT_32 DIGIT = 4;
static constexpr hint::UINT_64 BASE = 10000;
Integer()
{
data = DataVec();
}
Integer(const Integer &input)
{
if (this != &input)
{
data = input.data;
}
}
Integer(Integer &&input) noexcept
{
if (this != &input)
{
data = std::move(input.data);
}
}
Integer &operator=(const Integer &input)
{
if (this != &input)
{
data = input.data;
}
return *this;
}
Integer &operator=(Integer &&input) noexcept
{
if (this != &input)
{
data = std::move(input.data);
}
return *this;
}
Integer &operator=(const std::string &input)
{
string_in(input);
return *this;
}
Integer &operator=(const char input[])
{
string_in(input);
return *this;
}
DataType first_num() const
{
if (length() == 0)
{
return 0;
}
return data[length() - 1];
}
void change_sign(bool is_neg)
{
data.change_sign(is_neg);
}
bool is_neg() const
{
return data.sign();
}
SizeType length() const
{
return data.length();
}
SizeType length_base10() const
{
size_t len = data.length();
if (len == 0)
{
return 1;
}
return (len - 1) * DIGIT + std::ceil(std::log10(first_num() + 1));
}
void string_in(const std::string &str){string_in(str.data(),str.size());}
/*[r101a-IO1] 反向移植 saffah_codex_6a_agg3 #129699 <https://duck.ac/submission/129699>： */
__attribute__((target("avx2"))) void string_in(const char *str,size_t str_size)
{
if(str_size==0){data=DataVec();return;}
const size_t begin=str[0]=='-'?1u:0u;
const size_t limbs=(str_size-begin+3)/4;
data=DataVec(limbs);
data.change_length(limbs);
size_t pos=str_size,i=0;
const __m256i ascii0=_mm256_set1_epi8('0');
const __m256i tens=_mm256_set1_epi16(0x010a);
const __m256i hundreds=_mm256_set1_epi32(0x00010064);
const __m256i reverse=_mm256_setr_epi32(7,6,5,4,3,2,1,0);
while(pos-begin>=32){
const __m256i bytes=_mm256_sub_epi8(_mm256_loadu_si256((const __m256i*)(str+pos-32)),ascii0);
const __m256i pairs=_mm256_maddubs_epi16(bytes,tens);
const __m256i limbs8=_mm256_madd_epi16(pairs,hundreds);
_mm256_storeu_si256((__m256i*)(data.data()+i),_mm256_permutevar8x32_epi32(limbs8,reverse));
pos-=32;i+=8;
}
while(pos-begin>=4){
const char *p=str+pos-4;
data[i++]=(unsigned)(p[0]-'0')*1000u+(unsigned)(p[1]-'0')*100u+(unsigned)(p[2]-'0')*10u+(unsigned)(p[3]-'0')*1u;
pos-=4;
}
if(pos>begin){
unsigned v=0;
for(size_t k=begin;k<pos;++k)v=v*10u+(unsigned)(str[k]-'0');
data[i]=v;
}
change_sign(begin!=0);
data.set_true_len();
}
std::string to_string() const
{
size_t pos=length();
if(!pos)return "0";
std::string result=is_neg()?"-":"";
result+=std::to_string(first_num());
const size_t prefix=result.size();
result.resize(prefix+4*(pos-1));
char *dst=result.data()+prefix;
--pos;
const __m128i hundred_recip=_mm_set1_epi16(5243),hundred=_mm_set1_epi16(100);
const __m128i ten_recip=_mm_set1_epi16(205),ten=_mm_set1_epi16(10);
const __m128i ascii=_mm_set1_epi16(0x3030);
while(pos>=8){
const __m128i high=_mm_shuffle_epi32(_mm_loadu_si128((const __m128i*)(data.data()+pos-4)),0x1b);
const __m128i low=_mm_shuffle_epi32(_mm_loadu_si128((const __m128i*)(data.data()+pos-8)),0x1b);
const __m128i value=_mm_packus_epi32(high,low);
const __m128i q=_mm_srli_epi16(_mm_mulhi_epu16(value,hundred_recip),3);
const __m128i r=_mm_sub_epi16(value,_mm_mullo_epi16(q,hundred));
const __m128i qt=_mm_srli_epi16(_mm_mullo_epi16(q,ten_recip),11);
const __m128i rt=_mm_srli_epi16(_mm_mullo_epi16(r,ten_recip),11);
const __m128i qo=_mm_sub_epi16(q,_mm_mullo_epi16(qt,ten));
const __m128i ro=_mm_sub_epi16(r,_mm_mullo_epi16(rt,ten));
const __m128i qp=_mm_add_epi16(_mm_or_si128(qt,_mm_slli_epi16(qo,8)),ascii);
const __m128i rp=_mm_add_epi16(_mm_or_si128(rt,_mm_slli_epi16(ro,8)),ascii);
_mm_storeu_si128((__m128i*)dst,_mm_unpacklo_epi16(qp,rp));
_mm_storeu_si128((__m128i*)(dst+16),_mm_unpackhi_epi16(qp,rp));
pos-=8;dst+=32;
}
while(pos){
const unsigned v=data[--pos];
std::memcpy(dst,DECIMAL_LU.digits[v],4);dst+=4;
}
return result;
}
hint::INT_32 abs_compare(const Integer &input) const
{
size_t len1 = length(), len2 = input.length();
return hint_arithm::abs_compare(data.data(), input.data.data(), len1, len2);
}
Integer operator/(const Integer &input) const
{
Integer result;
size_t len1 = length(), len2 = input.length();
if (abs_compare(input) < 0)
{
return result;
}
if (len2 == 0)
{
throw("Can't divide by zero\n");
}
auto ptr1 = data.data();
auto ptr2 = input.data.data();
hint_arithm::abs_div<BASE>(ptr1, ptr2, result.data, len1, len2, false);
result.data.set_true_len();
result.change_sign(is_neg() != input.is_neg());
return result;
}
};
__attribute__((target("avx2"))) static const char* token_end(const char*p,const char*end){
const __m256i space=_mm256_set1_epi8(32);
while(end-p>=32){
unsigned mask=(unsigned)_mm256_movemask_epi8(_mm256_cmpgt_epi8(_mm256_loadu_si256((const __m256i*)p),space));
if(mask!=~0u)return p+__builtin_ctz(~mask);
p+=32;
}
while(p<end&&(unsigned char)*p>32)++p;
return p;
}
int main()
{
std::ios::sync_with_stdio(false);std::cin.tie(nullptr);
size_t cap=1u<<18,used=0;
char*raw=(char*)std::malloc(cap);
if(!raw)return 1;
while(true){
if(used==cap){size_t next=cap*2;char*q=(char*)std::realloc(raw,next);if(!q){std::free(raw);return 1;}raw=q;cap=next;}
const size_t got=std::fread(raw+used,1,cap-used,stdin);used+=got;
if(!got)break;
}
const char*p=raw,*end=raw+used;
Integer a,b;
while(p<end&&(unsigned char)*p<=32)++p;
const char*q=token_end(p,end);a.string_in(p,q-p);p=q;
while(p<end&&(unsigned char)*p<=32)++p;
q=token_end(p,end);b.string_in(p,q-p);
std::free(raw);
std::cout<<(a/b).to_string()<<'\n';
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 4.201 ms | 1 MB + 140 KB | Accepted | Score: 100 | 显示更多 |