src/real/transcendental/cl_R_atan2.cc

   1 // atan().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_real.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_float.h"
  13 #include "cl_F_tran.h"
  14 #include "cl_N.h"
  15 #include "cl_R.h"
  16
  17 const cl_R atan (const cl_R& x, const cl_R& y)
  18 {
  19 // Methode:
  20 // y=0 -> bei x>0: 0 als Ergebnis,
  21 //        bei x<0: pi als Ergebnis.
  22 //        bei x=0: Error.
  23 // x=0 -> bei y>0: pi/2 als Ergebnis.
  24 //        bei y<0: -pi/2 als Ergebnis.
  25 //        bei y=0: Error.
  26 // Falls x und y beide rational: beide in Floats umwandeln.
  27 // 0 <= |y| <= x  ->  atan(y/x)
  28 // 0 <= |x| <= y  ->  pi/2 - atan(x/y)
  29 // 0 <= |x| <= -y  ->  -pi/2 - atan(x/y)
  30 // 0 <= |y| <= -x  ->  für y>=0: pi + atan(y/x), für y<0: -pi + atan(y/x)
  31
  32         if (eq(y,0)) {
  33                 // y=0 (exakt)
  34                 if (zerop(x)) // x=0 -> Error
  35                         { cl_error_division_by_0(); }
  36                 if (minusp(x)) // x<0 -> pi in Default-Float-Genauigkeit
  37                         { return cl_pi(); }
  38                 return 0; // x>0 -> 0
  39         }
  40         elif (eq(x,0)) {
  41                 // x=0 (exakt)
  42                 if (zerop(y)) // y=0 -> Error
  43                         { cl_error_division_by_0(); }
  44                 if (minusp(y)) // y<0 -> -pi/2
  45                         { return - scale_float(cl_pi(),-1); }
  46                 return scale_float(cl_pi(),-1); // y>0 -> pi/2
  47         } else {
  48                 Mutable(cl_R,x);
  49                 Mutable(cl_R,y);
  50                 // Check special case of rational numbers:
  51                 if (rationalp(x))
  52                         if (rationalp(y)) {
  53                                 // x,y in Floats umwandeln:
  54                                 x = cl_float(The(cl_RA)(x));
  55                                 y = cl_float(The(cl_RA)(y));
  56                         }
  57                 // x,y nicht exakt =0, x/y und y/x werden Floats sein.
  58                 if (abs(x) >= abs(y)) {
  59                         // |x| >= |y|
  60                         var cl_F z = atanx(The(cl_F)(y/x));
  61                         // Division war erfolgreich, also x/=0.
  62                         if (minusp(x))
  63                                 // x<0 -> pi bzw. -pi addieren:
  64                                 if (!minusp(y))
  65                                         // y>=0 -> atan(y/x) + pi
  66                                         return z + cl_pi(z);
  67                                 else
  68                                         // y<0 -> atan(y/x) - pi
  69                                         return z - cl_pi(z);
  70                         else
  71                                 return z;
  72                 } else {
  73                         // |x| < |y|
  74                         var cl_F z = atanx(The(cl_F)(x/y));
  75                         // von pi/2 bzw. -pi/2 subtrahieren:
  76                         if (!minusp(y))
  77                                 // y>=0 -> pi/2 - atan(x/y)
  78                                 return scale_float(cl_pi(z),-1) - z;
  79                         else
  80                                 // y<0 -> -pi/2 - atan(x/y)
  81                                 return - scale_float(cl_pi(z),-1) - z;
  82                 }
  83         }
  84 }