src/complex/transcendental/cl_C_acos.cc

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