Groundhog Chasing Death

题目：
https://ac.nowcoder.com/acm/contest/5674/E
calculate $\prod\limits_{i=a}^{b}\prod\limits_{j=c}^{d}gcd(x^i,y^j)$modulo $998244353$

思路：
$\begin{aligned} &\prod\limits_{i=a}^{b}\prod\limits_{j=c}^{d}gcd(x^i,y^j)\\ =&\prod\limits_{i=a}^{b}\prod\limits_{j=c}^{d}\prod\limits_{k=1}^{n}p_k^{min(a_k\cdot i,b_k\cdot j)}\quad p_k为x和y共同的素因子,a_k为x中p_k的幂,同理b_k\\ =&\prod\limits_{k=1}^{n}\prod\limits_{i=a}^{b}\prod\limits_{j=c}^{d}p_k^{min(a_k\cdot i,b_k\cdot j)} \end{aligned}$
然后分类讨论：
枚举 $k$， $i$，当 $j\le\lfloor\frac{a_k\cdot i}{b_k}\rfloor$时， $min(a_k\cdot i,b_k\cdot j)=b_k\cdot j$，当 $j\ge\lfloor\frac{a_k\cdot i}{b_k}\rfloor+1$时， $min(a_k\cdot i,b_k\cdot j)=a_k\cdot i$。
当然还要对 $\lfloor\frac{a_k\cdot i}{b_k}\rfloor$进行分类讨论，取其中一种情况说明。
当 $c\le\lfloor\frac{a_k\cdot i}{b_k}\rfloor\le d$时
$\begin{aligned} &\prod\limits_{k=1}^{n}\prod\limits_{i=a}^{b}\prod\limits_{j=c}^{d}p_k^{min(a_k\cdot i,b_k\cdot j)}\\ =&\prod\limits_{k=1}^{n}\prod\limits_{i=a}^{b}p_k^{b_k\cdot \frac{(c+\lfloor\frac{a_k\cdot i}{b_k}\rfloor)\cdot (\lfloor\frac{a_k\cdot i}{b_k}\rfloor-c+1)}{2}+a_k\cdot i\cdot (d-\lfloor\frac{a_k\cdot i}{b_k}\rfloor+1)} \end{aligned}$
另外，可以预处理出 $\forall i\in[1,100000],p_k^i,p_k^{i*100000}$的值。则欧拉降幂后 $p_k^x=p_k^{100000\cdot a+b}$
$a\le100000,b\le100000$。这样可以省掉一个快速幂。

代码：

#include<bits/stdc++.h>
#pragma GCC optimize(2)
#define ll long long
#define mod 998244353
using namespace std;
struct node
{ll p,num;node(int x,int y):p(x),num(y) {}
};
struct node1
{ll p,num,num1;
}w2[200];
struct node2
{ll sum[100001],sum1[100001];
}w3[200];
ll a,b,c,d,x,y,ans=1;
set<ll>w,w1;
set<ll>::iterator it;
map<ll,ll>ma,ma1;
int si,k=0;
int cmp(node a,node b)
{return a.num<b.num;
}
int main()
{cin>>a>>b>>c>>d>>x>>y;ll tmp=x;for(ll i=2; i*i<=tmp; i++){if(x%i==0){ll num=0;while(x%i==0){num++;x/=i;}w.insert(i),ma[i]=num;}}if(x!=1)w.insert(x),ma[x]=1;tmp=y;for(ll i=2; i*i<=tmp; i++){if(y%i==0){ll num=0;while(y%i==0){num++;y/=i;}w1.insert(i),ma1[i]=num;}}if(y!=1)w1.insert(y),ma1[y]=1;for(it=w.begin(); it!=w.end(); it++)if(w1.find((*it))!=w1.end())w2[k].p=(*it),w2[k].num=ma[(*it)],w2[k++].num1=ma1[(*it)];for(int i=0;i<k;i++){w3[i].sum[0]=w3[i].sum1[0]=1;for(int j=1;j<=100000;j++)w3[i].sum[j]=w3[i].sum[j-1]*w2[i].p%mod;w3[i].sum1[1]=w3[i].sum[100000];for(int j=2;j<=100000;j++)w3[i].sum1[j]=w3[i].sum1[j-1]*w3[i].sum1[1]%mod;}for(ll i=a;i<=b;i++){for(int j=0;j<k;j++){ll locate=w2[j].num*i/w2[j].num1;if(locate<c)ans=ans*w3[j].sum[(w2[j].num*i*(d-c+1)%(mod-1))%100000]%mod*w3[j].sum1[(w2[j].num*i*(d-c+1)%(mod-1))/100000]%mod;else if(locate>=c&&locate<=d)ans=ans*w3[j].sum[(w2[j].num1*((c+locate)*(locate-c+1)/2)%(mod-1))%100000]%mod*w3[j].sum1[(w2[j].num1*((c+locate)*(locate-c+1)/2)%(mod-1))/100000]%mod*w3[j].sum[(w2[j].num*i*(d-locate)%(mod-1))%100000]%mod*w3[j].sum1[(w2[j].num*i*(d-locate)%(mod-1))/100000]%mod;elseans=ans*w3[j].sum[(w2[j].num1*((c+d)*(d-c+1)/2)%(mod-1))%100000]%mod*w3[j].sum1[(w2[j].num1*((c+d)*(d-c+1)/2)%(mod-1))/100000]%mod;}}cout<<ans;return 0;
}