/
genprime.f90
77 lines (71 loc) · 1.62 KB
/
genprime.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
program main
implicit none
integer, parameter :: i4b = selected_int_kind(9)
real times, timet, time
integer(kind=i4b) :: strt, stp, x, z
character(len=20) :: buffer
! NOTE: getarg is an extension (not strictly F90/95) but is supported in most
! compilers.
call getarg(1,buffer)
read(buffer,*) strt
call getarg(2,buffer)
read(buffer,*) stp
do x = strt, stp, strt
call cpu_time(times)
Z = genprime(x)
call cpu_time(timet)
time = timet - times
print '(a6, i8, a11, f10.5, a19, i10, a)', &
'Found ',x,' primes in ',time,' seconds (last was ',z,')'
end do
contains
function genprime(mx) result (res)
integer(kind=i4b), intent(in) :: mx
integer(kind=i4b) :: res, cnt, cur
cur = 0
cnt = 0
do while(cnt.lt.mx)
if(isprime(cur)) then
cnt = cnt + 1
end if
cur = cur + 1
end do
res = cur - 1
end function genprime
function isprime(x) result (res)
integer(kind=i4b), intent(in) :: x
integer(kind=i4b) :: lim, y
logical :: res
if(x.lt.2) then
res = .false.
return
end if
if(x.eq.2) then
res = .true.
return
end if
if(modulo(x,2).eq.0) then
res = .false.
return
end if
if(x.lt.9) then
res = .true.
return
end if
if(modulo(x + 1, 6) .ne. 0) then
if(modulo(x - 1, 6) .ne. 0) then
res = .false.
return
end if
end if
lim = int(sqrt(real(x)) + 1.0)
do y = 3, lim, 2
if(modulo(x,y).eq.0) then
res = .false.
return
end if
end do
res = .true.
return
end function isprime
end program main