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

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

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

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

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

Program jiayi

  Implicit Real *8(A, C-H, O-S, U-Z)
  Integer, Parameter :: nk = 10000
  Dimension :: ns(nk), na(nk), nt(nk), prime(20)
  Integer :: i, j, b, it, l

  Do i = 1, nk
    na(i) = 0
    ns(i) = 0
    nt(i) = 0
  End Do

  na(1) = 1
  ns(1) = 2

  n = 1
  Do j = 1, 5000
    n = n + 1

    Do i = 2, nk
      nd = mod(na(i-1), n)*10 + nt(i)
      nt(i) = nd/n
      na(i) = mod(nd, n)
    End Do

    Do i = 1, nk
      na(i) = nt(i)
      ns(i) = ns(i) + nt(i)
    End Do

    Do i = nk, 2, -1
      ns(i-1) = ns(i-1) + ns(i)/10
      ns(i) = mod(ns(i), 10)
    End Do

  End Do

  Print *, 'e=2.'
  Write (*, '(8x,77i1)')(ns(i), i=2, 77)
  Write (*, '(80i1)')(ns(i), i=78, nk)

  b = 1

loop_utmost: Do i = 2, nk - 10
    tmp = 0
    If (ns(i)/=0) Then
loop_tmp: Do j = i, i + 9 !loop 10 times to get the 10 digit number, and save it as "tmp"
        tmp = tmp*10 + ns(j)
      End Do loop_tmp

      l = sqrt(tmp+1) !choose the squareroot of the number itself as the biggest divisor to save memory

loop_check: Do it = 2, l !To check whether tmp can be devided by any numbers in the range or not
        If (mod(tmp,real(it))==0) Then
          Cycle loop_utmost

        Else If (it<=l) Then
          Cycle loop_check
        Else
          b = b + 1
          If (b<=20) Then
            prime(b) = tmp
            Cycle loop_utmost
          Else
            Exit
          End If
        End If
      End Do loop_check
    End If
  End Do loop_utmost

  Print *, 'The ten digital primes in e:'
  Write (*, '(i5,f15.1)')(prime(i), i=1, b)

End Program jiayi

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] syntaxhighlighter_viewsource syntaxhighlighter_copycode

Program jiayi
  Implicit None !// 推荐用 Implicit None
  Integer, Parameter :: nk = 10000
  Integer :: ns(nk), na(nk), nt(nk) , n , nd
  Real(Kind=8) :: prime(20) , tmp !// tmp 要定义，很重要
  Integer :: i, j, b, it , l

  Do i = 1, nk
    na(i) = 0
    ns(i) = 0
    nt(i) = 0
  End Do

  na(1) = 1
  ns(1) = 2

  n = 1
  Do j = 1, 5000
    n = n + 1

    Do i = 2, nk
      nd = mod(na(i-1), n)*10 + nt(i)
      nt(i) = nd/n
      na(i) = mod(nd, n)
    End Do

    Do i = 1, nk
      na(i) = nt(i)
      ns(i) = ns(i) + nt(i)
    End Do

    Do i = nk, 2, -1
      ns(i-1) = ns(i-1) + ns(i)/10
      ns(i) = mod(ns(i), 10)
    End Do

  End Do

  Print *, 'e=2.'
  Write (*, '(8x,77i1)')(ns(i), i=2, 77)
  Write (*, '(80i1)')(ns(i), i=78, nk)

  b = 0

loop_utmost: &
  Do i = 2, nk - 10
    tmp = 0.0D0
    If (ns(i)==0) Cycle loop_utmost
    Do j = i, i + 9 !loop 10 times to get the 10 digit number, and save it as "tmp"
      tmp = tmp * 10.0D0 + ns(j)
    End Do
    l = sqrt(tmp+1) !choose the squareroot of the number itself as the biggest divisor to save memory
    Do it = 2, l !To check whether tmp can be devided by any numbers in the range or not
      If (abs(mod(tmp,real(it)))<0.01D0) Cycle loop_utmost !// real 类型不要做相等判断
    End Do !// 此处逻辑上有问题，修改了 Do 和 If 代码
    b = b + 1
    If (b>20) Exit
    prime(b) = tmp
    Cycle loop_utmost
  End Do loop_utmost

  Print *, 'The ten digital primes in e:'
  Write (*, '(f15.1)')(prime(i), i=1, min(b,20)) !// 此处改为 min(b,20)

End Program jiayi