src/complex/transcendental/cl_C_expt_C.cc

   1 // expt().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/complex.h"
   8
   9 // Implementation.
  10
  11 #include "complex/cl_C.h"
  12 #include "cln/real.h"
  13 #include "real/cl_R.h"
  14 #include "cln/rational.h"
  15 #include "rational/cl_RA.h"
  16 #include "integer/cl_I.h"
  17 #include "base/cl_N.h"
  18
  19 namespace cln {
  20
  21 // Methode:
  22 // Falls y rational:
  23 //   Falls y Integer:
  24 //     Falls y=0: Ergebnis 1,
  25 //       [Nach CLTL folgendermaßen:
  26 //         x reell:
  27 //           x rational -> Fixnum 1
  28 //           x Float -> (float 1 x)
  29 //         x komplex:
  30 //           x komplex rational -> Fixnum 1
  31 //           sonst: #C(1.0 0.0) im Float-Format des Real- bzw. Imaginärteils von x
  32 //       ]
  33 //     Falls x rational oder komplex rational oder |y| klein:
  34 //       x^|y| durch wiederholtes Quadrieren und Multiplizieren und evtl.
  35 //       Kehrwert-Bilden ermitteln.
  36 //     Sonst wie bei 'y Float'.
  37 //   Falls y Ratio m/n:
  38 //     Es gilt (expt x m/n) = (expt (expt x 1/n) m).
  39 //     Falls x in Q(i) liegt (also rational oder komplex rational ist):
  40 //       Sollte x^(m/n) in Q(i) liegen, so auch eine n-te Wurzel x^(1/n)
  41 //       (und bei n=2 oder n=4 damit auch alle n-ten Wurzeln x^(1/n) ).
  42 //       Falls x rational >=0: n-te Wurzel aus x nehmen. Ist sie rational,
  43 //         deren m-te Potenz als Ergebnis.
  44 //       Falls x rational <=0 oder komplex rational und n Zweierpotenz:
  45 //         n-te Wurzel aus x nehmen (mehrfaches sqrt). Ist sie rational oder
  46 //         komplex rational, deren m-te Potenz als Ergebnis.
  47 //         [Beliebige n betrachten!??]
  48 //     Falls n Zweierpotenz und |m|,n klein: n-te Wurzel aus x nehmen
  49 //       (mehrfaches sqrt), davon die m-te Potenz durch wiederholtes
  50 //       Quadrieren und Multiplizieren und evtl. Kehrwert-Bilden.
  51 //     Sonst wie bei 'y Float'.
  52 // Falls y Float oder komplex:
  53 //   Falls (zerop x):
  54 //     Falls Realteil von y >0 :
  55 //       liefere 0.0 falls x und y reell, #C(0.0 0.0) sonst.
  56 //     Sonst Error.
  57 //   Falls y=0.0:
  58 //     liefere 1.0 falls x und y reell, #C(1.0 0.0) sonst.
  59 //   Sonst: (exp (* (log x) y))
  60 // Das Ergebnis liegt in Q(i), falls x in Q(i) liegt und 4y ein Integer ist.??
  61 // Genauigkeit erhöhen, log2(|y|) Bits mehr??
  62 // Bei x oder y rational und der andere Long-Float: bitte kein Single-Float!??
  63
  64 // Liefert x^0.
  65 inline const cl_N expt_0 (const cl_N& x)
  66 {
  67 #ifdef STRICT_CLTL
  68         // y=0 -> 1 im Format von x.
  69         if (realp(x)) {
  70                 DeclareType(cl_R,x);
  71                 if (rationalp(x)) {
  72                         DeclareType(cl_RA,x);
  73                         return 1;
  74                 } else {
  75                         DeclareType(cl_F,x);
  76                         return cl_float(1,x);
  77                 }
  78         } else {
  79                 DeclareType(cl_C,x);
  80                 var cl_R f = contagion(realpart(x),imagpart(x));
  81                 if (rationalp(f)) {
  82                         DeclareType(cl_RA,f);
  83                         return 1;
  84                 } else {
  85                         DeclareType(cl_F,f);
  86                         // #C(1.0 0.0)
  87                         return complex_C(cl_float(1,f),cl_float(0,f));
  88                 }
  89         }
  90 #else
  91         // Exponent exakt 0 -> Ergebnis exakt 1
  92         cl_unused x;
  93         return 1;
  94 #endif
  95 }
  96
  97 inline const cl_R contagion (const cl_N& x)
  98 {
  99         if (realp(x)) {
 100                 DeclareType(cl_R,x);
 101                 return x;
 102         } else {
 103                 DeclareType(cl_C,x);
 104                 return contagion(realpart(x),imagpart(x));
 105         }
 106 }
 107
 108 const cl_N expt (const cl_N& x, const cl_N& y)
 109 {
 110         if (realp(y)) {
 111             DeclareType(cl_R,y);
 112             if (rationalp(y)) {
 113                 DeclareType(cl_RA,y);
 114                 // y rational
 115                 if (integerp(y)) {
 116                         DeclareType(cl_I,y);
 117                         // y Integer
 118                         if (eq(y,0))
 119                                 return expt_0(x); // Liefere 1
 120                         if (fixnump(y)) // |y| klein ?
 121                                 return expt(x,y); // exakt ausrechnen
 122                         if (realp(x)) {
 123                                 DeclareType(cl_R,x);
 124                                 if (rationalp(x)) {
 125                                         DeclareType(cl_RA,x);
 126                                         return expt(x,y); // exakt ausrechnen
 127                                 }
 128                         } else {
 129                                 DeclareType(cl_C,x);
 130                                 if (rationalp(realpart(x)) && rationalp(imagpart(x)))
 131                                         return expt(x,y); // exakt ausrechnen
 132                         }
 133                         // x nicht exakt und |y| groß
 134                 } else {
 135                         DeclareType(cl_RT,y);
 136                         // y Ratio
 137                         var const cl_I& m = numerator(y);
 138                         var const cl_I& n = denominator(y);
 139                         if (realp(x)) {
 140                                 DeclareType(cl_R,x);
 141                                 if (rationalp(x)) {
 142                                         DeclareType(cl_RA,x);
 143                                         if (minusp(x))
 144                                                 goto complex_rational;
 145                                         // x rational >=0
 146                                         var cl_RA w;
 147                                         if (rootp(x,n,&w)) // Wurzel rational?
 148                                                 return expt(w,m);
 149                                 }
 150                         } else {
 151                                 DeclareType(cl_C,x);
 152                                 if (rationalp(realpart(x)) && rationalp(imagpart(x)))
 153                                         goto complex_rational;
 154                         }
 155                         if (false) {
 156                                 complex_rational:
 157                                 // x in Q(i)
 158                                 var uintC k = power2p(n);
 159                                 if (k) {
 160                                         // n Zweierpotenz = 2^(k-1). n>1, also k>1
 161                                         Mutable(cl_N,x);
 162                                         k--;
 163                                         do { x = sqrt(x); }
 164                                            while (--k > 0);
 165                                         return expt(x,m);
 166                                 }
 167                         }
 168                         if (fixnump(m) && fixnump(n)) { // |m| und n klein?
 169                                 var uintV _n = FN_to_UV(n);
 170                                 if ((_n & (_n-1)) == 0) { // n Zweierpotenz?
 171                                         Mutable(cl_N,x);
 172                                         until ((_n = _n >> 1) == 0)
 173                                               { x = sqrt(x); }
 174                                         return expt(x,m);
 175                                 }
 176                         }
 177                 }
 178             }
 179         }
 180         // allgemeiner Fall (z.B. y Float oder komplex):
 181         if (zerop(x)) { // x=0.0 ?
 182                 if (zerop(y))  // y=0.0?
 183                         return expt_0(x); // Liefere 1
 184                 if (rationalp(realpart(y))) // Realteil von y >0 exakt.
 185                         return 0;
 186                 if (!plusp(realpart(y))) // Realteil von y <=0 ?
 187                         throw division_by_0_exception();
 188                 else {
 189                         var cl_R f = contagion(contagion(x),contagion(y));
 190                         // ein Float, da sonst x = Fixnum 0 gewesen wäre
 191                     {   DeclareType(cl_F,f);
 192                         var cl_F f0 = cl_float(0,f);
 193                         return complex_C(f0,f0); // #C(0.0 0.0)
 194                     }
 195                 }
 196         }
 197         return exp(log(x)*y);
 198 }
 199
 200 }  // namespace cln