标题:c语言写的完整程序看不懂请老师逐行解释一下,谢谢!
取消只看楼主
ysr2857
Rank: 6Rank: 6
等 级:贵宾
威 望:28
帖 子:767
专家分:65
注 册:2020-2-10
结帖率:100%
已结贴  问题点数:20 回复次数:9 
c语言写的完整程序看不懂请老师逐行解释一下,谢谢!
下面是网上找到的用vc编程的fft快速高精度乘法程序,我看不懂,请老师看看,逐行解释一下,或翻译成vb可调用程序,谢谢!
#include<bits/stdc++.h>
using namespace std;
//complex是stl自带的定义复数的容器
typedef complex<double> cp;
#define N 2097153
//pie表示圆周率π
const double pie=acos(-1);
int n;
cp a[N],b[N];
int rev[N],ans[N];
char s1[N],s2[N];
//读入优化
int read(){
    int sum=0,f=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){sum=(sum<<3)+(sum<<1)+ch-'0';ch=getchar();}
    return sum*f;
}
//初始化每个位置最终到达的位置
{
    int len=1<<k;
    for(int i=0;i<len;i++)
    rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
}
//a表示要操作的系数,n表示序列长度
//若flag为1,则表示FFT,为-1则为IFFT(需要求倒数)
void fft(cp *a,int n,int flag){
    for(int i=0;i<n;i++)
    {
     //i小于rev[i]时才交换,防止同一个元素交换两次,回到它原来的位置。
      if(i<rev[i])swap(a[i],a[rev[i]]);
    }
    for(int h=1;h<n;h*=2)//h是准备合并序列的长度的二分之一
    {
    cp wn=exp(cp(0,flag*pie/h));//求单位根w_n^1
     for(int j=0;j<n;j+=h*2)//j表示合并到了哪一位
     {
      cp w(1,0);
       for(int k=j;k<j+h;k++)//只扫左半部分,得到右半部分的答案
       {
         cp x=a[k];
         cp y=w*a[k+h];
         a[k]=x+y;  //这两步是蝴蝶变换
         a[k+h]=x-y;
         w*=wn; //求w_n^k
       }
     }
     }
     //判断是否是FFT还是IFFT
     if(flag==-1)
     for(int i=0;i<n;i++)
     a[i]/=n;
}
int main(){
    n=read();
    scanf("%s%s",s1,s2);
    //读入的数的每一位看成多项式的一项,保存在复数的实部
    for(int i=0;i<n;i++)a[i]=(double)(s1[n-i-1]-'0');
    for(int i=0;i<n;i++)b[i]=(double)(s2[n-i-1]-'0');
    //k表示转化成二进制的位数
    int k=1,s=2;
     while((1<<k)<2*n-1)k++,s<<=1;
    init(k);
    //FFT 把a的系数表示转化为点值表示
    fft(a,s,1);
    //FFT 把b的系数表示转化为点值表示
    fft(b,s,1);
    //FFT 两个多项式的点值表示相乘
    for(int i=0;i<s;i++)
    a[i]*=b[i];
    //IFFT 把这个点值表示转化为系数表示
    fft(a,s,-1);
    //保存答案的每一位(注意进位)
    for(int i=0;i<s;i++)
    {
    //取实数四舍五入,此时虚数部分应当为0或由于浮点误差接近0
    ans[i]+=(int)(a[i].real()+0.5);
    ans[i+1]+=ans[i]/10;
    ans[i]%=10;
    }
    while(!ans[s]&&s>-1)s--;
    if(s==-1)printf("0");
    else
    for(int i=s;i>=0;i--)
    printf("%d",ans[i]);
    return 0;
}
搜索更多相关主题的帖子: 表示 i++ for sum int 
2020-03-29 17:58
ysr2857
Rank: 6Rank: 6
等 级:贵宾
威 望:28
帖 子:767
专家分:65
注 册:2020-2-10
得分:0 
请老师看一下这个程序能运行吗?是否真的完整?
2020-03-29 22:59
ysr2857
Rank: 6Rank: 6
等 级:贵宾
威 望:28
帖 子:767
专家分:65
注 册:2020-2-10
得分:0 
回复 3楼 wmf2014
谢谢您的指导!
vb是可以表示复数的,用到两个或三个数组,写成c=a() “+”b() “i”的形式。
fft的原理我也不懂,转化过程中把带小数点的数值四舍五入都取整数,好像是可以保证整数都正确都准确无误的。
还有一个法是数论变换法,不用复数,直接得到整数。
你发的这个链接我慢慢看看,对我来说很有用,谢谢!
2020-03-30 00:11
ysr2857
Rank: 6Rank: 6
等 级:贵宾
威 望:28
帖 子:767
专家分:65
注 册:2020-2-10
得分:0 
回复 3楼 wmf2014
我发到vb里那篇文章中有个分治法乘法,我试了一下不能用,不知道哪里有错?您有空再看看吧,谢谢!
2020-03-30 00:15
ysr2857
Rank: 6Rank: 6
等 级:贵宾
威 望:28
帖 子:767
专家分:65
注 册:2020-2-10
得分:0 
回复 6楼 wmf2014
嗯嗯好,是两回事。真不懂,有人说是可以翻译的,不过可能是猜测?谢谢指导!
2020-03-30 00:50
ysr2857
Rank: 6Rank: 6
等 级:贵宾
威 望:28
帖 子:767
专家分:65
注 册:2020-2-10
得分:0 
这个程序能运行吗?速度是否快?结果是否正确?
2020-03-30 08:09
ysr2857
Rank: 6Rank: 6
等 级:贵宾
威 望:28
帖 子:767
专家分:65
注 册:2020-2-10
得分:0 
回复 9楼 wmf2014
谢谢您!这个结果不对吧?不知道是哪儿有问题了。
2020-03-30 11:15
ysr2857
Rank: 6Rank: 6
等 级:贵宾
威 望:28
帖 子:767
专家分:65
注 册:2020-2-10
得分:0 
回复 9楼 wmf2014
又在网上搜索了一个程序,请您再看看是否完整?分两种,递归型:(105行,黏贴就乱了,整理一下)
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <math.h>
#include <conio.h>
#define N 150010const double pi = 3.141592653;
char s1[N>>1], s2[N>>1];
double rea[N], ina[N], reb[N], inb[N], Retmp[N], Intmp[N];
int ans[N>>1];
void FFT(double *reA, double *inA, int n, int flag)
{
    if(n == 1) return;    int k, u, i;
    double reWm = cos(2*pi/n), inWm = sin(2*pi/n);
    if(flag) inWm = -inWm;
    double reW = 1.0, inW = 0.0;
    /* 把下标为偶数的值按顺序放前面,下标为奇数的值按顺序放后面 */
    for(k = 1,u = 0; k < n; k += 2,u++)
{
        Retmp[u] = reA[k];
        Intmp[u] = inA[k];
    }
    for(k = 2; k < n; k += 2)
{
        reA[k/2] = reA[k];
        inA[k/2] = inA[k];
    }
    for(k = u,i = 0;
 k < n && i < u; k++, i++)
{
        reA[k] = Retmp[i];
        inA[k] = Intmp[i];
    }
    /* 递归处理 */
    FFT(reA, inA, n/2, flag);
    FFT(reA + n/2, inA + n/2, n/2, flag);
    for(k = 0; k < n/2; k++)
{
        int tag = k+n/2;
        double reT = reW * reA[tag] - inW * inA[tag];
        double inT = reW * inA[tag] + inW * reA[tag];
        double reU = reA[k], inU = inA[k];
        reA[k] = reU + reT;
        inA[k] = inU + inT;
        reA[tag] = reU - reT;
        inA[tag] = inU - inT;
        double rew_t = reW * reWm - inW * inWm;
         double inw_t = reW * inWm + inW * reWm;
         reW = rew_t;
        inW = inw_t;
    }
}
 int main()
{
#if 0
    freopen("in.txt","r",stdin);
#endif
    while(~scanf("%s%s", s1, s2))
{
        memset(ans, 0 , sizeof(ans));
        memset(rea, 0 , sizeof(rea));
        memset(ina, 0 , sizeof(ina));
        memset(reb, 0 , sizeof(reb));
        memset(inb, 0 , sizeof(inb));
        /* 计算长度为 2 的幂的长度len */
        int i, lent, len = 1, len1, len2;
        len1 = strlen(s1);
        len2 = strlen(s2);
        lent = (len1 > len2 ? len1 : len2);
        while(len < lent) len <<= 1;
        len <<= 1;
        /* 系数反转并添加 0 使长度凑成 2 的幂 */
        for(i = 0; i < len; i++)
{
            if(i < len1) rea[i] = (double)s1[len1-i-1] - '0';
            if(i < len2) reb[i] = (double)s2[len2-i-1] - '0';
            ina[i] = inb[i] = 0.0;
        }
        /* 分别把向量 a, 和向量 b 的系数表示转化为点值表示 */
        FFT(rea, ina, len, 0);
        FFT(reb, inb, len, 0);
        /* 点值相乘得到向量 c 的点值表示 */
        for(i = 0; i < len; i++)
{
            double rec = rea[i] * reb[i] - ina[i] * inb[i];
            double inc = rea[i] * inb[i] + ina[i] * reb[i];
           rea[i] = rec; ina[i] = inc;
        }
        /* 向量 c 的点值表示转化为系数表示 */
        FFT(rea, ina, len, 1);
        for(i = 0; i < len; i++)
{
            rea[i] /= len;
            ina[i] /= len;
        }
        /* 进位 */
        for(i = 0; i < len; i++)
            ans[i] = (int)(rea[i] + 0.5);
        for(i = 0; i < len; i++)
{
            ans[i+1] += ans[i] / 10;
            ans[i] %= 10;
       }
        int len_ans = len1 + len2 + 2;
        while(ans[len_ans] == 0 && len_ans > 0) len_ans--;
        for(i = len_ans; i >= 0; i--)
            printf("%d", ans[i]);
        printf("\n");
    }
   return 0;
}
————————————————
版权声明:本文为CSDN博主「o-pqy-o」的原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接及本声明。
原文链接:https://blog.
2020-03-31 09:03
ysr2857
Rank: 6Rank: 6
等 级:贵宾
威 望:28
帖 子:767
专家分:65
注 册:2020-2-10
得分:0 
迭代型:(共119行,乱了,再整理一下)
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <math.h>
#include <conio.h>
#define N 150010const double pi = 3.141592653;
char s1[N>>1], s2[N>>1];
double rea[N], ina[N], reb[N], inb[N];
int ans[N>>1];
 void Swap(double *x, double *y)
{
    double t = *x;
    *x = *y;
    *y = t;
}
 int Rev(int x, int len)
{
    int ans = 0;
    int i;
   for(i = 0;
 i < len; i++)
{
       ans <<= 1;
        ans |= (x & 1);
        x >>= 1;
    }
    return ans;
}
 void FFT(double *reA, double *inA, int n, bool flag)
{
   int s;
    double lgn = log((double)n) / log((double)2);
    int i;
    for(i = 0; i < n; i++)
{  
      int j = Rev(i, lgn);
        if(j > i)
{
            Swap(&reA[i], &reA[j]);
            Swap(&inA[i], &inA[j]);
        }
    }
   for(s = 1; s <= lgn; s++)
{   
     int m = (1<<s);
        double reWm = cos(2*pi/m), inWm = sin(2*pi/m);
        if(flag) inWm = -inWm;
        int k;
        for(k = 0; k < n; k += m)
{  
          double reW = 1.0, inW = 0.0;
            int j;
           for(j = 0; j < m / 2; j++)
{  
              int tag = k+j+m/2;
                double reT = reW * reA[tag] - inW * inA[tag];
                double inT = reW * inA[tag] + inW * reA[tag];
               double reU = reA[k+j], inU = inA[k+j];
                reA[k+j] = reU + reT;
                inA[k+j] = inU + inT;
                reA[tag] = reU - reT;
               inA[tag] = inU - inT;
               double rew_t = reW * reWm - inW * inWm;
                 double inw_t = reW * inWm + inW * reWm;
                reW = rew_t;
               inW = inw_t;
           }
        }
    }
    if(flag)
{
        for(i = 0;
 i < n; i++)
{
            reA[i] /= n;
            inA[i] /= n;
        }
   }
}
 int main()
{
#if 0
   freopen("in.txt","r",stdin);
#endif
   while(~scanf("%s%s", s1, s2))
{
       memset(ans, 0 , sizeof(ans));
        memset(rea, 0 , sizeof(rea));
        memset(ina, 0 , sizeof(ina));
        memset(reb, 0 , sizeof(reb));
        memset(inb, 0 , sizeof(inb));
        int i, lent, len = 1, len1, len2;
        len1 = strlen(s1);
        len2 = strlen(s2);
        lent = (len1 > len2 ? len1 : len2);
        while(len < lent) len <<= 1;
        len <<= 1;
        for(i = 0;
 i < len; i++)
{
           if(i < len1) rea[i] = (double)s1[len1-i-1] - '0';
            if(i < len2) reb[i] = (double)s2[len2-i-1] - '0';
            ina[i] = inb[i] = 0.0;
        }
        FFT(rea, ina, len, 0);
        FFT(reb, inb, len, 0);
        for(i = 0; i < len; i++)
{
           double rec = rea[i] * reb[i] - ina[i] * inb[i];
            double inc = rea[i] * inb[i] + ina[i] * reb[i];
            rea[i] = rec; ina[i] = inc;
        }
        FFT(rea, ina, len, 1);
        for(i = 0; i < len; i++)
           ans[i] = (int)(rea[i] + 0.4);
        for(i = 0; i < len; i++)
{
           ans[i+1] += ans[i] / 10;
            ans[i] %= 10;
        }
        int len_ans = len1 + len2 + 2;
        while(ans[len_ans] == 0 && len_ans > 0)
 len_ans--;
        for(i = len_ans; i >= 0; i--)
            printf("%d", ans[i]);
       printf("\n");
    }
    return 0;
}
————————————————
迭代型比递归型稍快一点。
2020-03-31 09:18
ysr2857
Rank: 6Rank: 6
等 级:贵宾
威 望:28
帖 子:767
专家分:65
注 册:2020-2-10
得分:0 
回复 9楼 wmf2014
接通知结贴,把积分送给你!
2020-04-01 08:12



参与讨论请移步原网站贴子:https://bbs.bccn.net/thread-500707-1-1.html




关于我们 | 广告合作 | 编程中国 | 清除Cookies | TOP | 手机版

编程中国 版权所有,并保留所有权利。
Powered by Discuz, Processed in 0.125454 second(s), 8 queries.
Copyright©2004-2024, BCCN.NET, All Rights Reserved