- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Hi

I use ifort compiler. in my code i need random integer number between 1<= i <=64.

i used random_number(x) and it is between 0<x<1. then for my purpose i use this transformation

i = floor(64*x) + 1

but when i use this transformation sometimes i = 65. after debug my code i find that when x = 0.9999

this mistake is happened.I guess the problem is due to rounding error. how to avoid from this mistake?

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

That seems like very unexpected behavior. I find this very concerning since I depend on floor working properly in my programs.

Can you provide a small reproducer? Also, what version of compiler?

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Indeed, you need to provide a reproducer. The only way for floor(64*x) to equal 64 is for x to be equal to 1, and RANDOM_NUMBER is supposed to never return 1; call random_number(x) is supposed to return x such that 0 <= x < 1. Here is a test program that does not use RANDOM_NUMBER that you can use to test your compiler:

program tstflr implicit none real x,dx integer i,ix ! dx=16.0*epsilon(x) do x = 1.0 -dx if(x == 1.0)stop if(floor(64*x) >= 64)then write(*,*)'When x = ',x,' floor(64*x) exceeds 63' stop end if write(*,'(1x,2ES15.8)')x,dx dx=0.5*dx end do end program

The output that I get from it:

9.99998093E-01 1.90734863E-06 9.99999046E-01 9.53674316E-07 9.99999523E-01 4.76837158E-07 9.99999762E-01 2.38418579E-07 9.99999881E-01 1.19209290E-07 9.99999940E-01 5.96046448E-08

This shows that x went through several values > 0.9999 before it became indistinguishable from 1.0 in single precision.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Laura S. wrote:

That seems like very unexpected behavior. I find this very concerning since I depend on floor working properly in my programs.

Can you provide a small reproducer? Also, what version of compiler?

my compiler version is 15.0.0 and 20140723

this is my code and if this mistake don't occur then in line 56 n1 must be printed as 0

but n1 printed as 2

program simple_ising implicit none integer, parameter :: n = 64 real, parameter :: t = 2.269 integer :: s( 1:n , 1:n ) real :: x, y real :: prob real :: finish, start integer :: i, j, ix, iy, r, l, u, d, mag, m1, mn, delta, en, n1 open(unit=1, file='a3') call cpu_time(start) mag = 0 en = 0 ! intial state is random call random_seed do i = 1, n, 1 do j = 1, n, 1 s(i,j) = -1 call random_number(x) if ( x .gt. 0.5 ) then s(i,j) = 1 end if end do end do mag = sum(s) print*,mag/real(n*n) n1=0 do i = 1, n*n, 1 do j = 1, n*n, 1 call random_number(x) ix = floor( n*x + 1.0 ) call random_number(y) iy = floor( n*y + 1.0) r = ix + 1 ; l = ix -1 ; u = iy + 1 ; d = iy - 1 if ( ix + 1 .gt. n ) r = 1 if ( ix - 1 .lt. 1 ) l = n if ( iy + 1 .gt. n ) u = 1 if ( iy - 1 .lt. 1 ) d = n if ( ix .gt. n .or. iy .gt. n) n1=n1+1 mn = s( r, iy) + s( l, iy) + s( ix, u) + s( ix, d ) delta = mn * s( ix, iy ) !delta=0 prob = exp((-2*delta)/t) !prob=1.1 call random_number(x) if ( x .lt. prob ) then s(ix,iy) = -s(ix,iy) en = en + 2*delta mag = mag + 2*s(ix,iy) end if end do write(1,*) i, mag/real(n*n), en/real(n*n) !pause end do call cpu_time(finish) print*,n1, (finish-start) end program simple_ising

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

On Line-48, instead of ix = floor( n*x + 1.0 ), use ix = floor( n*x) + 1. Similarly for iy on Line-48. The following test program illustrates the difference and the importance of truncation immediately after multiplying and before adding 1:

program tflr implicit none integer ix1,ix2 real x ! x=nearest(1.0,-1.0) ix1=floor(64*x+1) ix2=floor(64*x)+1 print *,ix1,ix2 end program

Note that the FLOOR function does not in general have the property

floor(A + B) = floor(A) + floor(B)

and, therefore, you have to write your program accordingly.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

mecej4 wrote:

On Line-48, instead of ix = floor( n*x + 1.0 ), use ix = floor( n*x) + 1. Similarly for iy on Line-48. The following test program illustrates the difference and the importance of truncation immediately after multiplying and before adding 1:

program tflr implicit none integer ix1,ix2 real x ! x=nearest(1.0,-1.0) ix1=floor(64*x+1) ix2=floor(64*x)+1 print *,ix1,ix2 end programNote that the FLOOR function does not in general have the property

floor(A + B) = floor(A) + floor(B)

and, therefore, you have to write your program accordingly.

yes. but in this case this property is true

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

mohammad b. wrote:

Quote:mecej4wrote:

On Line-48, instead of ix = floor( n*x + 1.0 ), use ix = floor( n*x) + 1. Similarly for iy on Line-48. Note that the FLOOR function does not in general have the property

floor(A + B) = floor(A) + floor(B)

and, therefore, you have to write your program accordingly.

yes. but in this case this property is true

No, it is not; if you can advance logical arguments to support that claim, I'll revise my opinion. On the other hand, making the changes (to replace floor(x*64.0+1.0) by floor(x*64)+1, and similarly with y) to your program in #4 and running shows that my explanation is correct (of how ix or iy exceeded the intended range). The result displayed for n1 is 0 after the changes, but you may see various non-zero values in your original program. If you wish to track down this issue more fully, you may set the random number seed to a predetermined value at the beginning of the program.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page