串珠素数问题

无心人

总也调不好程序 明明算法对 但通不过编译

无心人

import Primes readInt :: IO Int readInt = readLn parta p a n = fromIntegral \$ sum (take ((div n 3) + 1) (iterate (1000*) (10^(p-1)*a))) partb p b n = fromIntegral \$ sum (take ((div n 3) + (p - 1)) (iterate (1000*) (100^(2-p)*b))) partc p c n = fromIntegral \$ sum (take (div n 3) (iterate (1000*) (10^p*c))) isczPrime n = [(ns, a, b, c) | ns <- [4..n], let p = mod ns 3, p /= 0, let (aa, bb, cc) = if (p == 1) then ([1,3,7,9], [0..9],[0..9]) else ([1..9], [1, 3, 7, 9], [0..9]), a <- aa, b <- bb, c <- cc, (a /= c) || (b /= c), let pp = (parta p a ns) + (partb p b ns) + (partc p c ns), isPrime pp ] main = do putStrLn "Input Digits: " n <- readInt print (isczPrime n) 100内的结果 p3_100.txt (37.68 KB, 下载次数: 1)

2008-11-21 20:15 上传

点击文件名下载附件

[ 本帖最后由 无心人 于 2008-11-21 20:15 编辑 ]

gxqcn

其中， abcabc...abca 型的有536个， abcabc...abcab 型的有463个， 总和为 999 个。 无心人的结果似乎有些问题， 居然含有 (4,1,1,2)，而 1121 = 19 * 59 总数目也少了些。

gxqcn

1. //////////////////////////////////////////////////////////////////////////
2. //
3. // 目的：搜索串珠素数（形如abcabc...a(ab)的素数）
4. // 设计：郭先强 ( gxqcn@163.com; HugeCalc@Gmail.com )
5. // 日期：2008-11-21
6. //
7. // Web: http://www.emath.ac.cn/
8. // BBS: http://bbs.emath.ac.cn/
9. //
10. //////////////////////////////////////////////////////////////////////////
12. // Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:
13. // Win32 Debug: Debug Multithreaded DLL
14. // Win32 Release: Multithreaded DLL
16. #include <stdio.h>
18. #include "../../../HugeCalc_API/CppAPI/Include/HugeCalc.h"
19. #include "../../../HugeCalc_API/CppAPI/Include/HugeInt.h" // 10进制系统
21. #pragma message( "automatic link to .../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
22. #pragma comment( lib, "../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
24. #define MAX_DIGITS 100
26. int main(int argc, char* argv[])
27. {
28. printf( "Call %s\n\n", HugeCalc::GetVer());
30. if ( HC_LICENSE_NONE == HugeCalc::GetLicenseLevel() )
31. {
32. printf( "\n警告：您未通过 HugeCalc.dll 的许可认证！" \
33. "\n解决方案可选下列方案之一：" \
34. "\n 一、请将本程序移动到“/CopyrightByGuoXianqiang/[../]”目录下运行；" \
35. "\n 二、或请在 HugeCalc.chm 中进行注册（一劳永逸）。" );
36. }
37. else
38. {
40. UINT32 n=0, a, b, c, s, count=0;
41. CHugeInt stepA(1001), stepB, stepC, p;
43. HugeCalc::EnableTimer();
44. HugeCalc::ResetTimer();
46. while(++n)
47. {
48. /* abcabc...abca 型 */
49. if ( 3*n + 1 > MAX_DIGITS )
50. {
51. break;
52. }
54. ( stepC = stepA ).DecRShift( 2 );
55. for ( a=1, p=stepA; a<10; a+=2, ++p+=stepA )
56. {
57. for ( b=0, s=(n+1)*a; b<10; ++b, s-=n*9 )
58. {
59. for ( c=0; c<10; ++c, s+=n )
60. {
61. if ( 0!=(s%3) && p.IsPrime() )
62. {
63. printf( "%s\tNo.%u\t(digits,a,b,c)=(%u,%u,%u,%u)\t%s\n",
64. HugeCalc::GetTimerStr(), ++count, 3*n+1, a, b, c, (LPCTSTR)p );
65. }
66. p += stepC;
67. }
68. }
69. }
72. /* abcabc...abcab 型 */
73. if ( 3*n + 2 > MAX_DIGITS )
74. {
75. break;
76. }
78. stepB = stepA;
79. stepC.DecLShift( 1 );
80. stepA.DecLShift( 1 );
82. for ( a=1, (p=stepA)+=stepB; a<10; ++a )
83. {
84. for ( b=1, s=(n+1)*(a+b); b<10; b+=2, s-=n*8-2, ++p+=stepB )
85. {
86. for ( c=0; c<10; ++c, s+=n )
87. {
88. if ( 0!=(s%3) && p.IsPrime() )
89. {
90. printf( "%s\tNo.%u\t(digits,a,b,c)=(%u,%u,%u,%u)\t%s\n",
91. HugeCalc::GetTimerStr(), ++count, 3*n+2, a, b, c, (LPCTSTR)p );
92. }
93. p += stepC;
94. }
95. }
96. }
99. ++stepA.DecLShift( 2 );
100. }
101. }
103. system( "pause" );
104. return 0;
105. }

复制代码

运行 2.38 秒即可输出全部结果。

无心人

可能有模式错误 刚学这个 好多东西要问别人 呵呵 我再修改去

无心人

找到问题了 b的模式有问题 ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ 在原帖处修改了

gxqcn

看不到你的新结果， 现在可一致了？ 这个用 HugeCalc 运行肯定要快很多。

无心人

我对比你结果对了 但我这里在优化的过程里 又产生了个问题 正找人问呢

无心人

我上面已经把程序改了 就是partb后面 起始数字是100^(2-p)*b 原来是10^(2-p)*b 马上把结果帖上去 优化后的问题还要等人回答呢 ========================================= 一直在修改的优化一点的程序 加上被3整除后出现了类型推定错误 经人指点，做了点修改，搞定了 import Primes readInt :: IO Int readInt = readLn parta p a n = fromIntegral \$ sum (take ((div n 3) + 1) (iterate (1000*) (10^(p-1)* toInteger a))) partb p b n = fromIntegral \$ sum (take ((div n 3) + (p - 1)) (iterate (1000*) (100^(2-p)* toInteger b))) partc p c n = fromIntegral \$ sum (take (div n 3) (iterate (1000*) (10^p* toInteger c))) isczPrime n = [(ns, a, b, c, pp) | ns <- [4..n], let p = ns \`mod\` 3, p /= 0, let (aa, bb, cc) = if (p == 1) then ([1,3,7,9], [0..9],[0..9]) else ([1..9], [1, 3, 7, 9], [0..9]), a <- aa, b <- bb, c <- cc, (a /= c) || (b /= c), let l = ns \`div\` 3, ((l + 1) * a + (l + p - 1) * b + l * c) \`mod\` 3 /= 0, let pp = (parta p a ns) + (partb p b ns) + (partc p c ns), isPrime pp ] main = do putStrLn "Input Digits: " n <- readInt print (isczPrime n) [[i] 本帖最后由 无心人 于 2008-11-21 21:42 编辑 [/i]]

无心人

下面讨论通用形 最多10个不同数字 最大100位数字 最小等于不同数字个数加一