src/complex/algebraic/cl_C_sqrt.cc

   1 // sqrt().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/complex.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_C.h"
  13 #include "cln/real.h"
  14
  15 namespace cln {
  16
  17 const cl_N sqrt (const cl_N& x)
  18 {
  19 // Methode:
  20 // x reell -> Für x>=0 klar, für x<0: sqrt(-x)*i.
  21 // x=a+bi ->
  22 //   Bestimme r=abs(x)=sqrt(a*a+b*b).
  23 //   Falls a>=0: Setze c:=sqrt((r+a)/2), d:=(b/(2*c) falls c>0, c falls c=0).
  24 //   Falls a<0: Setze d:=sqrt((r-a)/2)*(1 falls b>=0, -1 falls b<0), c:=b/(2*d).
  25 //   Damit ist c>=0, 2*c*d=b, c*c=(r+a)/2, d*d=(r-a)/2, c*c-d*d=a, c*c+d*d=r,
  26 //   also c+di die gesuchte Wurzel.
  27         if (realp(x)) {
  28                 DeclareType(cl_R,x);
  29                 if (!minusp(x))
  30                         return sqrt(x);
  31                 else
  32                         return complex_C(0,sqrt(-x));
  33         } else {
  34                 DeclareType(cl_C,x);
  35                 var const cl_R& a = realpart(x);
  36                 var const cl_R& b = imagpart(x);
  37                 var cl_R r = cl_hypot(a,b); // r = (abs x)
  38                 if (!minusp(a)) {
  39                         // a>=0
  40                         var cl_R c = sqrt((r+a)/2);
  41                         var cl_R d = (!zerop(c) ? b/(2*c) : c);
  42                         return complex_C(c,d);
  43                 } else {
  44                         var cl_R d = sqrt((r-a)/2);
  45                         if (minusp(b))
  46                                 d = -d;
  47                         var cl_R c = b/(2*d);
  48                         return complex_C(c,d);
  49                 }
  50         }
  51 }
  52
  53 }  // namespace cln