技术宅的结界

 找回密码
 立即注册→加入我们

QQ登录

只需一步,快速开始

搜索
热搜: 下载 VB C 实现 编写
查看: 888|回复: 3
收起左侧

大加减乘除 c源码 (测试)

[复制链接]

25

主题

91

帖子

1209

积分

用户组: 版主

UID
1821
精华
6
威望
57 点
宅币
938 个
贡献
36 次
宅之契约
0 份
在线时间
210 小时
注册时间
2016-7-12
发表于 2016-8-3 15:42:58 | 显示全部楼层 |阅读模式

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

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

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

编译 命令行
详见帖 https://www.0xaa55.com/forum.php ... &page=1#pid6571
有bug请反馈 不擅长c代码 写的太难看
[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;    
}

hello_world.c

15.64 KB, 下载次数: 2

评分

参与人数 1威望 +10 宅币 +20 贡献 +10 收起 理由
0xAA55 + 10 + 20 + 10 加油!

查看全部评分

本帖被以下淘专辑推荐:

1010

主题

2246

帖子

5万

积分

用户组: 管理员

一只技术宅

UID
1
精华
200
威望
265 点
宅币
16901 个
贡献
33733 次
宅之契约
0 份
在线时间
1608 小时
注册时间
2014-1-26
发表于 2016-8-5 17:43:10 | 显示全部楼层
建议将add mul divi等函数重命名为complex_add complex_mul complex_div,然后将这些不导出的函数用static修饰。

1

主题

6

帖子

93

积分

用户组: 小·技术宅

UID
3349
精华
1
威望
4 点
宅币
74 个
贡献
0 次
宅之契约
0 份
在线时间
7 小时
注册时间
2018-1-14
发表于 2018-1-14 10:18:48 | 显示全部楼层
支持楼主。不擅长c代码都能写这么好!最主要是代码贴出来了很方便参考。
就是有些函数部分地方我看不懂,能加些注释就好了
000111000 000111000 000111000 000111000 000111000 UN

0

主题

41

帖子

45

积分

用户组: 初·技术宅

UID
3351
精华
0
威望
2 点
宅币
0 个
贡献
0 次
宅之契约
0 份
在线时间
0 小时
注册时间
2018-1-14
发表于 2018-1-14 12:51:48 | 显示全部楼层
支持楼主+1

本版积分规则

QQ|申请友链|Archiver|手机版|小黑屋|技术宅的结界 ( 滇ICP备16008837号|网站地图

GMT+8, 2018-12-16 07:36 , Processed in 0.104875 second(s), 40 queries , Gzip On.

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表