src/complex/transcendental/cl_C_acos.cc

   1 // acos().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/complex.h"
   8
   9
  10 // Implementation.
  11
  12 #include "complex/cl_C.h"
  13 #include "cln/real.h"
  14 #include "real/cl_R.h"
  15 #include "cln/rational.h"
  16 #include "rational/cl_RA.h"
  17 #include "cln/float.h"
  18
  19 namespace cln {
  20
  21 inline const cl_F pi (const cl_R& v)
  22 {
  23         if (rationalp(v))
  24                 return pi();
  25         else {
  26                 DeclareType(cl_F,v);
  27                 return pi(v);
  28         }
  29 }
  30
  31 const cl_N acos (const cl_N& z)
  32 {
  33 // Methode:
  34 // Wert und Branch Cuts nach der Formel CLTL2, S. 312:
  35 //   arccos(z) = log(z+i*sqrt(1-z^2))/i = pi/2 - arcsin(z)
  36 // Sei z=x+iy.
  37 // Falls y=0:
  38 //   Falls x rational:
  39 //     Bei x=1: Ergebnis 0.
  40 //     Bei x=1/2: Ergebnis pi/3.
  41 //     Bei x=0: Ergebnis pi/2.
  42 //     Bei x=-1/2: Ergebnis 2pi/3.
  43 //     Bei x=-1: Ergebnis pi.
  44 //     Sonst x in Float umwandeln.
  45 //   Falls x>1: Ergebnis i ln(x+sqrt(x^2-1)).
  46 // Sonst errechne u+iv = arsinh(-y+ix) wie oben, Ergebnis (pi/2-v)+iu.
  47
  48         cl_C_R u_v;
  49         if (realp(z)) {
  50                 DeclareType(cl_R,z);
  51                 // y=0
  52                 var const cl_R& x = z;
  53                 var cl_F xf;
  54                 if (rationalp(x)) {
  55                         DeclareType(cl_RA,x);
  56                         // x rational
  57                         if (integerp(x)) {
  58                                 DeclareType(cl_I,x);
  59                                 // x Integer
  60                                 if (eq(x,0)) // x=0 -> Ergebnis pi/2
  61                                         return scale_float(pi(),-1);
  62                                 if (eq(x,1)) // x=1 -> Ergebnis 0
  63                                         return 0;
  64                                 if (eq(x,-1)) // x=-1 -> Ergebnis pi
  65                                         return pi();
  66                                 xf = cl_float(x);
  67                         } else {
  68                                 DeclareType(cl_RT,x);
  69                                 // x Ratio
  70                                 if (eq(denominator(x),2)) { // Nenner = 2 ?
  71                                         if (eq(numerator(x),1)) // x=1/2 -> Ergebnis pi/3
  72                                                 return pi()/3;
  73                                         if (eq(numerator(x),-1)) // x=-1/2 -> Ergebnis 2pi/3
  74                                                 return scale_float(pi(),1)/3;
  75                                 }
  76                                 xf = cl_float(x);
  77                         }
  78                 } else {
  79                         DeclareType(cl_F,x);
  80                         xf = x;
  81                 }
  82                 // x Float
  83          {      var cl_F& x = xf;
  84                 if (cl_I(1) < x)
  85                         // x>1
  86                         return complex_C(0,ln(x+sqrt(square(x)-1)));
  87                 u_v = asinh(0,x);
  88         }} else {
  89                 DeclareType(cl_C,z);
  90                 u_v = asinh(-imagpart(z),realpart(z));
  91         }
  92         var cl_R& u = u_v.realpart;
  93         var cl_R& v = u_v.imagpart;
  94         var cl_F archimedes = pi(v); // pi im Float-Format von v
  95         return complex(scale_float(archimedes,-1)-v,u); // (pi/2-v)+iu
  96 }
  97
  98 }  // namespace cln