大数乘法c源码 -fft
本帖最后由 Ayala 于 2016-7-28 23:39 编辑用c写代码对我来说好痛苦!
double sin(double x);
double cos(double x);
double acos(double x);
double fmod(double x,double y);
double fabs(double x);
double floor(double x);
double ceil(double x);
typedef struct
{
double real;
double img;
}complex;
#ifdef cp
#undef cp
typedef complex cp;
#else
typedef complex cp;
#endif
#ifndef max
#define max(a,b) (((a) > (b)) ? (a) : (b))
#endif
#define N 1<<16
__int64 a;
__int64 b;
__int64 c;
double cc;
double PI;
const int LEN = 2, MOD = 100;
const double MODf = 100.0;
#ifdef x86_xp
#ifndef _ftol2
__int64 _ftol2(double x)
{
__int16 e,ee;
__int64 r;
__asm
{
fld qword ptr x
fstsw word ptr
mov ax,word ptr
or ah,0xc
mov word ptr ,ax
fldcw word ptr
fistp qword ptr
fldcw word ptr
mov eax,dword ptr
mov edx,dword ptr
}
}
#endif
#endif
void add(complex a,complex b,complex *c)
{
//_fpreset();
c->real =a.real + b.real;
c->img =a.img + b.img;
}
void sub(complex a,complex b,complex *c)
{
//_fpreset();
c->real=a.real-b.real;
c->img=a.img-b.img;
}
void mul(complex a,complex b,complex *c)
{
//_fpreset();
c->real =a.real*b.real- a.img*b.img;
c->img =a.real*b.img + a.img*b.real;
}
void divi(complex a,complex b,complex *c)
{
//_fpreset();
c->real =(a.real*b.real+a.img*b.img)/(b.real*b.real+b.img*b.img);
c->img =(a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img);
}
void bit_reverse_copy(cp a[], int n, cp b[])
{
int i, j, k, u, m;
for (u = 1, m = 0; u < n; u <<= 1, ++m);
for (i = 0; i < n; ++i)
{
j = i; k = 0;
for (u = 0; u < m; ++u, j >>= 1)
k = (k << 1) | (j & 1);
b = a;
//printf("b[%d]=%.2lf+%.2lfi\n",k,b.real,b.img);
}
}
void FFT(cp _x[], int n, int flag)
{
static cp x={0.0};
int i, j, k, kk, p, m;
double alpha;
cp wn,w,u,t,up,down;
//_fpreset();
bit_reverse_copy(_x, n, x);
alpha = 2.0 * acos(-1.0);
if (flag) alpha = -alpha;
for (i = 1, m = 0; i < n; i <<= 1, ++m);
for (i = 0, k = 2; i < m; ++i, k <<= 1)
{
wn.real = cos(alpha / k);
wn.img= sin(alpha / k);
for (j = 0; j < n; j += k)
{
w.real = 1.0;
w.img= 0.0;
kk = k >> 1;
for (p = 0; p < kk; ++p)
{
mul(w, x,&t);
add(x,t,&up);
sub(x,t,&down);
x = up;
x = down;
mul(w,wn,&u);
w = u;
}
}
}
memcpy(_x, x, sizeof(cp) * n);
}
void polynomial_multiply(__int64 a[], int na, __int64 b[], int nb, __int64 c[], int *nc)
{
static cp x, y,u;
int i, n;
//_fpreset();
i = max(na, nb) << 1;
for (n = 1; n < i; n <<= 1);
for (i = 0; i < na; ++i)
{
x.real = a;
x.img= 0;
//printf("x[%i]=%.2lf+%.2lfi\n",i,x.real,x.img);
}
for (; i < n; ++i)
{
x.real = 0;
x.img= 0;
}
FFT(x, n, 0);
for (i = 0; i < nb; ++i)
{
y.real = b;
y.img= 0;
//printf("y[%i]=%.2lf+%.2lfi\n",i,y.real,y.img);
}
for (; i < n; ++i)
{
y.real = 0;
y.img= 0;
}
FFT(y, n, 0);
for (i = 0; i < n; ++i)
{
//printf("x[%i]=%.2lf+%.2lfi\n",i,x.real,x.img);
//printf("y[%i]=%.2lf+%.2lfi\n",i,y.real,y.img);
mul(x,y,&u);
//printf("u[%i]=%.2lf+%.2lfi\n",i,u.real,u.img);
}
FFT(u, n, 1);
for (i = 0; i < n; ++i)
{
//printf("c[%d] .lf = %.53lf\n",i,u.real / n + 0.5);
c = _ftol2(u.real / n + 0.5);
//printf("c[%d]I64d = %I64d\n",i,c);
}
for (i = na + nb - 1; i > 1 && !c; --i);
*nc = i;
}
void convert(char *s, __int64 a[], int * n)
{
int len = strlen(s), i, j, k,r;
for (r = 0, i = len - LEN; i >= 0; i -= LEN)
{
for (j = k = 0; j < LEN; ++j)
{
k = k * 10 + (s - '0');
}
a = k;
//printf("@207 k = %d\n",k);
}
i += LEN;
if (i)
{
for (j = k = 0; j < i; ++j)
{
k = k * 10 + (s - '0');
}
a = k;
//printf("@217 k = %d\n",k);
}
*n=r;
//printf("@220 r = %d\n",r);
}
void print(__int64 a[], int n)
{
printf("%I64d", a[--n]);
while (n) printf("%02I64d", a[--n]);
puts("");
}
char buf;
#ifndef EOF
#define EOF (-1)
#endif
#ifndef false
#define false 0
#endif
#ifndef true
#define true 1
#endif
/*
void conc()
{
SetConsoleCP(936);
SetConsoleOutputCP(936);
SetConsoleMode(GetStdHandle(-10),0x7F);//STD_INPUT_HANDLE
}
*/
int mainCRTStartup()
{
int i,na, nb, nc,sign,e;
__int64 t1, t2;
double f1,f2;
//conc();
printf("******************************************************\n");
printf("* *\n");
printf("******************************************************\n");
while (scanf("%s", buf) != EOF)
{
sign = false;
if (buf == '-')
{
sign = !sign;
convert(buf + 1, a, &na);
}
else convert(buf, a, &na);
scanf("%s", buf);
if (buf == '-')
{
sign = !sign;
convert(buf + 1, b, &nb);
}
else convert(buf, b, &nb);
polynomial_multiply(a, na, b, nb, c, &nc);
t1 = 0;
for (i = 0; i < nc; ++i)
{
t2 = t1 + c;
c = t2 % MOD;
t1 = t2 / MOD;
}
for (; t1; t1 /= MOD) c = t1 % MOD;
if (sign) putchar('-');
print(c, nc);
printf("******************************************************\n");
}
return 0;
}
@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
咦，附件呢。。 ...
看样子我还得说明一下解压的办法，首先这后缀是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能不能打开我还不知道。。但我还是自己上传一个单压缩文件免得有人不会解压。。。
請問一下，除了加上這些preprocessor之外，還需要加上甚麼嗎? compile之後出現error
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <complex.h> 支持 ！！ 支持 ！！ 支持 ！！
页:
[1]