2

我希望我没有重复任何问题，但建议的主题没有提供任何类似的问题。我有一个检查数字是否是素数的功能。现在这是寻找素数的最慢方法。

subroutine is_prime_slow(num, stat)
        implicit none
        logical :: stat
        integer :: num
        integer :: i
        if ((num .le. 3) .and. (num .gt. 1)) then
            stat = .true.
            return
        end if

        ! write(*,*) 'limit = ',limit
        do i = 2,num - 1
            ! write(*,*) 'mod(',limit,i,') = ',mod(limit,i)
            if (mod(num,i) == 0) then
                stat = .false.
                return
            end if
        end do
        stat = .true.   
        return     
    end

现在假设我对其进行了一些改进。

subroutine is_prime_slow(num, stat)
    implicit none
    logical :: stat
    integer :: num
    integer :: i
    if ((num .le. 3) .and. (num .gt. 1)) then
        stat = .true.
        return
    end if
    ! IMPROVEMENT
    if ((mod(num,2) == 0) .or. (mod(num,3) == 0) .or. (mod(num,5) == 0) .or. (mod(num,7) == 0)) then
        stat = .false.
        return
    end if

    ! write(*,*) 'limit = ',limit
    do i = 2,num - 1
        ! write(*,*) 'mod(',limit,i,') = ',mod(limit,i)
        if (mod(num,i) == 0) then
            stat = .false.
            return
        end if
    end do
    stat = .true.   
    return     
end

现在第二个版本排除了一半以上的数字，例如所有可以被 2、3、5、7 整除的数字。当我使用 linux 'time' 程序定时执行时，'改进'版本的执行速度怎么可能一样慢？我错过了一些明显的技巧吗？

Searching the first 900000 numbers:
1st: 4m28sec
2nd  4m26sec
4

2 回答 2

7

无论如何，2、3、5 和 7 的倍数很快就会被原始算法拒绝，因此跳过它们根本不会提高性能。该算法花费大部分时间的地方是证明具有大素因数的数字是复合的。要从根本上提高性能，您应该使用更好的primality test，例如 Miller-Rabin。

一个更简单的改进是只测试高达sqrt(num)，而不是num-1. 如果这没有直接意义，请考虑一个合数的最小质因数可以有多大。此外，如果您正在寻找从 1 到 N 的素数，使用素数筛或仅通过您已经找到的素数测试可分性可能会更有效。

于 2013-09-22T18:01:44.860 回答
2

我最近刚刚编写了类似的代码;-)

! Algorithm taken from https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
subroutine eratosthenes_sieve(n, primes)
  implicit none
  integer,intent(in)    :: n
  integer,allocatable,intent(out)   :: primes(:)
  integer               :: i, j, maxPrime, stat
  logical               :: A(n)

  maxPrime = floor(sqrt(real(n)))
  A = .true.

  do i=2,maxPrime
    j = i*i
    do
      A(j) = .false.
      j = j + i ; if ( j .gt. n ) exit
    enddo
  enddo !i

  allocate( primes( count(A)-1 ), stat=stat )
  if ( stat /= 0 ) stop 'Cannot allocate memory!'

  j = 1
  do i=2,n ! Skip 1
    if ( .not. A(i) ) cycle
    primes( j ) = i
    j = j + 1 ; if ( j > size(primes) ) exit
  enddo
end subroutine

该算法为您提供不超过某个数字的所有素数，因此您可以轻松检查您的素数是否包括在内：

if ( any(number == prime) ) write(*,*) 'Prime found:',number
于 2013-09-22T18:01:37.493 回答