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[N+2];
  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)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",a[i]);
  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. }
Time limit exceeded #stdin #stdout 5s 313152KB
stdin
2^57885161
stdout
Standard output is empty