取e小数点后20个10位数的质数，求debug！

wjy8888840 · 发表于 2014-9-26 21:21:31

运行后发现第一部分求e的后10000位可以得到结果
第二部分就卡住了

我觉得可能是定义或者mod的问题，因为求的是10位数的质数，用tmp（10位数组）除数得到的余数不准确附了代码
求指明方向！！！！拜托了

[Fortran] 纯文本查看 复制代码

01Program jiayi
02 
03  Implicit Real *8(A, C-H, O-S, U-Z)
04  Integer, Parameter :: nk = 10000
05  Dimension :: ns(nk), na(nk), nt(nk), prime(20)
06  Integer :: i, j, b, it, l
07 
08  Do i = 1, nk
09    na(i) = 0
10    ns(i) = 0
11    nt(i) = 0
12  End Do
13 
14  na(1) = 1
15  ns(1) = 2
16 
17  n = 1
18  Do j = 1, 5000
19    n = n + 1
20 
21    Do i = 2, nk
22      nd = mod(na(i-1), n)*10 + nt(i)
23      nt(i) = nd/n
24      na(i) = mod(nd, n)
25    End Do
26 
27    Do i = 1, nk
28      na(i) = nt(i)
29      ns(i) = ns(i) + nt(i)
30    End Do
31 
32    Do i = nk, 2, -1
33      ns(i-1) = ns(i-1) + ns(i)/10
34      ns(i) = mod(ns(i), 10)
35    End Do
36 
37  End Do
38 
39  Print *, 'e=2.'
40  Write (*, '(8x,77i1)')(ns(i), i=2, 77)
41  Write (*, '(80i1)')(ns(i), i=78, nk)
42 
43  b = 1
44 
45loop_utmost: Do i = 2, nk - 10
46    tmp = 0
47    If (ns(i)/=0) Then
48loop_tmp: Do j = i, i + 9 !loop 10 times to get the 10 digit number, and save it as "tmp"
49        tmp = tmp*10 + ns(j)
50      End Do loop_tmp
51 
52      l = sqrt(tmp+1) !choose the squareroot of the number itself as the biggest divisor to save memory
53 
54loop_check: Do it = 2, l !To check whether tmp can be devided by any numbers in the range or not
55        If (mod(tmp,real(it))==0) Then
56          Cycle loop_utmost
57 
58        Else If (it<=l) Then
59          Cycle loop_check
60        Else
61          b = b + 1
62          If (b<=20) Then
63            prime(b) = tmp
64            Cycle loop_utmost
65          Else
66            Exit
67          End If
68        End If
69      End Do loop_check
70    End If
71  End Do loop_utmost
72 
73  Print *, 'The ten digital primes in e:'
74  Write (*, '(i5,f15.1)')(prime(i), i=1, b)
75 
76End Program jiayi

fcode · 发表于 2014-9-27 18:05:50

结果好像不止20个，还挺多

vvt · 发表于 2014-9-27 08:44:14

本帖最后由 vvt 于 2014-9-27 08:50 编辑

有多处问题，我已帮你修改并标出。
计算结果为

The ten digital primes in e:
7427466391.0
7413596629.0
6059563073.0
3490763233.0
2988075319.0
1573834187.0
7021540891.0
5408914993.0
6480016847.0
9920695517.0
1838606261.0
6062613313.0
3845830007.0
1692836819.0
4425056953.0
2505695369.0
5490598793.0
1782154249.0
8215424999.0
9229576351.0

[Fortran] 纯文本查看 复制代码

01Program jiayi
02  Implicit None !// 推荐用 Implicit None
03  Integer, Parameter :: nk = 10000
04  Integer :: ns(nk), na(nk), nt(nk) , n , nd
05  Real(Kind=8) :: prime(20) , tmp !// tmp 要定义，很重要
06  Integer :: i, j, b, it , l
07 
08  Do i = 1, nk
09    na(i) = 0
10    ns(i) = 0
11    nt(i) = 0
12  End Do
13 
14  na(1) = 1
15  ns(1) = 2
16 
17  n = 1
18  Do j = 1, 5000
19    n = n + 1
20 
21    Do i = 2, nk
22      nd = mod(na(i-1), n)*10 + nt(i)
23      nt(i) = nd/n
24      na(i) = mod(nd, n)
25    End Do
26 
27    Do i = 1, nk
28      na(i) = nt(i)
29      ns(i) = ns(i) + nt(i)
30    End Do
31 
32    Do i = nk, 2, -1
33      ns(i-1) = ns(i-1) + ns(i)/10
34      ns(i) = mod(ns(i), 10)
35    End Do
36 
37  End Do
38 
39  Print *, 'e=2.'
40  Write (*, '(8x,77i1)')(ns(i), i=2, 77)
41  Write (*, '(80i1)')(ns(i), i=78, nk)
42 
43  b = 0
44 
45loop_utmost: &
46  Do i = 2, nk - 10
47    tmp = 0.0D0
48    If (ns(i)==0) Cycle loop_utmost
49    Do j = i, i + 9 !loop 10 times to get the 10 digit number, and save it as "tmp"
50      tmp = tmp * 10.0D0 + ns(j)
51    End Do
52    l = sqrt(tmp+1) !choose the squareroot of the number itself as the biggest divisor to save memory
53    Do it = 2, l !To check whether tmp can be devided by any numbers in the range or not
54      If (abs(mod(tmp,real(it)))<0.01D0) Cycle loop_utmost !// real 类型不要做相等判断
55    End Do !// 此处逻辑上有问题，修改了 Do 和 If 代码
56    b = b + 1
57    If (b>20) Exit
58    prime(b) = tmp
59    Cycle loop_utmost
60  End Do loop_utmost
61 
62  Print *, 'The ten digital primes in e:'
63  Write (*, '(f15.1)')(prime(i), i=1, min(b,20)) !// 此处改为 min(b,20)
64 
65End Program jiayi

[求助] 取e小数点后20个10位数的质数，求debug！

评分

规矩勋章

美女勋章

热心勋章

星光勋章

新人勋章

贡献勋章

管理勋章

帅哥勋章

爱心勋章

元老勋章

水王勋章