fork download 
  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 LMULT(short a[],short b[],short c[],int k,double d[],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[2*N],e[2*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(y){
  31. 		if(y&1)LMULT(a,b,a,k,d,e,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. 	SUBn(a,c,a,k);
  39. 	end=clock();
  40. 	printf("%.3fs\n",(double)(end-start)/CLOCKS_PER_SEC);
  41. 	//for(i=N;i--;)printf("%04d")
  42. 	return 0;
  43. }
  44.  
  45. void SUBn(short a[],short b[],short c[],int n){
  46. 	int i;
  47. 	int t=0;
  48. 	for(i=0;i<n;i++){
  49. 		t+=a[i]-b[i]+base;
  50. 		c[i]=t%base;
  51. 		t/=base;
  52. 		t--;
  53. 	}
  54. 	return;
  55. }
  56.  
  57. int tib(int a,int k){
  58. 	int i=k,z=0;
  59. 	while(i--){
  60. 		z=(z<<1)+(a&1);
  61. 		a>>=1;
  62. 	}
  63. 	return z;
  64. }
  65.  
  66. void PRE(short a[],short b[],double d[],int k){
  67. 	int i,j;
  68. 	for(i=0;i<(1<<k);i++){
  69. 		j=tib(i,k);
  70. 		d[0+i*2]=(double)a[j];
  71. 		d[1+i*2]=(double)b[j];
  72. 	}
  73. 	return;
  74. }
  75.  
  76. void MID(double d[],double f[],int k){
  77. 	int i,j,k2;
  78. 	k2=1<<k;
  79. 	f[0+0*2]=d[0+0*2]*d[1+0*2];
  80. 	f[1+0*2]=0.0;
  81. 	for(i=1,j=k2-1;i<k2/2;i++,j--){
  82. 		f[0+i*2]=(d[0+i*2]*d[1+i*2]+d[0+j*2]*d[1+j*2])*0.5;
  83. 		f[1+i*2]=(-d[0+i*2]*d[0+i*2]+d[1+i*2]*d[1+i*2]+d[0+j*2]*d[0+j*2]-d[1+j*2]*d[1+j*2])*0.25;
  84. 	}
  85. 	f[0+k2]=d[0+k2]*d[1+k2];
  86. 	f[1+k2]=0.0;
  87. 	i=0;
  88. 	j=tib(i,k);
  89. 	d[0+j*2]=f[0+i*2];
  90. 	d[1+j*2]=-f[1+i*2];
  91. 	i=k2/2;
  92. 	j=tib(k2-i,k);
  93. 	d[0+j*2]=f[0+i*2];
  94. 	d[1+j*2]=-f[1+i*2];
  95. 	for(i=1;i<k2/2;){
  96. 		j=tib(i,k);
  97. 		d[0+j*2]=f[0+i*2];
  98. 		d[1+j*2]=-f[1+i*2];
  99. 		j=tib(k2-i,k);
  100. 		d[0+j*2]=f[0+i*2];
  101. 		d[1+j*2]=f[1+i*2];
  102. 		i++;
  103. 	}
  104. 	return;
  105. }
  106.  
  107. void INV(double d[],short c[],int k,long long t[]){
  108. 	int i,k2=1<<k;
  109. 	for(i=0;i<k2;i++)t[i]=(long long)round(d[0+i*2]/k2);
  110. 	for(i=0;i<k2;i++){
  111. 		c[i]=t[i]%base;
  112. 		t[i+1]+=t[i]/base;
  113. 	}
  114. 	return;
  115. }
  116.  
  117. void FFT(double d[],int k){
  118. 	int i,j,k5,k2=1<<k,cnt;
  119. 	double t[2],r,pr;
  120. 	for(i=0;i<k2*2;i+=4){
  121. 		t[0]=d[i+2];
  122. 		t[1]=d[i+3];
  123. 		d[i+2]=d[i+0]-t[0];
  124. 		d[i+3]=d[i+1]-t[1];
  125. 		d[i+0]+=t[0];
  126. 		d[i+1]+=t[1];
  127. 	}
  128. 	for(i=0;i<k2*2;i+=8){
  129. 		t[0]=d[i+4];
  130. 		t[1]=d[i+5];
  131. 		d[i+4]=d[i+0]-t[0];
  132. 		d[i+5]=d[i+1]-t[1];
  133. 		d[i+0]+=t[0];
  134. 		d[i+1]+=t[1];
  135. 		t[0]=       +d[i+7];
  136. 		t[1]=-d[i+6]       ;
  137. 		d[i+6]=d[i+2]-t[0];
  138. 		d[i+7]=d[i+3]-t[1];
  139. 		d[i+2]+=t[0];
  140. 		d[i+3]+=t[1];
  141. 	}
  142. 	for(j=2;j<k;j++){
  143. 		k5=1<<j;
  144. 		pr=PI/k5;
  145. 		for(i=0;i<k2;i+=k5){
  146. 			for(cnt=0;cnt<k5;cnt++){
  147. 				r=pr*cnt;
  148. 				t[0]= cos(r)*d[0+(i+k5)*2]+sin(r)*d[1+(i+k5)*2];
  149. 				t[1]=-sin(r)*d[0+(i+k5)*2]+cos(r)*d[1+(i+k5)*2];
  150. 				d[0+(i+k5)*2]=d[0+i*2]-t[0];
  151. 				d[1+(i+k5)*2]=d[1+i*2]-t[1];
  152. 				d[0+i*2]+=t[0];
  153. 				d[1+i*2]+=t[1];
  154. 				i++;
  155. 			}
  156. 		}
  157. 	}
  158. 	return;
  159. }
  160.  
  161. void PRE2(short a[],double d[],int k){
  162. 	int i,j;
  163. 	for(i=0;i<(1<<k)*2;i++){
  164. 		j=tib(i,k);
  165. 		d[0+i*2]=(double)a[0+j*2];
  166. 		d[1+i*2]=(double)a[1+j*2];
  167. 	}
  168. 	return;
  169. }
  170.  
  171.  
  172. void MID2(double d[],double e[],double f[],int k){
  173. 	int i,j,k2;
  174. 	double s[2],t[2],u[2],v[2],w[2],x[2],y[2],r,pr;
  175. 	k2=1<<k;
  176. 	f[0+0*2]=(d[0+0*2]+d[1+0*2])*(e[0+0*2]+e[1+0*2]);
  177. 	f[1+0*2]=0.0;
  178. 	f[0+k2 ]=(d[0+0*2]-d[1+0*2])*(e[0+0*2]-e[1+0*2]);
  179. 	f[1+k2 ]=0.0;
  180. 	pr=PI/(k2/2);
  181. 	for(i=1,j=k2/2-1;i<k2/4;i++,j--){
  182. 		r=i*pr;
  183. 		u[0]=0.5*(1+sin(r));
  184. 		u[1]=0.5*cos(r);
  185. 		s[0]=d[0+i*2]-d[0+j*2];
  186. 		s[1]=d[1+i*2]+d[1+j*2];
  187. 		t[0]=e[0+i*2]-e[0+j*2];
  188. 		t[1]=e[1+i*2]+e[1+j*2];
  189. 		v[0]=u[0]*s[0]-u[1]*s[1];
  190. 		v[1]=u[0]*s[1]+u[1]*s[0];
  191. 		w[0]=u[0]*t[0]-u[1]*t[1];
  192. 		w[1]=u[0]*t[1]+u[1]*t[0];
  193. 		x[0]=d[0+i*2]-v[0];
  194. 		x[1]=d[1+i*2]-v[1];
  195. 		y[0]=e[0+i*2]-w[0];
  196. 		y[1]=e[1+i*2]-w[1];
  197. 		f[0+i*2]=x[0]*y[0]-x[1]*y[1];
  198. 		f[1+i*2]=x[0]*y[1]+x[1]*y[0];
  199. 		x[0]=d[0+j*2]+v[0];
  200. 		x[1]=-(-d[1+j*2]+v[1]);
  201. 		y[0]=e[0+j*2]+w[0];
  202. 		y[1]=-(-e[1+j*2]+w[1]);
  203. 		f[0+j*2]=x[0]*y[0]-x[1]*y[1];
  204. 		f[1+j*2]=x[0]*y[1]+x[1]*y[0];
  205. 	}
  206. 	f[0+k2/2]=d[0+k2/2]*e[0+k2/2]-d[1+k2/2]*e[1+k2/2];
  207. 	f[1+k2/2]=-(d[0+k2/2]*e[1+k2/2]+d[1+k2/2]*e[0+k2/2]);
  208. 	i=0;
  209. 	j=tib(i,k);
  210. 	d[0+j*2]=f[0+i*2];
  211. 	d[1+j*2]=-f[1+i*2];
  212. 	i=k2/2;
  213. 	j=tib(k2-i,k);
  214. 	d[0+j*2]=f[0+i*2];
  215. 	d[1+j*2]=-f[1+i*2];
  216. 	for(i=1;i<k2/2;){
  217. 		j=tib(i,k);
  218. 		d[0+j*2]=f[0+i*2];
  219. 		d[1+j*2]=-f[1+i*2];
  220. 		j=tib(k2-i,k);
  221. 		d[0+j*2]=f[0+i*2];
  222. 		d[1+j*2]=f[1+i*2];
  223. 		i++;
  224. 	}
  225. 	return;
  226. }
  227.  
  228. void LMULT(short a[],short b[],short c[],int k,double d[],double f[],long long t[]){
  229. 	if(k<3)k=3;
  230. 	PRE(a,b,d,k);
  231. 	FFT(d,k);
  232. 	MID(d,f,k);
  233. 	FFT(d,k);
  234. 	INV(d,c,k,t);
  235. 	return;
  236. }
  237.  
  238. void LSQUARE(short a[],short c[],int k,double d[],double f[],long long t[]){
  239. 	if(k<3)k=3;
  240. 	PRE2(a,d,k-1);
  241. 	FFT(d,k-1);
  242. 	MID2(d,d,f,k);
  243. 	FFT(d,k);
  244. 	INV(d,c,k,t);
  245. 	return;
  246. }
Not running #stdin #stdout 0s 0KB
comments (?)
stdin 
2^57885161
stdout 
Standard output is empty
https://ideone.com/5U2gT0
language:
C (gcc 8.3)
created:
11 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!