3

I've implemented code to find the polar coordinates of a point in 2D space. if the point lies in the 1st or 2nd Quadrant, 0<=theta<=pi and if it lies in the 3rd or 4th Quadrant, -pi <= theta <= 0.

      module thetalib

      contains 
      real function comp_theta( x1, x2)
      implicit none

      real        , intent(in)    :: x1, x2
      real                        :: x1p, x2p
      real                        :: x1_c=0.0, x2_c=0.0
      real                        :: pi=4*atan(1.0)

      x1p = x1 - x1_c
      x2p = x2 - x2_c

!  - Patch
      !if ( x1p == 0 .and. x2p /= 0 ) then
      !   comp_theta = sign(pi/2.0, x2p)
      !else
      !   comp_theta = atan ( x2p / x1p )
      !endif

      comp_theta = atan( x2p / x1p)

      if ( x1p >= 0.0 .and. x2p >= 0.0 ) then
         comp_theta = comp_theta
      elseif ( x1p < 0 .and. x2p >= 0.0 ) then
         comp_theta = pi + comp_theta
      elseif( x1p < 0.0 .and. x2p < 0.0 ) then
         comp_theta = -1* (pi - comp_theta)
      elseif ( x1p >= 0.0 .and. x2p < 0.0 ) then
         comp_theta = comp_theta
      endif

      return
      end function comp_theta

      end module thetalib

      program main

      use thetalib

      implicit none

!     Quadrant 1
      print *, "(0.00, 1.00): ", comp_theta(0.00, 1.00)
      print *, "(1.00, 0.00): ", comp_theta(1.00, 0.00)
      print *, "(1.00, 1.00): ", comp_theta(1.00, 1.00)

!     Quadrant 2
      print *, "(-1.00, 1.00): ", comp_theta(-1.00, 1.00)
      print *, "(-1.00, 0.00): ", comp_theta(-1.00, 0.00)

!     Quadrant 3
      print *, "(-1.00, -1.00): ", comp_theta(-1.00, -1.00)


!     Quadrant 4
      print *, "(0.00, -1.00): ", comp_theta(0.00, -1.00)
      print *, "(1.00, -1.00): ", comp_theta(1.00, -1.00)

      end program main

In the function thetalib::comp_theta, when there is a division by zero and the numerator is +ve, fortran evaluates it to be -infinity and when the numerator is -ve, it evaluates it to be +infinity ( see output )

 (0.00, 1.00):   -1.570796    
 (1.00, 0.00):   0.0000000E+00
 (1.00, 1.00):   0.7853982    
 (-1.00, 1.00):    2.356194    
 (-1.00, 0.00):    3.141593    
 (-1.00, -1.00):   -2.356194    
 (0.00, -1.00):    1.570796    
 (1.00, -1.00):  -0.7853982

This baffled me. I've also implemented the patch you see to work around it. And to investigate it further, I setup a small test:

  program main

  implicit none

  real          :: x1, x2

  x1 = 0.0 - 0.0 ! Reflecting the x1p - 0.0
  x2 = 1.0

  write(*,*) "x2/x1=", x2/x1

  x2 = -1.0
  write(*,*) "x2/x1=", x2/x1

  end program main

This evaluates to:

 x2/x1=       Infinity
 x2/x1=      -Infinity

My fortran version:

$ ifort --version
ifort (IFORT) 19.0.1.144 20181018
Copyright (C) 1985-2018 Intel Corporation.  All rights reserved.

And I have three questions:

  1. Why there are signed infinite values?
  2. How are the signs determined?
  3. Why does infinity take the signs shown in outputs for both thetalib::comp_theta and the test program?
newkid
  • 1,368
  • 1
  • 11
  • 27
  • 2
    Won't it be easier to use the `atan2(y,x)` function than to spin you own? – Steve Jan 24 '19 at 18:29
  • @Steve yeah definitely for this particular case but I'd like to know why the signs don't match. – newkid Jan 24 '19 at 18:34
  • 1
    Are you asking why there are signed infinite values, why each is the result in the short test, or how to avoid them? – francescalus Jan 24 '19 at 20:24
  • @francescalus yes, (a) why there are signed infinite values? (b) how the signs are determined - both in the short test case as well as `thetalib::comp_theta`?How to avoid would follow from the answers to (a) and (b). Will edit the question to reflect this. – newkid Jan 24 '19 at 20:27
  • @IAmNerd2000 I try replicating it with the test program in the line `x1 = 0.0 - 0.0` but with different results. Limits was the first thing which crossed my mind too but I wasn't able to reproduce it. – newkid Jan 24 '19 at 21:03
  • You can also try checking if the number equals itself. if not then. it is infinite. EX: `if ( x2x1 .eq. x2x1)` then GOOD number. if not then infinity. It also may be that the value holding x1 is calculated by the computer where all bits in the number are set to 1 (-infinity) and simply dividing 1/(allbits =1) => changes the sign bit of the value to -infinity. same for -1/(allbits=1) changes sign bit to +Infinity. – IAmNerd2000 Jan 24 '19 at 21:16

2 Answers2

1

You can also try checking if the number equals itself. if not then. it is infinite.

EX: if ( x2x1 .eq. x2x1) then GOOD number. if not then infinity.

It also may be that the value holding x1 is calculated by the computer where all bits in the number are set to 1 (-infinity) and doing a bitwise divide you get the following:

which, is actually a subtraction operation where (0....001 - 01...111) = -Infinity and (0....001 - 11.....111) = +Infinity I would look up bitwise divide and see the info on that. More is done, but I don't need to explain the details.

IAmNerd2000
  • 761
  • 1
  • 4
  • 12
  • This explains a part of the answer i.e. the output in the test program which is already very nice. Thanks! Any ideas on what happens in the function `comp_theta` ? – newkid Jan 24 '19 at 21:48
  • I notice that you have the following `x1p < 0` in your elseif statement when it should be `x1p < 0.0` otherwise it is evaluating x1p as an integer – IAmNerd2000 Jan 24 '19 at 23:57
1

That there are signed infinite values follows from the compiler supporting IEEE arithmetic with the real type.

For motivation, consider real non-zero numerator and denominator. If these are both of the same sign then the quotient is a real (finite) positive number. If they are of opposite sign the quotient is a real (finite) negative number.

Consider the limit 1/x as x tends to zero from below. For any strictly negative value of x the value is negative. For continuity considerations the limit can be taken to be negative infinity.

So, when the numerator is non-zero, the quotient will be positive infinity if the numerator and denominator are of the same sign, and negative if of opposite sign. Recall also, that the zero denominator may be signed.

If you wish to examine the number, to see whether it is finite you can use the procedure IEEE_IS_FINITE of the intrinsic module ieee_arithmetic. Further, that module has the procedure IEEE_CLASS which provides useful information about its argument. Among other things:

  • whether it is a positive or negative normal number;
  • whether it is a positive or negative infinite value;
  • whether it is a positive or negative zero.
francescalus
  • 30,576
  • 16
  • 61
  • 96