// duck.ac 1004 -- exact product of two 1,000,000-digit decimals.
// base 10^8 limbs (125000 each), transform length 2^18, three small NTT primes
// (40370177 / 41680897 / 42729473, all = 1 mod 2^18), AVX2 Shoup NTT in 64-bit
// lanes with fully lazy range handling, Garner CRT, base-1e8 carry, decimal out.
// Compiles with plain -O2: all intrinsic code sits in target("avx2") functions.
typedef unsigned int u32;
typedef unsigned long long u64;
#include <immintrin.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
static inline u32 mpw(u32 a, u32 e, u32 p) {
u64 r = 1, b = a;
while (e) { if (e & 1) r = r * b % p; b = b * b % p; e >>= 1; }
return (u32)r;
}
static const int NL = 125000;
static const int NN = 1 << 18;
static const int NC = 2 * NL - 1;
static const u32 P1 = 40370177u, P2 = 41680897u, P3 = 42729473u;
// ---------------- twiddle tables + transforms ----------------
struct Ntt {
u32 P, G;
u32 *W, *WS, *Wi, *WSi;
int offs[32];
u32 NINV, NINVS;
int prepared;
void prepare(int n) {
if (prepared == n) return;
prepared = n;
if (!W) {
static u32 *pool = 0;
if (!pool) pool = (u32 *)_mm_malloc(4 * (size_t)(1 << 18) * 4 + 256, 64);
W = pool; WS = pool + (1 << 18); Wi = pool + 2 * (1 << 18); WSi = pool + 3 * (1 << 18);
}
u32 root = mpw(G, (P - 1) / n, P), rooti = mpw(root, P - 2, P);
const int top = n >> 1;
const int base = n - 1 - top;
offs[__builtin_ctz(n)] = base;
{
u32 r4 = (u32)((u64)root * root % P * root % P * root % P);
u32 ri4 = (u32)((u64)rooti * rooti % P * rooti % P * rooti % P);
u64 c0 = 1, c1 = root, c2 = (u64)root * root % P, c3 = (u64)c2 * root % P;
u64 d0 = 1, d1 = rooti, d2 = (u64)rooti * rooti % P, d3 = (u64)d2 * rooti % P;
u32 *w = W + base, *ws = WS + base, *wi = Wi + base, *wsi = WSi + base;
int j = 0;
for (; j + 4 <= top; j += 4) {
w[j] = (u32)c0; ws[j] = (u32)((c0 << 32) / P);
w[j+1] = (u32)c1; ws[j+1] = (u32)((c1 << 32) / P);
w[j+2] = (u32)c2; ws[j+2] = (u32)((c2 << 32) / P);
w[j+3] = (u32)c3; ws[j+3] = (u32)((c3 << 32) / P);
wi[j] = (u32)d0; wsi[j] = (u32)((d0 << 32) / P);
wi[j+1] = (u32)d1; wsi[j+1] = (u32)((d1 << 32) / P);
wi[j+2] = (u32)d2; wsi[j+2] = (u32)((d2 << 32) / P);
wi[j+3] = (u32)d3; wsi[j+3] = (u32)((d3 << 32) / P);
c0 = c0 * r4 % P; c1 = c1 * r4 % P; c2 = c2 * r4 % P; c3 = c3 * r4 % P;
d0 = d0 * ri4 % P; d1 = d1 * ri4 % P; d2 = d2 * ri4 % P; d3 = d3 * ri4 % P;
}
for (; j < top; j++) {
w[j] = (u32)c0; ws[j] = (u32)((c0 << 32) / P);
wi[j] = (u32)d0; wsi[j] = (u32)((d0 << 32) / P);
c0 = c0 * root % P; d0 = d0 * rooti % P;
}
}
int off = base;
for (int len = n >> 1; len >= 2; len >>= 1) {
int half = len >> 1;
off -= half;
offs[__builtin_ctz(len)] = off;
const u32 *sw = W + off + half, *sws = WS + off + half;
const u32 *si = Wi + off + half, *sis = WSi + off + half;
for (int j = 0; j < half; j++) {
W[off + j] = sw[2 * j]; WS[off + j] = sws[2 * j];
Wi[off + j] = si[2 * j]; WSi[off + j] = sis[2 * j];
}
}
u32 ninv = mpw((u32)n, P - 2, P);
NINV = ninv; NINVS = (u32)(((u64)ninv << 32) / P);
}
__attribute__((target("avx2")))
void dif(u64 *a, int n) const {
const __m256i Pv = _mm256_set1_epi64x((long long)P);
const __m256i P2 = _mm256_set1_epi64x(2 * (long long)P);
const __m256i P2m1 = _mm256_set1_epi64x(2 * (long long)P - 1);
for (int len = n; len >= 2; len >>= 1) {
int half = len >> 1;
const u32 *w = W + offs[__builtin_ctz(len)], *s = WS + offs[__builtin_ctz(len)];
if (len == 2) {
for (int i = 0; i < n; i += 2) {
u64 u = a[i], v = a[i + 1];
u64 sum = u + v;
sum -= (2 * (u64)P) & (u64)(-(long long)(sum >= 2 * (u64)P));
a[i] = sum; a[i + 1] = u + 2 * (u64)P - v;
}
continue;
}
for (int i = 0; i < n; i += len) {
u64 *lo = a + i, *hi = a + i + half;
int j = 0;
if (half >= 8) {
for (; j + 4 <= half; j += 4) {
__m256i u = _mm256_loadu_si256((const __m256i *)(lo + j));
__m256i v = _mm256_loadu_si256((const __m256i *)(hi + j));
__m256i wv = _mm256_cvtepu32_epi64(_mm_loadu_si128((const __m128i *)(w + j)));
__m256i sv = _mm256_cvtepu32_epi64(_mm_loadu_si128((const __m128i *)(s + j)));
__m256i sum = _mm256_add_epi64(u, v);
__m256i cmp = _mm256_cmpgt_epi32(sum, P2m1);
sum = _mm256_sub_epi64(sum, _mm256_and_si256(cmp, P2));
__m256i d = _mm256_sub_epi64(_mm256_add_epi64(u, P2), v);
__m256i r = _mm256_sub_epi64(_mm256_mul_epu32(d, wv),
_mm256_mul_epu32(_mm256_srli_epi64(_mm256_mul_epu32(d, sv), 32), Pv));
_mm256_storeu_si256((__m256i *)(lo + j), sum);
_mm256_storeu_si256((__m256i *)(hi + j), r);
}
}
for (; j < half; j++) {
u64 u = lo[j], v = hi[j];
u64 sum = u + v;
sum -= (2 * (u64)P) & (u64)(-(long long)(sum >= 2 * (u64)P));
u64 d = u + 2 * (u64)P - v;
u64 q = (u64)(u32)d * s[j] >> 32;
lo[j] = sum; hi[j] = (u64)(u32)d * w[j] - q * P;
}
}
}
}
__attribute__((target("avx2")))
void dit(u64 *a, int n) const {
const __m256i Pv = _mm256_set1_epi64x((long long)P);
const __m256i P2 = _mm256_set1_epi64x(2 * (long long)P);
for (int len = 2; len <= n; len <<= 1) {
int half = len >> 1;
const u32 *w = Wi + offs[__builtin_ctz(len)], *s = WSi + offs[__builtin_ctz(len)];
if (len == 2) {
for (int i = 0; i < n; i += 2) {
u64 u = a[i], v = a[i + 1];
a[i] = u + v; a[i + 1] = u + 2 * (u64)P - v;
}
continue;
}
for (int i = 0; i < n; i += len) {
u64 *lo = a + i, *hi = a + i + half;
int j = 0;
if (half >= 8) {
for (; j + 4 <= half; j += 4) {
__m256i u = _mm256_loadu_si256((const __m256i *)(lo + j));
__m256i v = _mm256_loadu_si256((const __m256i *)(hi + j));
__m256i wv = _mm256_cvtepu32_epi64(_mm_loadu_si128((const __m128i *)(w + j)));
__m256i sv = _mm256_cvtepu32_epi64(_mm_loadu_si128((const __m128i *)(s + j)));
__m256i r = _mm256_sub_epi64(_mm256_mul_epu32(v, wv),
_mm256_mul_epu32(_mm256_srli_epi64(_mm256_mul_epu32(v, sv), 32), Pv));
_mm256_storeu_si256((__m256i *)(lo + j), _mm256_add_epi64(u, r));
_mm256_storeu_si256((__m256i *)(hi + j), _mm256_sub_epi64(_mm256_add_epi64(u, P2), r));
}
}
for (; j < half; j++) {
u64 u = lo[j], v = hi[j];
u64 q = (u64)(u32)v * s[j] >> 32;
u64 r = (u64)(u32)v * w[j] - q * P;
lo[j] = u + r; hi[j] = u + 2 * (u64)P - r;
}
}
}
const __m256i nv = _mm256_set1_epi64x((long long)NINV);
const __m256i ns = _mm256_set1_epi64x((long long)NINVS);
int i = 0;
for (; i + 4 <= n; i += 4) {
__m256i x = _mm256_loadu_si256((const __m256i *)(a + i));
__m256i r = _mm256_sub_epi64(_mm256_mul_epu32(x, nv),
_mm256_mul_epu32(_mm256_srli_epi64(_mm256_mul_epu32(x, ns), 32), Pv));
_mm256_storeu_si256((__m256i *)(a + i), r);
}
for (; i < n; i++) {
u64 q = (u64)(u32)a[i] * NINVS >> 32;
a[i] = (u64)(u32)a[i] * NINV - q * P;
}
}
};
// ---------------- helpers ----------------
__attribute__((target("avx2")))
static inline __m256d to4d(__m256i x) {
__m256i p = _mm256_permutevar8x32_epi32(x, _mm256_setr_epi32(0, 2, 4, 6, 1, 3, 5, 7));
return _mm256_cvtepi32_pd(_mm256_castsi256_si128(p));
}
__attribute__((target("avx2")))
static inline __m256i to4q(__m256d d) {
return _mm256_cvtepu32_epi64(_mm256_cvttpd_epi32(d));
}
__attribute__((target("avx2")))
static inline __m256i red32(__m256i x, __m256i Pv, __m256d invP, __m256i Pm1) {
__m256d dx = to4d(x);
__m256i q = to4q(_mm256_mul_pd(dx, invP));
__m256i r = _mm256_sub_epi64(x, _mm256_mul_epu32(q, Pv));
__m256i m = _mm256_cmpgt_epi64(r, Pm1);
return _mm256_sub_epi64(r, _mm256_and_si256(m, Pv));
}
// ---------------- driver ----------------
static u32 *resv[3];
static u64 la[NL], lb[NL];
static u64 coef[NC + 8];
static u32 t2a[NC + 8];
static u64 *gW, *gB2;
__attribute__((target("avx2")))
static void one_prime(const Ntt &T, u32 P, u32 *out) {
const __m256i Pv = _mm256_set1_epi64x((long long)P);
const __m256i Pm1 = _mm256_sub_epi64(Pv, _mm256_set1_epi64x(1));
const __m256d invP = _mm256_set1_pd(1.0 / (double)P);
for (int i = 0; i < NN; i++) gW[i] = 0;
for (int i = NL; i < NN; i++) gB2[i] = 0;
for (int i = 0; i < NL; i++) {
u64 x = la[i]; if (x >= P) { x -= P; if (x >= P) x -= P; }
u64 y = lb[i]; if (y >= P) { y -= P; if (y >= P) y -= P; }
gW[i] = x; gB2[i] = y;
}
T.dif(gW, NN);
T.dif(gB2, NN);
for (int i = 0; i < NN; i += 4) {
__m256i x = _mm256_loadu_si256((const __m256i *)(gW + i));
__m256i y = _mm256_loadu_si256((const __m256i *)(gB2 + i));
__m256i xr = red32(x, Pv, invP, Pm1), yr = red32(y, Pv, invP, Pm1);
__m256d dx = to4d(xr), dy = to4d(yr);
__m256i q = to4q(_mm256_mul_pd(_mm256_mul_pd(dx, dy), invP));
__m256i pr = _mm256_mul_epu32(xr, yr);
__m256i r = _mm256_sub_epi64(pr, _mm256_mul_epu32(q, Pv));
__m256i neg = _mm256_cmpgt_epi64(_mm256_setzero_si256(), r);
r = _mm256_add_epi64(r, _mm256_and_si256(neg, Pv));
__m256i mm = _mm256_cmpgt_epi64(r, Pm1);
r = _mm256_sub_epi64(r, _mm256_and_si256(mm, Pv));
_mm256_storeu_si256((__m256i *)(gW + i), r);
}
T.dit(gW, NN);
for (int i = 0; i < NC; i += 4) {
__m256i x = _mm256_loadu_si256((const __m256i *)(gW + i));
__m256i r = red32(x, Pv, invP, Pm1);
__m256i pk = _mm256_permutevar8x32_epi32(r, _mm256_setr_epi32(0, 2, 4, 6, 1, 3, 5, 7));
_mm_storeu_si128((__m128i *)(out + i), _mm256_castsi256_si128(pk));
}
}
__attribute__((target("avx2")))
static void garner_all(void) {
const u32 inv1 = mpw(P1 % P2, P2 - 2, P2);
const u64 M2v = (u64)P1 * P2;
const u32 inv2 = mpw((u32)(M2v % P3), P3 - 2, P3);
const __m256i P2v = _mm256_set1_epi64x((long long)P2);
const __m256i P3v = _mm256_set1_epi64x((long long)P3);
const __m256i P1v = _mm256_set1_epi64x((long long)P1);
const __m256i P2m1 = _mm256_sub_epi64(P2v, _mm256_set1_epi64x(1));
const __m256i P3m1 = _mm256_sub_epi64(P3v, _mm256_set1_epi64x(1));
const __m256d invP2 = _mm256_set1_pd(1.0 / (double)P2);
const __m256d invP3 = _mm256_set1_pd(1.0 / (double)P3);
const __m256i inv1v = _mm256_set1_epi64x((long long)inv1);
const __m256i inv2v = _mm256_set1_epi64x((long long)inv2);
for (int i = 0; i < NC; i += 4) {
__m256i r1 = _mm256_cvtepu32_epi64(_mm_loadu_si128((const __m128i *)(resv[0] + i)));
__m256i r2 = _mm256_cvtepu32_epi64(_mm_loadu_si128((const __m128i *)(resv[1] + i)));
__m256i r3 = _mm256_cvtepu32_epi64(_mm_loadu_si128((const __m128i *)(resv[2] + i)));
__m256i d = _mm256_sub_epi64(_mm256_add_epi64(r2, P2v), r1);
__m256i mm = _mm256_cmpgt_epi64(d, P2m1);
d = _mm256_sub_epi64(d, _mm256_and_si256(mm, P2v));
__m256d dd = to4d(d);
__m256i q = to4q(_mm256_mul_pd(_mm256_mul_pd(dd, _mm256_set1_pd((double)inv1)), invP2));
__m256i t1 = _mm256_sub_epi64(_mm256_mul_epu32(d, inv1v), _mm256_mul_epu32(q, P2v));
__m256i ng = _mm256_cmpgt_epi64(_mm256_setzero_si256(), t1);
t1 = _mm256_add_epi64(t1, _mm256_and_si256(ng, P2v));
__m256i mg = _mm256_cmpgt_epi64(t1, P2m1);
t1 = _mm256_sub_epi64(t1, _mm256_and_si256(mg, P2v));
__m256i y = _mm256_add_epi64(r1, _mm256_mul_epu32(t1, P1v));
_mm256_storeu_si256((__m256i *)(coef + i), y);
u64 yv[4], y3[4];
_mm256_storeu_si256((__m256i *)yv, y);
for (int k = 0; k < 4; k++) {
u64 yy = yv[k];
u64 qq = (u64)((double)yy * (1.0 / (double)P3));
u64 rr = yy - qq * (u64)P3;
while (rr >= P3) rr -= P3;
y3[k] = rr;
}
__m256i y3v = _mm256_loadu_si256((const __m256i *)y3);
__m256i d3 = _mm256_sub_epi64(_mm256_add_epi64(r3, P3v), y3v);
__m256i m3 = _mm256_cmpgt_epi64(d3, P3m1);
d3 = _mm256_sub_epi64(d3, _mm256_and_si256(m3, P3v));
__m256d d3d = to4d(d3);
__m256i q3 = to4q(_mm256_mul_pd(_mm256_mul_pd(d3d, _mm256_set1_pd((double)inv2)), invP3));
__m256i t2 = _mm256_sub_epi64(_mm256_mul_epu32(d3, inv2v), _mm256_mul_epu32(q3, P3v));
__m256i n3 = _mm256_cmpgt_epi64(_mm256_setzero_si256(), t2);
t2 = _mm256_add_epi64(t2, _mm256_and_si256(n3, P3v));
__m256i m4 = _mm256_cmpgt_epi64(t2, P3m1);
t2 = _mm256_sub_epi64(t2, _mm256_and_si256(m4, P3v));
u64 t2v[4];
_mm256_storeu_si256((__m256i *)t2v, t2);
t2a[i] = (u32)t2v[0]; t2a[i + 1] = (u32)t2v[1];
t2a[i + 2] = (u32)t2v[2]; t2a[i + 3] = (u32)t2v[3];
}
}
static char obuf[2100000];
static u32 dig[NC + 16];
int main() {
gW = (u64 *)_mm_malloc(8 * (size_t)NN + 256, 64);
gB2 = (u64 *)_mm_malloc(8 * (size_t)NN + 256, 64);
for (int i = 0; i < 3; i++) resv[i] = (u32 *)malloc(4 * (size_t)(NC + 16));
static char ibuf[2 << 20];
size_t len = fread(ibuf, 1, sizeof(ibuf) - 1, stdin);
ibuf[len] = 0;
size_t p = 0;
while (p < len && (ibuf[p] < '0' || ibuf[p] > '9')) p++;
size_t pa = p; while (p < len && ibuf[p] >= '0' && ibuf[p] <= '9') p++;
size_t ea = p;
while (p < len && (ibuf[p] < '0' || ibuf[p] > '9')) p++;
size_t pb = p; while (p < len && ibuf[p] >= '0' && ibuf[p] <= '9') p++;
size_t eb = p;
for (int side = 0; side < 2; side++) {
const char *s = ibuf + (side ? pb : pa);
int L = (int)(side ? eb - pb : ea - pa);
u64 *out = side ? lb : la;
int n = 0, i = L;
while (i > 0) {
int j = i - 8; if (j < 0) j = 0;
u64 v = 0;
for (int k = j; k < i; k++) v = v * 10 + (u64)(s[k] - '0');
out[n++] = v;
i = j;
}
}
static Ntt T1, T2, T3;
T1.P = P1; T1.G = 3; T2.P = P2; T2.G = 5; T3.P = P3; T3.G = 3;
T1.W = T1.WS = T1.Wi = T1.WSi = 0; T1.prepared = 0;
T2.W = T2.WS = T2.Wi = T2.WSi = 0; T2.prepared = 0;
T3.W = T3.WS = T3.Wi = T3.WSi = 0; T3.prepared = 0;
T1.prepare(NN);
one_prime(T1, P1, resv[0]);
T2.prepare(NN);
one_prime(T2, P2, resv[1]);
T3.prepare(NN);
one_prime(T3, P3, resv[2]);
garner_all();
const u64 B = 100000000ull;
const u64 M2v = (u64)P1 * P2;
const u64 m1 = M2v / B, m0 = M2v % B;
int nd = 0;
u64 carry = 0;
for (int i = 0; i < NC; i++) {
u64 v = coef[i] + m0 * (u64)t2a[i] + carry;
dig[nd] = (u32)(v % B);
carry = v / B + m1 * (u64)t2a[i];
nd++;
}
while (carry) { dig[nd++] = (u32)(carry % B); carry /= B; }
while (nd > 1 && dig[nd - 1] == 0) nd--;
char *out = obuf;
int o = sprintf(out, "%u", dig[nd - 1]);
for (int i = nd - 2; i >= 0; i--) o += sprintf(out + o, "%08u", dig[i]);
out[o++] = '\n';
fwrite(out, 1, o, stdout);
return 0;
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 75.817 ms | 22 MB + 332 KB | Accepted | Score: 100 | 显示更多 |