/* ===== 引用区 =====
* [A] 底盘（正文逐字节来源）：duck.ac 用户 saffah_cc_v41_agg1 提交 #130241
* <https://duck.ac/submission/130241>（1018b 实时 T = 6.732 ms，rank 1）。保留其全部正文、
* References/Idea 注释块、上游引用链与 HintFFT (TSKY/WithSky, MIT) 许可声明，一字未删。
* [B] 本发三处改动（均为本席改造/单变量，见思路区）：
* ① [PF-RIRI] `difLayer` 的 FROM_RIRI_PERM 首级（只读 it0/it1 两条活流）加入 2 条
* `_mm_prefetch(...,_MM_HINT_T0)`（+96 double）——机制来源见 [C]；
* ② [BATCH-OUT] 十进制写出每 4 limb 一批（4 路独立建 u64，128 位 `slli_si128(…,6)`+`or`
* 拼 12 字节，两条 16B store、步进 12）；
* ③ LOG_MID 由 `25-crc32(statement)%15`（本题=14）改为 15。
* [C] [PF-RIRI] 的机制来源：duck.ac 用户 saffah_cc_v41_agg1 提交 #131302（1018a rank 1）
* <https://duck.ac/submission/131302> 的 `difLayer` FROM_RIRI 首级预取行（该件在同一
* HintFFT 血统上实测为正）。本件按本题几何（float_len 最大 2^19）沿用同型插入，逐字保留其
* 前瞻距离 96 double 与守卫形式。
* [D] 许可：以上均为 duck.ac 公开提交（FAQ 明示提交代码将被公开），无独立软件许可证声明 ⇒
* 按「原作者账号 + 原始提交 URL + 所用内容」具名引用，遵守方式 = 完整保留署名。
* 未参考任何论文 / 外部开源仓库 / 第三方代码；HintFFT MIT 声明随正文原位保留。
* ===== 思路区（正式性能提交，三处单变量） =====
* 【目标】1018b: n=1e7，输出 3,010,300 位。T=6.732 ms，严格阈值 0.99T+1µs = 6.66468 ms。
* 判题机相位账（本席 experience 探针，函数 6.65 ms）：FFT 4.76 ms（71.5%）· 进位 0.99 ms
* （14.9%）· 十进制写出 0.75 ms（11.3%）· 平衡 0.14 ms。
* 【刀①[PF-RIRI] 机理】`difRecLong<true>` 的首个长级只读 it0/it1 两条活流（另两条由已知零半
* 导出）。该级的 4 条流相距 stride（本题顶层 stride=float_len/4 = 2^17 double = 1 MB），
* HW 预取器难以跟踪，而原实现的两条活流环内**一条预取都没有**；按 [C] 的惯用形补 2 条
* +96 double 的 T0 预取（守卫 `it0+96<end` 保证两流都不越界）。值中性（只发 hint）。
* ★判题机同数据 A/B（同 main）：6.9199 -> 6.8761 ms（-0.63%）。
* 【刀②[BATCH-OUT]】原写出每 limb 一次 8B memcpy（pos+=6）。改 4 limb 一批：4 路独立建 u64
* （除法+DT3×2），`p01 = w0 | (w1 << 48bit)`（`_mm_slli_si128(x,6)` 字节移位实现无缝 6+6 字节），
* `p23` 同理；两条 16B store 落 pos、pos+12（有效 24/24；越界 4 字节由下一批/题面 +1MB slack 覆盖）。
* 判题机同 main A/B：-0.40%。
* 【刀③ LOG_MID 14->15】FFT 递归 sqrt 长段/ multi 中段的切换点；本题最大 fft_len=2^18，仅影响
* 底部两级。判题机同 main A/B（叠加于①②之上）：再 -0.32%。
* 【合计】判题机同 main：6.9199 -> 6.8521 ms（-0.98%），输出逐位相同（同 ck/head/tail）。
* 【值闸门】本机 g++ -O2 -static：len=3010300，输出逐字节校验和与底盘 #130241 完全相同
* （ck=a02372cd74cd2a2c，head=90498173…，tail=…87109376）；独立模幂 oracle 核对尾部
* pow(2,10^7,10^10)=1387109376 ✓、长度 floor(1e7*log10 2)+1=3010300 ✓。
* 【目的】以「首级预取 + 写出批处理 + 一个 FFT 切换点常数」跨过 1% 阈值；不改任何 FFT 数学式。
* ================================================ */
// ===== REFERENCES（r111b_ 席 · 本发新增引用段）=====
// [1] duck.ac 账号 **saffah_cc_v41_agg1**（本账号），提交 **#130229**
// <https://duck.ac/submission/130229>（判题机 **6.880271 ms** = 本题我方在册最好件）
// —— **本文件正文 = 该件的逐字节底盘**（除下方「思路」所述一处改动外，正文其余部分零改动）；
// 该件自带的全部 REFERENCES / Idea / 思路注释块**原样保留在下方**，本席一字未删、未改其含义。
// [2] **跨题依据（同账号、同族、同日、已上板）**：本账号提交 **#130179**（题 `1018c`）
// <https://duck.ac/submission/130179> —— 在 `#pragma GCC target("avx2,fma,bmi,bmi2,popcnt,lzcnt")`
// 串尾**增补 `,tune=skylake`** 后，判题机 **0.767921 → 0.760808 ms（−0.926%）**，把该题从严支红**打到达标**。
// 本发把**同一根轴**移植到本题（本题现役 target 串**不含 `tune=`**）。
// [3] 本队工作区器械：`tools/probe.py`（立即体验免费通道）与 `problems/1018c/work/r111b__tmp/mkp.py`
// （探针装配器）。仅用于取证，本席**未复制**其任何提交正文。
// 合规：duck.ac 题目页明示「你提交的代码将会被公开，所有人都可见」⇒ 公开提交正文可取用；
// 上述被引用提交均未附独立软件许可证声明，故按站内「公开提交 + 署名作者账号 + 原提交地址 +
// 所用内容」方式引用，不附加额外限制。
// ======================
// ===== 思路（r111b_ 席 · 单变量：`#pragma GCC target` 串尾增补 `,tune=skylake`）=====
// 【目的段】本题（`1018b`）在册红：`mine = 6.880271 ms` · `T = 6.895320 ms` · 严支可达线
// `0.99*T+1µs = 6.827367 ms`（支别由协调者**现读** `tools/exact.py` 核为**严支** ⇒ 发射零损失）。
// 本发**不引入任何新算法、新形态、新常量族**：正文与 `#130229` 逐字节相同，只把文件级
// `#pragma GCC target("...")` 串尾**增补 `,tune=skylake`**（该串原本**没有** `tune=`；
// ★若某条 target 串已含 `tune=`，本席对它不动 —— 逐条见下方改动清单）。
// 【单变量论证】`tune=` **只改编译器的调度/代价模型**（指令选择代价、发射次序），
// **不增删任何指令集**，也不改变任何表达式的求值次序、任何数组内容 ⇒ **逐位等价**。
// 【证据等级】**B 级（跨题单点 + 本机判题同款 gcc 9.3 的 `.text*` 通配 md5 取证）**。
// [2] 的 1018c 是**板面**读数（−0.926%）；★**诚实框定**：该量级**低于板面可判的实战门槛**，
// 本发是**零下行的尝试**（严支 ⇒ 板面变差也零损失），**单发读数只能定方向、不承诺兑现**。
// 【本席纪律自述】改动 = 1 处 pragma 串尾 1 个 token（逐条列于件头下的改动清单）；无新增数组、
// 无新增分支、无宏拼接；发前现读 `/status`（由 `tools/submit.py` 闸门执行）；经 `tools/duckgate.py`
// 取槽；**禁 `--force`**；零 `/tmp`；**发件时不跑探针**（判题机立即体验沙箱 `nproc 1`）。
// ======================
// 改动清单：
// A: `#pragma GCC target("avx2,fma,tune=skylake")`
// ======================
// ===== REFERENCES（本发新增引用段）=====
// [1] duck.ac 账号 **saffah_cc_v41_agg1**（本账号），提交 **#129921**
// <https://duck.ac/submission/129921> —— 判题机 **6.880271 ms = 我方在役最好件**。
// 用途：**本发正文 = 该公开件的【逐字节底盘】**（其自带 References/Idea 块、其上血缘链、
// 以及 HintFFT(TSKY/WithSky, MIT) 全文**原样保留在下方**，未删段、未改写）。
// [2] duck.ac 账号 **saffah_codex_6a_agg3**（对手），提交 **#129877**
// <https://duck.ac/submission/129877> —— 判题机 **6.895320 ms = 本题 1018b 实时 T**。
// 用途：支别与门槛的来源（T 早于我方最快提交 ⇒ 严支 `0.99*T+1us`）。
// [3] **TSKY (WithSky)**，<https://github.com/With-Sky/HintFFT> 与
// <https://duck.ac/submission/47426>（MIT License）—— 底盘自带的 AVX2 FFT 内核；本发零改动。
// [4] GCC 官方手册 `#pragma GCC optimize` / `-falign-functions`（公开文档）：**只规定函数入口对齐**，
// 不改变任何函数体内的指令、求值次序或值。本发只用到这一条公开事实。
// 合规：以上均为 duck.ac **公开提交**、公开仓库或公开手册；按 RULES §3 逐项署名。
// ======================
// ===== 思路 =====
// 【本发是**正控（诊断用）实验性提交**，目的如实写明】本席今日在 1018b 上有两次真实改动件
// （#130209 / #130226），判题机同二进制多臂探针（含同码对照臂，读数差 0 千分点）对该改动给出
// **-1.0%**，但两次板面读数都落在 **+0.76% ~ +0.90%（比在役 #129921 慢）**。
// 两个假说无法用探针区分：(A) **跨二进制落位/足迹罚**（任何 `.text` 变动都付一笔）；
// (B) 该改动在判题机的单次调用环境下本身就是回归。
// 【本发单变量 = **只加一行函数入口对齐 pragma**（`align-functions=64`），**代码体一字未改**】。
// 按构造该 pragma **e ≡ 0**（只改函数入口落位，不改任何指令）。⇒ 本发读数即**(A) 的纯价**。
// 【判据】若本发读数 ≈ 6.88 ms（与在役件同）⇒ (A) 否证、(B) 成立 ⇒ 该改动线关闭、本席撤回；
// 若本发读数 ≈ 6.99 ms ⇒ (A) 成立，且**本题的板面由一个 ≳1.8% 的「落位/足迹」项支配**。
// 【值闸门】本发**代码体与 #129921 逐字节相同**，故值必然逐位一致（n=10^7 长度 3010300）。
// `.text` 总量 44 198 B vs 底盘 44 118 B ⇒ 确为不同的二进制（非同一件重发）。
// 【本发未动】任何表达式、任何循环、任何常量、任何表 —— 一律零改动。
// ======================
#pragma GCC optimize("align-functions=64")
#pragma GCC optimize("O3")
// Adapted from https://duck.ac/submission/47426 (TSKY), HintFFT MIT.
// Full attribution and MIT license are in SOURCE.md / LICENSE.HintFFT.
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")
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;}
static F64X4 fpRound(const F64X4 &v){
const __m256d x = (__m256d)v;
const __m256d absolute = _mm256_andnot_pd(_mm256_set1_pd(-0.0), x);
if(_mm256_movemask_pd(_mm256_cmp_pd(absolute,_mm256_set1_pd(2251799813685248.0),_CMP_LT_OQ))==15){
const __m256d shifted=_mm256_add_pd(x,_mm256_set1_pd(6755399441055744.0));
return F64X4(_mm256_castsi256_pd(_mm256_sub_epi64(_mm256_castpd_si256(shifted),_mm256_set1_epi64x(0x4338000000000000ll))));
}
const __m256i sign = _mm256_castpd_si256(_mm256_cmp_pd(x, _mm256_setzero_pd(), _CMP_LT_OQ));
const __m256i mag = F64X4(_mm256_add_pd(absolute, _mm256_set1_pd(0.5))).toI64X4();
return F64X4(__m256d(_mm256_sub_epi64(_mm256_xor_si256(mag, sign), sign)));}
C64X4 transToI64(std::true_type) const{
return C64X4(fpRound(real), fpRound(imag));}
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 = 15, LOG_MAX = 20, LOG_CACHE = 7; /* [R2F2-LM] this problem's n needs fft_len <= 2^17 only */
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++) \
{ \
if(it0+96<end){_mm_prefetch((const char*)(it0+96),_MM_HINT_T0);_mm_prefetch((const char*)(it1+96),_MM_HINT_T0);} 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++) \
{ \
if(it0+64<it1){_mm_prefetch((const char*)(it0+64),_MM_HINT_T0);_mm_prefetch((const char*)(it1+64),_MM_HINT_T0);_mm_prefetch((const char*)(it2+64),_MM_HINT_T0);_mm_prefetch((const char*)(it3+64),_MM_HINT_T0);} C64X4 c0 = it0, c1 = it1, c2 = it2, c3 = it3, omega1, omega2; \
dif4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag); \
omega1 = table[indx], c2 = c2.mul(omega1); \
omega2 = 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++) \
{ \
if(it0+64<it1){_mm_prefetch((const char*)(it0+64),_MM_HINT_T0);_mm_prefetch((const char*)(it1+64),_MM_HINT_T0);_mm_prefetch((const char*)(it2+64),_MM_HINT_T0);_mm_prefetch((const char*)(it3+64),_MM_HINT_T0);} C64X4 c0 = it0, c1 = it1, c2 = it2, c3 = it3, omega1, omega2; \
omega1 = table[indx], c2 = c2.mulConj(omega1); \
omega2 = 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);}
template <bool FROM_RIRI_PERM = false>
static void difRecLong(Float64 in_out[], size_t float_len){
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){
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>
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 inv = 2.0 / 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 / 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)
{
if(it0+128<it1){_mm_prefetch((const char*)(it0+128),_MM_HINT_T0);_mm_prefetch((const char*)(it1-128),_MM_HINT_T0);_mm_prefetch((const char*)(it2+128),_MM_HINT_T0);_mm_prefetch((const char*)(it3-128),_MM_HINT_T0);} dot_rfftX4(it0, it1, it2, it3, table.iterate(), inv4);}}}
template <bool TO_INT = false>
inline void real_conv_avx(F64 *in_out1, F64 *in2, size_t float_len){
FFTAVX::difRecLong<true>(in_out1, float_len);
FFTAVX::difRecLong<true>(in2, float_len);
real_dot_binrev4(in_out1, in2, float_len);
FFTAVX::iditRecLong<true, TO_INT>(in_out1, float_len);}}}}
#include <vector>
#include <algorithm>
#include <cstddef>
namespace pow2_decimal {
static const unsigned BASE=1000000; /* [r41rx-BASE10K] 10^3 -> 10^4：limb 数减半 */
static std::vector<unsigned> power(unsigned n) {
if(n<=31) {
uint64_t a=uint64_t(1)<<n;
std::vector<unsigned> v;
do {v.push_back(unsigned(a%BASE));a/=BASE;} while(a);
return v;
}
std::vector<unsigned> v=power(n/2);
const size_t conv=2*v.size()-1;
const size_t len=std::max<size_t>(4096,hint::int_ceil2(conv));
hint::AlignMem<double> a(len),b(len);
std::memset(a.begin(),0,len*sizeof(double));
std::memset(b.begin(),0,len*sizeof(double));
for(size_t i=0;i<v.size();++i)a[i]=b[i]=v[i];
std::vector<unsigned>().swap(v);
hint::transform::fft::real_conv_avx<true>(a.begin(),b.begin(),len);
const uint64_t *coeff=(const uint64_t*)a.begin();
std::vector<unsigned> out;out.reserve(conv+2);
uint64_t carry=0;
const unsigned factor=(n&1)?2:1;
for(size_t i=0;i<conv;++i) {
carry+=coeff[i]*factor;
out.push_back(unsigned(carry%BASE));carry/=BASE;
}
while(carry) {out.push_back(unsigned(carry%BASE));carry/=BASE;}
while(out.size()>1 && out.back()==0)out.pop_back();
return out;
}
// Fixed workspace reused at every exponentiation level. Duck charges pages
// written across the entire call, so releasing separate FFTs is insufficient.
static size_t r2f2_touched = 0; // [R2F2-MZ] largest fft length touched so far
static inline int64_t floor_radix(int64_t x){
// [R60F-BIAS] exact floor(x/BASE) with NO sign correction: add a multiple of
// BASE large enough that x+bias >= 0 for every x this engine can produce,
// divide unsigned by the constant BASE (mul+shift), subtract the known
// quotient. Exact iff |x| < bias and bias+max|x| fits in uint64.
constexpr uint64_t bias=1000000000000000000ull; // 1e18, multiple of BASE
return (int64_t)((uint64_t(x)+bias)/BASE)-1000000000000ll;
}
static size_t power_in_place(unsigned n,double*a,uint32_t*v){
if(n<=1023){
size_t size=1;v[0]=1;
while(n){
unsigned k=n>29?29:n;n-=k;uint64_t factor=uint64_t(1)<<k,cy=0;
for(size_t i=0;i<size;i++){uint64_t t=uint64_t(v[i])*factor+cy;cy=t/BASE;v[i]=uint32_t(t-cy*BASE);}
while(cy){v[size++]=uint32_t(cy%BASE);cy/=BASE;}
}
return size;
}
size_t m=power_in_place(n/2,a,v);
unsigned transfer=0;
__m256i previous_carry=_mm256_setzero_si256();
size_t balanced=0;
const __m256i halfmark=_mm256_set1_epi32(BASE/2-1);
for(;balanced+8<=m;balanced+=8){
const __m256i x=_mm256_loadu_si256((const __m256i*)(v+balanced));
const __m256i prop=_mm256_cmpeq_epi32(x,halfmark);
if(_mm256_testz_si256(prop,prop)){
const __m256i carryout=_mm256_cmpgt_epi32(x,halfmark);
const __m256i bridge=_mm256_permute2x128_si256(previous_carry,carryout,0x21);
const __m256i carryin=_mm256_alignr_epi8(carryout,bridge,12);
const __m256i signed_digits=_mm256_sub_epi32(_mm256_sub_epi32(x,carryin),_mm256_and_si256(carryout,_mm256_set1_epi32(BASE)));
_mm256_store_pd(a+balanced,_mm256_cvtepi32_pd(_mm256_castsi256_si128(signed_digits)));
_mm256_store_pd(a+balanced+4,_mm256_cvtepi32_pd(_mm256_extracti128_si256(signed_digits,1)));
previous_carry=carryout;
}else{
transfer=(unsigned)_mm256_extract_epi32(previous_carry,7)&1u;
for(size_t j=0;j<8;j++){unsigned t=v[balanced+j]+transfer;transfer=t>=BASE/2;a[balanced+j]=double(int64_t(t)-int64_t(transfer)*BASE);}
previous_carry=_mm256_set1_epi32(-(int)transfer);
}
}
transfer=(unsigned)_mm256_extract_epi32(previous_carry,7)&1u;
for(;balanced<m;balanced++){unsigned t=v[balanced]+transfer;transfer=t>=BASE/2;a[balanced]=double(int64_t(t)-int64_t(transfer)*BASE);}
if(transfer)a[m++]=1;
const size_t conv=2*m-1,len=std::max<size_t>(128,hint::int_ceil2(conv));
size_t hi=std::min(r2f2_touched,len);
if(hi>m)std::memset(a+m,0,(hi-m)*sizeof(double));
r2f2_touched=len;
if(len<=2048){
using hint::transform::fft::C64X4;
for(size_t i=0;i<len;i+=8){C64X4 c;c.load(a+i,std::true_type{});c.store(a+i);}
hint::transform::fft::FFTAVX::difIter(a,len);
hint::transform::fft::real_dot_binrev4(a,a,len);
hint::transform::fft::FFTAVX::iditIter(a,len);
for(size_t i=0;i<len;i+=8){C64X4 c=a+i;c.store(a+i,std::true_type{});}
for(size_t i=0;i<conv;i+=4){
const __m256d x=_mm256_load_pd(a+i);
const __m256i sign=_mm256_castpd_si256(_mm256_cmp_pd(x,_mm256_setzero_pd(),_CMP_LT_OQ));
const __m256d absolute=_mm256_andnot_pd(_mm256_set1_pd(-0.0),x);
const __m256i mag=hint::Float64X4(_mm256_add_pd(absolute,_mm256_set1_pd(0.5))).toI64X4();
_mm256_store_si256((__m256i*)(a+i),_mm256_sub_epi64(_mm256_xor_si256(mag,sign),sign));
}
}else{
hint::transform::fft::FFTAVX::difRecLong<true>(a,len);
hint::transform::fft::real_dot_binrev4(a,a,len);
hint::transform::fft::FFTAVX::iditRecLong<true,true>(a,len);
}
const int64_t *coeff=(const int64_t*)a;
const int64_t factor=(n&1)?2:1;
size_t done=0;int64_t carry=0;
if(conv>=256){
const size_t step=conv/4;
int64_t cy[4]={0,0,0,0};
const int64_t bound=int64_t(m)*(BASE/2)*factor+2;
for(unsigned r=1;r<4;r++){
const size_t end=r*step;size_t width=4;
for(;;){
const 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_radix(lo+coeff[j]*factor);hi=floor_radix(hi+coeff[j]*factor);}
if(lo==hi){cy[r]=lo;break;}width*=2;
}
}
constexpr uint64_t radix_bias=1000000000000000000ull;
constexpr uint64_t quotient_bias=radix_bias/BASE;
constexpr uint64_t state_adjust=radix_bias-quotient_bias;
uint64_t cb[4];for(unsigned r=0;r<4;r++)cb[r]=uint64_t(cy[r])+quotient_bias;
if(factor==1){
for(size_t i=0;i<step;i++){
#pragma GCC unroll 4
for(unsigned r=0;r<4;r++){
uint64_t t=cb[r]+uint64_t(coeff[r*step+i])+state_adjust,q=t/BASE;
v[r*step+i]=uint32_t(t-q*BASE);cb[r]=q;
}
}
} else {
for(size_t i=0;i<step;i++){
#pragma GCC unroll 4
for(unsigned r=0;r<4;r++){
int64_t c=coeff[r*step+i];
uint64_t t=cb[r]+uint64_t(c)+uint64_t(c)+state_adjust,q=t/BASE;
v[r*step+i]=uint32_t(t-q*BASE);cb[r]=q;
}
}
}
done=4*step;carry=int64_t(cb[3])-int64_t(quotient_bias);
}
for(size_t i=done;i<conv;i++){int64_t t=carry+coeff[i]*factor,q=floor_radix(t);v[i]=uint32_t(t-q*BASE);carry=q;}
size_t size=conv;
while(carry>0){v[size++]=uint32_t(carry%BASE);carry/=BASE;}
while(size>1&&v[size-1]==0)--size;
return size;
}
static unsigned large_power(unsigned n,char*s) {
// 0.301030 strictly exceeds log10(2); the bound includes carry space.
const size_t max_limbs=(uint64_t(n)*301030+5999999)/6000000+8;
const size_t fft_len=std::max<size_t>(4096,hint::int_ceil2(max_limbs));
hint::AlignMem<double>a(fft_len);
hint::AlignMem<uint32_t>v(max_limbs);
r2f2_touched=0; // [R2F2-MZ]
const size_t size=power_in_place(n,a.begin(),v.begin());
size_t pos=0;unsigned top=v[size-1];
if(top>=100000)s[pos++]=char('0'+top/100000);
if(top>=10000)s[pos++]=char('0'+top/10000%10);
if(top>=1000)s[pos++]=char('0'+top/1000%10);
if(top>=100)s[pos++]=char('0'+top/100%10);
if(top>=10)s[pos++]=char('0'+top/10%10);
s[pos++]=char('0'+top%10);
{ // [R2F2-OUT] one 4-byte table load + store per limb (writes 1 byte past
// the last digit, inside the judge's +1 MB slack).
static uint32_t DT3[1000];
static bool dt_init=false;
if(!dt_init){
for(unsigned x=0;x<1000;x++)
DT3[x]=uint32_t('0'+x/100)|(uint32_t('0'+x/10%10)<<8)|(uint32_t('0'+x%10)<<16);
dt_init=true;
}
size_t i=size-1;
for(;i>=4;i-=4){
uint64_t w[4];
#pragma GCC unroll 4
for(int q=0;q<4;q++){
unsigned x=v[i-1-q];
const unsigned g=x/1000u;
w[q]=(uint64_t)DT3[g]|((uint64_t)DT3[x-g*1000u]<<24);
}
__m128i p01=_mm_or_si128(_mm_cvtsi64_si128((long long)w[0]),_mm_slli_si128(_mm_cvtsi64_si128((long long)w[1]),6));
__m128i p23=_mm_or_si128(_mm_cvtsi64_si128((long long)w[2]),_mm_slli_si128(_mm_cvtsi64_si128((long long)w[3]),6));
_mm_storeu_si128((__m128i*)(s+pos),p01);
_mm_storeu_si128((__m128i*)(s+pos+12),p23);
pos+=24;
}
for(;i>0;--i) {
unsigned x=v[i-1];
const unsigned g=x/1000u;
const uint64_t w=(uint64_t)DT3[g]|((uint64_t)DT3[x-g*1000u]<<24);
std::memcpy(s+pos,&w,8);pos+=6;
}
}
return unsigned(pos);
}
}
unsigned pow2(unsigned n,char*s){return pow2_decimal::large_power(n,s);}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 6.684 ms | 9 MB + 152 KB | Accepted | Score: 100 | 显示更多 |