Fortran Coder

标题: 取e小数点后20个10位数的质数,求debug! [打印本页]

作者: wjy8888840    时间: 2014-9-26 21:21
标题: 取e小数点后20个10位数的质数,求debug!
运行后发现 第一部分求e的后10000位可以得到结果
第二部分就卡住了

我觉得可能是定义或者mod的问题,因为求的是10位数的质数,用tmp(10位数组)除数得到的余数不准确附了代码
求指明方向!!!!拜托了
[Fortran] 纯文本查看 复制代码
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
本帖最后由 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] 纯文本查看 复制代码
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



作者: fcode    时间: 2014-9-27 18:05
结果好像不止20个,还挺多




欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2