Day5
$\texttt{T1}$
简要题意
有两个正整数 \(a<b\le 10^9\)，给出 \(\dfrac{a}{b}\) 的小数点后 \(19\) 位，要求还原 \(a,b\)，保证有解。
solution
一个科技：\(\texttt{Stern-Brocot tree}(SBT)\)，可以参考这个博客学习。
先给出 $O(n)$ 找的代码
#include<bits/stdc++.h>
#define LL long long
#define LD long double
#define int LL
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
using namespace std;
const LD eps=1e-19;
int A,B;
#define Iabs(x) ((x)>0?(x):-(x))
#define C (x-(LD)(m)/n)
void find(LD x,int a=0,int b=1,int c=1,int d=0)
{
int m=a+c,n=b+d;
if(Iabs(C)<eps) return A=m,B=n,void();
if(x<(LD)(m)/n) return find(x,a,b,m,n);
return find(x,m,n,c,d);
}
signed main()
{
ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);LD x;cin>>x;find(x);cout<<A<<" "<<B;
return 0;
}
注意到博客中的一句话：称从 \(SBT\) 上任意一个节点 \(\dfrac{a}{b}\) 沿着相同方向跳到需要换方向的第一个点为一次跳跃，那么该点到根节点至多跳跃 \(O(\log\max(a,b))\) 次。
于是跳的过程中二分优化一下，复杂度 \(O(\log^2 V)\)。
$\texttt{code}$
#include<bits/stdc++.h>
#define LL long long
#define LD long double
#define int LL
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
using namespace std;
const LD eps=1e-19;
int A,B;
#define Iabs(x) ((x)>0?(x):-(x))
#define C (x-(LD)(m)/n)
void find(LD x,int a=0,int b=1,int c=1,int d=0)
{
int m=a+c,n=b+d,l=1,r=1e9,mid,ans=-1;
if(Iabs(C)<eps) return A=m,B=n,void();
if(x<(LD)(m)/n)
{
while(l<=r)
{
int mid=(l+r)>>1;m=c+mid*a,n=d+mid*b;
if(C>=-eps) ans=mid,r=mid-1;else l=mid+1;
}ans--;
return find(x,a,b,c+ans*a,d+ans*b);
}
while(l<=r)
{
int mid=(l+r)>>1;m=a+mid*c,n=b+mid*d;
if(C<=eps) ans=mid,r=mid-1;else l=mid+1;
}ans--;
return find(x,a+ans*c,b+ans*d,c,d);
}
signed main()
{
ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);LD x;cin>>x;find(x);cout<<A<<" "<<B;
return 0;
}
$\texttt{T3}$
简要题意
有 \(1\le T\le 10^6\) 次询问，每次询问正整数 \(x,p\)，\(p\) 为素数，令 \(n=xp^2\le 10^{18}\)，问是否存在三个正整数 \(a,b,c\)，满足 \(ab+bc+ca=n\)。有的话给出构造，否则输出 \(-1\) 。
solution
打表注意到只有 \(n=4,18\) 是无解的。
打表
namespace DB
{
const int N=1e5;
struct node{int x,y,z;}db[N+5];
inline void init()
{
for(int i=1;i<=N;i++)
{
bool ok=1;
for(int a=1;a*a<=i;a++)
{
for(int b=1;b*b<=i;b++)
{
int t=i-a*b;if(t<=0) continue;
if(t%(a+b)==0){db[i]={a,b,t/(a+b)};ok=0;break;}
}
if(!ok) break;
}
}
}
}using DB::db;
这里尽管理论复杂度是 \(O(N^2)\) 的，但是实际跑 \(10^5\) 大约 \(0.2s\) 不到。我大概跑了 \(N=10^6\) 验证只有那两个无解。
下面做法不依赖与 \(n=xp^2\)
考虑钦定其中一个数 \(c\) 为 \(r\)，两边同时加 \(r^2\) 得：\((r+a)(r+b)=n+r^2\)
考虑找一个素数 \(q\mid n+r^2\)，发现 \(\lceil\) 是否存在 \(r\) 满足条件 \(\rfloor\) 这是一个二次剩余问题。我们不妨假设我们已经会快速求二次剩余了（后面其实不需要）。
此时，我们还需满足 \(a,b>0\)，即 \(q>r,\dfrac{n+r^2}{q}>r\)。首先 \(q>r\) 是容易满足的。
而 \(\dfrac{n+r^2}{q}>r\Leftarrow \dfrac{n+r^2}{q}>q\Leftarrow n+r^2>q^2\Leftarrow q^2<n\)。
由于我们预处理了 \(n\le 10^5\) 的情况，于是我们只能猜测在 \(q\le 300\) 以内能找到满足 \(-n\) 是模 \(q\) 的二次剩余，设素数为 \(q\)，二次剩余为 \(r\)。
此时 \(a=q-r,b=\dfrac{n+r^2}{q}-r\)，于是我们找到了解。
再来解决上面的两个遗留问题。
- 快速求二次剩余，这里我考场上做得很不光明正大，是贺板子的，后面不想改了。优秀做法：发现只要对 \(300\) 以内素数求二次剩余，预处理即可。
- 猜测在 \(q\le 300\) 以内能找到满足 \(-n\) 是模 \(q\) 的二次剩余。随机一个数是模 \(p\) 的二次剩余的概率大概为 \(\dfrac{1}{2}\)，于是大胆猜想在 \(300\) 内都没找到的概率很小很小。如果有数竞强神证证就好了。
于是我就通过过了这道题。复杂度 \(O(T)\)，常数 \(60\) 左右，实际在随机情况下跑得飞快，出题人很懒，也是随机造数据。
考场代码
#include<bits/stdc++.h>
#define P pair<int,int>
#define fi first
#define se second
#define LL long long
#define int LL
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
using namespace std;
namespace IO
{
const int _Pu=21e7+5,_d=32;
char buf[_Pu],obuf[_Pu],*p1=buf+_Pu,*p2=buf+_Pu,*p3=obuf,*p4=obuf+_Pu-_d;
inline void fin()
{
memmove(buf,p1,p2-p1);
int rlen=fread(buf+(p2-p1),1,p1-buf,stdin);
if(p1-rlen>buf) buf[p2-p1+rlen]=EOF;p1=buf;
}
inline void fout(){fwrite(obuf,p3-obuf,1,stdout),p3=obuf;}
inline int rd()
{
if(p1+_d>p2) fin();int isne=0,x=0;
for(;!isdigit(*p1);++p1) isne=(*p1=='-');x=(*p1++-'0');
for(;isdigit(*p1);++p1) x=x*10+(*p1-'0');
if(isne) x=-x;return x;
}
inline void wr(int x,char end='\n')
{
if(!x) return *p3++='0',*p3++=end,void();
if(x<0) *p3++='-',x=-x;
char sta[20],*top=sta;
do{*top++=(x%10)+'0';x/=10;}while(x);
do{*p3++=*--top;}while(top!=sta);(*p3++)=end;
}
}//快读（贺的板子），注意 _Pu 一定要开大，否则会少读/少输
LL T,x,p,n;
const LL un[]={1,2,4,6,10,18,22,30,42,58,70,78,102,130,190,210,330,462};//事实上只有 4,18 是 xp^2 型
namespace DB
{
const int N=1e5;
struct node{int x,y,z;}db[N+5];
inline void init()
{
for(int i=1;i<=N;i++)
{
bool ok=1;
for(int a=1;a*a<=i;a++)
{
for(int b=1;b*b<=i;b++)
{
int t=i-a*b;if(t<=0) continue;
if(t%(a+b)==0){db[i]={a,b,t/(a+b)};ok=0;break;}
}
if(!ok) break;
}
}
}
}using DB::db;//打表
namespace er
{
const int pr[]={3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293};
int n,mod,w;
struct comp
{
int x,y;
inline comp operator *(comp X)
{
return {(1ll*x*X.x%mod+1ll*y*X.y%mod*w%mod)%mod,
(1ll*x*X.y%mod+1ll*y*X.x%mod)%mod};
}
};
inline int ksm(int x,int p)
{
int s=1;
for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);
return s;
}
inline comp ksm1(comp x,int p)
{
comp s={1,0};
for(;p;(p&1)&&(s=s*x,1),x=x*x,p>>=1);
return s;
}
inline int Find()
{
int x;
while(1)
{
x=rand()%mod;w=((LL)x*x%mod-n+mod)%mod;
if(ksm(w,(mod-1)>>1)==mod-1) break;
}
return ksm1({x,1},(mod+1)>>1).x;
}
inline P fer(int a)
{
for(int i=0;i<60;i++)
{
mod=pr[i];n=mod-a%mod;if(!n||n==mod) continue;
if(ksm(n,(mod-1)>>1)==mod-1) continue;
int ans=Find();return {mod,ans};
}
return {-1,-1};
}
}using er::fer;//不光明正大的地方，你们可以把这里改成预处理300以内二次剩余
signed main()
{
// fr(construct)
T=IO::rd();DB::init();
while(T--)
{
x=IO::rd(),p=IO::rd();n=x*p*p;bool ok=1;
for(int i=0;i<18;i++) if(n==un[i]){ok=0;break;}//自己这里判断改一下
if(!ok){IO::wr(-1);continue;}
if(n&1){IO::wr(1,' ');IO::wr(1,' ');IO::wr((n-1)/2);continue;}//防止挂分特判
if(n%4==0){IO::wr(2,' ');IO::wr(2,' ');IO::wr((n-4)/4);continue;}//防止挂分特判
if(n<=1e5){IO::wr(db[n].x,' ');IO::wr(db[n].y,' ');IO::wr(db[n].z);continue;}
P t=fer(n);int q=t.fi,r=t.se;r=min(r,q-r);int Q=(n+r*r)/q;
int a=q-r,b=Q-r;
IO::wr(a,' ');IO::wr(b,' ');IO::wr(r);
}
return IO::fout(),0;
}
后记
考场上写完这题拍了好久就摆烂了，发现最后 \(\texttt{rk14}\)，是这几次最好，很开心。也有很多大佬切了这题，卷王，\(Linshey\) 两位 \(\texttt{AK}\)，珂怕。
正解是个很傻逼的结论，参考 oeis，实际上就是把 \(\texttt{un}\) 数组放进去搜。发现我这乱搞还是很优秀的，可以看做正解。感受那股劲！
Day6
$\texttt{T3}$
solution
首先，随机化、退火、贪心找规律等乱搞我们在考场上拿了 \([50,60]\) 的分，还是很高的。
考虑观察到一件事情：我们只需对于所有素数 \(p\)，确定 \(a_p\) 即可确定所有 \(a_i\)。
由于 \(a_i\) 完全积性，我们想到 \(\texttt{Jacobi}\) 符号，即考虑 \(\left(\dfrac{i}{3}\right)\)，我们令 \(a_3=-1,a_p=\left(\dfrac{p}{3}\right)(p\neq 3)\)，这样得到的 \(a\) 数列就有一个为 \(3\) 的循环节（不严谨，但你可以这样理解）。
于是这样的构造就比较优秀，实际 \(S=7,85pts\)，其中 \(a_3\gets 1\) 经测试不会使答案更优秀。
考虑在此基础上进行调整，对于前缀和的绝对值过大的情况，通过修改之前质数处的取值以控制 \(S\) 范围，可以做到 \(S = 5\)。
这句话是 \(Renshey\) 题解里的说法。我具体来说：
inline bool chk(){for(int i=1;i<=n-70;i++) if(abs(s[i+70]-s[i])>6) return 0;return 1;}
while(1)
{
sol();int s=0,w=-1;
for(int i=1;i<=n;i++)
{
s+=a[i];
if(abs(s)>4/*or5*/){w=i;break;}
}
if(w==-1) break;
for(int i=w;i;i--) if(!d[i]&&a[i]==sig(::s[w]))
{
a[i]=-a[i];sol();
if(chk()) break;a[i]=-a[i];
}
}
至于为啥我也不知道，反正这样调是对的。而且能做到让 \(S=4\)，如果要 \(S=5\) 把代码中 \(4\) 改成 \(5\) 即可。其中 \(4\) 跑了 \(12s\)，\(5\) 跑了不到 \(1s\)。
构造的代码
#include<bits/stdc++.h>
#define LL long long
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
using namespace std;
const int n=1e6;
int a[n+5],s[n+5],d[n+5],pr[n+5];
#define sig(x) ((x)>0?1:-1)
inline void sol(){for(int i=2;i<=n;i++){if(d[i]) a[i]=a[i/d[i]]*a[d[i]];s[i]=s[i-1]+a[i];}}
inline bool chk(){for(int i=1;i<=n-70;i++) if(abs(s[i+70]-s[i])>6) return 0;return 1;}
inline int f (int S)
{
if (S >= 100000) return ceil(1.0 * (1100000 - S) / 100000);
if (S >= 10000) return ceil(1.0 * (190000 - S) / 9000);
if (S >= 1000) return ceil(1.0 * (28000 - S) / 900);
if (S >= 500) return ceil(1.0 * (2500 - S) / 50);
if (S >= 100) return ceil(1.0 * (4700 - 3 * S) / 80);
if (S >= 17) return ceil(100 * pow(S, -0.5) + 45);
if (S >= 7) return ceil(100 * pow(log2(S) / 2, -0.5));
if (S == 6) return 92;
return 100;
}
inline void checkans()
{
int s=0,mx=0;
for(int i=1;i<=n;i++)
{
if(abs(a[i])!=1) return cerr<<"WA!a["<<i<<"]"<<" is not equal 1 or -1!",void();
s+=a[i],mx=max(mx,abs(s));
}
for(int i=1;i<=n;i++) for(int j=i;j<=n;j+=i) if(a[j]!=a[i]*a[j/i]) return cerr<<"WA!a["<<i<<"]*a["<<j/i<<"] is not equal a["<<j<<"]!",void();
int fen=f(mx);
if(fen==100) cerr<<"AC!";
else cerr<<"points = "<<fen<<"!";
}//假设你把答案数组放到a数组，编号从1开始，调用checkans()片段能帮助你 判断答案是否符合条件 and 查看你获得的分数
int main()
{
freopen("occultism.out","w",stdout);
for(int i=2;i<=n;i++) for(int j=2*i;j<=n;j+=i) if(!d[j]) d[j]=i;a[1]=s[1]=1;
for(int i=2;i<=n;i++)
{
if(!d[i]) a[i]=i%3==1?1:-1,pr[++pr[0]]=i;
else a[i]=a[i/d[i]]*a[d[i]];s[i]=s[i-1]+a[i];
}
while(1)
{
sol();int s=0,w=-1;
for(int i=1;i<=n;i++)
{
s+=a[i];
if(abs(s)>4/*or5*/){w=i;break;}
}
if(w==-1) break;
for(int i=w;i;i--) if(!d[i]&&a[i]==sig(::s[w]))
{
a[i]=-a[i];sol();
if(chk()) break;a[i]=-a[i];
}
}
checkans();
for(int i=1;i<=n;i++) (a[i]<0)&&(putchar('-')),putchar('1'),putchar(' ');
return 0;
}
后记
考场上乱搞我是 \(59pts\)，当时想过素数模 \(3\) 了，但是看打的表模 \(3\) 没有规律就没尝试了。还是应该勇于尝试的。