src/real/transcendental/cl_R_atan2.cc

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