src/complex/transcendental/cl_C_expt_C.cc

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