fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include <time.h>
  5. #define PI 3.1415926535897932384626433832795
  6. #define base 10000
  7. #define K 23
  8. #define N (1<<K)
  9. void LMULT3(short a[],short b[],short c[],int k,double d[],double e[],double f[],long long t[]);
  10. void LSQUARE(short a[],short c[],int k,double d[],double f[],long long t[]);
  11. void SUBn(short a[],short b[],short c[],int n);
  12.  
  13. int main(void){
  14.     clock_t start,end;
  15. 	static short a[N],b[N],c[N];
  16. 	static long long t[N];
  17. 	static double d[N],e[N],f[N];
  18. 	int i,k,x,y;
  19. 	double r;
  20. 	printf("x^yの形で入力\n");
  21. 	scanf("%d^%d",&x,&y);
  22. 	start=clock();
  23. 	r=log10(x);
  24. 	k=(int)ceil(log2(ceil(2*r/4)));
  25. 	a[0]=1;
  26. 	b[0]=x%base;
  27. 	b[1]=x/base;
  28. 	c[0]=1;
  29. 	i=1;
  30. 	while(1){
  31. 		if(y&1)LMULT3(a,b,a,k,d,e,f,t);
  32. 		y>>=1;
  33. 		if(!y)break;
  34. 		LSQUARE(b,b,k,d,e,t);
  35. 		i*=2;
  36. 		k=(int)ceil(log2(ceil(2*i*r/4)));
  37. 	}
  38. 	if(!k)k++;
  39. 	SUBn(a,c,a,k);
  40. 	end=clock();
  41. 	printf("%.3fs\n",(double)(end-start)/CLOCKS_PER_SEC);
  42. 	for(i=4;i--;)printf("%04d",a[i]);
  43. 	return 0;
  44. }
  45.  
  46. void SUBn(short a[],short b[],short c[],int n){
  47. 	int i;
  48. 	int t=0;
  49. 	for(i=0;i<n;i++){
  50. 		t+=a[i]-b[i]+base;
  51. 		c[i]=t%base;
  52. 		t/=base;
  53. 		t--;
  54. 	}
  55. 	return;
  56. }
  57.  
  58. int tib(int a,int k){
  59. 	int i=k,z=0;
  60. 	while(i--){
  61. 		z=(z<<1)+(a&1);
  62. 		a>>=1;
  63. 	}
  64. 	return z;
  65. }
  66.  
  67. void INV2(double d[],short c[],int k,long long t[]){
  68. 	int i,k2=1<<k;
  69. 	for(i=(1<<k)/2;i--;)d[1+i*2]*=-1;
  70. 	for(i=0;i<k2;i++)t[i]=(long long)(d[i]/(k2/2)+0.5);
  71. 	for(i=0;i<k2;i++){
  72. 		c[i]=t[i]%base;
  73. 		t[i+1]+=t[i]/base;
  74. 	}
  75. 	return;
  76. }
  77.  
  78. void FFT(double d[],int k){
  79. 	int i,j,k5,k2=1<<k,cnt;
  80. 	double t[2],r,pr;
  81. 	for(i=0;i<k2*2;i+=4){
  82. 		t[0]=d[i+2];
  83. 		t[1]=d[i+3];
  84. 		d[i+2]=d[i+0]-t[0];
  85. 		d[i+3]=d[i+1]-t[1];
  86. 		d[i+0]+=t[0];
  87. 		d[i+1]+=t[1];
  88. 	}
  89. 	for(i=0;i<k2*2;i+=8){
  90. 		t[0]=d[i+4];
  91. 		t[1]=d[i+5];
  92. 		d[i+4]=d[i+0]-t[0];
  93. 		d[i+5]=d[i+1]-t[1];
  94. 		d[i+0]+=t[0];
  95. 		d[i+1]+=t[1];
  96. 		t[0]=       +d[i+7];
  97. 		t[1]=-d[i+6]       ;
  98. 		d[i+6]=d[i+2]-t[0];
  99. 		d[i+7]=d[i+3]-t[1];
  100. 		d[i+2]+=t[0];
  101. 		d[i+3]+=t[1];
  102. 	}
  103. 	for(j=2;j<k;j++){
  104. 		k5=1<<j;
  105. 		pr=PI/k5;
  106. 		for(i=0;i<k2;i+=k5){
  107. 			for(cnt=0;cnt<k5;cnt++){
  108. 				r=pr*cnt;
  109. 				t[0]= cos(r)*d[0+(i+k5)*2]+sin(r)*d[1+(i+k5)*2];
  110. 				t[1]=-sin(r)*d[0+(i+k5)*2]+cos(r)*d[1+(i+k5)*2];
  111. 				d[0+(i+k5)*2]=d[0+i*2]-t[0];
  112. 				d[1+(i+k5)*2]=d[1+i*2]-t[1];
  113. 				d[0+i*2]+=t[0];
  114. 				d[1+i*2]+=t[1];
  115. 				i++;
  116. 			}
  117. 		}
  118. 	}
  119. 	return;
  120. }
  121.  
  122. void PRE2(short a[],double d[],int k){
  123. 	int i,j;
  124. 	for(i=0;i<(1<<k)*2;i++){
  125. 		j=tib(i,k);
  126. 		d[0+i*2]=(double)a[0+j*2];
  127. 		d[1+i*2]=(double)a[1+j*2];
  128. 	}
  129. 	return;
  130. }
  131.  
  132. void MID3(double d[],double e[],double f[],int k){
  133. 	int i,j,k2;
  134. 	double s[4],t[2],u[2],r,pr;
  135. 	k2=1<<k;
  136. 	f[0+0*2]=d[0+0*2]*e[0+0*2]+d[1+0*2]*e[1+0*2];
  137. 	f[1+0*2]=d[0+0*2]*e[1+0*2]+d[1+0*2]*e[0+0*2];
  138. 	pr=-PI/(k2/4);
  139. 	for(i=1,j=k2/2-1;i<k2/4;i++,j--){
  140. 		r=i*pr;
  141. 		u[0]=cos(r)+1;
  142. 		u[1]=sin(r);
  143. 		s[0]=d[0+i*2]-d[0+j*2];
  144. 		s[1]=d[1+i*2]+d[1+j*2];
  145. 		s[2]=e[0+i*2]-e[0+j*2];
  146. 		s[3]=e[1+i*2]+e[1+j*2];
  147. 		t[0]=-s[0]*s[3]-s[1]*s[2];
  148. 		t[1]=+s[0]*s[2]-s[1]*s[3];
  149. 		f[0+i*2]=-0.25*(+t[0]*u[1]+t[1]*u[0])+d[0+i*2]*e[0+i*2]-d[1+i*2]*e[1+i*2];
  150. 		f[1+i*2]=+0.25*(-t[1]*u[1]+t[0]*u[0])+d[0+i*2]*e[1+i*2]+d[1+i*2]*e[0+i*2];
  151. 		f[0+j*2]=-0.25*(+t[0]*u[1]+t[1]*u[0])+d[0+j*2]*e[0+j*2]-d[1+j*2]*e[1+j*2];
  152. 		f[1+j*2]=-0.25*(-t[1]*u[1]+t[0]*u[0])+d[0+j*2]*e[1+j*2]+d[1+j*2]*e[0+j*2];
  153. 	}
  154. 	f[0+k2/2]=d[0+k2/2]*e[0+k2/2]-d[1+k2/2]*e[1+k2/2];
  155. 	f[1+k2/2]=(d[0+k2/2]*e[1+k2/2]+d[1+k2/2]*e[0+k2/2]);
  156. 	for(i=0;i<k2/2;i++){
  157. 		j=tib(i,k-1);
  158. 		d[0+j*2]=f[0+i*2];
  159. 		d[1+j*2]=-f[1+i*2];
  160. 	}
  161. 	return;
  162. }
  163.  
  164. void LMULT3(short a[],short b[],short c[],int k,double d[],double e[],double f[],long long t[]){
  165. 	if(k<3)k=3;
  166. 	PRE2(a,d,k-1);
  167. 	FFT(d,k-1);
  168. 	PRE2(b,e,k-1);
  169. 	FFT(e,k-1);
  170. 	MID3(d,e,f,k);
  171. 	FFT(d,k-1);
  172. 	INV2(d,c,k,t);
  173. 	return;
  174. }
  175.  
  176. void LSQUARE(short a[],short c[],int k,double d[],double f[],long long t[]){
  177. 	if(k<3)k=3;
  178. 	PRE2(a,d,k-1);
  179. 	FFT(d,k-1);
  180. 	MID3(d,d,f,k);
  181. 	FFT(d,k-1);
  182. 	INV2(d,c,k,t);
  183. 	return;
  184. }
Time limit exceeded #stdin #stdout 5s 313152KB
comments (?)
stdin
copy 
2^57885161
stdout
copy 
Standard output is empty
https://ideone.com/9jNwyP
language:
C (gcc 8.3)
created:
12 years ago
visibility:
secret

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!