結果
| 問題 | No.3333 Consecutive Power Sum (Large) |
| コンテスト | |
| ユーザー |
 Vi24E
|
| 提出日時 | 2026-08-04 02:38:31 |
| 言語 | C++17 (gcc 15.2.0 + boost 1.90.0) |
| 結果 |
AC
|
| 実行時間 | 47 ms / 10,000 ms |
| + 180µs | |
| コード長 | 18,433 bytes |
| 記録 | |
| コンパイル時間 | 2,969 ms |
| コンパイル使用メモリ | 292,464 KB |
| 実行使用メモリ | 20,596 KB |
| 最終ジャッジ日時 | 2026-08-04 02:38:39 |
| 合計ジャッジ時間 | 7,916 ms |
|
ジャッジサーバーID （参考情報） |
judge3_0 / judge1_0 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 3 |
| other | AC * 63 |
ソースコード
// This code was generated by AI.
#include <bits/stdc++.h>
using namespace std;
using u64 = uint64_t;
using u128 = __uint128_t;
static u128 parse_u128(const string& s){ u128 x=0; for(char c:s) x=x*10+(c-'0'); return x; }
static u64 isqrt_u128(u128 n){
u64 x=(u64)sqrt((long double)n);
while((u128)x*x>n) --x;
while((u128)(x+1)*(x+1)<=n) ++x;
return x;
}
struct Mont128 {
u128 mod, ninv, r1, r2;
static inline u128 mul_hi(u128 x,u128 y){
u64 x0=(u64)x, x1=(u64)(x>>64);
u64 y0=(u64)y, y1=(u64)(y>>64);
u128 z11=(u128)x1*y1;
u128 z10=(u128)x1*y0;
u128 z01=(u128)x0*y1;
u128 z00=(u128)x0*y0;
u128 mid=(z00>>64)+(u64)z10+(u64)z01;
return z11+(z10>>64)+(z01>>64)+(mid>>64);
}
explicit Mont128(u128 n):mod(n){
u128 x=n;
for(int bits=1;bits<128;bits<<=1) x*=2-x*n;
ninv=-x;
r1=(-n)%n;
r2=r1;
for(int i=0;i<128;i++){
r2<<=1;
if(r2>=n) r2-=n;
}
}
inline u128 mul(u128 x,u128 y) const {
u128 lo=x*y;
u128 m=lo*ninv;
u128 t=mul_hi(x,y)+mul_hi(m,mod)+(lo!=0);
return t>=mod?t-mod:t;
}
inline u128 init(u128 x) const { return mul(x%mod,r2); }
inline u128 norm(u128 x) const { return mul(x,1); }
inline u128 add(u128 x,u128 y) const { u128 z=x+y; if(z>=mod) z-=mod; return z; }
inline u128 pow_mont(u128 x,u128 e) const {
u128 a=r1;
while(e){ if(e&1) a=mul(a,x); x=mul(x,x); e>>=1; }
return a;
}
inline u128 pow(u128 x,u128 e) const { return norm(pow_mont(init(x),e)); }
};
static u128 gcd128(u128 a,u128 b){ while(b){u128 r=a%b;a=b;b=r;}return a; }
static bool is_prime128(u128 n){
if(n<2) return false;
static const uint32_t small[]={2,3,5,7,11,13,17,19,23,29,31,37};
for(uint32_t p:small){ if(n==p) return true; if(n%p==0) return false; }
Mont128 md(n);
u128 d=n-1; int s=0; while(!(d&1)){d>>=1;++s;}
auto witness=[&](u128 a){
if(a%n==0) return false;
u128 x=md.pow_mont(md.init(a),d);
if(x==md.r1 || x==md.mod-md.r1) return false;
for(int r=1;r<s;r++){
x=md.mul(x,x);
if(x==md.mod-md.r1) return false;
}
return true;
};
// Deterministic over n < 3.317e24, hence over the problem range.
static const uint32_t bases[]={2,3,5,7,11,13,17,19,23,29,31,37,41};
for(uint32_t a:bases) if(witness(a)) return false;
return true;
}
static inline u64 mul_mod64(u64 a,u64 b,u64 m){ return (u128)a*b%m; }
static u64 pow_mod64(u64 a,u64 e,u64 m){u64 r=1;while(e){if(e&1)r=mul_mod64(r,a,m);a=mul_mod64(a,a,m);e>>=1;}return r;}
static u64 rho64(u64 n){
if(n%2==0)return 2;
auto rng=[]()->u64{ static u64 x=0x9e3779b97f4a7c15ULL ^ (u64)chrono::high_resolution_clock::now().time_since_epoch().count(); x^=x<<7;x^=x>>9;return x;};
while(true){
u64 y=rng()%(n-1)+1,c=rng()%(n-1)+1,m=128;
u64 g=1,r=1,q=1,x=0,ys=0;
auto f=[&](u64 v){return (mul_mod64(v,v,n)+c)%n;};
while(g==1){
x=y;for(u64 i=0;i<r;i++)y=f(y);
for(u64 k=0;k<r&&g==1;k+=m){
ys=y;u64 lim=min(m,r-k);
for(u64 i=0;i<lim;i++){y=f(y);u64 d=x>y?x-y:y-x;q=mul_mod64(q,d,n);}
g=std::gcd(q,n);
}
r<<=1;
}
if(g==n){do{ys=f(ys);u64 d=x>ys?x-ys:ys-x;g=std::gcd(d,n);}while(g==1);}
if(g!=n)return g;
}
}
static u64 tonelli(uint32_t n,uint32_t p){
if(p==2) return n&1;
if(n==0) return 0;
if(p%4==3) return pow_mod64(n,(p+1)/4,p);
uint32_t q=p-1,s=0;while((q&1)==0){q>>=1;++s;}
uint32_t z=2;while(pow_mod64(z,(p-1)/2,p)!=p-1)++z;
u64 c=pow_mod64(z,q,p),x=pow_mod64(n,(q+1)/2,p),t=pow_mod64(n,q,p);uint32_t m=s;
while(t!=1){
uint32_t i=1;u64 tt=mul_mod64(t,t,p);while(tt!=1){tt=mul_mod64(tt,tt,p);++i;}
u64 b=pow_mod64(c,1ULL<<(m-i-1),p);
x=mul_mod64(x,b,p);c=mul_mod64(b,b,p);t=mul_mod64(t,c,p);m=i;
}
return x;
}
struct FBPrime{uint32_t p,r1,r2;uint16_t lg;};
struct Atom{u64 x,q;};
struct Relation{vector<u64> parity; array<Atom,2> a; uint8_t cnt;};
struct Partial{vector<u64> parity; Atom a;};
static vector<uint32_t> primes_upto(int B){
vector<bool> isp(B+1,true);isp[0]=isp[1]=false;
vector<uint32_t> ps;
for(int i=2;i<=B;i++)if(isp[i]){ps.push_back(i);if((int64_t)i*i<=B)for(int j=i*i;j<=B;j+=i)isp[j]=false;}
return ps;
}
static u128 quadratic_sieve_factor(u128 n, int B=7000, int EXTRA=48, u64 MAX_BASE=3500000){
// Parameters tuned for 65--80 bit inputs.
const int BLOCK=1<<15;
vector<uint32_t> plist=primes_upto(B);
// A small Knuth--Schroeppel-style multiplier search.
static const uint32_t kcand[]={1,3,5,7,11,13,15,17,19,21,23,29,31,33,35,37,39,41,43,47};
uint32_t bestk=1; long double bestscore=-1e100L;
for(uint32_t k:kcand){
u128 g=gcd128(n,k);if(g>1&&g<n)return g;
long double sc=-0.5L*log((long double)k);
for(uint32_t p:plist){if(p>100)break;u64 a=(u64)((n%p)*(k%p)%p);if(a==0)sc+=log((long double)p)/(p-1);else if(p==2 || pow_mod64(a,(p-1)/2,p)==1)sc+=2*log((long double)p)/(p-1);}
if(sc>bestscore){bestscore=sc;bestk=k;}
}
u128 kn=n*bestk;
u64 m=isqrt_u128(kn);if((u128)m*m<kn)++m;
if((u128)m*m==kn){u128 g=gcd128(m,n);if(g>1&&g<n)return g; bestk=1;kn=n;m=isqrt_u128(n);if((u128)m*m<n)++m;}
vector<FBPrime> fb;fb.reserve(plist.size()/2+8);
for(uint32_t p:plist){
u64 np=(u64)(n%p);if(np==0)return p;
u64 a=(u64)((kn%p));
if(p==2){
uint32_t y=a&1, mm=m&1;uint32_t r=(y+2-mm)&1;
fb.push_back({p,r,r,64});
}else{
u64 ls=pow_mod64(a,(p-1)/2,p);
if(ls!=1 && a!=0)continue;
uint32_t y=(uint32_t)tonelli((uint32_t)a,p), mm=m%p;
uint32_t r1=(y+p-mm)%p;
uint32_t r2=((y? p-y:0)+p-mm)%p;
fb.push_back({p,r1,r2,(uint16_t)llround(log2((long double)p)*64)});
}
}
const int F=fb.size(), PW=(F+63)>>6;
vector<Relation> rels;rels.reserve(F+EXTRA+8);
unordered_map<u64,Partial> partial;partial.reserve((F+EXTRA)*8);
vector<uint16_t> score(BLOCK);
const u64 LPBOUND=(u64)B*B;
for(u64 base=0; rels.size()<(size_t)(F+EXTRA); base+=BLOCK){
fill(score.begin(),score.end(),0);
for(const auto& z:fb){
uint32_t p=z.p;
uint32_t bm=base%p;
uint32_t s1=(z.r1+p-bm)%p;
for(uint32_t i=s1;i<BLOCK;i+=p) score[i]+=z.lg;
if(z.r2!=z.r1){uint32_t s2=(z.r2+p-bm)%p;for(uint32_t i=s2;i<BLOCK;i+=p)score[i]+=z.lg;}
}
for(uint32_t i=0;i<BLOCK && rels.size()<(size_t)(F+EXTRA);++i){
u64 x=base+i;
u128 qq=(u128)(m+x)*(m+x)-kn;
if(qq==0 || qq>numeric_limits<u64>::max()) continue;
u64 q=(u64)qq;
int bits=64-__builtin_clzll(q);
if((int)score[i]+14*64 < bits*64) continue;
u64 rem=q;
vector<u64> par(PW);
for(int j=0;j<F;j++){
uint32_t p=fb[j].p;int e=0;
while(rem%p==0){rem/=p;e^=1;}
if(e)par[j>>6]^=1ULL<<(j&63);
if(rem==1)break;
}
if(rem==1){
Relation r; r.parity=move(par);r.a[0]={x,q};r.cnt=1;rels.push_back(move(r));
}else if(rem<=LPBOUND){
auto it=partial.find(rem);
if(it==partial.end()) partial.emplace(rem,Partial{move(par),{x,q}});
else{
Relation r;r.parity=move(par);for(int w=0;w<PW;w++)r.parity[w]^=it->second.parity[w];
r.a[0]=it->second.a;r.a[1]={x,q};r.cnt=2;rels.push_back(move(r));partial.erase(it);
}
}
}
// q stays in 64 bits while base is below roughly 4e6 for worst-case n.
if(base>MAX_BASE && rels.size()<(size_t)(F+EXTRA)) break;
}
if(rels.size()<= (size_t)F) return 0;
const int R=rels.size(), CW=(R+63)>>6;
struct BasisRow{vector<u64> p,c;bool used=false;};
vector<BasisRow> basis(F);
vector<vector<u64>> deps;deps.reserve(EXTRA);
for(int i=0;i<R;i++){
vector<u64> p=rels[i].parity,c(CW);c[i>>6]|=1ULL<<(i&63);
bool inserted=false;
for(int col=F-1;col>=0;--col) if((p[col>>6]>>(col&63))&1ULL){
if(!basis[col].used){basis[col].used=true;basis[col].p=move(p);basis[col].c=move(c);inserted=true;break;}
for(int w=0;w<PW;w++)p[w]^=basis[col].p[w];
for(int w=0;w<CW;w++)c[w]^=basis[col].c[w];
}
if(!inserted){bool zero=true;for(u64 v:p)if(v){zero=false;break;}if(zero)deps.push_back(move(c));}
}
Mont128 md(n);
for(auto &dep:deps){
vector<int> cnt(F);
unordered_map<u64,int> extra;extra.reserve(16);
u128 X=md.r1;
for(int ri=0;ri<R;ri++)if((dep[ri>>6]>>(ri&63))&1ULL){
for(int ai=0;ai<rels[ri].cnt;ai++){
auto [x,q]=rels[ri].a[ai];
X=md.mul(X,md.init((u128)(m+x)%n));
u64 rem=q;
for(int j=0;j<F;j++){
uint32_t p=fb[j].p;
while(rem%p==0){rem/=p;cnt[j]++;}
if(rem==1)break;
}
if(rem>1)extra[rem]++;
}
}
u128 Y=md.r1;
for(int j=0;j<F;j++) if(cnt[j]){
if(cnt[j]&1){Y=0;break;}
Y=md.mul(Y,md.pow_mont(md.init(fb[j].p),cnt[j]/2));
}
if(Y==0)continue;
for(auto [p,e]:extra){if(e&1){Y=0;break;}Y=md.mul(Y,md.pow_mont(md.init(p),e/2));}
if(Y==0)continue;
u128 xn=md.norm(X),yn=md.norm(Y);
u128 d=gcd128(xn>yn?xn-yn:yn-xn,n);
if(d>1&&d<n)return d;
u128 s=xn+yn;if(s>=n)s-=n;
d=gcd128(s,n);if(d>1&&d<n)return d;
}
return 0;
}
static u128 rho128(u128 n){
if(n%2==0)return 2;
auto rng=[]()->u64{static u64 x=0x243f6a8885a308d3ULL^(u64)chrono::high_resolution_clock::now().time_since_epoch().count();x^=x<<7;x^=x>>9;return x;};
Mont128 md(n);
auto rnd=[&](){return (((u128)rng()<<64)|rng())%n;};
while(true){
u128 y=md.init(rnd()), c=md.init(rnd()+1), x=0, ys=0;
u64 r=1, m=256;u128 g=1;
auto f=[&](u128 v){return md.add(md.mul(v,v),c);};
while(g==1){
x=y;for(u64 i=0;i<r;i++)y=f(y);
for(u64 k=0;k<r&&g==1;k+=m){
ys=y;u128 q=md.r1;u64 lim=min(m,r-k);
for(u64 i=0;i<lim;i++){y=f(y);u128 d=x>y?x-y:y-x;q=md.mul(q,d);}
g=gcd128(q,n);
}
r<<=1;
}
if(g==n){do{ys=f(ys);u128 d=x>ys?x-ys:ys-x;g=gcd128(d,n);}while(g==1);}
if(g!=n)return g;
}
}
static void factor_rec(u128 n, vector<u128>& out){
if(n==1)return;
static const uint32_t small[]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97};
for(uint32_t p:small)if(n%p==0){out.push_back(p);factor_rec(n/p,out);return;}
if(is_prime128(n)){out.push_back(n);return;}
u64 s=isqrt_u128(n);if((u128)s*s==n){factor_rec(s,out);factor_rec(s,out);return;}
u128 d=0;
if(n<=numeric_limits<u64>::max())d=rho64((u64)n);
else {
d=quadratic_sieve_factor(n,7000,48,3500000);
if(d==0||d==n)d=quadratic_sieve_factor(n,10000,64,5000000);
if(d==0||d==n)d=rho128(n);
}
factor_rec(d,out);factor_rec(n/d,out);
}
static inline u128 pow_cap(u128 a,int e,u128 cap){u128 r=1;for(int i=0;i<e;i++){if(a && r>cap/a)return cap+1;r*=a;}return r;}
static u128 power_sum_small(int e,u128 l,u128 w){
if(e==1)return l*w+w*(w-1)/2;
if(e==2)return l*w*(l+w-1)+w*(w-1)/2*(2*w-1)/3;
if(e==3)return w*(2*l+w-1)/2*l*(l+w-1)+w*w*(w-1)/2*(2*l+w-1)/2;
if(e==4)return l*w*(l+w-1)*(l*(l+w-1)+w*(w-1))+w*(w-1)/2*(2*w-1)/3*(3*w*(w-1)-1)/5;
if(e==5)return w*(2*l+w-1)/2*(l*l*l*l+(2*l+w)*(w-1)/2*(2*l*(3*l+2*w-1)+2*w*(w-1)-1)/3);
abort();
}
struct Sol {
int e;
u64 l, r;
};
static vector<u64> divisors_limited_fast(const map<u128,int>& fac, u64 limit){
vector<u64> ds;
ds.reserve(1 << 12);
ds.push_back(1);
for(auto [pp,e]:fac){
if(pp > limit) continue;
u64 p=(u64)pp;
size_t old=ds.size();
// Existing divisors are exactly those with exponent 0 for p.
for(size_t i=0;i<old;i++){
u64 x=ds[i];
for(int k=1;k<=e;k++){
if(x > limit/p) break;
x*=p;
ds.push_back(x);
}
}
}
return ds;
}
// All E=1 widths are < 2^41. Three 14-bit stable passes suffice.
static void radix_sort_42(vector<u64>& a){
if(a.size()<2) return;
constexpr unsigned B=14;
constexpr unsigned R=1u<<B;
constexpr u64 MASK=R-1;
vector<u64> b(a.size());
array<uint32_t,R> cnt{};
for(unsigned sh=0;sh<42;sh+=B){
cnt.fill(0);
for(u64 x:a) ++cnt[(x>>sh)&MASK];
uint32_t s=0;
for(unsigned i=0;i<R;i++){
uint32_t c=cnt[i]; cnt[i]=s; s+=c;
}
for(u64 x:a) b[cnt[(x>>sh)&MASK]++]=x;
a.swap(b);
}
}
static inline bool square_u128(u128 x, u64& root){
root=isqrt_u128(x);
return (u128)root*root==x;
}
static inline u128 cube_poly(u64 x, u128 a){
return (u128)x*x*x + a*x;
}
static u64 cubic_root_for(u128 A, u128 a){
u64 x=(u64)cbrtl((long double)A);
if(x==0) x=1;
// Integer Newton from an upper estimate.
for(int it=0;it<8;it++){
u128 x2=(u128)x*x;
u64 nx=(u64)((2*x2*x + A)/(3*x2+a));
if(nx>=x) break;
x=nx;
}
while(cube_poly(x,a)>A) --x;
while(cube_poly(x+1,a)<=A) ++x;
return x;
}
static inline void append_u64(string& out, u64 x){
char buf[24];
auto [p,ec]=to_chars(buf,buf+24,x);
out.append(buf,p);
}
static inline void append_fixed18(string& out, u64 x){
char buf[18];
for(int i=17;i>=0;--i){buf[i]=char('0'+x%10);x/=10;}
out.append(buf,18);
}
static inline void append_u128(string& out, u128 x){
constexpr u64 BASE=1000000000000000000ULL;
if(x<BASE){append_u64(out,(u64)x);return;}
u64 hi=(u64)(x/BASE);
u64 lo=(u64)(x-(u128)hi*BASE);
append_u64(out,hi);
append_fixed18(out,lo);
}
static inline void append_line(string& out,int e,u128 l,u128 r){
append_u64(out,(u64)e); out.push_back(' ');
append_u128(out,l); out.push_back(' ');
append_u128(out,r); out.push_back('\n');
}
int main(){
ios::sync_with_stdio(false);
cin.tie(nullptr);
string ss;
if(!(cin>>ss)) return 0;
u128 N=parse_u128(ss);
vector<u128> pf;
factor_rec(N,pf);
sort(pf.begin(),pf.end());
map<u128,int> fn;
for(u128 p:pf) ++fn[p];
// E=1. Store only the width; L is strictly decreasing as width increases.
map<u128,int> f2=fn;
++f2[2];
u128 twoN=2*N;
u64 limit1=isqrt_u128(twoN-1);
vector<u64> widths1=divisors_limited_fast(f2,limit1);
int v2total=1;
for(u128 t=N;(t&1)==0;t>>=1) ++v2total;
size_t wr=0;
for(u64 w:widths1){
// w and 2N/w must have opposite parity.
bool ok=(w&1) || (v2total<64 && __builtin_ctzll(w)==v2total);
if(ok) widths1[wr++]=w;
}
widths1.resize(wr);
radix_sort_42(widths1);
vector<Sol> rest;
rest.reserve(64);
static constexpr int D[5]={0,0,6,2,30};
for(int e=2;e<=4;e++){
map<u128,int> fd=fn;
int d=D[e];
for(int p=2;p*p<=d;p++) while(d%p==0){++fd[p];d/=p;}
if(d>1) ++fd[d];
u64 hi=1;
while(power_sum_small(e,1,hi)<=N) hi<<=1;
u64 lo=0;
while(lo+1<hi){
u64 mid=(lo+hi)>>1;
if(power_sum_small(e,1,mid)<=N) lo=mid;
else hi=mid;
}
vector<u64> widths=divisors_limited_fast(fd,lo);
if(e==2){
for(u64 w:widths){
u128 q=12*N/w;
u128 w2=(u128)w*w;
if(q+1<w2) continue;
u128 t=q-w2+1;
if(t%3) continue;
u64 x;
if(!square_u128(t/3,x)) continue;
if(x<w+1 || ((x^w)&1)==0) continue;
u64 l=(x-w+1)>>1;
rest.push_back({2,l,l+w-1});
}
}else if(e==3){
for(u64 w:widths){
u128 A=8*N/w;
u128 a=(u128)w*w-1;
u64 x=cubic_root_for(A,a);
if(cube_poly(x,a)!=A) continue;
if(x<w+1 || ((x^w)&1)==0) continue;
u64 l=(x-w+1)>>1;
rest.push_back({3,l,l+w-1});
}
}else{
for(u64 w:widths){
u128 w2=(u128)w*w;
u128 T=3*w2*w2-5*w2+2+60*N/w;
u128 disc=240*T;
u64 s;
if(!square_u128(disc,s)) continue;
u128 a=w2-1;
if((u128)s<30*a) continue;
u128 dy=(u128)s-30*a;
if(dy%30) continue;
u64 x;
if(!square_u128(dy/30,x)) continue;
if(x<w+1 || ((x^w)&1)==0) continue;
u64 l=(x-w+1)>>1;
rest.push_back({4,l,l+w-1});
}
}
}
int maxe=0;
for(u128 t=N;t>=2;t>>=1) ++maxe;
for(int e=5;e<=maxe;e++){
vector<u128> pw(1,0);
for(u64 i=1;;i++){
u128 v=pow_cap(i,e,N);
if(v>N) break;
pw.push_back(v);
}
size_t r=1;
u128 sum=0;
for(size_t l=1;l<pw.size();l++){
while(r<pw.size() && sum+pw[r]<=N){sum+=pw[r];++r;}
if(sum==N) rest.push_back({e,(u64)l,(u64)(r-1)});
if(r>l) sum-=pw[l];
}
}
sort(rest.begin(),rest.end(),[](const Sol& a,const Sol& b){
if(a.e!=b.e) return a.e<b.e;
if(a.l!=b.l) return a.l<b.l;
return a.r<b.r;
});
const size_t total=widths1.size()+rest.size();
string out;
out.reserve(total*58+32);
append_u64(out,total);
out.push_back('\n');
// widths descending => L ascending.
for(auto it=widths1.rbegin();it!=widths1.rend();++it){
u64 w=*it;
u128 z=twoN/w;
u128 l=(z-w+1)>>1;
append_line(out,1,l,l+w-1);
}
for(const Sol& s:rest) append_line(out,s.e,s.l,s.r);
fwrite(out.data(),1,out.size(),stdout);
return 0;
}