大数乘法c源码 -fft

Ayala

用c写代码对我来说好痛苦!

1. 
   
2. 
   
3. double sin(double x);
   
4. double cos(double x);
   
5. 
   
6. double acos(double x);
   
7. double fmod(double x,double y);
   
8. double fabs(double x);
   
9. double floor(double x);
   
10. double ceil(double x);
    
11. 
    
12. 
    
13. 
    
14. typedef struct
    
15. {
    
16.         double real;
    
17.         double img;
    
18. }complex;
    
19. 
    
20. #ifdef cp
    
21. #undef cp
    
22. typedef complex cp;
    
23. #else
    
24. typedef complex cp;
    
25. #endif
    
26. 
    
27. #ifndef max
    
28. #define max(a,b) (((a) > (b)) ? (a) : (b))
    
29. #endif
    
30. 
    
31. #define N 1<<16
    
32. 
    
33. 
    
34. __int64 a[N];
    
35. __int64 b[N];
    
36. __int64 c[N<<1];
    
37. double cc[N<<1];
    
38. double PI;
    
39. 
    
40. const int LEN = 2, MOD = 100;
    
41. const double MODf = 100.0;
    
42. 
    
43. 
    
44. #ifdef x86_xp
    
45. #ifndef _ftol2
    
46. 
    
47. __int64 _ftol2(double x)
    
48. {
    
49.   __int16 e,ee;
    
50.   __int64 r;
    
51.   __asm
    
52.   {
    
53.     fld qword ptr x
    
54.    
    
55.     fstsw word ptr [ee]
    
56.     mov ax,word ptr [ee]
    
57.     or ah,0xc
    
58.     mov word ptr [e],ax
    
59.    
    
60.     fldcw word ptr [e]
    
61.     fistp qword ptr [r]
    
62.     fldcw word ptr [ee]
    
63.    
    
64.     mov eax,dword ptr [r]
    
65.     mov edx,dword ptr [r+4]
    
66.   }
    
67. }
    
68. #endif
    
69. #endif
    
70. 
    
71. void   add(complex   a,complex   b,complex   *c)      
    
72. {
    
73.   //_fpreset();
    
74.         c->real        =a.real        +        b.real;      
    
75.         c->img        =a.img        +        b.img;      
    
76. }
    
77. 
    
78. void   sub(complex   a,complex   b,complex   *c)      
    
79. {
    
80.   //_fpreset();
    
81.         c->real=a.real-b.real;      
    
82.         c->img=a.img-b.img;      
    
83. }
    
84. 
    
85. void   mul(complex   a,complex   b,complex   *c)      
    
86. {      
    
87.   //_fpreset();
    
88.         c->real        =a.real*b.real  -   a.img*b.img;      
    
89.         c->img        =a.real*b.img   +   a.img*b.real;      
    
90. }      
    
91. 
    
92. 
    
93. void   divi(complex   a,complex   b,complex   *c)      
    
94. {      
    
95.   //_fpreset();
    
96.         c->real        =(a.real*b.real+a.img*b.img)/(b.real*b.real+b.img*b.img);      
    
97.         c->img        =(a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img);
    
98. }       
    
99. 
    
100.    
     
101. void bit_reverse_copy(cp a[], int n, cp b[])   
     
102. {   
     
103.     int i, j, k, u, m;   
     
104.     for (u = 1, m = 0; u < n; u <<= 1, ++m);   
     
105.     for (i = 0; i < n; ++i)   
     
106.     {   
     
107.         j = i; k = 0;   
     
108.         for (u = 0; u < m; ++u, j >>= 1)   
     
109.             k = (k << 1) | (j & 1);   
     
110.         b[k] = a[i];   
     
111.         //printf("b[%d]=%.2lf+%.2lfi\n",k,b[k].real,b[k].img);
     
112.     }   
     
113. }
     
114. 
     
115. 
     
116. 
     
117. void FFT(cp _x[], int n, int flag)   
     
118. {   
     
119.     static cp x[N << 1]={0.0};   
     
120.     int i, j, k, kk, p, m;   
     
121.     double alpha;   
     
122.    
     
123.     cp wn,w,u,t,up,down;
     
124.     //_fpreset();
     
125.    
     
126.     bit_reverse_copy(_x, n, x);
     
127.    
     
128.     alpha = 2.0 * acos(-1.0);
     
129.     if (flag) alpha = -alpha;   
     
130.    
     
131.     for (i = 1, m = 0; i < n; i <<= 1, ++m);   
     
132.    
     
133.     for (i = 0, k = 2; i < m; ++i, k <<= 1)   
     
134.     {   
     
135.         wn.real = cos(alpha / k);
     
136.         wn.img  = sin(alpha / k);
     
137.         for (j = 0; j < n; j += k)   
     
138.         {   
     
139.             w.real = 1.0;  
     
140.             w.img  = 0.0;
     
141.             kk = k >> 1;   
     
142.             for (p = 0; p < kk; ++p)   
     
143.             {      
     
144.                 mul(w, x[j + p + kk],&t);
     
145.                 add(x[j + p],t,&up);
     
146.                 sub(x[j + p],t,&down);
     
147.                
     
148.                 x[j + p]      = up;   
     
149.                 x[j + p + kk] = down;   
     
150.                
     
151.                 mul(w,wn,&u);
     
152.                 w = u;   
     
153.             }   
     
154.         }   
     
155.     }   
     
156.     memcpy(_x, x, sizeof(cp) * n);   
     
157. }   
     
158. 
     
159. 
     
160. void polynomial_multiply(__int64 a[], int na, __int64 b[], int nb, __int64 c[], int *nc)   
     
161. {   
     
162.     static cp x[N << 1], y[N << 1],u[N << 1];
     
163.     int i, n;
     
164.     //_fpreset();
     
165.     i = max(na, nb) << 1;   
     
166.     for (n = 1; n < i; n <<= 1);
     
167.    
     
168.     for (i = 0; i < na; ++i)
     
169.     {
     
170.       x[i].real = a[i];  
     
171.       x[i].img  = 0;
     
172.       //printf("x[%i]=%.2lf+%.2lfi\n",i,x[i].real,x[i].img);
     
173.     }
     
174.    
     
175.     for (; i < n; ++i)
     
176.     {
     
177.       x[i].real = 0;   
     
178.       x[i].img  = 0;
     
179.     }
     
180.    
     
181.     FFT(x, n, 0);   
     
182.    
     
183.    
     
184.     for (i = 0; i < nb; ++i)
     
185.     {
     
186.       y[i].real = b[i];   
     
187.       y[i].img  = 0;
     
188.       //printf("y[%i]=%.2lf+%.2lfi\n",i,y[i].real,y[i].img);
     
189.     }
     
190.    
     
191.     for (; i < n; ++i)
     
192.     {
     
193.       y[i].real = 0;   
     
194.       y[i].img  = 0;
     
195.     }
     
196.    
     
197.     FFT(y, n, 0);   
     
198.    
     
199.     for (i = 0; i < n; ++i)
     
200.     {
     
201.       //printf("x[%i]=%.2lf+%.2lfi\n",i,x[i].real,x[i].img);
     
202.       //printf("y[%i]=%.2lf+%.2lfi\n",i,y[i].real,y[i].img);
     
203.       mul(x[i],y[i],&u[i]);   
     
204.       //printf("u[%i]=%.2lf+%.2lfi\n",i,u[i].real,u[i].img);
     
205.     }
     
206.    
     
207.     FFT(u, n, 1);   
     
208.    
     
209.     for (i = 0; i < n; ++i)     
     
210.     {   
     
211.         //printf("c[%d] .lf = %.53lf\n",i,u[i].real / n + 0.5);
     
212.         c[i] = _ftol2(u[i].real / n + 0.5);  
     
213.         //printf("c[%d]  I64d = %I64d\n",i,c[i]);
     
214.     }   
     
215.     for (i = na + nb - 1; i > 1 && !c[i - 1]; --i);
     
216.     *nc = i;
     
217. }   
     
218. 
     
219. 
     
220. void convert(char *s, __int64 a[], int * n)   
     
221. {   
     
222.     int len = strlen(s), i, j, k,r;
     
223.     for (r = 0, i = len - LEN; i >= 0; i -= LEN)   
     
224.     {   
     
225.         for (j = k = 0; j < LEN; ++j)  
     
226.         {  
     
227.             k = k * 10 + (s[i + j] - '0');   
     
228.         }
     
229.         a[r++] = k;   
     
230.         //printf("@207 k = %d\n",k);
     
231.     }   
     
232.     i += LEN;   
     
233.     if (i)   
     
234.     {   
     
235.         for (j = k = 0; j < i; ++j)
     
236.         {
     
237.             k = k * 10 + (s[j] - '0');   
     
238.         }
     
239.         a[r++] = k;   
     
240.         //printf("@217 k = %d\n",k);
     
241.     }   
     
242.     *n=r;
     
243.     //printf("@220 r = %d\n",r);
     
244. }   
     
245.    
     
246. void print(__int64 a[], int n)   
     
247. {   
     
248.     printf("%I64d", a[--n]);   
     
249.     while (n) printf("%02I64d", a[--n]);   
     
250.     puts("");   
     
251. }   
     
252. 
     
253. 
     
254. char buf[N + 10];   
     
255. 
     
256. #ifndef EOF
     
257.   #define EOF (-1)
     
258. #endif
     
259. 
     
260. #ifndef false
     
261.   #define false 0
     
262. #endif
     
263. 
     
264. #ifndef true
     
265.   #define true 1
     
266. #endif
     
267. 
     
268. 
     
269. /*
     
270. void conc()
     
271. {
     
272.   
     
273.   SetConsoleCP(936);
     
274.   SetConsoleOutputCP(936);
     
275.   SetConsoleMode(GetStdHandle(-10),0x7F);//STD_INPUT_HANDLE
     
276. }
     
277. */
     
278. 
     
279. int mainCRTStartup()   
     
280. {   
     
281.     int i,na, nb, nc,sign,e;   
     
282.     __int64 t1, t2;   
     
283.     double f1,f2;
     
284.         
     
285.     //conc();
     
286.     printf("******************************************************\n");
     
287.     printf("*                                                    *\n");
     
288.     printf("******************************************************\n");
     
289.     while (scanf("%s", buf) != EOF)   
     
290.     {   
     
291.         sign = false;   
     
292.         if (buf[0] == '-')   
     
293.         {   
     
294.             sign = !sign;     
     
295.             convert(buf + 1, a, &na);   
     
296.         }   
     
297.         else convert(buf, a, &na);   
     
298.         scanf("%s", buf);   
     
299.         if (buf[0] == '-')   
     
300.         {   
     
301.             sign = !sign;   
     
302.             convert(buf + 1, b, &nb);   
     
303.         }   
     
304.         else convert(buf, b, &nb);   
     
305.         
     
306.    
     
307.           polynomial_multiply(a, na, b, nb, c, &nc);
     
308.            t1 = 0;   
     
309.             for (i = 0; i < nc; ++i)   
     
310.             {   
     
311.               t2 = t1 + c[i];   
     
312.               c[i] = t2 % MOD;   
     
313.               t1 = t2 / MOD;   
     
314.             }   
     
315.             for (; t1; t1 /= MOD) c[nc++] = t1 % MOD;   
     
316.             if (sign) putchar('-');   
     
317.             print(c, nc);
     
318.         
     
319.         printf("******************************************************\n");
     
320.     }   
     
321.     return 0;   
     
322. }
     
323. 
     
324.    

复制代码

@echo off
:re
cls
echo /*********************************************/
echo /     构建命令                                /
echo /*********************************************/
.\tools\x86\cl.exe .\src\hello_world.c /I".\inc\crt" /D"x86_xp" /I".\inc\api" /Fp:strict /Fa"Debug\hello_world.asm" /Fo"Debug\hello_world.obj" /c /MD

echo /*********************************************/
echo /                                             /
echo /*********************************************/


.\tools\x86\link.exe .\Debug\hello_world.obj /LIBPATH:".\lib\wxp\i386" /LIBPATH:".\lib\Crt\i386"  /OUT:"Debug\hello_world.exe" /NOLOGO /SUBSYSTEM:CONSOLE /MACHINE:X86 "kernel32.lib"
echo /*********************************************/
echo /                                             /
echo /*********************************************/


pause
goto re

0xAA55

咦，附件呢。。

Ayala

...
_hello_world.7z.001.zip (1000 KB, 下载次数: 10)

2016-7-29 19:51 上传

点击文件名下载附件

_hello_world.7z.002.zip (1000 KB, 下载次数: 3)

2016-7-29 19:51 上传

点击文件名下载附件

_hello_world.7z.003.zip (1000 KB, 下载次数: 3)

2016-7-29 19:51 上传

点击文件名下载附件

_hello_world.7z.004.zip (909.34 KB, 下载次数: 3)

2016-7-29 19:51 上传

点击文件名下载附件

0xAA55

看样子我还得说明一下解压的办法，首先这后缀是zip但它不是zip压缩包，因此改后缀，去掉所有的“.zip”。
那么4个文件，应该是这样的文件名：
_hello_world.7z.001
_hello_world.7z.002
_hello_world.7z.003
_hello_world.7z.004
然后选中第一个文件，也就是_hello_world.7z.001，右键用7z打开就行了。
RAR能不能打开我还不知道。。但我还是自己上传一个单压缩文件免得有人不会解压。。。
_hello_world.7z (3.81 MB, 下载次数: 34)

2016-7-30 15:51 上传

点击文件名下载附件

elephantcaptain

請問一下，除了加上這些preprocessor之外，還需要加上甚麼嗎? compile之後出現error
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <complex.h>

jasonchen

支持    ！！

jasonchen

支持    ！！

jasonchen

支持    ！！