src/complex/transcendental/cl_C_log2.cc

   1 // log().
   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 "base/cl_N.h"
  16
  17 namespace cln {
  18
  19 const cl_N log (const cl_N& a, const cl_N& b)
  20 {
  21 // Methode:
  22 // (log a b) =
  23 //   falls b reell, >0:
  24 //     (complex (/ (log (abs a)) (log b)) (/ (phase a) (log b))), genauer:
  25 //     falls a reell, >0: bekannt
  26 //     falls (= a 0): Error
  27 //     sonst: (phase a) errechnen, ein Float.
  28 //            b (falls rational) ins selbe Float-Format umwandeln,
  29 //            Imaginärteil := (/ (phase a) (log dieses_b)).
  30 //            Falls a rational: (log (abs a) b).
  31 //            Falls a komplex mit rationalem Real- und Imaginärteil,
  32 //              Betragsquadrat  (expt (abs a) 2)  exakt ausrechnen als
  33 //              (+ (expt (realpart a) 2) (expt (imagpart a) 2)).
  34 //              Setze  Realteil := (/ (log Betragsquadrat b) 2).
  35 //              [Eventuell wird hierbei (log b) ein zweites Mal ausgerechnet,
  36 //               aber dies sowieso nur in Single-Precision.]
  37 //            Sonst bilde (abs a), ein Float, und (log (abs a)), ein Float,
  38 //              wandle b (falls rational) ins selbe Float-Format um,
  39 //              setze  Realteil := (/ (log (abs a)) (log dieses_b)).
  40 //   sonst: (/ (log a) (log b))
  41         if (realp(b)) {
  42             DeclareType(cl_R,b);
  43             if (plusp(b)) {
  44                 // b ist reell und >0
  45                 if (realp(a)) {
  46                         DeclareType(cl_R,a);
  47                         if (plusp(a))
  48                                 // a und b sind beide reell und >0
  49                                 return log(a,b);
  50                 }
  51                 // b ist reell und >0, a aber nicht.
  52
  53                 // Imaginärteil (/ (phase a) (log b)) errechnen:
  54                 var cl_F im;
  55                 {
  56                         var cl_R angle = phase(a);
  57                         if (eq(angle,0)) // = Fixnum 0 <==> (= a 0) -> Error
  58                                 { throw division_by_0_exception(); }
  59                  {      DeclareType(cl_F,angle);
  60                         var cl_F bf = cl_somefloat(b,angle); // (float b)
  61                         im = angle / ln(bf);
  62                 }}
  63
  64                 // Realteil (/ (log (abs a)) (log b)) errechnen:
  65                 var cl_R re;
  66                 if (realp(a)) {
  67                         DeclareType(cl_R,a);
  68                         if (rationalp(a)) {
  69                                 // a rational -> (log (abs a) b) errechnen:
  70                                 re = log(abs(a),b); // NB: (abs a) > 0
  71                                 goto re_ok;
  72                         }
  73                 } else {
  74                         DeclareType(cl_C,a);
  75                         if (rationalp(realpart(a)) && rationalp(imagpart(a))) {
  76                                 // a komplex mit rationalem Real- und Imaginärteil a1,a2
  77                                 var const cl_R& a1 = realpart(a);
  78                                 var const cl_R& a2 = imagpart(a);
  79                                 re = log(square(a1)+square(a2),b) / 2;
  80                                 goto re_ok;
  81                         }
  82                 }
  83                 // Keine Chance für rationalen Realteil.
  84                 {
  85                         var cl_F abs_a = The(cl_F)(abs(a));
  86                         var cl_F log_abs_a = ln(abs_a);
  87                         var cl_F bf = cl_somefloat(b,log_abs_a); // (float b)
  88                         re = log_abs_a / ln(bf);
  89                 }
  90                 re_ok:
  91
  92                 return complex_C(re,im);
  93             }
  94         }
  95
  96         // normaler komplexer Fall
  97         return log(a) / log(b);
  98 }
  99
 100 }  // namespace cln