注册 登录
编程论坛 VB6论坛

求助:如何把利用快速傅里叶变换的大数乘法变成vb程序?

ysr2857 发布于 2020-08-30 23:47, 6557 次点击
貌似一种新的乘法快速计算方法已经提交论文,理论上可以达到大数乘法的效率极限。
文章链接:多项式乘法到快速傅里叶变换;此文介绍的非常详细,极力推荐。
文章链接:使用快速傅里叶变换计算大整数乘法;快速傅里叶变换,使用算法设计思想中的分治法,降低傅里叶变换的时间复杂度到 O(N logN)。
傅里叶变换,算法的时间复杂度还是 O(N2)。关键在于:直接进行离散傅里叶变换的计算复杂度是 O(N2)。快速傅里叶变换可以计算出与直接计算相同的结果,但只需要 O(N logN) 的计算复杂度。 N logN 和 N2 之间的差别是巨大的。例如,当 N = 106 时,在一个每秒运算百万次的计算机上,粗略地说,它们之间就是占用 30 秒 CPU 时间和两星期 CPU 时间的差别。
快速傅里叶变换的要点如下:一个界长为 N 的离散傅里叶变换可以重新写成两个界长各为 N/2 的离散傅里叶变换之和。其中一个变换由原来 N 个点中的偶数点构成,另一个变换由奇数点构成。这个过程可以递归地进行下去,直到我们将全部数据细分为界长为 1 的变换。
什么是界长为 1 的傅里叶变换呢?它正是把一个输入值复制成它的一个输出值的恒等运算。要实现以上算法,最容易的情况是原始的 N 为 2 的整幂次项,如果数据集的界长不是 2 的幂次时,则可添上一些零值,直到 2 的下一幂次。在这个算法中,每递归一次需 N 阶运算,共需要 log N 次递归,所以快速傅里叶变换算法的时间复杂度是 O(N logN)。

关键是把各位数字看成了多项式系数?如何把多项式表示法转换为点值表示法?
系数的数字都不是相等的,不规则变化的,与单位圆上均匀分布的点值不同,半径相等?
如何利用对称性,周期性?

不懂,大大的不懂!原理都不懂,咋变成程序?
47 回复
#2
ysr28572020-09-05 21:51
好像明白了一点:
我们以两个 3 ( N = 3 ) 位数 a = 678 和 b = 432 的乘积来说明以上过程吧。
a = (678)'10 = 6x10^2 + 7x10^1 + 8x10^0
b = (432)'10 = 4x10^2 + 3x10^1 + 2x10^0
由此:
c = a x b = 10^4xc4 + 10^3xc3 + 10^2xc2 + 10^1xc1 + 10^0xc0
   = 10000x24 + 1000x46 + 100x65 + 10x38 + 1x16
   = 292896
如果按以上方法计算大整数的乘法,时间复杂度是 O(N^2)。
在我们的例子中,乘积 c = 292896,共 6 位数字,N 需要扩展到 2^3 = 8。那么,向量 {ai} 和向量 {bj} 如下所示:
{ a7, a6, a5, a4, a3, a2, a1, a0 } = { 0, 0, 0, 0, 0, 6, 7, 8 }
{ b7, b6, b5, b4, b3, b2, b1, b0 } = { 0, 0, 0, 0, 0, 4, 3, 2 }

ω^0~ω^7分别是1,0.7-0.7i,-i,-0.7-0.7i,-1,-0.7+0.7i,i,0.7+0.7i.
可见ω^1和ω^7是共轭的,ω^2和ω^6是共轭的,等,这就是对称性。
ω^0和ω^8是相等的,ω^4和ω^12是相等的,等,这就是周期性。

向量A0就是8,7,6依次乘以ω^0,再加起来得到21.
以此类推得到:
{ A7, A6, A5, A4, A3, A2, A1, A0 } = { 12.9+10.9i, 2+7i, 3.1-1.1i, 7, 3.1+1.1i, 2-7i, 12.9-10.9i, 21 }
{ B7, B6, B5, B4, B3, B2, B1, B0 } = { 4.1+6.1i, -2+3i, -0.1-1.9i, 3, -0.1+1.9i, -2-3i, 4.1-6.1i, 9 }

现在,将她们逐项相乘得到向量 {Ck},即 { C7, C6, C5, C4, C3, C2, C1, C0 }
= { -13.6+123.4i, -25-8i, -2.4-5.8i, 21, -2.4+5.8i, -25+8i, -13.6-123.4i, 189 }

这样,我们就使用离散傅里叶变换和逆变换计算出了向量 {ai} 和向量 {bj} 的卷积向量 {ck},如下所示:
{ c7, c6, c5, c4, c3, c2, c1, c0 } = { 0, 0, 0, 0, 24, 46, 65, 38, 16 }

这些两位数错位相加,就得到292896。

问题是如何利用对称性和周期性?使这个算法的时间复杂度由 O(N^2),变成O(N logN)?

各位老师,请给与指点!谢谢!

[此贴子已经被作者于2020-9-5 21:53编辑过]

#3
ysr28572020-09-05 22:03
注意到(可以观察到)向量Ai和向量Bj,都有共轭复数,对称出现的。
#4
ysr28572020-09-05 22:12
由向量Ai和向量Bj,得到向量Ck仍然需要相乘?虽然向量Ck也有共轭复数,且是对称出现的。
#5
ysr28572020-09-05 22:31
还可注意到:3*7=21,9*21=189.
所以,可能是对应项相乘的。

希望指点,帮忙编成vb程序!谢谢!
#6
ysr28572020-09-05 23:26
还可注意到:3*7=21,9*21=189.
所以,可能是对应项相乘的,就是由向量Ai和向量Bj,得到向量Ck是对应项相乘的。

希望指点,帮忙编成vb程序!谢谢!
#7
ysr28572020-09-07 21:48
问题点数成了0?
希望老师指点!欢迎探讨和沟通!
谢谢指点!如果程序应许,我可以给您加分!
#8
ysr28572020-09-07 23:20
如下程序怎么不运行?如何用呢?

'VB版快速傅里叶变换模块
'模块说明
'pr(),pi()为输入序列的实部与虚部,fr(),fi()分别是输出序列的实部与虚部。
'数组长度为n, n=2^k.
'L为0时做正变换,为1时做反变换。

Public Sub kkfft(ByRef pr() As Double, _
ByRef pi() As Double, _
ByRef fr() As Double, _
ByRef fi() As Double, _
ByVal n As Integer, _
ByVal k As Integer, ByVal L As Integer)

Dim it As Integer
Dim m As Integer
Dim iis As Integer
Dim i As Integer
Dim j As Integer
Dim nv As Integer
Dim l0 As Integer
Dim p, q, s, vr, vi, poddr, poddi As Double
For it = 0 To n - 1
m = it
iis = 0
For i = 0 To k - 1
j = Int(m / 2)
iis = 2 * iis + (m - 2 * j)
m = j
Next
fr(it) = pr(iis)
fi(it) = pi(iis)
Next
pr(0) = 1#
pi(0) = 0#
p = 6.283185306 / (1# * n)
pr(1) = Cos(p)
pi(1) = -Sin(p)
If L <> 0 Then pi(1) = -pi(1)
For i = 2 To n - 1
p = pr(i - 1) * pr(1)
q = pi(i - 1) * pi(1)
s = (pr(i - 1) + pi(i - 1)) * (pr(1) + pi(1))
pr(i) = p - q
pi(i) = s - p - q
Next
For it = 0 To n - 2 Step 2
vr = fr(it)
vi = fi(it)
fr(it) = vr + fr(it + 1)
fi(it) = vi + fi(it + 1)
fr(it + 1) = vr - fr(it + 1)
fi(it + 1) = vi - fi(it + 1)
Next
m = n / 2
nv = 2
For l0 = k - 2 To 0 Step -1
m = m / 2: nv = 2 * nv
For it = 0 To (m - 1) * nv Step nv
For j = 0 To (nv / 2) - 1
p = pr(m * j) * fr(it + j + nv / 2)
q = pi(m * j) * fi(it + j + nv / 2)
s = pr(m * j) + pi(m * j)
s = s * (fr(it + j + nv / 2) + fi(it + j + nv / 2))
poddr = p - q
poddi = s - p - q
fr(it + j + nv / 2) = fr(it + j) - podd
fi(it + j + nv / 2) = fi(it + j) - poddi
fr(it + j) = fr(it + j) + poddr
fi(it + j) = fi(it + j) + poddi
Next
Next
Next
If L <> 0 Then
For i = 0 To n - 1
fr(i) = fr(i) / (1# * n)
fi(i) = fi(i) / (1# * n)
Next
End If
End Sub

#9
ysr28572020-09-10 10:03
如下是VB版快速离散傅立叶变换,能不能用?怎么用?

Public Function fft(ByRef Data() As Double) As Double()
    ReDim ffft(128, 2) As Double
    Dim length As Integer
    length = UBound(Data, 1) + 1
'    Dim numArray(length - 1, 2) As Double
      
    Dim index As Integer
    Dim num5 As Integer
    Dim num6 As Integer
    Dim num7 As Integer
    Dim num10 As Integer
    Dim num3 As Integer
    Dim num2 As Integer
    Dim num11 As Integer
    Dim num9 As Integer
    num9 = length
      
    Dim num8 As Integer
    num8 = CInt(Math.Log(CDbl(num9)) / Math.Log(2#))
      
    Dim numArray2(128) As Double
    Dim numArray3(128) As Double
    Dim numArray4(128) As Double
    Dim numArray5(128) As Double
    For index = 0 To num9 - 1
        numArray2(index) = Data(index)
        numArray3(index) = 0#
    Next
    Dim a As Double
    Dim num14 As Double
      
   num14 = 6.28318530717959 / CDbl(num9)
    index = 0
    While index < (num9 \ 2)
        numArray4(index) = Math.Sin(a)
        numArray5(index) = Math.Cos(a)
        a = a + num14
        index = index + 1
    Wend
    num7 = num9
    num3 = 1
    For num2 = 1 To num8
        num7 = num7 / 2
        num6 = 0
        For num11 = 1 To num3
            num10 = 0
            index = num6
            While index <= ((num7 + num6) - 1)
                num5 = index + num7
                a = numArray2(index) - numArray2(num5)
                num14 = numArray3(index) - numArray3(num5)
                numArray2(index) = numArray2(index) + numArray2(num5)
                numArray3(index) = numArray3(index) + numArray3(num5)
                If num10 = 0 Then
                    numArray2(num5) = a
                    numArray3(num5) = num14
                Else
                    numArray2(num5) = (a * numArray5(num10)) + (num14 * numArray4(num10))
                    numArray3(num5) = (num14 * numArray5(num10)) - (a * numArray4(num10))
                End If
                num10 = num10 + num3
                index = index + 1
           Wend
            num6 = (num6 + num7) + num7
        Next
        num3 = num3 + num3
    Next
    num5 = num9 \ 2
    For index = 1 To (num9 - 1)
        num6 = num9
        If num5 < index Then
            Dim num12 As Double
              
          num12 = numArray2(index)
            numArray2(index) = numArray2(num5)
            numArray2(num5) = num12
            num12 = numArray3(index)
            numArray3(index) = numArray3(num5)
            numArray3(num5) = num12
        End If
        num6 = num6 / 2
        Do While num5 >= num6
            num5 = num5 - num6
            num6 = num6 / 2
            If num5 = 0 Then
                Exit Do
            End If
        Loop
        num5 = num5 + num6
    Next
    For index = 0 To num9 - 1
        numArray(index, 0) = numArray2(index)
        numArray(index, 1) = numArray3(index)
        numArray(index, 2) = ((numArray2(index)) ^ 2# + (numArray3(index)) ^ 2#) ^ 0.5
    Next
   fft = numArray
End Function
默认了数组为128




#10
ysr28572020-09-10 10:05
原来的说明是:输出结果1维是实部,2维是虚部。
#11
ysr28572020-09-11 23:02
如下是vc版FFT模板,能不能用?如何用?

void fft(cp *a,int n,int inv)
{
    int bit=0;
    while ((1<<bit)<n)bit++;
    fo(i,0,n-1)
    {
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
        if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交换两次(就是没交换)
    }
    for (int mid=1;mid<n;mid*=2)//mid是准备合并序列的长度的二分之一
    {
        cp temp(cos(pi/mid),inv*sin(pi/mid));//单位根,pi的系数2已经约掉了
        for (int i=0;i<n;i+=mid*2)//mid*2是准备合并序列的长度,i是合并到了哪一位
        {
            cp omega(1,0);
            for (int j=0;j<mid;j++,omega*=temp)//只扫左半部分,得到右半部分的答案
            {
                cp x=a[i+j],y=omega*a[i+j+mid];
                a[i+j]=x+y,a[i+j+mid]=x-y;//这个就是蝴蝶变换什么的
            }
        }
    }
}
#12
ysr28572020-09-13 20:42
9楼的可调用程序如何调用?老提示参数类型不对或缺少定义?
#13
ysr28572020-09-14 08:47
希望各位老师和好朋友继续给予指导!非常感谢!
#14
cwa99582020-09-14 08:59
建议:
你把相关的数据代码放上,出错的提示,运算的数据(数组)。
没有相关数据,别人怎么测试,
傅立叶变换,听起来就是高大上的,因为别人没有用到,就不会去造数据的。
#15
ysr28572020-09-14 09:32
没有相关数据,我就是输入一个数值计算,结果是数据(数组)类型不符?我是这样用的:
x=fft("123456789")
#16
ysr28572020-09-14 09:34
非常感谢老师的指导和帮助!欢迎批评!
不知道怎么弄了,该怎么用?
#17
ysr28572020-09-14 09:45
代码如下:
Public Function fft(ByRef Data() As Double) As Double()
    ReDim ffft(128, 2) As Double
    Dim length As Integer
    length = UBound(Data, 1) + 1
'    Dim numArray(length - 1, 2) As Double
      
    Dim index As Integer
    Dim num5 As Integer
    Dim num6 As Integer
    Dim num7 As Integer
    Dim num10 As Integer
    Dim num3 As Integer
    Dim num2 As Integer
    Dim num11 As Integer
    Dim num9 As Integer
    num9 = length
      
    Dim num8 As Integer
    num8 = CInt(Math.Log(CDbl(num9)) / Math.Log(2#))
      
    Dim numArray2(128) As Double
    Dim numArray3(128) As Double
    Dim numArray4(128) As Double
    Dim numArray5(128) As Double
    For index = 0 To num9 - 1
        numArray2(index) = Data(index)
        numArray3(index) = 0#
    Next
    Dim a As Double
    Dim num14 As Double
      
   num14 = 6.28318530717959 / CDbl(num9)
    index = 0
    While index < (num9 \ 2)
        numArray4(index) = Math.Sin(a)
        numArray5(index) = Math.Cos(a)
        a = a + num14
        index = index + 1
    Wend
    num7 = num9
    num3 = 1
    For num2 = 1 To num8
        num7 = num7 / 2
        num6 = 0
        For num11 = 1 To num3
            num10 = 0
            index = num6
            While index <= ((num7 + num6) - 1)
                num5 = index + num7
                a = numArray2(index) - numArray2(num5)
                num14 = numArray3(index) - numArray3(num5)
                numArray2(index) = numArray2(index) + numArray2(num5)
                numArray3(index) = numArray3(index) + numArray3(num5)
                If num10 = 0 Then
                    numArray2(num5) = a
                    numArray3(num5) = num14
                Else
                    numArray2(num5) = (a * numArray5(num10)) + (num14 * numArray4(num10))
                    numArray3(num5) = (num14 * numArray5(num10)) - (a * numArray4(num10))
                End If
                num10 = num10 + num3
                index = index + 1
           Wend
            num6 = (num6 + num7) + num7
        Next
        num3 = num3 + num3
    Next
    num5 = num9 \ 2
    For index = 1 To (num9 - 1)
        num6 = num9
        If num5 < index Then
            Dim num12 As Double
              
          num12 = numArray2(index)
            numArray2(index) = numArray2(num5)
            numArray2(num5) = num12
            num12 = numArray3(index)
            numArray3(index) = numArray3(num5)
            numArray3(num5) = num12
        End If
        num6 = num6 / 2
        Do While num5 >= num6
            num5 = num5 - num6
            num6 = num6 / 2
            If num5 = 0 Then
                Exit Do
            End If
        Loop
        num5 = num5 + num6
    Next
    For index = 0 To num9 - 1
        numArray(index, 0) = numArray2(index)
        numArray(index, 1) = numArray3(index)
        numArray(index, 2) = ((numArray2(index)) ^ 2# + (numArray3(index)) ^ 2#) ^ 0.5
    Next
   fft = numArray
End Function

Private Sub Command1_Click()
Dim x
x = fft("123456789")
Text1 = x
End Sub

Private Sub Command2_Click()
Text1 = ""
End Sub

#18
ysr28572020-09-14 09:51
如下图片就是启动程序的结果:
只有本站会员才能查看附件,请 登录
#19
cwa99582020-09-14 16:23
fft()函数你定义的传递数据是数组,你把一个字符串传递过去,肯定出错的啊。

Function fft(ByRef Data() As Double) As Double()
返回的就是一个数,哪来的有虚部,实部数据呢?
到底是要算什么?
#20
ysr28572020-09-14 22:32
谢谢指导和帮助!
这个可调用程序是复制人家的,不知道这个数组如何调用?如何输入?
实部和虚部是人家的说明和提示。不知道,不明白程序的原理。
#21
风吹过b2020-09-15 21:35
稍微搜索了一下 大数乘法,
搜到了就是 五种算法:
1、模拟小学乘法:最简单的乘法竖式手算的累加型;
2、分治乘法:最简单的是Karatsuba乘法,一般化以后有Toom-Cook乘法;
3、快速傅里叶变换FFT:(为了避免精度问题,可以改用快速数论变换FNTT),时间复杂度O(N lgN lglgN)。具体可参照Schönhage–Strassen algorithm;
4、中国剩余定理:把每个数分解到一些互素的模上,然后每个同余方程对应乘起来就行;
5、Furer’s algorithm:在渐进意义上FNTT还快的算法。不过好像不太实用,本文就不作介绍了。大家可以参考维基百科Fürer’s algorithm

然后看了一介绍,Java 中,采用的是
用二重循环直接相乘
采用 Karatsuba algorithm
采用 Toom-Cook multiplication

然后看了一下,竟然是数论。好吧,我根本看不懂。
按java的算法类型,建议你照着 Karatsuba 乘法 去写吧!数位不够时,可以使用递归就是了。

或者你去翻译:GMP(The GNU MP Bignum Library)。这里面有你所需要的算法,不过是 C++ 的。

附:找到的 Karatsuba 算法, C++ 的,稍微琢磨了一下,估计要先写一个大数加减法的。主要是看不太懂。
程序代码:
/**

 * Karatsuba乘法

 */
public static long karatsuba(long num1, long num2){
    //递归终止条件
    if(num1 < 10 || num2 < 10) return num1 * num2;

    // 计算拆分长度
    int size1 = String.valueOf(num1).length();
    int size2 = String.valueOf(num2).length();
    int halfN = Math.max(size1, size2) / 2;

    /* 拆分为a, b, c, d */
    long a = Long.valueOf(String.valueOf(num1).substring(0, size1 - halfN));
    long b = Long.valueOf(String.valueOf(num1).substring(size1 - halfN));
    long c = Long.valueOf(String.valueOf(num2).substring(0, size2 - halfN));
    long d = Long.valueOf(String.valueOf(num2).substring(size2 - halfN));

    // 计算z2, z0, z1, 此处的乘法使用递归
    long z2 = karatsuba(a, c);
    long z0 = karatsuba(b, d);
    long z1 = karatsuba((a + b), (c + d)) - z0 - z2;

    return (long)(z2 * Math.pow(10, (2*halfN)) + z1 * Math.pow(10, halfN) + z0);
}




--------------------
最后,我决定我再不关注这个贴子了,看不懂。

[此贴子已经被作者于2020-9-15 21:36编辑过]

#22
ysr28572020-09-15 23:06
谢谢关注和指导!这些东西很重要,我好好学习一下,感觉还是快速傅里叶变换或数论变换最快,都需要学习,我也是不懂,谢谢您!
#23
ysr28572020-10-09 17:47
ω是如何计算的呢?
为了计算离散傅里叶变换,我们令:
ω = e^(-2πi/N )= e^(-2πi/8) = e^(-πi/4) = cos(-π/4) + i sin(-π/4) = √2 / 2 - i √2 / 2 ≈ 0.7-0.7i
而ω^0~ω^8都以此为基础,以此类推,如ω^0=e^0=1,
ω^1=e^(-πi/4) = cos(-π/4) + i sin(-π/4) = √2 / 2 - i √2 / 2 ≈ 0.7-0.7i,
ω^2=e^(-πi/2)=0+i*(-1)=-i,
……。

[此贴子已经被作者于2020-10-9 23:21编辑过]

#24
ysr28572020-10-09 17:49
ω是如何计算的呢?
为了计算离散傅里叶变换,如当N=8时,我们令:
ω = e^(-2πi/N )= e^(-2πi/8) = e^(-πi/4) = cos(-π/4) + i sin(-π/4) = √2 / 2 - i √2 / 2 ≈ 0.7-0.7i
而ω^0~ω^8都以此为基础,以此类推,如ω^0=e^0=1,
ω^1=e^(-πi/4) = cos(-π/4) + i sin(-π/4) = √2 / 2 - i √2 / 2 ≈ 0.7-0.7i,
ω^2=e^(-πi/2)=0+i*(-1)=-i,
……。

[此贴子已经被作者于2020-10-9 23:22编辑过]

#25
ysr28572020-10-09 18:14
总计进行了两次快速傅里叶变换(被乘数和乘数各一次),得到:
{ A7, A6, A5, A4, A3, A2, A1, A0 } = { 12.9+10.9i, 2+7i, 3.1-1.1i, 7, 3.1+1.1i, 2-7i, 12.9-10.9i, 21 }
{ B7, B6, B5, B4, B3, B2, B1, B0 } = { 4.1+6.1i, -2+3i, -0.1-1.9i, 3, -0.1+1.9i, -2-3i, 4.1-6.1i, 9 }

一次对应项相乘,得到向量 {Ck},即 { C7, C6, C5, C4, C3, C2, C1, C0 }
= { -13.6+123.4i, -25-8i, -2.4-5.8i, 21, -2.4+5.8i, -25+8i, -13.6-123.4i, 189 }

一次逆变换计算出了向量 {ai} 和向量 {bj} 的卷积向量 {ck},如下所示:
{ c7, c6, c5, c4, c3, c2, c1, c0 } = { 0, 0, 0, 0, 24, 46, 65, 38, 16 }

最后错位相加得到结果。
#26
ysr28572020-10-23 18:50
例如如下向量:
{ A7, A6, A5, A4, A3, A2, A1, A0 } = { 12.9+10.9i, 2+7i, 3.1-1.1i, 7, 3.1+1.1i, 2-7i, 12.9-10.9i, 21 }
显然以7为对称中心,对称项都是共轭复数,只有7和21是不对称的,如何快速得到?
先得到7,还是先得到 3.1+1.1i?然后由对称性和其它算法(而不再进行乘法)推演出来其他项?
#27
ysr28572020-12-24 21:04
改造后的蝶形运算程序和结果:
Private Sub Command1_Click()
Dim xr() As Double, a As String
a = Trim(Text1)
ReDim xr(0 To Len(a) - 1)
For i1 = 0 To Len(a) - 1
xr(i1) = Mid(a, i1 + 1, 1)
  Next
Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double
Dim xi()
n = Len(a) '求数组大小,其值必须是2的幂
m = 0
l = 2
pi = 3.14159265358979
Do
 l = l + l
 m = m + 1
 Loop Until l > n
 n = l / 2
ReDim xi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  wr = 1
  wi = 0
  w1r = Cos(t)
  w1i = -Sin(t)
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * wr - xi(q) * wi
   ti = xr(q) * wi + xi(q) * wr
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   
   p = p + le
Loop Until p > n - 1

wr = wr * w1r - wi * w1i
wi = wr * w1i - wi * w1r
r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
   xr(i) = Sqr(xr(i) ^ 2 + xi(i) ^ 2)
   Text2 = Text2 & "  " & xr(i)
   Next
End Sub

Private Sub Command2_Click()
Text1 = ""
Text2 = ""
End Sub
运算结果:
Text1=00000678,Text2=21  7  9  5  21  7  9  5,
Text1=00000432,Text2=  9  3  1  5  9  3  1  5.
#28
ysr28572020-12-24 23:59
程序出现了错误,怪不得了呀,修改了一下,如下为代码和计算结果:
Private Sub Command1_Click()
Dim xr() As Double, a As String
a = Trim(Text1)
ReDim xr(0 To Len(a) - 1)
For i1 = 0 To Len(a) - 1
xr(i1) = Mid(a, i1 + 1, 1)
  Next
Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double
Dim xi()
n = Len(a) '求数组大小,其值必须是2的幂
m = 0
l = 2
pi = 3.14159265358979
Do
 l = l + l
 m = m + 1
 Loop Until l > n
 n = l / 2
ReDim xi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  wr = 1
  wi = 0
  t = pi / le1
  w1r = Cos(t)
  w1i = -Sin(t)
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * wr - xi(q) * wi
   ti = xr(q) * wi + xi(q) * wr
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   p = p + le
Loop Until p > n - 1

wr = wr * w1r - wi * w1i
wi = wr * w1i - wi * w1r
r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
   xr(i) = Sqr(xr(i) ^ 2 + xi(i) ^ 2)
   Text2 = Text2 & "  " & xr(i)
   Next
End Sub

Private Sub Command2_Click()
Text1 = ""
Text2 = ""
End Sub
计算结果:
Text1=00000678,Text2= 21  5.19615242270664  2.60761871483253  2.80919177949488  21  5.19615242270664  2.60761871483253  2.80919177949488.
#29
ysr28572020-12-25 00:04
改造后输出复数的蝶形运算程序代码及结果:
Private Sub Command1_Click()
Dim xr() As Double, a As String
a = Trim(Text1)
ReDim xr(0 To Len(a) - 1)
For i1 = 0 To Len(a) - 1
xr(i1) = Mid(a, i1 + 1, 1)
  Next
Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double
Dim xi()
n = Len(a) '求数组大小,其值必须是2的幂
m = 0
l = 2
pi = 3.14159265358979
Do
 l = l + l
 m = m + 1
 Loop Until l > n
 n = l / 2
ReDim xi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  wr = 1
  wi = 0
  t = pi / le1
  w1r = Cos(t)
  w1i = -Sin(t)
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * wr - xi(q) * wi
   ti = xr(q) * wi + xi(q) * wr
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   p = p + le
Loop Until p > n - 1

wr = wr * w1r - wi * w1i
wi = wr * w1i - wi * w1r
r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
   
   Text2 = Text2 & "  " & xr(i) & "+" & xi(i) & "i"
   Next
End Sub

Private Sub Command2_Click()
Text1 = ""
Text2 = ""
End Sub
计算结果:
Text1=00000678,Text2= 21+0i  -4.24264068711929+3i  -1.31801948466055+-2.25i  -1.68198051533947+2.25i  -21+0i  4.24264068711929+-3i  1.31801948466055+2.25i  1.68198051533947+-2.25i
#30
ysr28572020-12-25 00:12
跟前面的例子比较如下:
 21+0i  -4.24264068711929+3i  -1.31801948466055+-2.25i  -1.68198051533947+2.25i  -21+0i  4.24264068711929+-3i  1.31801948466055+2.25i  1.68198051533947+-2.25i

 12.9+10.9i, 2+7i, 3.1-1.1i, 7, 3.1+1.1i, 2-7i, 12.9-10.9i, 21

不一样,咋回事呢?希望老师指点!谢谢!
#31
ysr28572020-12-25 00:19
Text1=00000432,Text2=
9+0i  -2.82842712474619+2i  -0.146446609406727+-0.25i  -1.12132034355964+1.5i  -9+0i  2.82842712474619+-2i  0.146446609406727+0.25i  1.12132034355964+-1.5i
和前面的比较:
4.1+6.1i, -2+3i, -0.1-1.9i, 3, -0.1+1.9i, -2-3i, 4.1-6.1i, 9

不一样,显然有错误,咋回事呢?
#32
ysr28572020-12-27 00:57
改造后的代码及运行结果:
 Private Sub Command1_Click()
Dim xr() As Double, a As String
a = Trim(Text1)
ReDim xr(0 To Len(a) - 1)
For i1 = 0 To Len(a) - 1
xr(i1) = Mid(a, i1 + 1, 1)
  Next
Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double
Dim xi()
n = Len(a) '求数组大小,其值必须是2的幂
m = 0
l = 2
pi = 3.14159265358979
Do
 l = l + l
 m = m + 1
 Loop Until l > n
 n = l / 2
ReDim xi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  wr = 1
  wi = 0
  t = pi / le1
  w1r = -Cos(t)
  w1i = Sin(t)
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * wr - xi(q) * wi
   ti = xr(q) * wi + xi(q) * wr
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   p = p + le
Loop Until p > n - 1

wr = wr * w1r - wi * w1i
wi = wr * w1i + wi * w1r
r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
   xr(i) = Sqr(xr(i) ^ 2 + xi(i) ^ 2)
   Text2 = Text2 & "  " & xr(i) & "+" & xi(i) & "i"
   Next
End Sub

Private Sub Command2_Click()
Text1 = ""
Text2 = ""
End Sub

 输入:Text1 =00070608,结果:
     21+0i  5.19615242270664+3i  8.91503909477451+-1.9142135623731i  12.2115949757928+9.49264068711931i  
7+0i  5.19615242270661+-2.99999999999997i  5.62840405338267+1.9142135623731i  12.2115949757928+-9.49264068711934i
与正确值比较:12.9+10.9i, 2+7i, 3.1-1.1i, 7, 3.1+1.1i, 2-7i, 12.9-10.9i, 21也是不对。
输入:Text1 =00030402,结果:
     9+0i  3.46410161513776+2i  2.30993212835938+1.9142135623731i  8.1410633171952+6.3284271247462i  
3+0i  3.46410161513775+-1.99999999999999i  5.08144347718889+-1.9142135623731i  8.14106331719521+-6.32842712474621i
与正确值比较:4.1+6.1i, -2+3i, -0.1-1.9i, 3, -0.1+1.9i, -2-3i, 4.1-6.1i, 9也是不对。
#33
ysr28572021-01-03 19:42
Private Sub Command1_Click()
Dim xr() As Double, a As String
a = Trim(Text1)
b = Trim(Text3)
ReDim xr(0 To Len(a) - 1): ReDim yr(0 To Len(b) - 1): ReDim zr(0 To Len(b) - 1)
For i1 = 0 To Len(a) - 1
xr(i1) = Mid(a, i1 + 1, 1)
yr(i1) = Mid(b, i1 + 1, 1)

  Next
Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double, tr1 As Double
Dim xi(): Dim yi(): Dim zi()
n = Len(a) '求数组大小,其值必须是2的幂
m = 0
l = 2
pi = 3.14159265358979
Do
 l = l + l
 m = m + 1
 Loop Until l > n
 n = l / 2
ReDim xi(n - 1): ReDim yi(n - 1): ReDim zi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  wr = 1
  wi = 0
  t = pi / le1
  w1r = Cos(t)
  w1i = -Sin(t)
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * wr - xi(q) * wi
   ti = xr(q) * wi + xi(q) * wr
   
   tr1 = yr(q) * wr - yi(q) * wi
   ti1 = yr(q) * wi + yi(q) * wr
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   
    yr(q) = yr(p) - tr1
   yi(q) = yi(p) - ti1
   yr(p) = yr(p) + tr1
   yi(p) = yi(p) + ti1
   
   p = p + le
Loop Until p > n - 1

wr = wr * w1r - wi * w1i
wi = wr * w1i + wi * w1r
r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
   zr(i) = xr(i) * yr(i) - xi(i) * yi(i): zi(i) = xr(i) * yi(i) + xi(i) * yr(i)
   Text2 = Text2 & "  " & zr(i) & "+" & zi(i) & "i"
   Next
End Sub

Private Sub Command2_Click()
Text1 = ""
Text2 = ""
Text3 = ""
End Sub


 输入:Text1 =80607000,Text2=20403000,结果:
 189+0i  48.1200576850889+-33.8492424049177i  -9.10907287525383+0.845301243345486i  8.07394476689+-4.83093657874164i
  21+0i  -5.62005768508875+4.15075759508251i  -6.76592712474619+-6.46840725563929i  24.7688582392567+7.92894805108087i
这个与正确值不同,不知道逆变换能不能出来正确结果?
#34
ysr28572021-01-03 19:43
正确结果是:将她们逐项相乘得到向量 {Ck},即 { C7, C6, C5, C4, C3, C2, C1, C0 }
= { -13.6+123.4i, -25-8i, -2.4-5.8i, 21, -2.4+5.8i, -25+8i, -13.6-123.4i, 189 }
#35
ysr28572021-01-04 11:49
Private Sub Command1_Click()
Dim xr() As Double, a As String
a = Trim(Text1)
b = Trim(Text3)
ReDim xr(0 To Len(a) - 1): ReDim yr(0 To Len(b) - 1): ReDim zr(0 To Len(b) - 1)
For i1 = 0 To Len(a) - 1
xr(i1) = Mid(a, i1 + 1, 1)
yr(i1) = Mid(b, i1 + 1, 1)

  Next
Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double, tr1 As Double
Dim xi(): Dim yi(): Dim zi()
n = Len(a) '求数组大小,其值必须是2的幂
m = 0
l = 2
pi = 3.14159265358979
Do
l = l + l
m = m + 1
Loop Until l > n
n = l / 2
ReDim xi(n - 1): ReDim yi(n - 1): ReDim zi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  wr = 1
  wi = 0
  t = pi / le1
  w1r = Cos(t)
  w1i = -Sin(t)
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * wr - xi(q) * wi
   ti = xr(q) * wi + xi(q) * wr
   
   tr1 = yr(q) * wr - yi(q) * wi
   ti1 = yr(q) * wi + yi(q) * wr
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   
    yr(q) = yr(p) - tr1
   yi(q) = yi(p) - ti1
   yr(p) = yr(p) + tr1
   yi(p) = yi(p) + ti1
   
   p = p + le
Loop Until p > n - 1

wr = wr * w1r - wi * w1i
wi = wr * w1i + wi * w1r
r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
   zr(i) = xr(i) * yr(i) - xi(i) * yi(i): zi(i) = xr(i) * yi(i) + xi(i) * yr(i)
   If InStr(zr(i), ".") = 0 Then
   zr(i) = zr(i)
   Else
   a1 = Left(zr(i), InStr(zr(i), ".") - 1)
   b1 = Mid(zr(i), InStr(zr(i), "."), 3)
   zr(i) = a1 & b1
   End If
   If InStr(zi(i), ".") = 0 Then
   zi(i) = zi(i)
   Else
   a2 = Left(zi(i), InStr(zi(i), ".") - 1)
   b2 = Mid(zi(i), InStr(zi(i), "."), 3)
   zi(i) = a2 & b2
   End If
   s = zr(i) & "/" & s
   s1 = zi(i) & "/" & s1
   Next
   Text2 = s1
End Sub

Private Sub Command2_Click()
Text1 = ""
Text2 = ""
Text3 = ""
End Sub

输入:Text1 =80607000,Text2=20403000,结果:
实部: 24.76/-6.76/-5.62/21/8.07/-9.10/48.12/189/
虚部:  7.92/-6.46/4.15/0/-4.83/0.84/-33.84/0/
与正确值比较:-13.6+123.4i, -25-8i, -2.4-5.8i, 21, -2.4+5.8i, -25+8i, -13.6-123.4i, 189也是不对。
#36
ysr28572021-01-06 21:38
Private Sub Command1_Click()
Dim xr() As Double, a As String
a = Trim(Text1)
b = Trim(Text3)
ReDim xr(0 To Len(a) - 1): ReDim yr(0 To Len(b) - 1): ReDim zr(0 To Len(b) - 1)
For i1 = 0 To Len(a) - 1
xr(i1) = Mid(a, i1 + 1, 1)
yr(i1) = Mid(b, i1 + 1, 1)

  Next
Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double, tr1 As Double
Dim xi(): Dim yi(): Dim zi()
n = Len(a) '求数组大小,其值必须是2的幂
m = 0
l = 2
pi = 3.14159265358979
Do
l = l + l
m = m + 1
Loop Until l > n
n = l / 2
ReDim xi(n - 1): ReDim yi(n - 1): ReDim zi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  wr = 1
  wi = 0
  t = pi / le1
  w1r = Cos(t)
  w1i = -Sin(t)
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * wr - xi(q) * wi
   ti = xr(q) * wi + xi(q) * wr
   
   tr1 = yr(q) * wr - yi(q) * wi
   ti1 = yr(q) * wi + yi(q) * wr
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   
    yr(q) = yr(p) - tr1
   yi(q) = yi(p) - ti1
   yr(p) = yr(p) + tr1
   yi(p) = yi(p) + ti1
   
   p = p + le
Loop Until p > n - 1

wr = wr * w1r - wi * w1i
wi = wr * w1i + wi * w1r
r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
   zr(i) = xr(i) * yr(i) - xi(i) * yi(i): zi(i) = xr(i) * yi(i) + xi(i) * yr(i)
   If InStr(zr(i), ".") = 0 Then
   zr(i) = zr(i)
   Else
   a1 = Left(zr(i), InStr(zr(i), ".") - 1)
   b1 = Mid(zr(i), InStr(zr(i), "."), 3)
   zr(i) = a1 & b1
   End If
   If InStr(zi(i), ".") = 0 Then
   zi(i) = zi(i)
   Else
   a2 = Left(zi(i), InStr(zi(i), ".") - 1)
   b2 = Mid(zi(i), InStr(zi(i), "."), 3)
   zi(i) = a2 & b2
   End If
   s = "/" & zr(i) & s
   s1 = "/" & zi(i) & s1
   Next
   Dim c(), d()
s2 = Split(s, "/")
s3 = Split(s1, "/")
   j = UBound(s2)
  For k = 1 To j
      n1 = n1 + 1
       ReDim Preserve c(0 To n1 - 1)
       ReDim Preserve d(0 To n1 - 1)
      c(n1 - 1) = s2(n1): d(n1 - 1) = s3(n1)
    Next
   Text2 = s
End Sub

Private Sub Command2_Click()
Text1 = ""
Text2 = ""
Text3 = ""
End Sub


输入:Text1 =80607000,Text2=20403000,结果:
 /24.76/-6.76/-5.62/21/8.07/-9.10/48.12/189
  /7.92/-6.46/4.15/0/-4.83/0.84/-33.84/0
与正确值比较:-13.6+123.4i, -25-8i, -2.4-5.8i, 21, -2.4+5.8i, -25+8i, -13.6-123.4i, 189也是不对。
#37
ysr28572021-01-06 23:02
Private Sub Command1_Click()
Dim xr() As Double, a As String
a = Trim(Text1)
b = Trim(Text3)
ReDim xr(0 To Len(a) - 1): ReDim yr(0 To Len(b) - 1): ReDim zr(0 To Len(b) - 1)
For i1 = 0 To Len(a) - 1
xr(i1) = Mid(a, i1 + 1, 1)
yr(i1) = Mid(b, i1 + 1, 1)

  Next
Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double, tr1 As Double
Dim xi(): Dim yi(): Dim zi()
n = Len(a) '求数组大小,其值必须是2的幂
m = 0
l = 2
pi = 3.14159265358979
Do
l = l + l
m = m + 1
Loop Until l > n
n = l / 2
ReDim xi(n - 1): ReDim yi(n - 1): ReDim zi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  wr = 1
  wi = 0
  t = pi / le1
  w1r = Cos(t)
  w1i = -Sin(t)
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * wr - xi(q) * wi
   ti = xr(q) * wi + xi(q) * wr
   
   tr1 = yr(q) * wr - yi(q) * wi
   ti1 = yr(q) * wi + yi(q) * wr
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   
    yr(q) = yr(p) - tr1
   yi(q) = yi(p) - ti1
   yr(p) = yr(p) + tr1
   yi(p) = yi(p) + ti1
   
   p = p + le
Loop Until p > n - 1

wr = wr * w1r - wi * w1i
wi = wr * w1i + wi * w1r
r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
   zr(i) = xr(i) * yr(i) - xi(i) * yi(i): zi(i) = xr(i) * yi(i) + xi(i) * yr(i)
   If InStr(zr(i), ".") = 0 Then
   zr(i) = zr(i)
   Else
   a1 = Left(zr(i), InStr(zr(i), ".") - 1)
   b1 = Mid(zr(i), InStr(zr(i), "."), 3)
   zr(i) = a1 & b1
   End If
   If InStr(zi(i), ".") = 0 Then
   zi(i) = zi(i)
   Else
   a2 = Left(zi(i), InStr(zi(i), ".") - 1)
   b2 = Mid(zi(i), InStr(zi(i), "."), 3)
   zi(i) = a2 & b2
   End If
   s = "/" & zr(i) & s
   s1 = "/" & zi(i) & s1
   Next
  s2 = nifft(Trim(s), Trim(s1))
   Text2 = s2
End Sub

Private Sub Command2_Click()
Text1 = ""
Text2 = ""
Text3 = ""
End Sub

Private Function nifft(sa As String, sb As String) As String
Dim l As Long, le As Long, le1 As Long, j As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double, tr1 As Double
Dim xi(): Dim yi(): Dim zi()
Dim xr(), yr(), zr()
s2 = Split(sa, "/")
s3 = Split(sb, "/")
   j = UBound(s2)
   n = j
  For k = 1 To j
      n1 = n1 + 1
       ReDim Preserve xr(0 To n1 - 1)
       ReDim Preserve yr(0 To n1 - 1)
      xr(n1 - 1) = s2(n1): yr(n1 - 1) = s3(n1)
    Next
   

 ReDim zr(0 To j - 1)

m = 0
l = 2
pi = 3.14159265358979
Do
l = l + l
m = m + 1
Loop Until l > n
n = l / 2
ReDim xi(n - 1): ReDim yi(n - 1): ReDim zi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  wr = 1
  wi = 0
  t = pi / le1
  w1r = Cos(t)
  w1i = -Sin(t)
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * wr - xi(q) * wi
   ti = xr(q) * wi + xi(q) * wr
   
   tr1 = yr(q) * wr - yi(q) * wi
   ti1 = yr(q) * wi + yi(q) * wr
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   
    yr(q) = yr(p) - tr1
   yi(q) = yi(p) - ti1
   yr(p) = yr(p) + tr1
   yi(p) = yi(p) + ti1
   
   p = p + le
Loop Until p > n - 1

wr = wr * w1r - wi * w1i
wi = wr * w1i + wi * w1r
r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
zr(i) = xr(i) + yr(i)
   If InStr(zr(i), ".") = 0 Then
   zr(i) = zr(i)
   Else
   a1 = Left(zr(i), InStr(zr(i), ".") - 1)
   b1 = Mid(zr(i), InStr(zr(i), "."), 3)
   zr(i) = a1 & b1
   End If
   
   s1 = Int(Val(zr(i) + 0.5))
   s = "/" & s1 & s
   Next
   
   
nifft = s

End Function

输入:Text1 =80607000,Text2=20403000,结果:
/48/30/38/-159/43/-31/54/237
与正确值比较:0, 0, 0, 24, 46, 65, 38, 16也是不对。
#38
ysr28572021-02-08 23:00
程序结果不对,可能是蝶形运算程序不对,请参考下面图片中的内容:(都是vc程序咱不懂)
只有本站会员才能查看附件,请 登录
只有本站会员才能查看附件,请 登录
只有本站会员才能查看附件,请 登录
#39
ysr28572021-02-26 09:48
这个结果对,但不节约时间:text1=80607000时,结果为:  21+0i  12.9497474683058+-10.9497474683058i  2.00000000000001+-7.00000000000002i  3.05025253169415+1.05025253169416i  7+1.61554255216634E-14i  3.0502525316942+-1.05025253169419i  1.99999999999997+6.99999999999994i  12.9497474683057+10.9497474683059i
代码如下:

Private Sub Command1_Click() '蝶形运算程序
Dim xr() As Double, a As String
 a = Trim(Text1)
 ReDim xr(0 To Len(a) - 1)
 For i1 = 0 To Len(a) - 1
 xr(i1) = Mid(a, i1 + 1, 1)
   Next
 Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
 Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
 Dim pi As Double, t As Double
 Dim xi()
 Dim yr(), yi()

 n = Len(a) '求数组大小,其值必须是2的幂
m = 0
 l = 2
 pi = 3.14159265358979
 Do
 l = l + l
 m = m + 1
 Loop Until l > n
 n = l / 2
 ReDim xi(n - 1)
 ReDim yr(n - 1): ReDim yi(n - 1)

 l = 1
 For i1 = 0 To Len(a) - 1
  tr1 = xr(0) * Cos(0): ti1 = xr(0) * Sin(0)
  tr2 = xr(4) * Cos((-2 * pi / 8) * i1): ti2 = xr(4) * Sin((-2 * pi / 8) * i1)
  tr3 = xr(2) * Cos((-2 * pi / 8) * i1 * 2): ti3 = xr(2) * Sin((-2 * pi / 8) * i1 * 2)
  yr(i1) = tr1 + tr2 + tr3: yi(i1) = ti1 + ti2 + ti3
  Next
  
 Do
   le = 2 ^ l
   le1 = le / 2
   wr = 1
   wi = 0
   t = pi / le1
   w1r = Cos(t)
   w1i = -Sin(t)
   Print l
   r = 0
 Do
   p = r
   Do
    q = p + le1
   
    tr = xr(q) * wr - xi(q) * wi
    ti = xr(q) * wi + xi(q) * wr
   
    xr(q) = xr(p) - tr
    xi(q) = xi(p) - ti
    xr(p) = xr(p) + tr
    xi(p) = xi(p) + ti
    Print p, q
   
    Print xr(p); xi(p); le1, xr(q); xi(q); le1
   
   
    p = p + le
 Loop Until p > n - 2



 wr = wr * w1r - wi * w1i
 wi = wr * w1i + wi * w1r
 r = r + 1
 Loop Until r > le1 - 1
 l = l + 1
 Loop Until l > m

 For i = 0 To n - 1 '仅输出模
Print xr(i), xi(i)
    Text2 = Text2 & "  " & yr(i) & "+" & yi(i) & "i"
    Next

 End Sub

 Private Sub Command2_Click()
 Text1 = ""
 Text2 = ""
 End Sub
#40
ysr28572021-02-28 08:52
程序弄不对,改进了一下,如下是结果和代码以及原理图片:
Text1=80607000,结果为: 8+16i  -1.29610059419054+-0.372221061679249i  3.05025253169417+8.70710678118656i  1.532843272421+4.59431073134172i  8+0i  1+-2.21998012670182i  8+-2.60660171779821i  8+-1.52862418649973i.
代码如下:

Private Sub Command1_Click() '蝶形运算程序
Dim xr() As Double, a As String
 a = Trim(Text1)
 ReDim xr(0 To Len(a) - 1)
 For i1 = 0 To Len(a) - 1
 xr(i1) = Mid(a, i1 + 1, 1)
   Next
 Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
 Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
 Dim pi As Double, t As Double
 Dim xi()
 n = Len(a) '求数组大小,其值必须是2的幂
m = 0
 l = 2
 pi = 3.14159265358979
 Do
 l = l + l
 m = m + 1
 Loop Until l > n
 n = l / 2
 ReDim xi(n - 1)

 l = 1
 Do
   le = 2 ^ l
   le1 = le / 2
   wr = 1
   wi = 0
   t = pi / le1
   w1r = Cos(t)
   w1i = -Sin(t)
   Print l
   r = 0
 Do
   p = r
   Do
    q = p + le1
   
    tr = xr(q) * wr - xi(q) * wi
    ti = xr(q) * wi + xi(q) * wr
    m1 = (p + Len(a) / 2) Mod Len(a)
    xr(p) = xr(p) + xr(m1) * Sin((-2 * pi / Len(a)) * r * 2 ^ (m - le1))
  xr(q) = xr(p) - xr(m1) * Sin((-2 * pi / Len(a)) * r * 2 ^ (m - le1))
    xi(p) = xr(p) + xr(m1) * Cos((-2 * pi / Len(a)) * r * 2 ^ (m - le1))
    xi(q) = xr(p) - xr(m1) * Cos((-2 * pi / Len(a)) * r * 2 ^ (m - le1))
   
    Print p, q
   
    Print xr(p); xi(p); le1, xr(q); xi(q); le1
   
   
    p = p + le
 Loop Until p > n - 2

 wr = wr * w1r - wi * w1i
 wi = wr * w1i + wi * w1r
 r = r + 1
 Loop Until r > le1 - 1
 l = l + 1
 Loop Until l > m

 For i = 0 To n - 1 '仅输出模
Print xr(i), xi(i)
    Text2 = Text2 & "  " & xr(i) & "+" & xi(i) & "i"
    Next

 End Sub

 Private Sub Command2_Click()
 Text1 = ""
 Text2 = ""
 End Sub

原理图片:
只有本站会员才能查看附件,请 登录
只有本站会员才能查看附件,请 登录
只有本站会员才能查看附件,请 登录
#41
ysr28572021-02-28 09:22
图片中的-j应该是-i吧?网上的东西胡写八道,怎么就弄对了?
#42
ysr28572021-03-02 00:22
改了一下程序,仍然不对,下面是结果和代码:
 输入:Text1=80607000,结果: 21+0i  14.467156727579+-8.67878402655563i  6.94974746830584+-4.94974746830583i  10.6787840265556+-0.467156727579004i  7+0i  1.532843272421+-3.32121597344437i  -2.94974746830584+4.94974746830583i  5.32121597344435+12.467156727579i.
代码如下:

Private Sub Command1_Click() '蝶形运算程序
Dim xr() As Double, a As String
 a = Trim(Text1)
 ReDim xr(0 To Len(a) - 1)
 For i1 = 0 To Len(a) - 1
 xr(i1) = Mid(a, i1 + 1, 1)
   Next
 Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
 Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
 Dim pi As Double, t As Double
 Dim xi()
 n = Len(a) '求数组大小,其值必须是2的幂
m = 0
 l = 2
 pi = 3.14159265358979
 Do
 l = l + l
 m = m + 1
 Loop Until l > n
 n = l / 2
 ReDim xi(n - 1)

 l = 1
 Do
   le = 2 ^ l
   le1 = le / 2
   wr = 1
   wi = 0
   t = pi / le1
   w1r = Cos(t)
   w1i = -Sin(t)
   Print l
   r = 0
 Do
   p = r
   Do
    q = p + le1
   
    tr = xr(q) * Cos((-2 * pi / 2 ^ le1) * r)
    ti = xr(q) * Sin((-2 * pi / 2 ^ le1) * r)
   
    xr(q) = xr(p) - tr
    xi(q) = xi(p) - ti
    xr(p) = xr(p) + tr
    xi(p) = xi(p) + ti
    Print p, q
   
    Print xr(p); xi(p); r, xr(q); xi(q); r
   
   
    p = p + le
 Loop Until p > n - 2

 wr = wr * w1r - wi * w1i
 wi = wr * w1i + wi * w1r
 r = r + 1
 Loop Until r > le1 - 1
 l = l + 1
 Loop Until l > m

 For i = 0 To n - 1 '仅输出模
Print xr(i), xi(i)
    Text2 = Text2 & "  " & xr(i) & "+" & xi(i) & "i"
    Next

 End Sub

 Private Sub Command2_Click()
 Text1 = ""
 Text2 = ""
 End Sub
#43
ysr28572021-03-02 00:38
又改了一下程序,仍然不对,下面是结果和代码:
 输入:Text1=80607000,结果:  21+0i  8.86549696282261+-1.36563225411292i  8.46715672757901+-2.67878402655563i  19.8202872861178+-3.88899163113719i  7+0i  -4.86549696282261+1.36563225411288i  -4.46715672757901+2.67878402655563i  8.17971271388218+3.88899163113723i代码如下:

Private Sub Command1_Click() '蝶形运算程序
Dim xr() As Double, a As String
a = Trim(Text1)
ReDim xr(0 To Len(a) - 1)
For i1 = 0 To Len(a) - 1
xr(i1) = Mid(a, i1 + 1, 1)
  Next
Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double
Dim xi()
n = Len(a) '求数组大小,其值必须是2的幂
m = 0
l = 2
pi = 3.14159265358979
Do
l = l + l
m = m + 1
Loop Until l > n
n = l / 2
ReDim xi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  wr = 1
  wi = 0
  t = pi / le1
  w1r = Cos(t)
  w1i = -Sin(t)
  Print l
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * Cos((-2 * pi / 2 ^ le1) * r * 2 ^ (m - le1))
   ti = xr(q) * Sin((-2 * pi / 2 ^ le1) * r * 2 ^ (m - le1))
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   Print p, q
   
   Print xr(p); xi(p); r, xr(q); xi(q); r
   
   
   p = p + le
Loop Until p > n - 2


r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
Print xr(i), xi(i)
   Text2 = Text2 & "  " & xr(i) & "+" & xi(i) & "i"
   Next

End Sub

Private Sub Command2_Click()
Text1 = ""
Text2 = ""
End Sub

#44
ysr28572021-03-02 00:54
哈哈哈!这回可能蒙对了,太激动了!!
 如下是结果和程序代码:
 输入:Text1=80607000,结果: 21+0i  12.9497474683058+-10.9497474683058i  2.00000000000001+-7i  3.05025253169417+1.05025253169416i  7+0i  3.05025253169417+-1.05025253169417i  1.99999999999999+7i  12.9497474683058+10.9497474683058i.

代码如下:

Private Sub Command1_Click() '蝶形运算程序
Dim xr() As Double, a As String
a = Trim(Text1)
ReDim xr(0 To Len(a) - 1)
For i1 = 0 To Len(a) - 1
xr(i1) = Mid(a, i1 + 1, 1)
  Next
Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double
Dim xi()
n = Len(a) '求数组大小,其值必须是2的幂
m = 0
l = 2
pi = 3.14159265358979
Do
l = l + l
m = m + 1
Loop Until l > n
n = l / 2
ReDim xi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  Print l
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * Cos((-2 * pi / 2 ^ l) * r)
   ti = xr(q) * Sin((-2 * pi / 2 ^ l) * r)
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   Print p, q
   
   Print xr(p); xi(p); r, xr(q); xi(q); r
   
   
   p = p + le
Loop Until p > n - 2


r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
Print xr(i), xi(i)
   Text2 = Text2 & "  " & xr(i) & "+" & xi(i) & "i"
   Next

End Sub

Private Sub Command2_Click()
Text1 = ""
Text2 = ""
End Sub


明天再说!!

[此贴子已经被作者于2021-3-2 01:11编辑过]

#45
ysr28572021-03-02 06:45
倒序程序修改了一下就正确了:
 输入:Text1=00000678,结果:80607000.
输入:Text1=00000432,结果:20403000.

代码如下:
Private Sub Command1_Click()
  Dim x_() As Double, a As String
  a = Trim(Text1)
  ReDim x_(1 To Len(a))
  For i1 = 1 To Len(a)
  x_(i1) = Mid(a, Len(a) - i1 + 1, 1)
    Next
  Dim n As Integer, i As Long, j As Long, mn As Long, lh As Long, t As Double, k As Long
  '位序倒置
n = Len(a) '求数组大小,其值必须是2的幂
lh = n / 2
  j = n / 2
  For i = 1 To n - 2


  Debug.Print i, j
  k = lh '下面是向右进位算法
Do
  If k > j Then Exit Do '高位是1吗
j = j - k '是的,高位置0
  k = k / 2 '准备次高位的权
Loop Until k = 0 '次高位的权若非0,则检查新的次高位
j = j + k '非则若最高位是0,则置1
  Text2 = Text2 & x_(j + 1)
  Next
  Text2 = x_(1) & x_(1 + Len(a) / 2) & Text2
  
  End Sub

  Private Sub Command2_Click()
  Text1 = ""
  Text2 = ""
  End Sub
#46
ysr28572021-03-02 22:17
又不对了困难了,啊啊!
输入:Text1 =80607000,Text2=20403000,结果:
/29/65/34/24/10/2/11/16
与正确值比较:0, 0, 0, 0, 24, 46, 65, 38, 16也是不对。

Private Sub Command1_Click()
Dim xr() As Double, a As String
a = Trim(Text1)
b = Trim(Text3)
sb1 = Len(a) + Len(b)
sb2 = Log(sb1) / Log(2)
If InStr(sb2, ".") = 0 Then
sb2 = sb2
Else
sb2 = Int(sb2) + 1
End If
sb = 2 ^ sb2
Print sb
ReDim xr(0 To Len(a) - 1): ReDim yr(0 To Len(b) - 1): ReDim zr(0 To Len(b) - 1)
For i1 = 0 To Len(a) - 1
xr(i1) = Mid(a, i1 + 1, 1)
yr(i1) = Mid(b, i1 + 1, 1)

  Next
Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double, tr1 As Double
Dim xi(): Dim yi(): Dim zi()
n = Len(a) '求数组大小,其值必须是2的幂
m = 0
l = 2
pi = 3.14159265358979
Do
l = l + l
m = m + 1
Loop Until l > n
n = l / 2
ReDim xi(n - 1): ReDim yi(n - 1): ReDim zi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
 
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * Cos((-2 * pi / 2 ^ l) * r)
   ti = xr(q) * Sin((-2 * pi / 2 ^ l) * r)
   
   tr1 = yr(q) * Cos((-2 * pi / 2 ^ l) * r)
   ti1 = yr(q) * Sin((-2 * pi / 2 ^ l) * r)
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   
    yr(q) = yr(p) - tr1
   yi(q) = yi(p) - ti1
   yr(p) = yr(p) + tr1
   yi(p) = yi(p) + ti1
   
   p = p + le
Loop Until p > n - 1


r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
   zr(i) = xr(i) * yr(i) - xi(i) * yi(i): zi(i) = xr(i) * yi(i) + xi(i) * yr(i)
   If InStr(zr(i), ".") = 0 Then
   zr(i) = zr(i)
   Else
   a1 = Left(zr(i), InStr(zr(i), ".") - 1)
   b1 = Mid(zr(i), InStr(zr(i), "."), 3)
   zr(i) = a1 & b1
   End If
   If InStr(zi(i), ".") = 0 Then
   zi(i) = zi(i)
   Else
   a2 = Left(zi(i), InStr(zi(i), ".") - 1)
   b2 = Mid(zi(i), InStr(zi(i), "."), 3)
   zi(i) = a2 & b2
   End If
   s = s & "/" & zr(i)
   s1 = s1 & "/" & zi(i)
   Next
  s2 = nifft(dxcx1(Trim(s)), dxcx1(Trim(s1)))
   Text2 = s2
End Sub

Private Sub Command2_Click()
Text1 = ""
Text2 = ""
Text3 = ""
Form1.Cls
End Sub

Private Function dxcx0(sa As String, sb As String) As String

Dim x_() As Double, a As String
 a = Trim(sa)
 ReDim x_(1 To sb)
 For i1 = 1 To sb
 x_(i1) = Mid(a, sb - i1 + 1, 1)
   Next
 Dim n As Integer, i As Long, j As Long, mn As Long, lh As Long, t As Double, k As Long
 '位序倒置
n = sb '求数组大小,其值必须是2的幂
lh = n / 2
 j = n / 2
 For i = 1 To n - 2


 Debug.Print i, j
 k = lh '下面是向右进位算法
Do
 If k > j Then Exit Do '高位是1吗
j = j - k '是的,高位置0
 k = k / 2 '准备次高位的权
Loop Until k = 0 '次高位的权若非0,则检查新的次高位
j = j + k '非则若最高位是0,则置1
 s = s & x_(j + 1)
 Next
 dxcx0 = x_(1) & x_(1 + sb / 2) & s
 

End Function

Private Function nifft(sa As String, sb As String) As String
Dim l As Long, le As Long, le1 As Long, j As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double, tr1 As Double
Dim xi(): Dim yi(): Dim zi()
Dim xr(), yr(), zr()
s2 = Split(sa, "/")
s3 = Split(sb, "/")
   j = UBound(s2)
   n = j
  For k = 1 To j
      n1 = n1 + 1
       ReDim Preserve xr(0 To n1 - 1)
       ReDim Preserve yr(0 To n1 - 1)
      xr(n1 - 1) = s2(n1): yr(n1 - 1) = s3(n1)
    Next
   

ReDim zr(0 To n - 1)

m = 0
l = 2
pi = 3.14159265358979
Do
l = l + l
m = m + 1
Loop Until l > n
n = l / 2
ReDim xi(n - 1): ReDim yi(n - 1): ReDim zi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * Cos((-2 * pi / 2 ^ l) * r)
   ti = xr(q) * Sin((-2 * pi / 2 ^ l) * r)
   
   tr1 = yr(q) * Cos((-2 * pi / 2 ^ l) * r)
   ti1 = yr(q) * Sin((-2 * pi / 2 ^ l) * r)
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   
    yr(q) = yr(p) - tr1
   yi(q) = yi(p) - ti1
   yr(p) = yr(p) + tr1
   yi(p) = yi(p) + ti1
   
   p = p + le
Loop Until p > n - 1


r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
zr(i) = (xr(i) - yi(i)) / n
   If InStr(zr(i), ".") = 0 Then
   zr(i) = zr(i)
   Else
   a1 = Left(zr(i), InStr(zr(i), ".") - 1)
   b1 = Mid(zr(i), InStr(zr(i), "."), 3)
   zr(i) = a1 & b1
   End If
   
   s1 = Int(Val(zr(i) + 0.5))
   s = "/" & s1 & s
   Next
   
   
nifft = s

End Function


Private Function dxcx1(sa As String) As String

Dim x_() As Double, a As String
 a = Trim(sa)
   

s2 = Split(sa, "/")

   j = UBound(s2)
   sb = j
   
    ReDim x_(1 To sb)
  For k = 1 To j
      n1 = n1 + 1
       ReDim Preserve x_(1 To n1)
      
      x_(n1) = s2(n1)
    Next
 Dim n As Integer, i As Long, mn As Long, lh As Long, t As Double
 '位序倒置
n = sb '求数组大小,其值必须是2的幂
lh = n / 2
 j = n / 2
 For i = 1 To n - 2


 Debug.Print i, j
 k = lh '下面是向右进位算法
Do
 If k > j Then Exit Do '高位是1吗
j = j - k '是的,高位置0
 k = k / 2 '准备次高位的权
Loop Until k = 0 '次高位的权若非0,则检查新的次高位
j = j + k '非则若最高位是0,则置1
 s = s & "/" & x_(j + 1)
 Next
 dxcx1 = "/" & x_(1) & "/" & x_(1 + sb / 2) & s
 
 End Function

#47
ysr28572021-03-02 22:22
0, 0, 0, 24, 46, 65, 38, 16依次错位相加得到292896,678*432=292896.
所以,这样才是正确的。

[此贴子已经被作者于2021-3-2 22:23编辑过]

#48
ysr28572021-03-03 10:39
改了一下(改了逆变换的符号)还是不对了困难了,啊啊!
输入:Text1 =80607000,Text2=20403000,结果:
/11/2/10/24/34/65/29/16
与正确值比较:0, 0, 0, 0, 24, 46, 65, 38, 16也是不对。

Private Sub Command1_Click()
Dim xr() As Double, a As String
a = Trim(Text1)
b = Trim(Text3)
sb1 = Len(a) + Len(b)
sb2 = Log(sb1) / Log(2)
If InStr(sb2, ".") = 0 Then
sb2 = sb2
Else
sb2 = Int(sb2) + 1
End If
sb = 2 ^ sb2
Print sb
ReDim xr(0 To Len(a) - 1): ReDim yr(0 To Len(b) - 1): ReDim zr(0 To Len(b) - 1)
For i1 = 0 To Len(a) - 1
xr(i1) = Mid(a, i1 + 1, 1)
yr(i1) = Mid(b, i1 + 1, 1)

  Next
Dim l As Long, le As Long, le1 As Long, n As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double, tr1 As Double
Dim xi(): Dim yi(): Dim zi()
n = Len(a) '求数组大小,其值必须是2的幂
m = 0
l = 2
pi = 3.14159265358979
Do
l = l + l
m = m + 1
Loop Until l > n
n = l / 2
ReDim xi(n - 1): ReDim yi(n - 1): ReDim zi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
 
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * Cos((-2 * pi / 2 ^ l) * r)
   ti = xr(q) * Sin((-2 * pi / 2 ^ l) * r)
   
   tr1 = yr(q) * Cos((-2 * pi / 2 ^ l) * r)
   ti1 = yr(q) * Sin((-2 * pi / 2 ^ l) * r)
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   
    yr(q) = yr(p) - tr1
   yi(q) = yi(p) - ti1
   yr(p) = yr(p) + tr1
   yi(p) = yi(p) + ti1
   
   p = p + le
Loop Until p > n - 1


r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
   zr(i) = xr(i) * yr(i) - xi(i) * yi(i): zi(i) = xr(i) * yi(i) + xi(i) * yr(i)
   If InStr(zr(i), ".") = 0 Then
   zr(i) = zr(i)
   Else
   a1 = Left(zr(i), InStr(zr(i), ".") - 1)
   b1 = Mid(zr(i), InStr(zr(i), "."), 3)
   zr(i) = a1 & b1
   End If
   If InStr(zi(i), ".") = 0 Then
   zi(i) = zi(i)
   Else
   a2 = Left(zi(i), InStr(zi(i), ".") - 1)
   b2 = Mid(zi(i), InStr(zi(i), "."), 3)
   zi(i) = a2 & b2
   End If
   s = s & "/" & zr(i)
   s1 = s1 & "/" & zi(i)
   Next
  s2 = nifft(dxcx1(Trim(s)), dxcx1(Trim(s1)))
   Text2 = s2
End Sub

Private Sub Command2_Click()
Text1 = ""
Text2 = ""
Text3 = ""
Form1.Cls
End Sub

Private Function dxcx0(sa As String, sb As String) As String

Dim x_() As Double, a As String
 a = Trim(sa)
 ReDim x_(1 To sb)
 For i1 = 1 To sb
 x_(i1) = Mid(a, sb - i1 + 1, 1)
   Next
 Dim n As Integer, i As Long, j As Long, mn As Long, lh As Long, t As Double, k As Long
 '位序倒置
n = sb '求数组大小,其值必须是2的幂
lh = n / 2
 j = n / 2
 For i = 1 To n - 2


 Debug.Print i, j
 k = lh '下面是向右进位算法
Do
 If k > j Then Exit Do '高位是1吗
j = j - k '是的,高位置0
 k = k / 2 '准备次高位的权
Loop Until k = 0 '次高位的权若非0,则检查新的次高位
j = j + k '非则若最高位是0,则置1
 s = s & x_(j + 1)
 Next
 dxcx0 = x_(1) & x_(1 + sb / 2) & s
 

End Function

Private Function nifft(sa As String, sb As String) As String
Dim l As Long, le As Long, le1 As Long, j As Long, r As Long, p As Long, q As Long, m As Byte
Dim wr As Double, w1 As Double, wlr As Double, wl1 As Double, tr As Double, t1 As Double
Dim pi As Double, t As Double, tr1 As Double
Dim xi(): Dim yi(): Dim zi()
Dim xr(), yr(), zr()
s2 = Split(sa, "/")
s3 = Split(sb, "/")
   j = UBound(s2)
   n = j
  For k = 1 To j
      n1 = n1 + 1
       ReDim Preserve xr(0 To n1 - 1)
       ReDim Preserve yr(0 To n1 - 1)
      xr(n1 - 1) = s2(n1): yr(n1 - 1) = s3(n1)
    Next
   

ReDim zr(0 To n - 1)

m = 0
l = 2
pi = 3.14159265358979
Do
l = l + l
m = m + 1
Loop Until l > n
n = l / 2
ReDim xi(n - 1): ReDim yi(n - 1): ReDim zi(n - 1)

l = 1
Do
  le = 2 ^ l
  le1 = le / 2
  
  r = 0
Do
  p = r
  Do
   q = p + le1
   
   tr = xr(q) * Cos((2 * pi / 2 ^ l) * r)
   ti = xr(q) * Sin((2 * pi / 2 ^ l) * r)
   
   tr1 = yr(q) * Cos((2 * pi / 2 ^ l) * r)
   ti1 = yr(q) * Sin((2 * pi / 2 ^ l) * r)
   
   xr(q) = xr(p) - tr
   xi(q) = xi(p) - ti
   xr(p) = xr(p) + tr
   xi(p) = xi(p) + ti
   
    yr(q) = yr(p) - tr1
   yi(q) = yi(p) - ti1
   yr(p) = yr(p) + tr1
   yi(p) = yi(p) + ti1
   
   p = p + le
Loop Until p > n - 1


r = r + 1
Loop Until r > le1 - 1
l = l + 1
Loop Until l > m

For i = 0 To n - 1 '仅输出模
zr(i) = (xr(i) - yi(i)) / n
   If InStr(zr(i), ".") = 0 Then
   zr(i) = zr(i)
   Else
   a1 = Left(zr(i), InStr(zr(i), ".") - 1)
   b1 = Mid(zr(i), InStr(zr(i), "."), 3)
   zr(i) = a1 & b1
   End If
   
   s1 = Int(Val(zr(i) + 0.5))
   s = "/" & s1 & s
   Next
   
   
nifft = s

End Function


Private Function dxcx1(sa As String) As String

Dim x_() As Double, a As String
 a = Trim(sa)
   

s2 = Split(sa, "/")

   j = UBound(s2)
   sb = j
   
    ReDim x_(1 To sb)
  For k = 1 To j
      n1 = n1 + 1
       ReDim Preserve x_(1 To n1)
      
      x_(n1) = s2(n1)
    Next
 Dim n As Integer, i As Long, mn As Long, lh As Long, t As Double
 '位序倒置
n = sb '求数组大小,其值必须是2的幂
lh = n / 2
 j = n / 2
 For i = 1 To n - 2


 Debug.Print i, j
 k = lh '下面是向右进位算法
Do
 If k > j Then Exit Do '高位是1吗
j = j - k '是的,高位置0
 k = k / 2 '准备次高位的权
Loop Until k = 0 '次高位的权若非0,则检查新的次高位
j = j + k '非则若最高位是0,则置1
 s = s & "/" & x_(j + 1)
 Next
 dxcx1 = "/" & x_(1) & "/" & x_(1 + sb / 2) & s
 
 End Function

1