大加减乘除 c源码 （测试）
本帖最后由 Ayala 于 2016-12-3 18:58 编辑编译 命令行
详见帖 https://www.0xaa55.com/forum.php?mod=viewthread&tid=1885&page=1#pid6571
有bug请反馈 不擅长c代码 写的太难看
/*
#include <Windows.h>
#include <WinBase.h>
#include <malloc.h>
*/
#ifndef Uint8B
typedef unsigned __int64 Uint8B;
#endif
#pragma pack(8)
/*
#if !defined(_LARGE_INTEGER) && !defined(LARGE_INTEGER)
typedef union _LARGE_INTEGER
{
struct
{
__int32 LowPart;
__int32 HighPart;
}u;
__int64 QuadPart;
}LARGE_INTEGER;
#endif
typedef struct
{
struct //
{
__int32 n;
__int32 max;
}len;
__int64*databuffer; //
void* filehandle;//
__int64 mod;
}polynomial;
typedef polynomial* ppolynomial;
*/
typedef struct
{
double real;
double img;
}complex;
typedef complex* pcomplex;
#pragma pack()
#ifdef cp
#undef cp
typedef complex cp;
#else
typedef complex cp;
#endif
#ifndef max
#define max(a,b) (((a) > (b)) ? (a) : (b))
#endif
#ifndef min
#define min(a,b) (((a) < (b)) ? (a) : (b))
#endif
#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
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);
void print(__int64 x[],int nx);
void convertd(__int64 a[],int na,__int64 b[],int *nb);
void trim(__int64 x[],int *nx);
__int64 polynomial_cmp(__int64 a[],int na,__int64 b[],int nb);
void polynomial_add(__int64 a[],int na,__int64 b[],int nb,__int64 c[],int *nc);
void polynomial_sub(__int64 a[],int na,__int64 b[],int nb,__int64 c[],int *nc);
void polynomial_mul(__int64 a[],int na,__int64 b[],int nb,__int64 c[],int *nc);
void polynomial_div(__int64 a[],int na,__int64 b[],int nb,__int64 c[],int *nc,__int64 d[],int *nd);
#define N 1<<16
__int64 a;
__int64 b;
__int64 c;
__int64 d;
__int64 t;
double PI;
#define _MOD 100
#define fmt "%02I64d"
#define LEN2
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);
}
__int64 IsTrue(__int64 a[],int na)
{ int i;__int64 t=0;
if (na)
{
for (i=0;i<na;i++) t |= a;
}
return t;
}
__int64 polynomial_cmp(__int64 a[],int na,__int64 b[],int nb)
{
static __int64 pe;int te;__int64 t;
polynomial_sub(a,na,b,nb,pe,&te);
if (0==te) t|=-1;
else t=(__int64)IsTrue(pe,te);
return t;
}
void polynomial_add(__int64 a[], int na, __int64 b[], int nb, __int64 c[], int *nc)
{
int i,j,k;__int64* pD;
if (na>nb)
{
k=nb;
j=na;
pD=a;
}
else
{
k=na;
j=nb;
pD=b;
}
*nc=j;
for (i=0;i<k;i++) c=a+b;
for (i=k;i<j;i++) c=pD;
//convertd(c,*nc,c,nc)
}
void polynomial_sub(__int64 a[], int na, __int64 b[], int nb, __int64 c[], int *nc)
{
int i,j=na,k=nb;
if (j<k)
{
*nc=0;
}
else
{
for (i=0;i<k;i++) c =a-b;
for ( ;i<j;i++) c =a;
for (i=0;i<j-1;i++)
{
//printf("c[%d]=%I64d\n",i,c);
if (c<0)
{
c+=_MOD;
c-=1;
}
}
trim(c,&j);
////printf("i=%d\n",i);
if (c<0) *nc=0;
else *nc=max(j,1);
////printf("i=%d\n",i);
}
}
void polynomial_mul(__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);
#ifdef x86_xp
c = _ftol2(u.real / n + 0.5);
#else
c = (__int64)(u.real / n + 0.5);
#endif ////printf("c[%d]I64d = %I64d\n",i,c);
}
for (i = na + nb - 1; i > 1 && !c; --i);
*nc = i;
}
void polynomial_div(__int64 a[], int na, __int64 b[], int nb, __int64 c[], int *nc,__int64 d[],int *nd)
{
static __int64 pa,pb,pc;
int i,j,k,t;
int ta,tca,tb,tc,te,rc;
__int64 t0,t1,t2,tt;
__int64 *pca;
signed __int64 flag;
if (na<nb)
{
c=0;
*nc =1;
for (i=0;i<na;i++) d=a;
*nd = na;
}
else
{
for (i=0,ta=na;i<na;i++) pa=a;
for (i=0,tb=nb;i<nb;i++) pb=b;
for (i=0,*nc=0;i<na;i++) c&=0;
for (k=ta-tb,j=k,*nc=k+1;j>=0;j--)
{
pca=&pa;
tca=tb;
flag = polynomial_cmp(pca,tca,pb,tb);
if (flag>0)//pca>pb
{
t0=pca;
tt=pb;
t1=t0 / (tt+1);
t2=t0 / tt;
convertd(pca,tca,pca,&tca);
//printf("line:%d t1=%I64d,t2=%I64d,tt=%I64d\nt0=%I64d\t\n",__LINE__,t1,t2,tt,t0);
tt=(t1+t2+1)/2;
while (t1!=t2)
{
tt=(t1+t2+1)/2;
pc=tt;tc=1;
//printf("line:%d t1=%I64d,t2=%I64d,tt=%I64d\n",__LINE__,t1,t2,tt);
//printf("line:%d,tb=%d,pb=\t",__LINE__,tb);
//print(pb,tb);
polynomial_mul(pc,tc,pb,tb,pc,&tc);
convertd(pc,tc,pc,&tc);
//printf("line:%d,tca=%d,pca=\t",__LINE__,tca);
//print(pca,tca);
//printf("line:%d,tc=%d,pc=\t",__LINE__,tc);
//print(pc,tc);
flag = polynomial_cmp(pca,tca,pc,tc);
if (flag==0)
{
break;
}
else if (flag>0)
{
t1=tt;
}
else
{
if (tt-t1<=1)
{
tt=t1;
break;
}
t2=tt;
}
}
/*
printf("line:%d,tt=%I64d\n",__LINE__,tt);
*/
c=tt;
pc=tt;
tc=1;
polynomial_mul(pc,tc,pb,tb,pc,&tc);
convertd(pc,tc,pc,&tc);
/*
printf("line:%d,tca=%d,pca=\t",__LINE__,tca);
print(pca,tca);
printf("line:%d,tc=%d,pc=\t",__LINE__,tc);
print(pc,tc);
*/
polynomial_sub(pca,tca,pc,tc,pca,&tca);
convertd(pa,ta,pa,&ta);
trim(pa,&ta);
/*
printf("line:%d,tca=%d,pca=\t",__LINE__,tca);
print(pca,tca);
printf("line:%d,ta=%d,pa=\t",__LINE__,ta);
print(pa,ta);
*/
}
else if (flag==0)
{
c=1;
/* */
polynomial_sub(pca,tca,pb,tb,pca,&tca);
convertd(pa,ta,pa,&ta);
trim(pa,&ta);
}
else c=0;
if ((tca==tb)&&(ta>tb))
{
pa+=_MOD*pa;
pa&=0;
ta-=1;
j=max(j,1);
}
else if (tca>tb)
{
printf("line:%d Error\n",__LINE__);
system("pause");
}
//printf("*****************************************************\n");
}
convertd(pa,ta,pa,&ta);
for (i=0;i<ta;i++) d=pa;
*nd = max(i,1);
convertd(c,*nc,c,nc);
trim(c,nc);
//print(c,*nc);
//free(pheap);
}
}
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 trim(__int64 s[],int *ns)
{
int i,j;
j=*ns;
if (j>1)
{
while (j>1&&s==0) j--;
*ns=j;
}
}
void convertd(__int64 s[],int ns, __int64 a[], int * n)
{
int i,ii;__int64 t1,t2;
if (ns)
{
for (i=0,t1=0; i < ns; ++i)
{
t2 = t1 + s;
a = t2 % _MOD;
t1 = t2 / _MOD;
}
for (; t1; t1 /= _MOD) a = t1 % _MOD;
*n=i;
}
////printf("@262 r = %d\n",r);
}
void print(__int64 a[], int n)
{
if (n)
{
printf("%I64d", a[--n]);
while (n>0) printf(fmt, a[--n]);
}
else printf("Error\n");
puts("");
}
char buf;
#ifndef EOF
#define EOF (-1)
#endif
#ifndef false
#define false 0
#endif
#ifndef true
#define true 1
#endif
__int64 ddd (__int64 a[],int* na)
{
int i,k;__int64 t=0;
////printf("na=%d\n",*na);
if ((1>=a) && (1>=*na)) goto done;
for (i=0;i<*na,0==a;++i);
////printf("a[%d]=%I64d\n",i,a);
a--;
////printf("a[%d]=%I64d\n",i,a);
if (k=*na,0==a[--k]) *na-=1;
////printf("i=%d,na=%d\n",i,*na);
for (k=0;k<i;k++) a+=_MOD-1;
for (i=0;i<=*na;t |= a);
*na=max(*na,1);
//system("pause");
done:
return t;
}
int mainCRTStartup()
{
int i,na, nb, nc, nd, nt,sign,e;
__int64 t1, t2,M;double f1,f2;
long Ticks;
//conc();
printf("******************************************************\n");
printf("* *\n");
printf("******************************************************\n");
nc&=0;
while (scanf("%s", buf) != EOF)
{
if (buf== 'c')
{
system("cls");
continue;
}
else if (buf== 'q')
{
break;
}
else if ((nc<N) && ( \
(buf== '+') || \
(buf== '-') || \
(buf== '*') || \
(buf== '/') || \
(buf== '!') \
))
{
for (i=0;i<nc;i++) a=c;
na=nc;
}
else
{
convert(buf, a, &na);
scanf("%s", buf);
}
if (buf== '+')
{
scanf("%s", buf);
convert(buf, b, &nb);
polynomial_add(a, na, b, nb, c, &nc);
convertd(c,nc,c,&nc);
print(c, nc);
}
else if (buf== '-')
{
scanf("%s", buf);
convert(buf, b, &nb);
polynomial_sub(a, na, b, nb, c, &nc);
convertd(c,nc,c,&nc);
print(c, nc);
}
else if (buf== '*')
{
scanf("%s", buf);
convert(buf, b, &nb);
polynomial_mul(a, na, b, nb, c, &nc);
convertd(c,nc,c,&nc);
print(c, nc);
}
else if (buf== '/')
{
scanf("%s", buf);
convert(buf, b, &nb);
polynomial_div(a, na, b, nb, c, &nc,d,&nd);
convertd(c,nc,c,&nc);
convertd(d,nd,d,&nd);
print(c, nc);
print(d, nd);
}
else if (buf== '!')
{
if (na<=8)
{
t=a;
t=a;
nt=na;
while (ddd(t,&nt))
{
polynomial_mul(a, na, t, nt, c, &nc);
convertd(c,nc,a,&na);
}
convertd(c,nc,c,&nc);
print(c, nc);
}
}
else if (buf== 'c')
{
system("cls");
}
else if (buf== 'q')
{
break;
}
else
{
printf("+-*/!\n");
}
printf("******************************************************\n");
}
return 0;
}
建议将add mul divi等函数重命名为complex_add complex_mul complex_div，然后将这些不导出的函数用static修饰。 支持楼主。不擅长c代码都能写这么好！最主要是代码贴出来了很方便参考。
就是有些函数部分地方我看不懂，能加些注释就好了:) 支持楼主+1
页:
[1]