fork(2) download
copy 
  1. #include<iostream>
  2. #include<cstdlib>
  3. #include<cstring>
  4. #include<algorithm>
  5. #include<cmath>
  6. #include<stdio.h>
  7. #include<vector>
  8. #include<set>
  9. using namespace std;
  10. const int mod=1e9+7; 
  11. typedef long long LL; 
  12. int power(int a,int n)
  13. {
  14. 	if(n==0)
  15. 	   return 1;
  16. 	int ret=power(a,n/2);
  17. 	ret=1LL*ret*ret%mod;
  18. 	if(n%2)
  19. 	   ret=1LL*ret*a%mod;
  20. 	return ret;
  21. }
  22.  
  23. #define inf 0x3f3f3f3f
  24. #define mem(a,b) memset(a,b,sizeof(a))
  25. #define REP(i,a,b) for(int i=a; i<=b; ++i)
  26. #define PER(i,a,b) for(int i=a; i>=b; --i)
  27. #define MP make_pair
  28. #define PB push_back
  29. #define fi first
  30. #define se second
  31. typedef pair<int,int> pii;
  32. int ex[1101],fac[1101],f[2101],inv[1101];
  33. const double PI = acos(-1.0);
  34. const int Bin = 11;
  35. const int P=1e9+7;
  36. const int maxn = (1<<Bin);
  37. const int md=1e9+7;
  38. const int mmd =(1<<15);
  39. vector<int> vec[12];
  40. struct Comp {
  41.     double re, im;
  42.     Comp(): re(0), im(0) {}
  43.     Comp(const double _re, const double _im) { re = _re; im = _im; }
  44.     Comp operator+(const Comp &a) { return Comp(re+a.re, im+a.im); }
  45.     Comp operator-(const Comp &a) { return Comp(re-a.re, im-a.im); }
  46.     Comp operator*(const Comp &a) { return Comp(re*a.re-im*a.im, re*a.im+im*a.re); }
  47.     Comp operator*(const double &a) { return Comp(re*a, im*a); }
  48.     Comp operator/(const double &a) { return Comp(re/a, im/a); }
  49.     void operator*=(const Comp &a) { (*this) = (*this) * a; }
  50.     void operator+=(const Comp &a) { re += a.re; im += a.im; }
  51.     void operator*=(const double &a) { re*=a; im*=a; }
  52.     void operator/=(const double &a) { re/=a; im/=a; }
  53.     Comp conj() { return Comp(re, -im); }
  54. };
  55.  
  56. vector<Comp> ww[Bin + 2];
  57. void Init() {
  58.     for(int i = 1, l = 2; l <= maxn; ++i, l <<= 1) {
  59.        for(int k = 0; k < l/2; ++k)
  60. 	    	ww[i].push_back(Comp(cos(2*PI*k/l),sin(2*PI*k/l)));
  61.     }
  62. }
  63.  
  64. void FFT(Comp *a, int len, int type) {
  65.     int i, j, k, l, tt;
  66.     for(i = 1, j = len>>1; i < len-1; ++i) {
  67.         if(i < j) swap(a[i], a[j]);
  68.         k = len>>1;
  69.         while(j >= k) j -= k, k >>= 1;
  70.         j += k;
  71.     }
  72.     Comp var, u, v;
  73.     for(l = 2, tt = 1; l <= len; l <<= 1, ++tt) {
  74.         for(k = 0; k < len; k += l) {
  75.             for(i = k; i < k+l/2; ++i) {
  76.                 var = ww[tt][i-k];
  77.                 if(type == -1) var.im = -var.im;
  78.                 u = a[i], v = var * a[i+l/2];
  79.                 a[i] = u+v, a[i+l/2] = u-v;
  80.             }
  81.         }
  82.     }
  83.     if(type == -1) for(i = 0; i < len; ++i) a[i] /= len;
  84. }
  85. Comp A[maxn + 5], B[maxn + 5], C[maxn + 5], DD[maxn + 5];
  86. int aa[maxn*2],bb[maxn*2],cc[maxn*2],dd[maxn*2];
  87. void Conv(int *a, int *b, int *c, int len) {
  88.     for(int i = 0; i < len; ++i) {
  89.         A[i] = Comp(a[i]>>15, b[i]>>15);
  90.         B[i] = Comp(a[i]&(mmd-1), b[i]&(mmd-1));
  91.     }
  92.     FFT(A, len, 1); FFT(B, len, 1);
  93.     for(int i = 0; i < len; ++i) {
  94.         int j = (len-i) & (len-1);
  95.         C[i] = (A[i]*A[i] - (A[j]*A[j]).conj()) * Comp(0, -0.25)+(A[i]*B[i] - (A[j]*B[j]).conj()) * Comp(0.5, 0);
  96.         DD[i] = (B[i]*B[i] - (B[j]*B[j]).conj()) * Comp(0, -0.25);
  97.     }
  98.     FFT(C, len, -1); FFT(DD, len, -1);
  99.     for(int i = 0; i < len; ++i) {
  100.            c[i] = (((LL(C[i].re+0.5) %mod)<<30)%mod  + ((LL(C[i].im+0.5)%mod ) <<15)%mod+ LL(DD[i].re+0.5))%mod ;
  101. 	}
  102. }
  103. inline int read(){
  104.     int f=1,x=0;char ch=getchar();
  105.     while (ch<'0'||ch>'9'){if (ch=='-')f=-1;ch=getchar();}
  106.     while (ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
  107.     return f*x;
  108. }
  109.  
  110. inline int qpow(int x,int p){int ans=1;
  111.     for(;p;p>>=1,x=1ll*x*x%P)if(p&1)ans=1ll*ans*x%P;
  112.     return ans;
  113. } 
  114. int temp[maxn],a[maxn],b[maxn],c[maxn],d[maxn],tmp[maxn];
  115. inline void GetInv(int *a,int *b,int n){
  116.     if (n==1) return void(b[0]=qpow(a[0],P-2));
  117.     GetInv(a,b,n>>1);
  118.     for (int i=0;i<n;++i)tmp[i]=a[i],tmp[n+i]=0;
  119.     int k;for (k=1;k<=n;k<<=1);
  120.     for (int i=n;i<k;++i)b[i]=0;
  121.     Conv(tmp,b,tmp,k);
  122.     for(int i=n;i<k;i++)
  123.       tmp[i]=0;
  124.     Conv(tmp,b,tmp,k);
  125.     for (int i=0;i<n;++i)b[i]=(2*b[i]%mod-tmp[i]+mod)%mod,b[n+i]=0;
  126. }
  127. inline void Getln(int *a,int *b,int *c,int n){
  128.     for (int i=0;i<n;++i)b[i]=1ll*(i+1)*a[i+1]%P;
  129.     GetInv(a,c,n);
  130.     for (int i=0;i<n;++i)tmp[i]=c[i],tmp[n+i]=0;
  131.     int k;for (k=1;k<=n;k<<=1);
  132.     Conv(tmp,b,tmp,k);
  133.     for (int i=1;i<=n;++i)b[i]=1ll*tmp[i-1]*qpow(i,P-2)%P,b[i+n]=0;b[0]=0;
  134. }
  135.  
  136. inline void Getexp(int *a,int *b,int *c,int *d,int n){
  137.     if (n==1)return void(b[0]=1);
  138.     Getexp(a,b,c,d,n>>1);
  139.     Getln(b,c,d,n);
  140.     for (int i=0;i<n;++i)tmp[i]=(a[i]-c[i]+P)%P,tmp[n+i]=0;++tmp[0];
  141.     int k;for (k=1;k<=n;k<<=1);
  142.     Conv(tmp,b,tmp,k);
  143.     for (int i=0;i<n;++i)b[i]=tmp[i],b[i+n]=0;
  144. }
  145. void got(int n,int deep)
  146. {
  147. 	int i;
  148. 	if(n==1)
  149. 	{
  150. 	    vec[deep].push_back(1);
  151. 		vec[deep].push_back(1);
  152. 		return;
  153. 	}
  154. 	got(n>>1,deep+1);
  155. 	vec[deep].resize(n+1);
  156. 	int d=n>>1;
  157. 	for(i=0;i<=n;i++)
  158. 	{
  159. 		if(i==0)
  160. 		   ex[i]=1;
  161. 		else
  162. 		   ex[i]=1LL*ex[i-1]*d%mod; 
  163. 	}
  164. 	int len=1;
  165. 	while(len<=d+d)
  166. 	   len<<=1;
  167. 	for(i=0;i<len;i++)
  168. 	{
  169. 	    if(i<=(n>>1))
  170. 	    {
  171. 	       aa[i]=1LL*ex[(n>>1)-i]*inv[(n>>1)-i]%mod;
  172. 		   bb[i]=1LL*fac[i]*vec[deep+1][i]%mod;
  173. 		}
  174. 		else
  175. 		   bb[i]=aa[i]=0;	
  176. 	}
  177. 	Conv(aa,bb,cc,len);
  178. 	for(i=0;i<=d;i++)
  179. 		cc[i]=1LL*cc[i+(n>>1)]*inv[i]%mod;
  180. 	for(i=d+1;i<len;i++)
  181. 	    cc[i]=0;
  182. 	for(i=0;i<len;i++)
  183. 	{
  184. 		if(i<=(n>>1))
  185. 		   aa[i]=vec[deep+1][i];
  186. 		else
  187. 		   aa[i]=0;
  188. 	}
  189. 	Conv(aa,cc,bb,len);
  190. 	if(n&1)
  191. 	{
  192. 		for(i=0;i<=n;i++)
  193. 		{
  194. 			vec[deep][i]=1LL*bb[i]*n%mod;
  195. 			vec[deep][i]=(vec[deep][i]+bb[i-1])%mod;
  196. 		}
  197. 	}
  198. 	else
  199. 	{
  200. 		for(i=0;i<=n;i++)
  201. 		   vec[deep][i]=bb[i];
  202. 	}
  203. }
  204. int CC(int n,int m)
  205. {
  206. 	return 1LL*fac[n]*inv[m]%mod*inv[n-m]%mod;
  207. }
  208. long long S,T;
  209. int n,m;
  210. int main()
  211. {
  212. 	Init();
  213. 	scanf("%lld%lld%d%d",&S,&T,&n,&m);
  214. 	int i,j,s,p,q;
  215. 	for(i=0;i<=m-n+1;i++)
  216. 	{
  217. 		if(i==0)
  218. 		  fac[i]=1;
  219. 		else
  220. 		  fac[i]=1LL*fac[i-1]*i%mod;
  221. 	}
  222. 	inv[m-n+1]=power(fac[m-n+1],mod-2);
  223. 	for(i=m-n;i>=0;i--)
  224. 	   inv[i]=1LL*inv[i+1]*(i+1)%mod;
  225. 	int len=1;
  226. 	while(len<=m-n)
  227. 	   len<<=1;
  228.     got(m-n,0);
  229.     S=(S-(m-n)+mod)%mod;
  230.     T%=mod;
  231.     memset(a,0,sizeof(a));
  232.     memset(b,0,sizeof(b));
  233.     for(i=1;i<=m-n;i++)
  234.     	b[i]=inv[i+1];
  235. 	b[0]=1;
  236. 	GetInv(b,a,len);  
  237. 	for(i=0;i<=m-n+1;i++)
  238. 	{
  239. 		if(i==0)
  240. 		   ex[i]=1;
  241. 		else
  242. 		   ex[i]=1LL*ex[i-1]*(T+1)%mod;
  243. 	}
  244. 	for(i=0;i<=m-n;i++)
  245. 	{
  246. 	    b[i]=1LL*(ex[i+1]-1)*inv[i+1]%mod;
  247.         if(b[i]<0)
  248.            b[i]+=mod;
  249. 	}
  250. 	for(i=m-n+1;i<len<<1;i++)
  251. 		a[i]=b[i]=0;
  252. 	Conv(a,b,f,len<<1);
  253. 	for(i=0;i<=m-n;i++)
  254. 	{
  255. 		f[i]=1LL*f[i]*fac[i]%mod;
  256. 		if(f[i]<0)
  257. 		   f[i]+=mod;
  258. 	}
  259. 	memset(a,0,sizeof(a));
  260. 	memset(b,0,sizeof(b));
  261. 	memset(c,0,sizeof(c));
  262. 	int x,y;
  263. 	for(i=0;i<=m-n;i++)
  264. 	{
  265. 		a[i]=1LL*f[i]*inv[i]%mod;
  266. 		if(i%2)
  267. 		   a[i]=(mod-a[i])%mod;
  268. 		if(i==0)
  269. 		{
  270. 			x=power(a[i],mod-2);
  271. 			y=power(a[i],n);
  272. 		}
  273. 		a[i]=1LL*a[i]*x%mod;
  274. 	}
  275. 	Getln(a,b,c,len);
  276. 	for(i=0;i<=m-n;i++)
  277. 	   b[i]=1LL*b[i]*n%mod;
  278. 	for(i=m-n+1;i<len;i++)
  279. 	   b[i]=0;
  280. 	memset(a,0,sizeof(a));
  281. 	memset(c,0,sizeof(c));
  282. 	memset(d,0,sizeof(d));
  283. 	Getexp(b,a,c,d,len);
  284. 	for(i=0;i<=m-n;i++)
  285. 	   a[i]=1LL*a[i]*y%mod;
  286. 	int awp=1;
  287. 	for(i=0;i<=m-n;i++)
  288. 	{
  289. 	    b[i]=1LL*awp*inv[i]%mod;
  290. 	    awp=1LL*awp*S%mod;
  291. 	}
  292. 	for(i=m-n+1;i<len<<1;i++)
  293. 	   a[i]=b[i]=0;
  294. 	Conv(a,b,c,len<<1);
  295. 	int ans=0;
  296. 	for(i=0;i<=m-n;i++)
  297. 		ans=(ans+1LL*c[i]*vec[0][i]%mod*fac[i])%mod;
  298. 	ans=1LL*ans*inv[m-n]%mod;
  299. 	printf("%d\n",ans);
  300. 	return 0;
  301. }
  302.  
Success #stdin #stdout 0s 15512KB
comments (?)
stdin
copy 
1000000000000000000 100000 1000000000 1000001000
stdout
copy 
846395340
https://ideone.com/DifbHi
language:
C++ (gcc 8.3)
created:
5 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!