串珠素数问题

无心人

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

无心人

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. //////////////////////////////////////////////////////////////////////////
    
11. 
    
12. // Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:
    
13. //      Win32 Debug:    Debug Multithreaded DLL
    
14. //      Win32 Release:  Multithreaded DLL
    
15. 
    
16. #include <stdio.h>
    
17. 
    
18. #include "../../../HugeCalc_API/CppAPI/Include/HugeCalc.h"
    
19. #include "../../../HugeCalc_API/CppAPI/Include/HugeInt.h"  // 10进制系统
    
20. 
    
21. #pragma message( "automatic link to .../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
    
22. #pragma comment( lib, "../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
    
23. 
    
24. #define MAX_DIGITS    100
    
25. 
    
26. int main(int argc, char* argv[])
    
27. {
    
28.     printf( "Call %s\n\n", HugeCalc::GetVer());
    
29. 
    
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.     {
    
39. 
    
40.         UINT32 n=0, a, b, c, s, count=0;
    
41.         CHugeInt stepA(1001), stepB, stepC, p;
    
42. 
    
43.         HugeCalc::EnableTimer();
    
44.         HugeCalc::ResetTimer();
    
45. 
    
46.         while(++n)
    
47.         {
    
48.             /* abcabc...abca 型  */
    
49.             if ( 3*n + 1 > MAX_DIGITS )
    
50.             {
    
51.                 break;
    
52.             }
    
53. 
    
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.             }
    
70. 
    
71. 
    
72.             /* abcabc...abcab 型  */
    
73.             if ( 3*n + 2 > MAX_DIGITS )
    
74.             {
    
75.                 break;
    
76.             }
    
77. 
    
78.             stepB = stepA;
    
79.             stepC.DecLShift( 1 );
    
80.             stepA.DecLShift( 1 );
    
81. 
    
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.             }
    
97. 
    
98. 
    
99.             ++stepA.DecLShift( 2 );
    
100.         }
     
101.     }
     
102. 
     
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位数字
   最小等于不同数字个数加一