技术宅的结界

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

QQ登录

只需一步,快速开始

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

大数乘法c源码 -fft

[复制链接]

25

主题

81

帖子

1088

积分

用户组: 版主

UID
1821
精华
6
威望
57 点
宅币
832 个
贡献
31 次
宅之契约
0 份
在线时间
196 小时
注册时间
2016-7-12
发表于 2016-7-28 23:38:46 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 Ayala 于 2016-7-28 23:39 编辑

用c写代码对我来说好痛苦!
[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[N];
__int64 b[N];
__int64 c[N<<1];
double cc[N<<1];
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 [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

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);    
}    


void polynomial_multiply(__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);
        c[i] = _ftol2(u[i].real / n + 0.5);  
        //printf("c[%d]  I64d = %I64d\n",i,c[i]);
    }    
    for (i = na + nb - 1; i > 1 && !c[i - 1]; --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[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 print(__int64 a[], int n)    
{    
    printf("%I64d", a[--n]);    
    while (n) printf("%02I64d", a[--n]);    
    puts("");    
}    


char buf[N + 10];    

#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[0] == '-')    
        {    
            sign = !sign;     
            convert(buf + 1, a, &na);    
        }    
        else convert(buf, a, &na);    
        scanf("%s", buf);    
        if (buf[0] == '-')    
        {    
            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[i];    
              c[i] = t2 % MOD;    
              t1 = t2 / MOD;    
            }    
            for (; t1; t1 /= MOD) c[nc++] = 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

评分

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

查看全部评分

本帖被以下淘专辑推荐:

995

主题

2207

帖子

5万

积分

用户组: 管理员

一只技术宅

UID
1
精华
197
威望
261 点
宅币
16459 个
贡献
32323 次
宅之契约
0 份
在线时间
1565 小时
注册时间
2014-1-26
发表于 2016-7-29 17:21:35 | 显示全部楼层
咦,附件呢。。
回复

使用道具 举报

25

主题

81

帖子

1088

积分

用户组: 版主

UID
1821
精华
6
威望
57 点
宅币
832 个
贡献
31 次
宅之契约
0 份
在线时间
196 小时
注册时间
2016-7-12
 楼主| 发表于 2016-7-29 19:52:57 | 显示全部楼层
...
_hello_world.7z.001.zip (1000 KB, 下载次数: 10)

995

主题

2207

帖子

5万

积分

用户组: 管理员

一只技术宅

UID
1
精华
197
威望
261 点
宅币
16459 个
贡献
32323 次
宅之契约
0 份
在线时间
1565 小时
注册时间
2014-1-26
发表于 2016-7-30 15:41:49 | 显示全部楼层
看样子我还得说明一下解压的办法,首先这后缀是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, 下载次数: 11)

0

主题

1

帖子

13

积分

用户组: 初·技术宅

UID
1920
精华
0
威望
1 点
宅币
10 个
贡献
0 次
宅之契约
0 份
在线时间
1 小时
注册时间
2016-8-24
发表于 2016-8-25 13:53:00 | 显示全部楼层
請問一下,除了加上這些preprocessor之外,還需要加上甚麼嗎? compile之後出現error
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <complex.h>

0

主题

24

帖子

54

积分

用户组: 小·技术宅

UID
2067
精华
0
威望
1 点
宅币
28 个
贡献
0 次
宅之契约
0 份
在线时间
0 小时
注册时间
2016-11-17
发表于 2016-11-17 10:31:59 | 显示全部楼层
支持    !!

0

主题

24

帖子

54

积分

用户组: 小·技术宅

UID
2067
精华
0
威望
1 点
宅币
28 个
贡献
0 次
宅之契约
0 份
在线时间
0 小时
注册时间
2016-11-17
发表于 2016-11-17 10:32:20 | 显示全部楼层
支持    !!

0

主题

24

帖子

54

积分

用户组: 小·技术宅

UID
2067
精华
0
威望
1 点
宅币
28 个
贡献
0 次
宅之契约
0 份
在线时间
0 小时
注册时间
2016-11-17
发表于 2016-11-17 10:32:45 | 显示全部楼层
支持    !!

本版积分规则

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

GMT+8, 2018-9-19 01:26 , Processed in 0.122109 second(s), 44 queries , Gzip On.

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

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