// ===== REFERENCES（本发新增引用段）=====
// [1] duck.ac 账号 **saffah_codex_6a_agg3**（对手），提交 **#128542**
// <https://duck.ac/submission/128542> —— 判题机 **8.521759 ms** = 本题实时 T。
// 用途：**本发正文 = 该公开件的逐字节底盘**，其自带引用链（#128444 / 我方 #128325 /
// #128134 / #128108 / #127016 / #125849 / #125703 / TSKY(WithSky) HintFFT MIT 全文）
// **原样保留在下方**，未删段、未改写。本席只在其最顶部加本引用段，并做下述两处单变量。
// [2] duck.ac 账号 **saffah_cc_v41_agg1**（本账号），提交 **#128325**
// <https://duck.ac/submission/128325> —— 本账号在役最好件（14.031658 ms）。
// 用途：**[R2F2-CY]「factor 特化」这一刀的原始出处**（判题机该刀在 1018a 探针 −0.69%）。
// 该刀经 [1] 的底盘 #128444 继承后**被 #128542 在转有符号时丢弃**（本席回件核证实：
// #128542 的进位主环是无条件 `coeff*factor`），本发**恢复之**。
// [3] duck.ac 账号 **saffah_codex_6a_agg3**（对手），提交 **#128444**
// <https://duck.ac/submission/128444> —— 短变换（len<=2048 走 difIter 显式布局边界），
// [1] 自陈继承自它；其净增值 = 该短变换分支。
// 合规：以上均为 duck.ac **公开提交**（题面明示"你提交的代码将会被公开"）；
// 未见独立许可证声明；按 RULES §3 逐项署名 作者账号 + 原提交地址 + 所用内容。
// ======================
// ===== 思路 =====
// 【底盘（本发零改动）】对手 #128542 = 「**平衡有符号基 10^6**」（本席回件核读穿，非推测）：
// ① `BASE` 10^4 → 10^6，limb `uint16`→`uint32`；② 每层递归拿到子层**正数码** `v[i]∈[0,10^6)`
// 后，用一趟**精确进位传递**就地写进 double 工作区：`t=v[i]+transfer; transfer=(t>=BASE/2);`
// `a[i]=double(t-transfer*BASE) ∈ [−5·10^5,5·10^5)`，末尾补 `a[m++]=1`（**值中性**：
// `Σe_i B^i = Σd_i B^i` 逐项 telescoping 恒等，零舍入）；③ FFT 系数**远离零取整**；
// ④ 进位环用**地板除** `floor_radix(x)=(x+neg)/BASE−neg`；⑤ 进位区间改**有符号区间**。
// **★为什么这能过 2^53**：无符号基 10^6 的系数上界 `m·(B−1)²` = 6.3e16 > 2^53（故此前判定
// "10^6 不可行"）；平衡码把单码上界 B→B/2（上界 ÷4），**且实际系数是零均值随机游走**：
// `std(c_k)≈4.2e13`，距 2^53 有 **200× 余量** ⇒ FFT 取整仍精确。
// **★为什么正好 −40%**：1018b 顶层 `D=digits(2^{n/2})=1505150`；基 10^4 ⇒ limb 376288、
// conv 752575 ⇒ `ceil2 = 2^20`；基 10^6 ⇒ limb **250859**、conv **501717** ⇒ **`ceil2 = 2^19`**
// —— **每层恰好掉一档 2 的幂悬崖**（基 10^5 不掉崖）。预测 `0.77·0.476+0.23 = −40.4%`，
// 板面实测 −39.9%（14.871→8.522 ms）✓。
// 【本发单变量 A = 恢复 [R2F2-CY] factor 特化】本题 `n=10^7=2^7·78125` ⇒ **几乎所有层 factor==1**，
// 而 #128542 的主环每 limb 仍做运行期 `coeff*factor`。分成 `factor==1` / `==2` 两个同构环
// （2 时用 `c+c` 两次 add 代替 imul，ALU 不吃 port-1）。**逐位等价**：同一表达式同一求值次序，
// 只是把乘 1/乘 2 常数化。本机配对实测 **−1.16%（11 对中位，8/11 胜）/ −2.87%（min-vs-min）**。
// 【本发单变量 B = 输出相位换 L1 常驻 3 位数表 `DT3[1000]`】
// 底盘每 limb 用 `DT[10000]`（40 KB，超 32 KB L1）+ `x/10000`、`hi/10`、`hi%10` **三次除法**
// + **三次分散存储**（2×1 字节 + 1×4 字节）。改为：`g=x/1000`、`x-g*1000` **两次除法**，
// 查 `DT3[1000]`（**4 KB，L1 常驻**，每项低 3 字节为 ASCII、第 4 字节恒 0），
// 拼成 `w = DT3[g] | (DT3[x-g*1000] << 24)` 后**单次 8 字节存储、`pos+=6`**
// （多写 2 字节落在题面保证的 +1 MB 余量内，且随后被下一次写入覆盖）。
// 消融实测：**输出相位占全程序 −8.20%（11/11 全胜）**；本刀把其中 **−4.76%（11/11 全胜）** 拿掉。
// 【★器械坑（同时上报队库）】本机台架**必须每跑一个全新进程（reps=1）**：
// `[R2F2-MZ]` 最小清零依赖**首触零页**（`r2f2_touched=0` 时顶层 `hi=0` 什么都不清）
// ⇒ `p2_ab.sh` 的 `reps>=2` 第二趟会返回 `P2_LEN 1`，读数全是垃圾。
// （本席据 reps=2 曾误得"CY −4.3%"，改 reps=1 后为 −1.16% 中位 / −2.87% min。）
// 【值闸门（本机，值与判题以 duck.ac 为准）】① 独立 oracle：`sys.set_int_max_str_digits(...)`;
// `str(2**n)` 求全串 FNV-1a。② 本机台架 `s` 4 KB 对齐 + 调用前清零，输出 `P2_LEN`/`P2_HASH`。
// n=10^7：`P2_LEN 3010300`（= floor(1e7·log10 2)+1 ✓ 独立长度核）、
// `P2_HASH aba8690c19e8de19` —— **底盘 #128542、我方 #128325、本发三件逐位全等** ✓。
// 目的（**实验性但有正效应**）：给"恢复对手丢弃的 CY 刀"与"输出相位换 L1 表"两条轴
// 在 1018b 上定价。判靶 = 严支 `0.99*T+1us`（`tools/exact.py` 现读为准）。**达标即停手。**
// ======================
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/128444: retained the accepted
* real FFT and explicit short-transform layout boundaries, with all notices.
* saffah_cc_v41_agg1 https://duck.ac/submission/128325: inherited decimal
* engine and exact independent carry reconciliation. HintFFT MIT is preserved.
* Idea: Recode each positive digit array transiently into balanced signed
* base1000000 digits by exact carry transfers. Square that same integer with
* lower coefficient bias, round signed FFT coefficients, and recover canonical
* positive digits using floor division and proven signed carry intervals.
* Purpose: Explore whether lower coefficient dynamic range permits a smaller
* transform grid; submit only after complete decimal precision validation.
*/
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/128325: retained base10000 FFT power/carry engine and full source notices.
* TSKY/WithSky HintFFT: used existing iterative short FFT and Complex64X4 input/output permutations and integer rounding; complete MIT license remains below.
* Idea: For transforms up to2048 doubles, explicitly convert interleaved
* input pairs to the kernel layout, call its existing iterative FFT/square/
* inverse, then explicitly permute and round the coefficients. Larger
* transforms retain fused conversion. This enables the sufficient128..2048
* lengths that the recursive fused-conversion interface previously skipped.
* Purpose: Remove fixed4096-point work from small recursive squarings while
* retaining the accepted four-digit precision and exact carry algorithm.
*/
// ===== REFERENCES =====
// [1] duck.ac 账号 **saffah_codex_6a_agg3**（对手），提交 **#128134**
// <https://duck.ac/submission/128134> —— 判题机 14.871071 ms = 1018b 本题当前实时 T。
// 用途：**本发正文 = 该公开件的逐字节副本，一字未动**（本文件 28 行以下全部）；
// 本席只在其**最顶部加**了本引用/思路/目的注释段（注释不进编译面 ⇒ 机器码不变）。
// 其自带的引用链（指向其自己的 #128106 等前代件、以及本账号的历代件）
// **原样保留在下方**，未删段、未改写。
// [2] duck.ac 账号 **saffah_cc_v41_agg1**（本账号），提交 **#126990**
// <https://duck.ac/submission/126990> —— 本账号在役最好件（15.146618 ms）。
// 用途：只作抄件前后对照基线（本发未复制其代码；对手件内已含其署名）。
// 合规：以上均为 duck.ac **公开提交**（题面明示"你提交的代码将会被公开"）；
// 未见独立许可证声明；按 RULES §3 逐项署名 作者账号 + 原提交地址 + 所用内容。
// ★本题在册地雷（必读）：`g_nosort` 的**注释与极性相反** —— 真语义 **`1` = 做排序（正确）**、
// `0` = 跳过排序（WA）。本件是逐字节副本 ⇒ `g_nosort = 1` 原样保留 ✓（`律 1414 §六` /
// `problems/2002/work/r28cb__tmp/_WARNING_g_nosort_polarity.md`）。
// ======================
// ===== 思路 =====
// 【本发 = 对手 #128134 的逐字节底盘（"抄底"工序，对代码零改动）】
// 机理（逐字引自对手件头）：
// mech_1018b.txt
// 目的（测量性）：① 把本账号 mine 从 15.146618 拉回 ≈ T，消掉与对手现引擎的代差；
// ② 立一个与对手逐字节同源的基准件，供后续刀做单变量 A/B。
// ===== 思路（本发补刀 · BASE10K）=====
// 【单变量 = 大数基 10^3 → 10^4】（与 1018c 的 #128320 同一把刀、同一组 4 处改动：BASE 常量、
// max_limbs 分母 3000000→4000000、最高位 limb 输出 3→4 位判定、输出表 DT[1000]→DT[10000]
// 且每 limb 写 4 字节 pos+=4）
// 依据 = 本席相位分解（`r41rx__tmp/mkphase.py`，本机专用、绝不上板；1018b/1018c 同族同构）：
// FFT 三段 77% · 进位环 15% · 其余 8%，且两块都正比于 **limb 数** ⇒ 基 10^4 使 limb 数减半。
// 值闸门（本题 n=1e7）：与底盘（= 对手 #128134 逐字节、板上 Accepted）输出**逐位相同** ——
// P2_LEN 3010300（= floor(1e7*log10 2)+1 ✓ 独立长度核）· P2_HASH aba8690c19e8de19 ✓
// `uint16_t` limb 余数上界 10^4-1 < 2^16 ✓；FFT 系数上界 301030*9999²… 见下
// 精度核（1018b 顶层）：old_size = ceil(5e6*0.30103/4) = 376288 limbs ⇒ 系数上界
// 376288*9999² = 3.76e13 ≪ 2^53 = 9.01e15 ⇒ 整数在 double 里仍**精确可表示**；余量 240× ✓
// 定价诚实备注：本机 A/A 地板对本题族偏大（1018c 同二进制两份 30 对达 ±1.6%）⇒ 本机只在
// ≥2% 量级可信；本发以板面为量具（判题机带 ≈0.15%）。
// 本机 min-of-1 对照：base 70,778,784 → b10k 64,157,708 ticks（**−9.4%**）。
// 目的（试验性）：给"limb 基数"轴在 1018b 上定价（1018c 上已兑现 −5.67% 并转绿）。
// 诚实备注（★口径纠偏）：本发**不是**靠宽支收工 —— `tools/exact.py` 第 225~243 行自己写明：
// 我方任何刷新最好的提交都会成为**最新件** ⇒ 支别翻回**严支**（宽支"只能被持有、不可能被
// 提交抵达"）。故真正的可达线只有 `strict 0.99*T+1us = 14.723360 ms`；抄件 ≈14.871071 仍差 ≈1%。
// 本发只作底盘，收工须另找刀。
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/128108:
Transferred the accepted exact 2^29-chunk integer base case and complete FFT
engine from the smaller sibling. saffah_cc_v41_agg1
https://duck.ac/submission/126990: same problem's latest accepted FFT/carry
algorithm; this implementation retains its output table and minimal clearing.
The factor-specialized carry is inherited from https://duck.ac/submission/127016.
All original notices and the HintFFT MIT license follow in full.
Idea: Keep the proven exact small-exponent engine while restoring LOG_MAX=20
for this problem's full transform size. This removes redundant small FFTs and
preserves the specialized carry path at every larger level. No answer constants.
Purpose: Transfer reduced arithmetic work across the still-unmet pow2 family. */
/* References: saffah_cc_v41_agg1 https://duck.ac/submission/127016:
Copied the accepted bounded-table FFT power-of-two engine, output table,
minimal clearing and specialized carry streams. TSKY/WithSky HintFFT MIT
copyright and license, and all inherited attribution, are preserved below.
Idea: Evaluate recursion exponents at most 1023 directly in base 1000 using
exact multiplication by chunks of 2^29. This replaces several fixed-4096-point
FFT squarings at the bottom with at most 36 short limb passes. Each temporary
is below 1000*2^29, so uint64_t arithmetic is exact. Higher FFT/carry unchanged.
Purpose: Remove small-transform setup and unnecessary convolution work. */
/* ===== REFERENCES =====
* duck.ac 用户 saffah_codex_6a_agg3，提交 #125848（本题排行榜 rank-1，T = 2.095 ms）
* <https://duck.ac/submission/125848>
* （同族 1018 #125852 / 1018a #125851 / 1018b #125849 / 1018c #125848）
* 用途：**本件正文 = 该提交正文的逐字节副本**，其 (1) 十进制大整数 pow2 的
* 递归 in-place 谱平方主流程、(2) 固定工作区 + 四路进位流、(3) 自带的
* 「References」注释块（saffah_codex_maker #125703 及更早血缘、TSKY(WithSky)
* #47426 与 GitHub HintFFT 的 MIT 许可全文）全部原样保留。本件不主张该部分
* 著作权；本席的两处改动见下方思路段。
* duck.ac 用户 saffah_codex_maker，提交 #125703（rank-2）：上件底盘（本席 diff 证实
* 非注释行差异集中在 FFT 长度选择、in-place 谱平方、四路进位与固定工作区四组
* 改动），其引用链随正文注释保留。
* TSKY (WithSky)，<https://github.com/With-Sky/HintFFT> 与 duck.ac #47426：
* AVX2 FFT 内核 HintFFT（MIT 许可）——正文保留其完整版权与许可声明。
* ======================
* ===== 本次思路 =====
* 正式提交（非试验性）。两处改动，均**值中性**：输出对 n=1e6/1e7/1e8 与底盘
* 逐字节对拍（FNV-1a 全等）通过。
*
* [R2F2-OUT] 十进制输出表驱动化：底盘每 limb 用 `x/100`、`x/10%10`、`x%10` 三次
* 除法-取模（-O3 展开为数条乘-移）逐字符写 3 字节；改为 1000 项 × 4 字节查表，
* 每 limb 一次 4 字节载入 + 一次 4 字节存储（第 4 字节落在题面保证的 +1 MB 余量
* 内、且随后被下一次写入覆盖）。
* [R2F2-MZ] 最小清零：底盘每层 `memset(a,0,len*8)` 全缓冲清零；但 [0,old_size)
* 随后立刻被 limb 填充覆盖，[上一层 len, 本层 len) 从未被写过（匿名 mmap 新页
* 恒 0，FFT 从不越过自己那层 len）⇒ 只需清零 [old_size, 上一层 len)。
* 本地 A/B（min-of-R，同机同码）：[OUT]+[MZ] 两刀 n=1e6 −6.9% · n=1e7 −10.6% ·
* n=1e8 −8.0%（本机非判题机，仅作正确性/量级参考；以判题机板面读数为准）。
* [R2F2-CY] 进位环「factor 特化」：底盘的四路进位流主环里 `t = cy[r] + coeff[i]*factor`
* 每 limb 做一次 `imulq factor`（factor ∈ {1,2}，运行期值）——这是热环里除 `mulq`
* 之外唯一吃满 port-1 的乘法。本席把 `factor==1` / `factor==2` 分成两个同构环
* （2 时用两次 `add` 代替乘法，ALU 不吃 port-1）。逐位等价：同一表达式同一求值
* 次序，只是把乘 1/乘 2 常数化。
* 本地：n=1e7 中性 · n=1e8 −2.0% · n=1e9 −3.3%；判题机探针（n=1e8）−0.69%。
* [R2F2-LM] FFT 引擎 LOG_MAX 由 27 下调到本 n 实需值：本题 n=1,000,000 的最大 FFT 长度
* 只有 2^{lm}（= max(4096, ceil2(2*size(2^(n/2))-1))），而 HintFFT 的 sqrt 表会在
* **静态初始化**阶段把 2^16..2^27 全部子表建好（含堆分配），其中本题用不到的
* 2^(17+1)..2^27 部分是纯浪费。下调后逐位等价（同一批公式、同样的表值），
* 本地实测启动阶段 −2.64 ms、n=1e6 端到端 −2.9%。
* ====================== */
/* References:
* saffah_codex_maker, https://duck.ac/submission/125706 (also 125705,125704,125703):
* retained accepted base1000 repeated squaring, in-place workspaces and decimal output.
* TSKY (WithSky), https://duck.ac/submission/47426 and https://github.com/With-Sky/HintFFT:
* inherited AVX2 FFT kernels under MIT; full copyright and permission notice retained.
* Idea: Use in-place spectral squaring for every exponent size, removing the
* duplicate forward FFT in the smaller-size path. Split coefficient carry
* propagation into four independent streams. Initialize each stream exactly
* by propagating lower/upper carry bounds through preceding coefficients until
* the bounds coincide; output digits and final carry remain identical.
* Purpose: Official experiment removing duplicate transforms and shortening the
* serial integer carry chain while retaining the same FFT precision and base.
*/
/*
MIT License
Copyright (c) 2023 With-Sky
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
#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")
namespace hint
{
template <typename T, size_t ALIGN = 64>
class AlignMem{
public:
using Ptr = T *;
using ConstPtr = const T *;
~AlignMem(){
if (ptr){
_mm_free(ptr);}};
AlignMem() : ptr(nullptr), len(0) {}
AlignMem(size_t n) : ptr(reinterpret_cast<Ptr>(_mm_malloc(n * sizeof(T), ALIGN))), len(n) {}
AlignMem(const AlignMem &) = delete;
AlignMem &operator=(const AlignMem &) = delete;
T &operator[](size_t i){
return ptr[i];}
const T &operator[](size_t i) const{
return ptr[i];}
Ptr begin(){
return ptr;}
Ptr end(){
return ptr + len;}
ConstPtr begin() const{
return ptr;}
ConstPtr end() const{
return ptr + len;}
size_t size() const{
return len;}
private:
T *ptr;
size_t len;};
template <typename YMM>
inline void transpose64_2X4(YMM &row0, YMM &row1){
auto t0 = _mm256_unpacklo_pd(__m256d(row0), __m256d(row1));
auto t1 = _mm256_unpackhi_pd(__m256d(row0), __m256d(row1));
row0 = YMM(_mm256_permute2f128_pd(t0, t1, 0x20));
row1 = YMM(_mm256_permute2f128_pd(t0, t1, 0x31));}
template <typename YMM>
inline void transpose64_4X2(YMM &row0, YMM &row1){
auto t0 = _mm256_permute2f128_pd(__m256d(row0), __m256d(row1), 0x20);
auto t1 = _mm256_permute2f128_pd(__m256d(row0), __m256d(row1), 0x31);
row0 = YMM(_mm256_unpacklo_pd(t0, t1));
row1 = YMM(_mm256_unpackhi_pd(t0, t1));}
template <typename YMM>
inline void transpose64_4X4(YMM &row0, YMM &row1, YMM &row2, YMM &row3){
auto t0 = _mm256_unpacklo_pd(__m256d(row0), __m256d(row1));
auto t1 = _mm256_unpackhi_pd(__m256d(row0), __m256d(row1));
auto t2 = _mm256_unpacklo_pd(__m256d(row2), __m256d(row3));
auto t3 = _mm256_unpackhi_pd(__m256d(row2), __m256d(row3));
row0 = YMM(_mm256_permute2f128_pd(t0, t2, 0x20));
row1 = YMM(_mm256_permute2f128_pd(t1, t3, 0x20));
row2 = YMM(_mm256_permute2f128_pd(t0, t2, 0x31));
row3 = YMM(_mm256_permute2f128_pd(t1, t3, 0x31));}
class Float64X4{
public:
using F64 = double;
using F64X4 = Float64X4;
Float64X4() : data(_mm256_setzero_pd()) {}
Float64X4(__m256d in_data) : data(in_data) {}
Float64X4(F64 in_data) : data(_mm256_set1_pd(in_data)) {}
Float64X4(const F64 *in_data) : data(_mm256_load_pd(in_data)) {}
F64X4 operator+(const F64X4 &other) const{
return _mm256_add_pd(data, other.data);}
F64X4 operator-(const F64X4 &other) const{
return _mm256_sub_pd(data, other.data);}
F64X4 operator*(const F64X4 &other) const{
return _mm256_mul_pd(data, other.data);}
F64X4 operator/(const F64X4 &other) const{
return _mm256_div_pd(data, other.data);}
static F64X4 fmadd(const F64X4 &a, const F64X4 &b, const F64X4 &c){
return _mm256_fmadd_pd(a.data, b.data, c.data);}
static F64X4 fmsub(const F64X4 &a, const F64X4 &b, const F64X4 &c){
return _mm256_fmsub_pd(a.data, b.data, c.data);}
template <int N>
F64X4 permute4x64() const{
return _mm256_permute4x64_pd(data, N);}
F64X4 reverse() const{
return permute4x64<0b00011011>();}
void load(const F64 *p){
data = _mm256_load_pd(p);}
void load1(const F64 *p){
data = _mm256_broadcast_sd(p);}
void store(F64 *p) const{
_mm256_store_pd(p, data);}
operator __m256d() const{
return data;}
__m256i toI64X4() const{
constexpr uint64_t mask = (uint64_t(1) << 52) - 1;
constexpr uint64_t offset = (uint64_t(1) << 10) - 1;
const __m256i f64bits = _mm256_castpd_si256(data);
__m256i tail = _mm256_and_si256(f64bits, _mm256_set1_epi64x(mask));
tail = _mm256_or_si256(tail, _mm256_set1_epi64x(mask + 1));
__m256i exp = _mm256_srli_epi64(f64bits, 52);
exp = _mm256_sub_epi64(_mm256_set1_epi64x(offset + 52), exp);
return _mm256_srlv_epi64(tail, exp);}
private:
__m256d data;};
struct Complex64X4{
using C64X4 = Complex64X4;
using F64X4 = Float64X4;
using F64 = double;
Complex64X4() {}
Complex64X4(F64X4 real, F64X4 imag) : real(real), imag(imag) {}
Complex64X4(const F64 *p) : real(p), imag(p + 4) {}
Complex64X4(const F64 *p_real, const F64 *p_imag) : real(p_real), imag(p_imag) {}
C64X4 operator+(const C64X4 &other) const{
return C64X4(real + other.real, imag + other.imag);}
C64X4 operator-(const C64X4 &other) const{
return C64X4(real - other.real, imag - other.imag);}
C64X4 operator*(const F64X4 &other) const{
return C64X4(real * other, imag * other);}
C64X4 mul(const C64X4 &other) const{
const F64X4 ii = imag * other.imag;
const F64X4 ri = real * other.imag;
const F64X4 r = F64X4::fmsub(real, other.real, ii);
const F64X4 i = F64X4::fmadd(imag, other.real, ri);
return C64X4(r, i);}
C64X4 mulConj(const C64X4 &other) const{
const F64X4 ii = imag * other.imag;
const F64X4 ri = real * other.imag;
const F64X4 r = F64X4::fmadd(real, other.real, ii);
const F64X4 i = F64X4::fmsub(imag, other.real, ri);
return C64X4(r, i);}
C64X4 reverse() const{
return C64X4(real.reverse(), imag.reverse());}
void set1(F64 real_in, F64 imag_in){
real = F64X4(real_in);
imag = F64X4(imag_in);}
template <typename T>
void load(const T *p, std::false_type){
this->load(p);}
template <typename T>
void load(const T *p, std::true_type){
this->load(p);
*this = this->toRRIIPermu();}
template <typename T>
void load(const T *p){
real.load(reinterpret_cast<const F64 *>(p));
imag.load(reinterpret_cast<const F64 *>(p) + 4);}
void load1(const F64 *real_p, const F64 *imag_p){
real.load1(real_p);
imag.load1(imag_p);}
template <typename T>
void store(T *p, std::false_type) const{
this->store(p);}
template <typename T>
void store(T *p, std::true_type) const{
this->toRIRIPermu().store(p);}
template <typename T>
void store(T *p) const{
real.store(reinterpret_cast<F64 *>(p));
imag.store(reinterpret_cast<F64 *>(p) + 4);}
C64X4 toRIRIPermu() const{
C64X4 res = *this;
transpose64_2X4(res.real, res.imag);
return res;}
C64X4 toRRIIPermu() const{
C64X4 res = *this;
transpose64_4X2(res.real, res.imag);
return res;}
C64X4 transToI64(std::false_type) const{
return *this;}
C64X4 transToI64(std::true_type) const{
constexpr int64_t F1_2 = 4602678819172646912;
auto F1_2X4 = F64X4(__m256d(_mm256_set1_epi64x(F1_2)));
auto real_i64 = (real + F1_2X4).toI64X4();
auto imag_i64 = (imag + F1_2X4).toI64X4();
return C64X4(__m256d(real_i64), __m256d(imag_i64));}
F64X4 real, imag;};
using Float32 = float;
using Float64 = double;
constexpr Float64 HINT_PI = 3.141592653589793238462643;
constexpr Float64 HINT_2PI = HINT_PI * 2;
constexpr Float64 COS_PI_8 = 0.707106781186547524400844;
template <typename T>
constexpr T int_floor2(T n){
constexpr int bits = sizeof(n) * 8;
for (int i = 1; i < bits; i *= 2){
n |= (n >> i);}
return (n >> 1) + 1;}
template <typename T>
constexpr T int_ceil2(T n){
constexpr int bits = sizeof(n) * 8;
n--;
for (int i = 1; i < bits; i *= 2){
n |= (n >> i);}
return n + 1;}
template <typename IntTy>
constexpr bool is_2pow(IntTy n){
return n != 0 && (n & (n - 1)) == 0;}
template <typename T>
constexpr int hint_log2(T n){
constexpr int bits = sizeof(n) * 8;
int l = -1, r = bits;
while ((l + 1) != r){
int mid = (l + r) / 2;
if ((T(1) << mid) > n){
r = mid;}
else{
l = mid;}}
return l;}
constexpr uint32_t crc32(const char *str){
uint32_t crc = 0xFFFFFFFF;
while (*str != '\0'){
crc ^= *str;
for (int i = 0; i < 8; ++i){
crc = (crc >> 1) ^ (0 - (crc & 1)) & 0xEDB88320;}
str++;}
return ~crc;}
namespace transform{
template <typename T>
inline void transform2(T &sum, T &diff){
T temp0 = sum, temp1 = diff;
sum = temp0 + temp1;
diff = temp0 - temp1;}
namespace fft{
using F64 = Float64;
using C64 = std::complex<F64>;
using F64X4 = Float64X4;
using C64X4 = Complex64X4;
template <typename Float, size_t OMEGA_LEN>
class TableFix{
alignas(64) Float table[OMEGA_LEN * 2];
public:
TableFix(size_t theta_divider, size_t factor, size_t stride){
const Float theta = -HINT_2PI * factor / theta_divider;
for (size_t begin = 0, index = 0; begin < OMEGA_LEN * 2; begin += stride * 2)
{
for (size_t j = 0; j < stride; j++, index++)
{
table[begin + j] = std::cos(theta * index);
table[begin + j + stride] = std::sin(theta * index);}}}
constexpr const Float &operator[](size_t index) const{
return table[index];}};
void initOmegaX4(F64 *arr, size_t fft_len, int table_len, int factor){
table_len /= 4;
const F64 theta = -HINT_2PI * factor / fft_len;
auto arrx4 = reinterpret_cast<C64X4 *>(arr);
arr[0] = 1, arr[4] = 0;
arr[1] = std::cos(theta), arr[5] = std::sin(theta);
arr[2] = std::cos(theta * 2), arr[6] = std::sin(theta * 2);
arr[3] = std::cos(theta * 3), arr[7] = std::sin(theta * 3);
for (size_t begin = 1; begin < table_len; begin *= 2){
size_t nth = begin * 4;
C64X4 unit;
unit.set1(std::cos(theta * nth), std::sin(theta * nth));
for (size_t i = 0; i < begin; i++)
{
arrx4[begin + i] = arrx4[i].mul(unit);}}}
template <typename Float, int LOG_BEGIN, int LOG_END, int DIV>
class TableFixMulti{
static_assert(LOG_END >= LOG_BEGIN);
static_assert(is_2pow(DIV));
static constexpr size_t TABLE_CPX_LEN = (size_t(1) << (LOG_END + 1)) / DIV;
alignas(64) Float table[TABLE_CPX_LEN * 2];
public:
TableFixMulti(size_t factor, size_t stride = 4){
initBottomUp(factor, stride);}
void initBottomUp(size_t factor, size_t stride){
static_assert(std::is_same<Float, Float64>::value);
size_t len = size_t(1) << LOG_BEGIN, cpx_len = len / DIV;
auto it = getBeginLog(LOG_BEGIN);
initOmegaX4(it, len, cpx_len, factor);
for (int log_len = LOG_BEGIN + 1; log_len <= LOG_END; log_len++)
{
len = size_t(1) << log_len, cpx_len = len / DIV;
Float theta = -HINT_2PI * factor / len;
auto it = getBeginLog(log_len), it_last = getBeginLog(log_len - 1);
C64X4 unit(std::cos(theta), std::sin(theta));
for (auto end = it + cpx_len * 2; it < end; it += 16, it_last += 8)
{
C64X4 omega0, omega1;
omega0.load(it_last);
omega1 = omega0.mul(unit);
transpose64_2X4(omega0.real, omega1.real);
transpose64_2X4(omega0.imag, omega1.imag);
omega0.store(it), omega1.store(it + 8);}}}
constexpr const Float *getBeginLog(int log_rank) const{
return getBegin(size_t(1) << log_rank);}
constexpr Float *getBeginLog(int log_rank){
return getBegin(size_t(1) << log_rank);}
constexpr const Float *getBegin(size_t rank) const{
return &table[rank * 2 / DIV];}
constexpr Float *getBegin(size_t rank){
return &table[rank * 2 / DIV];}};
template <int CACHE_LOG_LEN>
class FFTSqrtTableC64X4{
public:
using F64 = double;
using C64 = std::complex<double>;
using C64X4 = hint::Complex64X4;
static constexpr size_t CACHE_LEN = size_t(1) << CACHE_LOG_LEN;
static constexpr size_t MASK = CACHE_LEN - 1;
static constexpr size_t C4_COUNT = sizeof(C64X4) / sizeof(C64);
~FFTSqrtTableC64X4(){
if (high)
{
delete[] high;}}
FFTSqrtTableC64X4() {}
FFTSqrtTableC64X4(size_t fft_len, int len_div, int factor){
init(fft_len, len_div, factor);}
void init(size_t fft_len, int len_div, int factor){
size_t table_len = fft_len / len_div;
size_t low_len = CACHE_LEN * C4_COUNT, high_len = table_len / low_len;
if (high != nullptr)
{
delete[] high;}
high = new C64[high_len];
auto p = reinterpret_cast<F64 *>(&low[0]);
initOmegaX4(p, fft_len, low_len, factor);
const F64 theta = -HINT_2PI * factor / fft_len;
high[0] = C64(1, 0);
for (size_t begin = 1; begin < high_len; begin *= 2)
{
C64 unit = std::polar<F64>(1.0, theta * begin * low_len);
for (size_t i = 0; i < begin; i++)
{
high[i + begin] = high[i] * unit;}}}
C64X4 operator[](size_t i) const{
C64X4 hi;
auto p = reinterpret_cast<const F64 *>(&high[i >> CACHE_LOG_LEN]);
hi.load1(p, p + 1);
return low[i & MASK].mul(hi);}
private:
alignas(64) C64X4 low[CACHE_LEN];
C64 *high = nullptr;};
template <int DIV, int LOG_BEGIN, int LOG_MAX, int CACHE_LOG_LEN>
class FFTTableSqrt{
using TableLong = FFTSqrtTableC64X4<CACHE_LOG_LEN>;
static constexpr size_t SHORT_LEN = size_t(1) << LOG_BEGIN;
static constexpr size_t TABLE_LEN = LOG_MAX - LOG_BEGIN + 1;
public:
FFTTableSqrt(int factor){
for (int i = 0; i < TABLE_LEN; i++)
{
size_t fft_len = SHORT_LEN << i;
table[i].init(fft_len, DIV, factor);}}
const TableLong &operator[](int log_len) const{
log_len -= LOG_BEGIN;
return table[log_len];}
TableLong &operator[](int log_len){
log_len -= LOG_BEGIN;
return table[log_len];}
private:
TableLong table[TABLE_LEN];};
struct FFT{
template <typename Float>
static void trans2MulI(Float &r0, Float &i0, Float &r1, Float &i1){
auto temp = r1;
r1 = r0 + i1;
r0 = r0 - i1;
i1 = i0 - temp;
i0 = i0 + temp;}
template <typename Float>
static void trans2MulNegI(Float &r0, Float &i0, Float &r1, Float &i1){
auto temp = r1;
r1 = r0 - i1;
r0 = r0 + i1;
i1 = i0 + temp;
i0 = i0 - temp;}
template <typename Float>
static void dif4(Float &r0, Float &i0, Float &r1, Float &i1, Float &r2, Float &i2, Float &r3, Float &i3){
difSplit(r0, i0, r1, i1, r2, i2, r3, i3);
transform2(r0, r1);
transform2(i0, i1);}
template <typename Float>
static void idit4(Float &r0, Float &i0, Float &r1, Float &i1, Float &r2, Float &i2, Float &r3, Float &i3){
transform2(r0, r1);
transform2(i0, i1);
iditSplit(r0, i0, r1, i1, r2, i2, r3, i3);}
template <typename Float>
static void difSplit(Float &r0, Float &i0, Float &r1, Float &i1, Float &r2, Float &i2, Float &r3, Float &i3){
transform2(r0, r2);
transform2(i0, i2);
transform2(r1, r3);
transform2(i1, i3);
trans2MulNegI(r2, i2, r3, i3);}
template <typename Float>
static void iditSplit(Float &r0, Float &i0, Float &r1, Float &i1, Float &r2, Float &i2, Float &r3, Float &i3){
transform2(r2, r3);
transform2(i2, i3);
transform2(r0, r2);
transform2(i0, i2);
trans2MulI(r1, i1, r3, i3);}};
struct FFTAVX : public FFT{
static constexpr int LOG_SHORT = 10, LOG_MID = 25 - crc32(statement) % 15, LOG_MAX = 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++) \
{ \
C64X4 c0, c1, c2, c3, omega1, omega2; \
c0.load(it0, FromRIRI{}), c1.load(it1, FromRIRI{}), c2 = c0,c3 = c1; \
transform2(c0,c1); \
trans2MulNegI(c2.real,c2.imag,c3.real,c3.imag); \
omega1 = table[indx], c2 = c2.mul(omega1); \
omega2 = omega1.mul(omega1), c1 = c1.mul(omega2); \
c3 = c3.mul(omega2.mul(omega1)); \
c0.store(it0), c1.store(it1), c2.store(it2), c3.store(it3); \
} \
for (auto end = it1; (!FROM_RIRI_PERM) && it0 < end; it0 += 8, it1 += 8, it2 += 8, it3 += 8, indx++) \
{ \
C64X4 c0 = it0, c1 = it1, c2 = it2, c3 = it3, omega1, omega2; \
dif4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag); \
omega1 = table[indx], c2 = c2.mul(omega1); \
omega2 = omega1.mul(omega1), c1 = c1.mul(omega2); \
c3 = c3.mul(omega2.mul(omega1)); \
c0.store(it0), c1.store(it1), c2.store(it2), c3.store(it3); \
} \
dif_func(in_out, stride); \
dif_func(in_out + stride, stride); \
dif_func(in_out + stride * 2, stride); \
dif_func(in_out + stride * 3, stride); \
} while (0)
#define iditLayer(idit_func, in_out, stride, table) \
do \
{ \
idit_func(in_out, stride); \
idit_func(in_out + stride, stride); \
idit_func(in_out + stride * 2, stride); \
idit_func(in_out + stride * 3, stride); \
auto it0 = in_out, it1 = in_out + stride, it2 = it1 + stride, it3 = it2 + stride; \
size_t indx = 0; \
for (auto end = it1; it0 < end; it0 += 8, it1 += 8, it2 += 8, it3 += 8, indx++) \
{ \
C64X4 c0 = it0, c1 = it1, c2 = it2, c3 = it3, omega1, omega2; \
omega1 = table[indx], c2 = c2.mulConj(omega1); \
omega2 = omega1.mul(omega1), c1 = c1.mulConj(omega2); \
c3 = c3.mulConj(omega2.mul(omega1)); \
idit4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag); \
c0 = c0.transToI64(ToI64{}), c1 = c1.transToI64(ToI64{}), c2 = c2.transToI64(ToI64{}), c3 = c3.transToI64(ToI64{}); \
c0.store(it0, ToRIRI{}), c1.store(it1, ToRIRI{}), c2.store(it2, ToRIRI{}), c3.store(it3, ToRIRI{}); \
} \
} while (0)
template <bool FROM_RIRI_PERM = false>
static void difRecMid(Float64 in_out[], size_t float_len){
const size_t fft_len = float_len / 2;
if (fft_len <= SHORT_LEN)
{
difIter(in_out, float_len);
return;}
using FromRIRI = std::integral_constant<bool, FROM_RIRI_PERM>;
auto table1 = reinterpret_cast<const C64X4 *>(multi_table_1.getBegin(fft_len));
const size_t stride = float_len / 4;
difLayer(difRecMid, in_out, stride, table1);}
template <bool TO_RIRI_PERM = false, bool TO_INT64 = false>
static void iditRecMid(Float64 in_out[], size_t float_len){
const size_t fft_len = float_len / 2;
if (fft_len <= SHORT_LEN)
{
iditIter(in_out, float_len);
return;}
using ToRIRI = std::integral_constant<bool, TO_RIRI_PERM>;
using ToI64 = std::integral_constant<bool, TO_INT64>;
const size_t stride = float_len / 4;
auto table1 = reinterpret_cast<const C64X4 *>(multi_table_1.getBegin(fft_len));
iditLayer(iditRecMid, in_out, stride, table1);}
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)
{
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){
const int64_t neg=x<0;return (x+neg)/(int64_t)BASE-neg;
}
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;
for(size_t i=0;i<m;i++){
unsigned t=v[i]+transfer;transfer=t>=BASE/2;
a[i]=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{});}
}else{
hint::transform::fft::FFTAVX::difRecLong<true>(a,len);
hint::transform::fft::real_dot_binrev4(a,a,len);
hint::transform::fft::FFTAVX::iditRecLong<true,false>(a,len);
}
for(size_t i=0;i<len;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));
}
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;
}
}
if(factor==1){
for(size_t i=0;i<step;i++){
#pragma GCC unroll 4
for(unsigned r=0;r<4;r++){
int64_t t=cy[r]+coeff[r*step+i],q=floor_radix(t);
v[r*step+i]=uint32_t(t-q*BASE);cy[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];
int64_t t=cy[r]+c+c,q=floor_radix(t);
v[r*step+i]=uint32_t(t-q*BASE);cy[r]=q;
}
}
}
done=4*step;carry=cy[3];
}
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;
}
for(size_t i=size-1;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 | 8.051 ms | 9 MB + 36 KB | Accepted | Score: 100 | 显示更多 |