大加减乘除 c源码（测试）

Ayala · 发表于 2016-8-3 15:42:58

欢迎访问技术宅的结界，请注册或者登录吧。

您需要登录才可以下载或查看，没有账号？立即注册→加入我们

×

本帖最后由 Ayala 于 2016-12-3 18:58 编辑

编译命令行
详见帖 https://www.0xaa55.com/forum.php ... &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 [ee]
mov ax,word ptr [ee]
or ah,0xc
mov word ptr [e],ax
fldcw word ptr [e]
fistp qword ptr [r]
fldcw word ptr [ee]
mov eax,dword ptr [r]
mov edx,dword ptr [r+4]
}
}
#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[N];
__int64 b[N];
__int64 c[N<<1];
__int64 d[N];
__int64 t[4];
double PI;
#define _MOD 100
#define fmt "%02I64d"
#define LEN 2
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[k] = a[i];
////printf("b[%d]=%.2lf+%.2lfi\n",k,b[k].real,b[k].img);
}
}
void FFT(cp _x[], int n, int flag)
{
static cp x[N << 1]={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[j + p + kk],&t);
add(x[j + p],t,&up);
sub(x[j + p],t,&down);
x[j + p] = up;
x[j + p + kk] = 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[i];
}
return t;
}
__int64 polynomial_cmp(__int64 a[],int na,__int64 b[],int nb)
{
static __int64 pe[N];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[i]=a[i]+b[i];
for (i=k;i<j;i++) c[i]=pD[i];
//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[i] =a[i]-b[i];
for ( ;i<j;i++) c[i] =a[i];
for (i=0;i<j-1;i++)
{
//printf("c[%d]=%I64d\n",i,c[i]);
if (c[i]<0)
{
c[i]+=_MOD;
c[i+1]-=1;
}
}
trim(c,&j);
////printf("i=%d\n",i);
if (c[i]<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[N << 1], y[N << 1],u[N << 1];
int i, n;
//_fpreset();
i = max(na, nb) << 1;
for (n = 1; n < i; n <<= 1);
for (i = 0; i < na; ++i)
{
x[i].real = a[i];
x[i].img = 0;
////printf("x[%i]=%.2lf+%.2lfi\n",i,x[i].real,x[i].img);
}
for (; i < n; ++i)
{
x[i].real = 0;
x[i].img = 0;
}
FFT(x, n, 0);
for (i = 0; i < nb; ++i)
{
y[i].real = b[i];
y[i].img = 0;
////printf("y[%i]=%.2lf+%.2lfi\n",i,y[i].real,y[i].img);
}
for (; i < n; ++i)
{
y[i].real = 0;
y[i].img = 0;
}
FFT(y, n, 0);
for (i = 0; i < n; ++i)
{
////printf("x[%i]=%.2lf+%.2lfi\n",i,x[i].real,x[i].img);
////printf("y[%i]=%.2lf+%.2lfi\n",i,y[i].real,y[i].img);
mul(x[i],y[i],&u[i]);
////printf("u[%i]=%.2lf+%.2lfi\n",i,u[i].real,u[i].img);
}
FFT(u, n, 1);
for (i = 0; i < n; ++i)
{
////printf("c[%d] .lf = %.53lf\n",i,u[i].real / n + 0.5);
#ifdef x86_xp
c[i] = _ftol2(u[i].real / n + 0.5);
#else
c[i] = (__int64)(u[i].real / n + 0.5);
#endif ////printf("c[%d] I64d = %I64d\n",i,c[i]);
}
for (i = na + nb - 1; i > 1 && !c[i - 1]; --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[N],pb[N],pc[N];
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]=0;
*nc =1;
for (i=0;i<na;i++) d[i]=a[i];
*nd = na;
}
else
{
for (i=0,ta=na;i<na;i++) pa[i]=a[i];
for (i=0,tb=nb;i<nb;i++) pb[i]=b[i];
for (i=0,*nc=0;i<na;i++) c[i]&=0;
for (k=ta-tb,j=k,*nc=k+1;j>=0;j--)
{
pca=&pa[j];
tca=tb;
flag = polynomial_cmp(pca,tca,pb,tb);
if (flag>0)//pca>pb
{
t0=pca[tca-1];
tt=pb[tb-1];
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[0]=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[j]=tt;
pc[0]=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[j]=1;
/* */
polynomial_sub(pca,tca,pb,tb,pca,&tca);
convertd(pa,ta,pa,&ta);
trim(pa,&ta);
}
else c[j]=0;
if ((tca==tb)&&(ta>tb))
{
pa[ta-2]+=_MOD*pa[ta-1];
pa[ta-1]&=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[i]=pa[i];
*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[i + j] - '0');
}
a[r++] = k;
////printf("@207 k = %d\n",k);
}
i += LEN;
if (i)
{
for (j = k = 0; j < i; ++j)
{
k = k * 10 + (s[j] - '0');
}
a[r++] = 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[j-1]==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[i];
a[i] = t2 % _MOD;
t1 = t2 / _MOD;
}
for (; t1; t1 /= _MOD) a[i++] = 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[N + 10];
#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[0]) && (1>=*na)) goto done;
for (i=0;i<*na,0==a[i];++i);
////printf("a[%d]=%I64d\n",i,a[i]);
a[i]--;
////printf("a[%d]=%I64d\n",i,a[i]);
if (k=*na,0==a[--k]) *na-=1;
////printf("i=%d,na=%d\n",i,*na);
for (k=0;k<i;k++) a[k]+=_MOD-1;
for (i=0;i<=*na;t |= a[i++]);
*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[0]== 'c')
{
system("cls");
continue;
}
else if (buf[0]== 'q')
{
break;
}
else if ((nc<N) && ( \
(buf[0]== '+') || \
(buf[0]== '-') || \
(buf[0]== '*') || \
(buf[0]== '/') || \
(buf[0]== '!') \
))
{
for (i=0;i<nc;i++) a[i]=c[i];
na=nc;
}
else
{
convert(buf, a, &na);
scanf("%s", buf);
}
if (buf[0]== '+')
{
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[0]== '-')
{
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[0]== '*')
{
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[0]== '/')
{
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[0]== '!')
{
if (na<=8)
{
t[0]=a[0];
t[1]=a[1];
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[0]== 'c')
{
system("cls");
}
else if (buf[0]== 'q')
{
break;
}
else
{
printf("+-*/!\n");
}
printf("******************************************************\n");
}
return 0;
}

复制代码

0xAA55 · 发表于 2016-8-5 17:43:10

建议将add mul divi等函数重命名为complex_add complex_mul complex_div，然后将这些不导出的函数用static修饰。

human_rights · 发表于 2018-1-14 10:18:48

支持楼主。不擅长c代码都能写这么好！最主要是代码贴出来了很方便参考。
就是有些函数部分地方我看不懂，能加些注释就好了

7KY6 · 发表于 2018-1-14 12:51:48

支持楼主+1

账号		自动登录	找回密码
密码			立即注册→加入我们

大加减乘除 c源码 （测试）

欢迎访问技术宅的结界，请注册或者登录吧。

本帖被以下淘专辑推荐:

大加减乘除 c源码（测试）