用VBA解欧拉计划题目（146）

grf1973

今天发个求助吧，折腾了几天没搞出来

第146题： 素数模式研究 '
使得n^2+1，n^2+3，n^2+7，n^2+9，n^2+13以及n^2+27为连续素数的最小的n是10
所有小于一百万的数中，满足这一条件的所有整数n的和为1242490。
在小于一亿五千万的数中，所有满足这一条件的正整数n的和是多少？

现在通过各类筛选、剪枝，已经做到<10^6，10秒左右得出结果1242490
但15*10^7花的时间显然远不止10秒的150倍。。。。。。
我需要一个更高效的判断质数的办法。

grf1973

贴上我的代码和分析思路。计算到100万没有问题。但要算到15亿，估计得好几个小时。

1. Dim PrimeArr, k&, lmt&, p&, xp As Double
   
2. 'p146: 素数模式研究 '
   
3. '使得n^2+1，n^2+3，n^2+7，n^2+9，n^2+13以及n^2+27均为连续素数的最小的n是10.所有小于一百万的数中，满足这一条件的所有整数n的和为1242490。'
   
4. '在小于一亿五千万的数中，所有满足这一条件的正整数n的和是多少？
   
5. '分析：1、显然n不能是3，7，13的倍数，而且n是偶数'
   
6. '      2、考察n2=n*n,
   
7.            '(1)n2 mod 3=1（如果n2 mod 3=0，那么n2+9 mod 3=0;如果n2 mod 3=2，那么n2+1 mod 3=0)-->n mod 3=1，2（没意义，因为已经判断过n不是3的倍数，那么n2 mod 3==1）
   
8.            '(2)n2 mod 5=0（如果n2 mod 5=1，那么n2+9 mod 5=0;同理如果n2 mod 5=2，3，4都会有类似问题）。由于n是偶数，所以n是10的倍数。
   
9.            '(3)n2 mod 7=2或3（理由同上)-->n mod 7=4 或 3 -->n2 mod 7=2
   
10. '      4、题目需连续质数，所以要判断中间的n^2+5，n^2+11，n^2+15，n^2+17，n^2+19，n^2+21，n^2+23，n^2+25必须是合数
    
11.            '(1)因为n2尾数是0，显然n^2+5，n^2+15，n^2+25是合数
    
12.            '(2)因为n2 mod 3=1，n^2+11，n^2+17，n^2+23是合数
    
13.            '(2)因为n2 mod 7=2，n^2+19 是合数
    
14.            '(3)因此只要判断n^2+21是否合数即可
    
15. Sub problem146()
    
16.     tm = Timer
    
17.     Nmax& = 10 ^ 6
    
18. '    Nmax& = 15 * 10 ^ 7
    
19.     Dim n#, n2#, kk%
    
20.     arr = Array(1, 3, 7, 9, 13, 27)
    
21.     PrimeArr = GetPrimeArray(Nmax)
    
22.     ReDim brr(1 To Nmax) As Byte
    
23.     For Each a In Array(3, 7, 13)  '先筛掉一遍3，7，13的倍数
    
24.         For j = a To Nmax Step a
    
25.             brr(j) = 1
    
26.         Next
    
27.     Next
    
28.    
    
29.     For n = 10 To Nmax Step 10
    
30.         If brr(n) = 0 Then
    
31.             If n Mod 7 <> 3 And n Mod 7 <> 4 Then GoTo nextn
    
32.             n2 = n * n: lmt = n + 2
    
33.             For kk = 0 To 5
    
34.                 If Is_Prime(n2 + arr(kk)) = False Then GoTo nextn
    
35.             Next
    
36.             If Is_Prime(n2 + 21) = True Then GoTo nextn
    
37. '            If Is_Prime(n2 + 19) = True Then GoTo nextn
    
38.             res = res + n
    
39.             Debug.Print n, res
    
40.         End If
    
41. nextn:    Next
    
42.     Debug.Print res, "共运行" & Timer - tm & "秒"
    
43. End Sub
    
44. 
    
45. Function Is_Prime(x#) As Boolean   '判断x是不是质数，直接用小于sqr(x)的质因数试除
    
46. '    lmt = Int(Sqr(x))
    
47.     For k = 1 To UBound(PrimeArr)
    
48.         p = PrimeArr(k)
    
49.         If p > lmt Then Exit For
    
50.         xp = x / p
    
51.         If xp = Int(xp) Then Is_Prime = False: Exit Function
    
52.     Next
    
53.     Is_Prime = True
    
54. End Function
    
55. 
    
56. Function GetPrimeArray(n&) 'by kagawa
    
57.     Dim i&, j&, k&, m&
    
58.     m = Int((n - 1) / 2): ReDim a(1 To m) As Byte
    
59.     For i = 1 To Sqr(n) \ 2
    
60.         If a(i) = 0 Then
    
61.             For j = i * 3 + 1 To m Step i * 2 + 1
    
62.                 a(j) = 1
    
63.             Next
    
64.         End If
    
65.     Next
    
66.    
    
67.     ReDim b&(1 To m): b(1) = 2:  k = 1   '当a(i)=0时，2i+1是质数
    
68.     For i = 1 To m
    
69.         If a(i) = 0 Then k = k + 1: b(k) = i * 2 + 1
    
70.     Next
    
71.     ReDim Preserve b&(1 To k)
    
72.     GetPrimeArray = b
    
73. End Function

复制代码

grf1973

产生质数数组的代码来自香川。

香川群子

一个大数位的奇数的素性测试，可以采用拉宾-米勒测试。

参考：
http://blog.csdn.net/u013948191/article/details/47270979

香川群子

首先，当数值达到亿的平方数量级时，VBA中，传统的Mod、int算法都已经没有足够的精度来计算了。

所以，只能采用素数的数学验证方法。

初步了解如下：
有费马定律证明： a^(p-1) Mod p = 1 时，有较大概率p为素数。
因此，需要利用大数位算法，进行大数位的幂取模简化算法。

grf1973

香川群子 发表于 2017-11-14 20:55
首先，当数值达到亿的平方数量级时，VBA中，传统的Mod、int算法都已经没有足够的精度来计算了。

所以， ...

是啊，搞死我了。我用CDEC转换后还是溢出。。。当质数大于2.8*10^14时，用快速幂取模很有可能溢出。不过在此之前速度确实快了不少。
初步运算结果：
N              累加N              累计时间
10                10                    2.5625
315410        315420           3.1875
927070        1242490         4.453125
2525870       3768360        7.890625
8146100       11914460      20.20313
16755190      28669650      39.54688
当N= 16781930时，在用tstprime= 2测试时，需要计算 281633174524903 * 281633174524903  mod P，溢出。

grf1973

又是我家儿子告诉我，用2, 3, 7, 61, 24251这五个质数测试，准确率非常高。

grf1973

作为一次失败的尝试，贴上代码，以求来者。。。。。

1. Dim PrimeArr, tstPrimeArr, tstN%, tstprime
   
2. Sub problem146()
   
3.     tm = Timer
   
4.     nmax = 15 * 10 ^ 7
   
5.     Dim n, n2, kk%
   
6.     arr = Array(1, 3, 7, 9, 13, 27)
   
7.     tstPrimeArr = Array(2, 3, 7, 61, 89)
   
8.     ReDim brr(1 To nmax) As Byte
   
9.     For Each a In Array(3, 7, 13)  '先筛掉一遍3，7，13的倍数
   
10.         For j = a To nmax Step a
    
11.             brr(j) = 1
    
12.         Next
    
13.     Next
    
14.     On Error GoTo ext
    
15.     For n = 10 To nmax Step 10
    
16.         If brr(n) = 0 Then
    
17.             If n Mod 7 <> 3 And n Mod 7 <> 4 Then GoTo nextn
    
18.             n2 = CDec(n) * CDec(n)
    
19.             For kk = 0 To 5
    
20.                 If Is_Prime(CDec(n2) + CDec(arr(kk))) = False Then GoTo nextn
    
21.             Next
    
22.             If Is_Prime(CDec(n2) + CDec(21)) = True Then GoTo nextn
    
23.             res = res + n
    
24.             Debug.Print n, res, Timer - tm
    
25.         End If
    
26. nextn:    Next
    
27.     Debug.Print res, "共运行" & Timer - tm & "秒"
    
28.     Exit Sub
    
29. ext:  Debug.Print "error in :N="; n, "tstprime="; tstprime, CDec(n2) + CDec(arr(kk))
    
30. End Sub
    
31. Public Function Is_Prime(n) As Boolean   '判断N是不是质数
    
32. 'Miller-Rabbin 素数测试算法
    
33. 'tstPrimeArr = Array(2, 3, 7, 61, 24251)   '据说用这五个数足够了
    
34.     For tstN = 0 To UBound(tstPrimeArr)
    
35.         tstprime = tstPrimeArr(tstN)
    
36.         If Pow(tstprime, n - 1, n) <> 1 Then Exit Function
    
37.     Next
    
38.     Is_Prime = True
    
39. End Function
    
40. 
    
41. Function Pow(a, b, c)   '快速幂a^b mod c  当aa超过2.8*10^14，aa*aa会超过dec精度出错！！
    
42.     ans = 1
    
43.     aa = a: xx = b
    
44.     aa = Modular(aa, c)  '预处理，使得a处于c的数据范围之下
    
45.     Do While xx > 0  '二进制的移位操作，相当于每次除以2，用二进制看，就是我们不断的遍历b的二进制位
    
46.         tmp = Modular(xx, 2)   '从后往前的二进制位
    
47.         If tmp = 1 Then ans = Modular(CDec(ans) * CDec(aa), c)  '如果b的二进制位不是0，那么我们的结果是要参与运算的
    
48.         aa = Modular(CDec(aa) * CDec(aa), c)    ' 不断的加倍
    
49.         xx = Int(xx / 2)   '二进制左移一位
    
50.     Loop
    
51.     Pow = ans
    
52. End Function
    
53. 
    
54. Function Modular(a, b)   '返回a mod b
    
55.     Dim nb
    
56.     If a < b Then
    
57.         Modular = a
    
58.     ElseIf a = b Then
    
59.         Modular = 0
    
60.     Else
    
61.         Modular = a - Int(a / b) * b
    
62.     End If
    
63. End Function

复制代码

grf1973

第7行：tstPrimeArr = Array(2, 3, 7, 61, 89) 可以改成
tstPrimeArr = Array(2, 3, 7, 61, 24251)  但结果没变，报错位置改变了。

香川群子

grf1973 发表于 2017-11-15 09:54
又是我家儿子告诉我，用2, 3, 7, 61, 24251这五个质数测试，准确率非常高。

'9585921133193329 是一个合数，但是能通过绝大多数的费马素性测试。
2，3，7，61，24251 这几个都能通过的。

…………当然，即使这一个通过了，后面的应该通不过，所以不影响程序结果。

…………但是，如果真要做可靠的素性测试，还是需要用随机的拉宾-米勒算法进行测试：

伪代码如下：
Function IsPrime_3&(n, Optional k& = 8) 'n为要测试的大数(奇数)、k为随机测试次数 默认取k=8 已经能达到十万级精度
    Dim a, i&, j&, m, n2, r, r2&

    n2 = x2y(n, 1) '大数算法得到 n-1 （偶数）

    m = n2 'n-1
    Do
        m = x_2(m): r2 = r2 + 1 '大数整除2算法提取 2^r2
    Loop Until Right(m, 1) Mod 2
    '上述方法、分解使得 N-1=2^R*M

    Randomize
    For j = 1 To k '通常随机测试8次 错误概率t=(1/4)^k < 1.525*10^-5
        r = Rnd * 16777216 '2^24
        a = x4y(x3y(n, r), 16777216) 'a=n*rnd '取得随机数A<N
        If PowMod2(a, m, n) <> 1 Then IsPrime_3 = 1: Exit Function '检查 A^M%N=1时继续 <>1时判断为合数结束函数过程
        For i = 0 To r2 '上述A^M%N=1测试通过后，继续进行下面2^i循环测试
            If PowMod2(a, x3y(m, Pow_2(2, i)), n) <> n2 Then IsPrime_3 = 1: Exit Function
             '检查 A^((2^i)*M)%N=N-1 时继续、<>n-1(即n2)时判断为合数结束函数过程
        Next
    Next
    '全部随机测试都通过时，即可极大概率地判断为素数。
End Function


其中的计算，在VBA中都需要进行【大数位计算】的自定义函数编写。