src/complex/transcendental/cl_C_expt_C.cc

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