|
#52
ysr28572021-03-14 14:06
VC版大数乘法程序（迭代型），不知道能不能运行？（我不懂VC程序，希望老师帮忙翻译成VB程序试试，谢谢！）
#include <stdio.h> #include <stdlib.h> #include <string.h> #include <stdbool.h> #include <math.h> #include <conio.h> #define N 150010 const 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; } |
Sub 倒序(X_() As Double)
Dim n As Integer, i As Long, j As Long, mn As Long, lh As Long, t As Double, k As Long
'位序倒置
n = UBound(X_) '求数组大小，其值必须是2的幂
lh = n / 2
j = n / 2
For i = 1 To n - 2
If i < j Then '倒序
t = X_(j)
X_(j) = X_(i)
X_(i) = t
End If
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
Next
End Sub
蝶形算法代码
Sub 蝶形算法(xr() As Double)
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 = UBound(xr) '求数组大小，其值必须是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)
Next
End Sub
检验的时候可以这样：
Sub 检验()
Dim y(63) As Double
For i = 0 To 64
y(i) = Sin(2 * 3.1415926 * i / 16)
Next
倒序 y()
蝶形算法 y()
'现在结果在y()中
End Sub
