src/complex/misc/cl_C_expt.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
  15 // Methode:
  16 // Für y>0:
  17 //   a:=x, b:=y.
  18 //   Solange b gerade, setze a:=a*a, b:=b/2. [a^b bleibt invariant, = x^y.]
  19 //   c:=a.
  20 //   Solange b:=floor(b/2) >0 ist,
  21 //     setze a:=a*a, und falls b ungerade, setze c:=a*c.
  22 //   Ergebnis c.
  23 // Für y=0: Ergebnis 1.
  24 // Für y<0: (/ (expt x (- y))).
  25
  26 // Assume y>0.
  27 inline const cl_N expt_pos (const cl_N& x, uintL y)
  28 {
  29         var cl_N a = x;
  30         var uintL b = y;
  31         while (!(b % 2)) { a = square(a); b = b >> 1; }
  32         var cl_N c = a;
  33         until (b == 1)
  34           { b = b >> 1;
  35             a = square(a);
  36             if (b % 2) { c = a * c; }
  37           }
  38         return c;
  39 }
  40
  41 const cl_N expt (const cl_N& x, sintL y)
  42 {
  43         if (realp(x)) {
  44                 DeclareType(cl_R,x);
  45                 // x reell -> schnellere Routine
  46                 return expt(x,y);
  47         }
  48         if (y==0) { return 1; } // y=0 -> Ergebnis 1
  49         var uintL abs_y = (y<0 ? (uintL)(-y) : y); // Betrag von y nehmen
  50         var cl_N z = expt_pos(x,abs_y); // (expt x (abs y))
  51         return (y<0 ? recip(z) : z); // evtl. noch Kehrwert nehmen
  52 }