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

Why does -Infinity not show as true in isnan - it always shows false. It is not a number?

I can fix it by checking if XFAC is zero, but that is a kludge?

Link Copied

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

+ and - infinity are numbers, they are just too large to represent.

A Not A Number is an invalid bit pattern, provided to cause a hardware trap (Signalling NAN), or to propigate NAN's without hardware trap (Quiet NAN).

https://en.wikipedia.org/wiki/IEEE_754

### Infinities

The infinities of the extended real number line can be represented in IEEE floating-point datatypes, just like ordinary floating-point values like 1, 1.5, etc. They are not error values in any way, though they are often (depends on the rounding) used as replacement values when there is an overflow. Upon a divide-by-zero exception, a positive or negative infinity is returned as an exact result. An infinity can also be introduced as a numeral (like C's "INFINITY" macro, or "∞" if the programming language allows that syntax).

IEEE 754 requires infinities to be handled in a reasonable way, such as

- (+∞) + (+7) = (+∞)
- (+∞) × (−2) = (−∞)
- (+∞) × 0 = NaN – there is no meaningful thing to do

A fast way to test for both is to test for exponent containing all 1's

Single precision exponent is 8-bits,

Double precision exponent is 12-bits

Quadruple precision is 15-bits

Use TRANSFER to convert (reinterpret) as the appropriate integer type, mask out everything except the exponent with IAND, then test for all 1's in the exponent field (msb is sign, next is exponent, finally mantissa).

This is much faster than using the IEEE_xxx functions.

Jim Dempsey

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

Jim

If I multiple -Infinity by 0 in Intel Fortran I get 0. So the last test is not possible.

```
if(XFAC .lt. 0.0001) then
write(*,*)" A lambda factor is zero or lower"
XFactor(IA,2) = 0.01
XFactor(IA,3) = 0.01
dum = log(XFAC)
dum1 = log10(XFAC)
ANaN(2) = isnan(dum*0)
ANaN(3) = isnan(dum1*0)
XX = 2143289344
print *, transfer(dum*0.0, 1.0)
print *, transfer(XX, 1.0)
343 Format(" Error NaN in line :: ",i6," Error :: ",f10.3," Logical ",L)
if(ANaN(2) .eq. .TRUE.) then
write(*,343)IA,XFAC,ANaN(2)
endif
```

I tried Transfer and got the following , the XX generates NaN, but dum or dum*0 does not. The only test that works is to check if less than some small value close to zero.

A bit to one side:

Getting a set of random number vectors that are independent is simple in Basic and very blasted hard in Fortran. The code to make sure each vector is different looks like this

```
!---------------------------------------------------------------------------------------------------------------------
!
! Sixth Gaussian Data Set
!
!---------------------------------------------------------------------------------------------------------------------
code = 6
allocate (RSIGMF(nTB), STAT = ALLOC_ERR)
if(ALLOC_ERR /= 0) then
write(*,110)code
end if
call get_seed(mul,6)
write(*,20)mul,mul**0.9645*sqrt(float(oldmul))*sqrt(float(oldmul))*sqrt(float(oldmul))
mul = mul**0.9645*sqrt(float(oldmul))*sqrt(float(oldmul))*sqrt(float(oldmul))
if(mul .gt. 400000) then
mul = sqrt(float(mul))
endif
if(mul .eq. oldmul) then
pause
stop " Not random F"
else
oldmul = sqrt(float(abs(mul)))
endif
call MKL_VSL_TEST(mul,rsigmF,amean,sigma,NTB)
```

1. take the time from the clock, use the hundreds and seconds to get a seed.

2. Then multiply every seed to a different power every time and use the sqrt of the last seed to generate a new mul, need to check it is not a large negative number of the sqrt or the whole mess goes birko. Remember Sgt Birko?

3. that then gives an independent set.

This problem has been a problem since I was trying random numbers in the late 1980s, Basic is so simple and works so well and Fortran is such a pain.

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

That algorithm is computationally expensive.

Can you try this for randomness:

```
module mod_rng
type rng_t
union
map
integer(8) :: seed
end map
map
integer(2) :: seedLow
integer(4) :: seedMid
integer(2) :: seedHigh
end map
end union
union
map
real :: temp
end map
map
integer(4) :: tempAsInteger
end map
end union
integer(4) :: bigPrime = 1000008259
integer(4) :: mask = Z'7FFFFF'
end type rng_t
type(rng_t) :: my_rng
contains
function rng() result(ret)
implicit none
real :: ret
logical, save :: doOnce = .true.
if(doOnce) then
call SYSTEM_CLOCK(my_rng%seed)
my_rng%seed = my_rng%seed + my_rng%seed + 1 ! assure odd number
doOnce = .false.
endif
! compute next seed
my_rng%seed = my_rng%seed * my_rng%bigPrime
my_rng%temp = 1.0
my_rng%tempAsInteger = ior(my_rng%tempAsInteger, iand(my_rng%seedMid, my_rng%mask))
ret = my_rng%temp - 1.0
end function rng
end module mod_rng
program Console26
use mod_rng
implicit none
integer :: i
do i=1, 10
print *, rng()
end do
end program Console26
```

That should get you 2^23-1 pseudo random numbers. With a little extra code we can get one more bit.

I haven't tested this (and it may be equivalent to one of the rng's in MKL.

But since you like playing with random number generators...

Jim Dempsey

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

I will note that if you are depending on IEEE exceptional value behavior, be sure to also specify /fp:strict.

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

Jim:

I do not enjoy statistics, I studied pure math and then did engineering to get a job, but stat is a pain, necessary but a pain.

I tried your program, it works quite well and generates a nice uniform distribution, the residuals compared between a normal PDF and your results is a strange Gaussian, a negative logisitic with some skew warp.

I will store it away and when I get a chance add it to the program.

I need Gaussian, but the blasted MKL is sensitive to the seed. I spent an enjoyable day playing with the MKL, their program needs a seed that is huge, in the order of several thousand. The manual from memory suggests 1 or 0 and this is a washout.

The problem is engineering often has limited data and the engineers just say Gaussian and get on with it.

If we add the Box Mueller algorithm, it will be Gaussian, I have some code somewhere.

In the end it was easiest to check for a zero input. I was wondering if I set a equal to 0ne/zero then I have infinity and I can compare the variable that is coming back as infinity ???

Thanks for the help.

John

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

>>I tried your program, it works quite well and generates a nice uniform distribution, the residuals compared between a normal PDF and your results is a strange Gaussian, a negative logisitic with some skew warp.

The provide algorithm was built for speed, but there is a gap in the distribution (that can be filled in with additional code). For 32-bit floating point you have

s eeeeeeee mmmmmmmmmmmmmmmmmmmmmmm

s = sign

e's = an 8-bit exponent (shift count+127)

m's = fraction sans the leading 1 bit (when number not exactly 0)

A floating point 1.0 will have e's of 127 and 0's for the m's (127 - 127 = 0 = no shift left or right).

The integer odd number seed multiplied by an integer prime number yields an odd number.

We extract the middle 23 bits which are pseudo random.

Now, I suspect the issue (peculiarity) you observed results from when the 23-bit pseudo random bit pattern has a leading 0, that when the hidden leading 1 bit of the 1.0 is removed, that when the result is normalized, that the 23-bit bit pattern will be shifted left (incrementing exponent each time) until a 1-bit is shifted out of the 23-bit bit pattern. This 1'bit becomes the hidden bit.

Now then, what is shifted in on the right side are 0 bits as opposed to random bit patterns. And this is what is the gap in the distribution.

To correct for this (or at least fill in the gap better), would be to fill in the shifted in 0's with that number of random bits.

To to this, after the -1.0, you store back into temp, extract the e's from tempAsInteger, rotate right 23 bits, subtract 127 to get the number of lsb 0-bits that require to be filled in with random bits. While this may sound complicated.

```
temp = temp - 1.0 ! instead of into ret
count=IAND(ISHIFT(tempAsInteger,-23), Z'FF') - 127 ! # trailing 0-bits
mask = ISHIFT(1,count) - 1
seed = seed * bigPrime
IOR(tempAsInteger,IAND(seedMid,mask))
ret = temp
```

Excepting for this last part, the fast rng algorithm is what I used on a PDP-8I back in the late 1960's.

Fast and reasonably good results.

Jim Demspey

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

Minor debugging issues with the improved code suggestion:

```
module mod_rng
type rng_t
union
map
integer(8) :: seed
end map
map
integer(2) :: seedLow
integer(4) :: seedMid
integer(2) :: seedHigh
end map
end union
union
map
real :: temp
end map
map
integer(4) :: tempAsInteger
end map
end union
integer(4) :: bigPrime = 1000008259
integer(4) :: mask = Z'7FFFFF'
end type rng_t
type(rng_t) :: my_rng
contains
function rng() result(ret)
implicit none
real :: ret
logical, save :: doOnce = .true.
integer :: count, mask
if(doOnce) then
call SYSTEM_CLOCK(my_rng%seed)
my_rng%seed = my_rng%seed + my_rng%seed + 1 ! assure odd number
doOnce = .false.
endif
! compute next seed
my_rng%seed = my_rng%seed * my_rng%bigPrime
my_rng%temp = 1.0
my_rng%tempAsInteger = ior(my_rng%tempAsInteger, iand(my_rng%seedMid, my_rng%mask))
my_rng%temp = my_rng%temp - 1.0
count = -(IAND(ISHFT(my_rng%tempAsInteger,-23), Z'FF') - 127) ! # trailing 0-bits
mask = ISHFT(1,count) - 1
my_rng%seed = my_rng%seed * my_rng%bigPrime
my_rng%tempAsInteger = IOR(my_rng%tempAsInteger,IAND(my_rng%seedMid,mask))
ret = my_rng%temp
end function rng
end module mod_rng
program Console26
use mod_rng
implicit none
integer :: i
do i=1, 10
print *, rng()
end do
end program Console26
```

Jim Dempsey

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

Jim:

It is a pretty good uniform generator,

These are the difference between exact and yours for 1 million entries in 0.1 bins, run three times.

Q1: Does an FFT show a pattern? My guess is yes.

Q2: Is 0.3 always low, or is it just not enough tests. My guess is not enough tests, but it would be a nice exercise for a class to discover.

Can you check the code please, it needed a few fixes, such as ISFHT and not ISHIFT so I assume you were just using pseudo code?

```
module mod_rng
type rng_t
union
map
integer(8) :: seed
end map
map
integer(2) :: seedLow
integer(4) :: seedMid
integer(2) :: seedHigh
end map
end union
union
map
real :: temp
end map
map
integer(4) :: tempAsInteger
end map
end union
integer(4) :: bigPrime = 1000008259
integer(4) :: mask = Z'7FFFFF'
end type rng_t
type(rng_t) :: my_rng
contains
function rng() result(ret)
implicit none
real :: ret
integer count
integer(1) mask
integer res
logical, save :: doOnce = .true.
if(doOnce) then
call SYSTEM_CLOCK(my_rng%seed)
my_rng%seed = my_rng%seed + my_rng%seed + 1 ! assure odd number
doOnce = .false.
endif
! compute next seed
my_rng%seed = my_rng%seed * my_rng%bigPrime
my_rng%temp = 1.0
my_rng%tempAsInteger = ior(my_rng%tempAsInteger, iand(my_rng%seedMid, my_rng%mask))
my_rng%temp = my_rng%temp - 1.0 ! instead of into ret
count=IAND(ISHFT(my_rng%tempAsInteger,-23), Z'FF') - 127 ! # trailing 0-bits
mask = ISHFT(1,count) - 1
my_rng%seed = my_rng%seed * my_rng%bigPrime
res = IOR(my_rng%tempAsInteger,IAND(my_rng%seedMid,mask))
ret = my_rng%temp
end function rng
end module mod_rng
program Console26
use mod_rng
implicit none
CHARACTER(1) input
integer :: i
open(2,file = "b.csv",status = "UNKNOWN")
do i=1, 1000000
write(2,100)I,rng()
100 Format(I7,",",F10.8)
end do
input = ' '
DO WHILE ((input .NE. 'n') .AND. (input .NE. 'y'))
WRITE (*, 110)
110 Format('Enter y or n: ',\)
READ (*, '(A)') input
END DO
end program Console26
```

We could add the Box Mueller method, this shows a sample, the author is not contactable, and nothing in the code can be copyrighted as it is copied from elsewhere, but it is a standard try to claim ownership.

```
!****************************************************
! 正規乱数生成（by ボックス＝ミューラー法）
! Normal random number generation (by box = Mueller method)
! 0 〜 20 の正規乱数をボックス＝ミューラー法により
!: Normal random numbers from 0 to 20 by the Box-Muller method
! 計算し、ヒストグラムを出力する。
! Compute and output a histogram.
!
! date name version
! 2018.12.03 mk-mode.com 1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!****************************************************
!
module const
implicit none
! SP: 単精度(4), DP: 倍精度(8)
! SP: single precision (4), DP: double precision (8)
integer, parameter :: SP = kind(1.0)
integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
integer(SP), parameter :: M = 20 ! 平均 ! average
real(DP), parameter :: S = 2.5_DP ! 標準偏差 ! standard deviation
integer(SP), parameter :: N = 100000 ! 発生させる乱数の個数 ! Number of random numbers to generate
real(DP), parameter :: PI = 4.0_DP * atan(1.0_DP) ! 円周率 ! Pi
real(DP), parameter :: SC = N / 100.0_DP ! ヒストグラム用スケール ! Histogram Scale
end module const
module box_muller
use const, only : SP, DP, S, M, PI
implicit none
private
public :: rnd_seed, rnd
contains
! 乱数の種の設定
! Set random number seed
subroutine rnd_seed
implicit none
integer(SP) :: seed_size, clock
integer(SP), allocatable :: seed(:)
call system_clock(clock)
call random_seed(size=seed_size)
allocate(seed(seed_size))
seed = clock
call random_seed(put=seed)
deallocate(seed)
end subroutine rnd_seed
! 正規乱数生成
! normal random number generator
! :param(out) integer(4) r(2)
subroutine rnd(r)
integer(SP), intent(out) :: r(2)
real(DP) :: r_u(2) ! [0, 1] の一様乱数 2 個 ! 2 uniform random numbers in [0, 1]
call random_number(r_u(1))
call random_number(r_u(2))
r(1) = int(S * sqrt(-2 * log(r_u(1))) * cos(2 * PI * r_u(2)) + M)
r(2) = int(S * sqrt(-2 * log(r_u(1))) * sin(2 * PI * r_u(2)) + M)
end subroutine rnd
end module box_muller
program rndnum_bm
use const, only : SP, N, M, SC
use box_muller
implicit none
integer(SP) :: i, j, r(2), hist(0:M * 2)
! 乱数の種の設定
! (正規乱数生成に使用する一様乱数の種)
! Set random number seed
! (uniform random number seed for normal random number generation)
call rnd_seed
! 正規乱数の生成
! generate normal random numbers
hist(:) = 0
do i = 1 , N
call rnd(r)
hist(r) = hist(r) + 1
end do
! 結果出力
! result output
do i = 0, M * 2
write (*, '(I3, ":", I5, " | ", A)') &
& i, hist(i), repeat("*", int(hist(i) / SC))
end do
end program rndnum_bm
```

Interestingly enough those guys were in part associated with the Princeton Statistical Group from WW2 fame, I came across a reference every human and their academic dog used as the standard reference for BM Method. I could not find a copy anywhere, so I went to Princeton and asked, a very nice human dug it out of an old box that had not been opened for 50 years and made a pdf. She said I was the first person to ask for it in her time at the place.

You are not supposed to reference stuff you have not read personally.

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

The code I posted earlier than last post, did have two errors (as it was not tested)

1) ISHIFT vs ISHFT

2) more importantly, the conversion of the exponent to shift count had the wrong sign. Resulting in the mask generation being shifted in the wrong direction.

The code posted on 10-08-2023 07:29 AM corrects for these errors.

Sorry for any inconvenience this may have caused.

Jim Dempsey

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