src/real/misc/cl_R_expt_I.cc

   1 // expt().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/real.h"
   8
   9
  10 // Implementation.
  11
  12 #include "real/cl_R.h"
  13 #include "cln/rational.h"
  14 #include "cln/integer.h"
  15 #include "integer/cl_I.h"
  16
  17 namespace cln {
  18
  19 // Methode:
  20 // Für y>0:
  21 //   a:=x, b:=y.
  22 //   Solange b gerade, setze a:=a*a, b:=b/2. [a^b bleibt invariant, = x^y.]
  23 //   c:=a.
  24 //   Solange b:=floor(b/2) >0 ist,
  25 //     setze a:=a*a, und falls b ungerade, setze c:=a*c.
  26 //   Ergebnis c.
  27 // Für y=0: Ergebnis 1.
  28 // Für y<0: (/ (expt x (- y))).
  29
  30 // Assume y>0.
  31 inline const cl_R expt_pos (const cl_R& x, const cl_I& y)
  32 {
  33         if (rationalp(x)) {
  34                 DeclareType(cl_RA,x);
  35                 return expt(x,y); // x rational -> schnellere Routine
  36         } else {
  37                 DeclareType(cl_F,x);
  38                 var cl_F a = x;
  39                 var cl_I b = y;
  40                 while (!oddp(b)) { a = square(a); b = b >> 1; }
  41                 var cl_F c = a;
  42                 until (eq(b,1))
  43                   { b = b >> 1;
  44                     a = square(a);
  45                     if (oddp(b)) { c = a * c; }
  46                   }
  47                 return c;
  48         }
  49 }
  50
  51 const cl_R expt (const cl_R& x, const cl_I& y)
  52 {
  53         if (eq(y,0)) { return 1; } // y=0 -> Ergebnis 1
  54         var bool y_negative = minusp(y);
  55         var cl_I abs_y = (y_negative ? -y : y); // Betrag von y nehmen
  56         var cl_R z = expt_pos(x,abs_y); // (expt x (abs y))
  57         return (y_negative ? recip(z) : z); // evtl. noch Kehrwert nehmen
  58 }
  59
  60 }  // namespace cln