src/complex/transcendental/cl_C_log2.cc

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