Ambiguity in FATAN2

forth

    Next

  • 1. Forth books on eBay
    I hope this post doesn't offend anyone, but if anyone is in need of reference books, I've put 13 books up for sale as a group on eBay. I'm not trying to get rich, just trying to get them to someone else interested in Forth. :-) Doug

Ambiguity in FATAN2

Postby David N. Williams » Sat, 14 Mar 2009 01:58:56 GMT

I think the ANS specification for FATAN2 is ambiguous:

-----------
12.6.2.1488 FATAN    "f-a-tan"       FLOATING EXT
     ( F: r1 -- r2 ) or ( r1 -- r2 )

     r2 is the principal radian angle whose tangent is r1.

12.6.2.1489 FATAN2   "f-a-tan-two"   FLOATING EXT

     ( F: r1 r2 -- r3 ) or ( r1 r2 -- r3 )

     r3 is the radian angle whose tangent is r1/r2.  An ambiguous
     condition exists if r1 and r2 are zero.

See: A.12.6.2.1489 FATAN2
-----------

Note that the FATAN2 spec is silent about whether the angle is
the principal angle, atan between -pi and +pi.  I wonder if this
is simply an oversight?

For reference, Open Math [1] and ISO C / single UNIX spec 3 [2]
specify the principal angle.

My questions:

1. What did the standards committee actually intend?

2. Are there any ANS-compatible systems that do not return the
    principal angle?

-- David


-----------

[1]  http://www.**--****.com/ #arctan

[2] The Single UNIX Specification Version 3:

        http://www.**--****.com/ 

     Select "Topic", then "Math Interfaces", then "atan2()":

        http://www.**--****.com/ 

Re: Ambiguity in FATAN2

Postby hansoft » Sat, 14 Mar 2009 02:18:47 GMT


[snip]
I really can't tell you. I took the specification from here:

    http://www.**--****.com/ 

I compared the results I got with those of gForth and it seemed fine.
I really can't tell you whether this is a vanilla implementation. Note
that it uses fatan as well for most situations.

Hans Bezemer



\ 4tH library - FATAN2 - Copyright 2008 J.L. Bezemer
\ You can redistribute this file and/or modify it under
\ the terms of the GNU General Public License

[UNDEFINED] fatan2 [IF]
[UNDEFINED] fatan [IF] include lib/asinacos.4th [THEN]

\ x == 0 and y == 0   invalid operation
\ x == 0 and y < 0    -(pi/2)
\ x == 0 and y > 0    pi/2
\ x > 0               arctan (y/x)
\ x < 0 and y < 0     arctan (y/x) - pi
\ otherwise           arctan (y/x) + pi

: fatan2                               ( sin[y] cos[x] -- rad)
  fdup f0= if                          \ if x equals 0
    fdrop fdup f0= if FE.INVALID ferror ! exit then
    f0< pi 2 s>f f/ if fnegate then    \ calculate the radian (equals
pi/2)
  else                                 \ if x doesn't equal zero
    fover fover f/ fatan               \ calculate arctan(y/x)
    fswap f0< if pi frot f0< if fnegate then f+ else fnip then
  then                                 \ adjust accordingly
;
[THEN]


Re: Ambiguity in FATAN2

Postby Andrew Haley » Sat, 14 Mar 2009 03:18:55 GMT










Yes.  I don't doubt that for a moment.

Andrew.

Re: Ambiguity in FATAN2

Postby C. G. Montgomery » Sat, 14 Mar 2009 10:39:29 GMT


The Fortran spec, which may be the original, definitely chooses the
principal angle.

Is it necessary to worry about negative zeros?  Fortran 2003 added the
stipulation that if x < 0 and y is a negative zero, the result should be
(an approximation to) -pi.

I doubt that the ambiguity was intentional.  It should be fixed.

cheers   cgm


Re: Ambiguity in FATAN2

Postby David N. Williams » Sat, 14 Mar 2009 21:54:28 GMT





Thanks!  I'll add 4tH to the list for eventual summary.

-- David

Re: Ambiguity in FATAN2

Postby David N. Williams » Sat, 14 Mar 2009 22:02:54 GMT





Good to know.


Actually, ISO C / single UNIX specifies a long list of special
cases for signed zero, NaN, and +/-Inf.


I think so, too.  My immediate problem is that I'd like to avoid
having to test in complex.fs whether the host system returns the
principal angle.

-- David

Re: Ambiguity in FATAN2

Postby Ed » Sun, 15 Mar 2009 17:06:28 GMT



SwiftForth 3.x:

-1e 1e fatan2 f. 5.4977871  ok
-1e 0e fatan2 f. 4.7123890  ok
-1e -1e fatan2 f. 3.9269908  ok




Re: Ambiguity in FATAN2

Postby David N. Williams » Mon, 16 Mar 2009 01:58:26 GMT




 >> ...
 >> My questions:
 >>
 >> 1. What did the standards committee actually intend?
 >>
 >> 2. Are there any ANS-compatible systems that do not return the
 >>     principal angle?
 >
 > SwiftForth 3.x:
 >
 > -1e 1e fatan2 f. 5.4977871  ok
 > -1e 0e fatan2 f. 4.7123890  ok
 > -1e -1e fatan2 f. 3.9269908  ok

A nonprincipal result!  Thanks!  But I guess it means I'll have
to keep the test for that in complex.fs. :-(

I ran your tests with pfe, gforth, and kforth on Mac OS X ppc
and Linux i686.  They all return the principal angle, e.g.:

gforth 0.7.0:

-1e 1e fatan2 f. -0.785398163397448  ok
-1e 0e fatan2 f. -1.5707963267949  ok
-1e -1e fatan2 f. -2.35619449019234  ok

Out of curiosity I tested the ISO C / Single UNIX 3 spec for the
special floating point numbers by copying and pasting the code
appended at the end of this message.

Pfe, gforth, and kforth agree with those specs on the same OS X
ppc and Linux intel systems, with the following exceptions:

* Gforth on my OS X ppc does not print the "-" in "-0e".  It
   does, however, do minus zero correctly internally.

* Neither gforth nor kforth have NAN or INF words.

I've turned these into Hayes-style, ttester tests here:

 http://www.**--****.com/ ~williams/archive/forth/utilities/#floatingpoint
 http://www.**--****.com/ ~williams/archive/forth/utilities/fatan2-test.fs

The test file conditionally defines NAN, INF, and -INF, so
gforth is satisfied.  It needs to be modified to run with
kforth.  I found it interesting that the tests work with
*exact* ttester floating point comparisons.

-- David


\ If y is +/-0 and x is < 0, +/-pi shall be returned.
  0e -1e fatan2 f.
-0e -1e fatan2 f.
\ If y is +/-0 and x is > 0, +/-0 shall be returned.
  0e  1e fatan2 f.
-0e  1e fatan2 f.
\ If y is < 0 and x is +/-0, -pi/2 shall be returned.
-1e  0e fatan2 f.
-1e -0e fatan2 f.
\ If y is > 0 and x is +/-0, pi/2 shall be returned.
  1e  0e fatan2 f.
  1e -0e fatan2 f.

\ Optional specs:

\ If either x or y is NaN, a NaN shall be returned.
NAN  1e fatan2 f.
  1e NAN fatan2 f.
NAN NAN fatan2 f.
\ If y is +/-0 and x is -0, +/-pi shall be returned.
  0e -0e fatan2 f.
-0e -0e fatan2 f.
\ If y is +/-0 and x is +0, +/-0 shall be returned.
  0e  0e fatan2 f.
-0e  0e fatan2 f.
\ For finite values of +/-y > 0, if x is -Inf, +/-pi shall be returned.
  1e -INF fatan2 f.
-1e -INF fatan2 f.
\ For finite values of +/-y > 0, if x is +Inf, +/-0 shall be returned.
  1e +INF fatan2 f.
-1e +INF fatan2 f.
\ For finite values of x, if y is +/-Inf, +/-pi/2 shall be returned.
+INF  1e fatan2 f.
+INF -1e fatan2 f.
+INF  0e fatan2 f.
+INF -0e fatan2 f.
-INF  1e fatan2 f.
-INF -1e fatan2 f.
-INF  0e fatan2 f.
-INF -0e fatan2 f.
\ If y is +/-Inf and x is -Inf, +/-3pi/4 shall be returned.
+INF -INF fatan2 f.
-INF -INF fatan2 f.
\ If y is +/-Inf and x is +Inf, +/-pi/4 shall be returned.
+INF +INF fatan2 f.
-INF +INF fatan2 f.

Re: Ambiguity in FATAN2

Postby mhx » Mon, 16 Mar 2009 02:39:51 GMT

David N. Williams" < XXXX@XXXXX.COM > writes Re: Ambiguity in FATAN2
[..]



I'm sure it will be considered a bug and fixed, when brought to their attention.

[..]


iForth 3.0.3:
FORTH> -1e 1e fatan2 +e. -7.85398163397448e-001 ok
FORTH> -1e 0e fatan2 +e. -1.57079632679490e+000 ok
FORTH> -1e -1e fatan2 +e. -2.35619449019234e+000 ok


[..]



Note that iForth has +NaN and -NaN, but not NaN.

FORTH> \ If y is +/-0 and x is < 0, +/-pi shall be returned. ok
FORTH> 0e -1e fatan2 f. 3.141593 ok
FORTH> -0e -1e fatan2 f. -3.141593 ok
FORTH> \ If y is +/-0 and x is > 0, +/-0 shall be returned. ok
FORTH> 0e 1e fatan2 f. 0.000000 ok
FORTH> -0e 1e fatan2 f. 0.000000 ok
FORTH> \ If y is < 0 and x is +/-0, -pi/2 shall be returned. ok
FORTH> -1e 0e fatan2 f. -1.570796 ok
FORTH> -1e -0e fatan2 f. -1.570796 ok
FORTH> \ If y is > 0 and x is +/-0, pi/2 shall be returned. ok
FORTH> 1e 0e fatan2 f. 1.570796 ok
FORTH> 1e -0e fatan2 f. 1.570796 ok
FORTH> ok
FORTH> \ Optional specs: ok
FORTH> ok
FORTH> \ If either x or y is NaN, a NaN shall be returned. ok
FORTH> -NAN 1e fatan2 f. -NAN ok
FORTH> 1e -NAN fatan2 f. -NAN ok
FORTH> -NAN -NAN fatan2 f. -NAN ok
FORTH> \ If y is +/-0 and x is -0, +/-pi shall be returned. ok
FORTH> 0e -0e fatan2 f. 3.141593 ok
FORTH> -0e -0e fatan2 f. -3.141593 ok
FORTH> \ If y is +/-0 and x is +0, +/-0 shall be returned. ok
FORTH> 0e 0e fatan2 f. 0.000000 ok
FORTH> -0e 0e fatan2 f. 0.000000 ok
FORTH> \ For finite values of +/-y > 0, if x is -Inf, +/-pi shall be returned. ok
FORTH> 1e -INF fatan2 f. 3.141593 ok
FORTH> -1e -INF fatan2 f. -3.141593 ok
FORTH> \ For finite values of +/-y > 0, if x is +Inf, +/-0 shall be returned. ok
FORTH> 1e +INF fatan2 f. 0.000000 ok
FORTH> -1e +INF fatan2 f. 0.000000 ok
FORTH> \ For finite values of x, if y is +/-Inf, +/-pi/2 shall be returned. ok
FORTH> +INF 1e fatan2 f. 1.570796 ok
FORTH> +INF -1e fatan2 f. 1.570796 ok
FORTH> +INF 0e fatan2 f. 1.570796 ok
FORTH> +INF -0e fatan2 f. 1.570796 ok
FORTH> -INF 1e fatan2 f. -1.570796 ok
FORTH> -INF -1e fatan2 f. -1.570796 ok
FORTH> -INF 0e fatan2 f. -1.570796 ok
FORTH> -INF -0e fatan2 f. -1.570796 ok
FORTH> \ If y is +/-Inf and x is -Inf, +/-3pi/4 shall be returned. ok
FORTH> +INF -INF fatan2 f. 2.356194 ok
FORTH> -INF -INF fatan2 f. -2.356194 ok
FORTH> \ If y is +/-Inf and x is +Inf, +/-pi/4 shall be returned. ok
FORTH> +INF +INF fatan2 f. 0.785398 ok
FORTH> -INF +INF fatan2 f. -0.785398 ok
FORTH>

-marcel


Re: Ambiguity in FATAN2

Postby hansoft » Mon, 16 Mar 2009 16:31:28 GMT

Well, here are the results for 4tH. 4tH neither knows infinity nor
NaN.

---8<---
include lib/ansfloat.4th
include lib/fatan2.4th

fclear
10 set-precision

-1 s>f  1 s>f fatan2 f. cr             \ -0.785398163397448  ok
-1 s>f  0 s>f fatan2 f. cr             \ -1.5707963267949  ok
-1 s>f -1 s>f fatan2 f. cr             \ -2.35619449019234  ok

\ If y is +/-0 and x is < 0, +/-pi shall be returned.
 0 s>f -1 s>f fatan2 f. cr
\ If y is +/-0 and x is > 0, +/-0 shall be returned.
 0 s>f  1 s>f fatan2 f. cr
\ If y is < 0 and x is +/-0, -pi/2 shall be returned.
-1 s>f  0 s>f fatan2 f. cr
\ If y is > 0 and x is +/-0, pi/2 shall be returned.
1 s>f   0 s>f fatan2 f. cr

\ Optional specs:

\ If y is +/-0 and x is +0, +/-0 shall be returned.
0 s>f 0 s>f fatan2 f. cr
---8<---

habe@linux-471m:~/4th.35d> 4th cxq demo/atan2tst.4th
-0.78539816
-1.57079632
-2.35619448
3.141592653
0.000000000
-1.57079632
1.570796326
0.000000000
habe@linux-471m:~/4th.35d>

Hans Bezemer

Re: Ambiguity in FATAN2

Postby Ed » Thu, 19 Mar 2009 11:15:57 GMT

avid N. Williams wrote:

When implementing FATAN2 a year ago, I wrote a test and
compared the results against several forths and MS-Fortran.
That's when I found discrepancies in the angle returned.

With some mathematical functions the result can be expressed
in alternate ways. IMO, a computer language standard should
define functions with a uniform result - else it just creates
complication for programmers. If one form isn't adequate,
then define two functions.

Below are cleaned-up versions of the test I used. It interrogates
all quadrants but has no provision for nan, inf etc. as these are
machine dependent. Included is a high-level FATAN2 which
works equally with floats on the data stack.

\ include FPOUT.F \ for F.R - refer http://dxforth.webhop.org/forth

0 [if] cr .( re-defining FATAN2 ...) cr

3.1415926535897932e fconstant PI

: FATAN2 ( y x -- r )
fdup f0< >r
fdup f0= if
fswap f<
if [ 0.5e pi f* ] fliteral
else [ -0.5e pi f* ] fliteral
then
else
f/ fatan
then
r> if
pi fover f0> if f- else f+ then
then ;

[then]

: .OUT ( x y -- ) cr
fover 1 10 f.r
fdup 1 10 f.r
fswap fatan2 8 20 f.r ;

cr cr
.( x y ATAN2)
cr

1e 0e .out
1e 1e .out
0e 1e .out
-1e 1e .out
cr

-1e 0e .out
-1e -1e .out
0e -1e .out
1e -1e .out
cr
\ 0e 0e .out .( *ambiguous condition* )
cr

real*8 x, y, z

write(*, 10)

write(*, 20)

x = 1
y = 0
z = atan2(y,x)
write(*, 30) x, y, z

x = 1
y = 1
z = atan2(y,x)
write(*, 30) x, y, z

x = 0
y = 1
z = atan2(y,x)
write(*, 30) x, y, z

x = -1
y = 1
z = atan2(y,x)
write(*, 30) x, y, z

write(*, 20)

x = -1
y = 0
z = atan2(y,x)
write(*, 30) x, y, z

x = -1
y = -1
z = atan2(y,x)
write(*, 30) x, y, z

x = 0
y = -1
z = atan2(y,x)
write(*, 30) x, y, z

x = 1
y = -1
z = atan2(y,x)
write(*, 30) x, y, z

write(*, 20)

C x = 0
C y = 0
C z = atan2(y,x)
C write(*, 30) x, y, z

10 format(' x y ATAN2')
20 format()
30 format(1x, 0pf10.1, 0pf10.1, 0pf20.8)

end

MS-Fortran 5.1:
x y ATAN2

1.0 .0 .00000000
1.0 1.0 .78539816
.0 1.0 1.57079633
-1.0 1.0 2.35619449

-1.0 .0 3.14159265
-1.0 -1.0 -2.35619449
.0 -1.0 -1.57079633
1.0 -1.0 -.78539816

SwiftForth 3:
x y ATAN2

1.0 0.0 0.00000000
1.0 1.0 0.78539816
0.0 1.0 1.57079630
-1.0 1.0 2.35619440

-1.0 0.0 3.14159260
-1.0 -1.0 3.92699080
0.0 -1.0 4.71238890
1.0 -1.0 5.49778710

VFX 4 (Ndp387.fth): \ seems to want x/y reversed
x y ATAN2

1.0 0.0 1.57079630
1.0 1.0 0.78539816
0.0 1.0 0.00000000
-1.0 1.0

Re: Ambiguity in FATAN2

Postby David N. Williams » Sat, 21 Mar 2009 05:50:24 GMT


 > [...]
 >
 > When implementing FATAN2 a year ago, I wrote a test and
 > compared the results against several forths and MS-Fortran.
 > That's when I found discrepancies in the angle returned.
 >
 > With some mathematical functions the result can be expressed
 > in alternate ways.  IMO, a computer language standard should
 > define functions with a uniform result - else it just creates
 > complication for programmers.  If one form isn't adequate,
 > then define two functions.

I certainly agree!

 > Below are cleaned-up versions of the test I used.  It interrogates
 > all quadrants but has no provision for nan, inf etc. as these are
 > machine dependent.  Included is a high-level FATAN2 which
 > works equally with floats on the data stack.
 >
 > [code and results deleted]

Thanks for this!  I've added your extra x-y values to
fatan2-test.fs, and I think that makes it fairly complete:

 http://www.**--****.com/ ~williams/archive/forth/utilities/fatan2-test.fs

So, with your results for VFX Forth and Win32Forth, plus what
was reported before, we have so far:

                      FATAN2 CONVENTIONS
   implementation    -pi to pi   0 to 2pi
   --------------------------------------
   pfe 0.33.70          yes         no
   gforth 0.7.0         yes         no
   kforth 1.4.0         yes         no
   iForth 3.0.3         yes         no
   4tH 3.5c             yes         no
   SwiftForth 3         no          yes
   VFX Forth  4         yes         no
   Win32Forth 6.13.00   yes         no


So far, nobody seems to think that the ANS 94 standard is
nonambiguous; and there is consensus among the responders that
FATAN2 should return the principal angle.

Please speak up if I've got something wrong!

-- David

Re: Ambiguity in FATAN2

Postby Krishna Myneni » Sat, 21 Mar 2009 08:25:47 GMT




Done.

---
krishna@ungee:~/kforth$ kforth fatan2-test

#ERRORS: 0
  ok
---

The kforth version is at

ftp://ccreweb.org/software/kforth/examples/

KM

Re: Ambiguity in FATAN2

Postby Bernd Paysan » Tue, 24 Mar 2009 19:15:24 GMT


 http://www.**--****.com/ ~williams/archive/forth/utilities/fatan2-test.fs

A suggestion: Define the constants pi/2 etc. in terms of pi:

 pi 4e f/ fconstant   pi/4
-pi 4e f/ fconstant  -pi/4
 pi 2e f/ fconstant   pi/2
-pi 2e f/ fconstant  -pi/2
 pi/2 3e f* fconstant  3pi/2

 pi/4 3e f* fconstant  3pi/4
-pi/4 3e f* fconstant -3pi/4

I'm having troubles in bigForth with your constants for these smaller
values; pi and -pi itself work (maybe not enough precision).

-- 
Bernd Paysan
"If you want it done right, you have to do it yourself"
 http://www.**--****.com/ ~paysan/

Re: Ambiguity in FATAN2

Postby David N. Williams » Wed, 25 Mar 2009 05:29:24 GMT


 > [...]
 > 
 http://www.**--****.com/ ~williams/archive/forth/utilities/fatan2-test.fs
 >
 > A suggestion: Define the constants pi/2 etc. in terms of pi:
 >
 >  pi 4e f/ fconstant   pi/4
 > -pi 4e f/ fconstant  -pi/4
 >  pi 2e f/ fconstant   pi/2
 > -pi 2e f/ fconstant  -pi/2
 >  pi/2 3e f* fconstant  3pi/2
 >
 >  pi/4 3e f* fconstant  3pi/4
 > -pi/4 3e f* fconstant -3pi/4
 >
 > I'm having troubles in bigForth with your constants for these smaller
 > values; pi and -pi itself work (maybe not enough precision).

Sounds like a good suggestion, avoiding most of the decimal to
binary conversions.  On IEEE 754 systems, floating point
multiplication by any power of two conserves accuracy exactly,
as long as there's no exponent overflow or underflow.
Multiplying by other integers should be good to at least 0.5
ulp.

But I'm curious what's going on.  Does the fact that pi and -pi
work for you mean that they're supplied by bigForth, so there's
no conversion from decimal via the conditional definitions in
the test file?

In any event, I've checked (with the help of hexfloat.fs and
Mathematica) that, on my IEEE 64-bit float system, the above
definitions and the 22-digit, rounded decimal definitions both
give the same raw fp values, correct to 53 "exactly rounded"
significand bits.

I've added your definitions, and kept the old ones as a
switchable option.

-- David

Similar Threads:

1.RfD: Ambiguity in FATAN2 version 1.2

2.RfD: Ambiguity in FATAN2

3.RfD: Ambiguity in FATAN2

On Mar 24, 10:41m,  XXXX@XXXXX.COM  (Anton Ertl)
wrote:
> "C. G. Montgomery" < XXXX@XXXXX.COM > writes:
>
> >Other computer languages which provide this function commonly specify the
> >"principal angle", the one between -pi and pi.
>
> >Proposal
>
> >The first sentence of the definition should be expanded to:
> > r3 is the radian angle, greater than -pi and less than or equal to pi,
> >whose tangent is r1/r2.
>
> The C standard specifies "[-pi,pi]", i.e., greater than or equal to
> -pi.
>
> And on Gforth (which uses C's atan2):
> -1e 0e fatan2 f.
>
> outputs
>
> -1.5707963267949
>
> and I guess that Gforth is not the only one.
>

W32F gives the same (assuming the precision is set to 14). Since pi
cannot be represented exactly shouldn't it be the closest
approximations of pi and -pi, allowing for rounding).

> So I suggest including -pi in the range.
>
> - anton

George Hubert

4.CfV: Ambiguity in FATAN2

5.RfD: Ambiguity in FATAN2 version 1.1

6. printf format ambiguity

7. Ambiguity in semantics of assignments?

8. pragma import ambiguity



Return to forth

 

Who is online

Users browsing this forum: No registered users and 96 guest