-Infinity

JohnNichols · ‎10-05-2023

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?

jimdempseyatthecove · ‎10-05-2023

+ 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

Further information on the concept of infinite: Infinity

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

JohnNichols · ‎10-06-2023

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.

jimdempseyatthecove · ‎10-06-2023

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

Steve_Lionel · ‎10-06-2023

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

JohnNichols · ‎10-07-2023

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

jimdempseyatthecove · ‎10-08-2023

>>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

jimdempseyatthecove · ‎10-08-2023

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

JohnNichols · ‎10-08-2023

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.

jimdempseyatthecove · ‎10-08-2023

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