Initial revision

[cln.git] / src / real / misc / cl_R_expt_I.cc

diff --git a/src/real/misc/cl_R_expt_I.cc b/src/real/misc/cl_R_expt_I.cc
new file mode 100644 (file)
index 0000000..93db843
--- /dev/null
+++ b/src/real/misc/cl_R_expt_I.cc
@@ -0,0 +1,56 @@
+// expt().
+
+// General includes.
+#include "cl_sysdep.h"
+
+// Specification.
+#include "cl_real.h"
+
+
+// Implementation.
+
+#include "cl_R.h"
+#include "cl_rational.h"
+#include "cl_integer.h"
+#include "cl_I.h"
+
+// Methode:
+// Für y>0:
+//   a:=x, b:=y.
+//   Solange b gerade, setze a:=a*a, b:=b/2. [a^b bleibt invariant, = x^y.]
+//   c:=a.
+//   Solange b:=floor(b/2) >0 ist,
+//     setze a:=a*a, und falls b ungerade, setze c:=a*c.
+//   Ergebnis c.
+// Für y=0: Ergebnis 1.
+// Für y<0: (/ (expt x (- y))).
+
+// Assume y>0.
+inline const cl_R expt_pos (const cl_R& x, const cl_I& y)
+{
+       if (rationalp(x)) {
+               DeclareType(cl_RA,x);
+               return expt(x,y); // x rational -> schnellere Routine
+       } else {
+               DeclareType(cl_F,x);
+               var cl_F a = x;
+               var cl_I b = y;
+               while (!oddp(b)) { a = square(a); b = b >> 1; }
+               var cl_F c = a;
+               until (eq(b,1))
+                 { b = b >> 1;
+                   a = square(a);
+                   if (oddp(b)) { c = a * c; }
+                 }
+               return c;
+       }
+}
+
+const cl_R expt (const cl_R& x, const cl_I& y)
+{
+       if (eq(y,0)) { return 1; } // y=0 -> Ergebnis 1
+       var cl_boolean y_negative = minusp(y);
+       var cl_I abs_y = (y_negative ? -y : y); // Betrag von y nehmen
+       var cl_R z = expt_pos(x,abs_y); // (expt x (abs y))
+       return (y_negative ? recip(z) : z); // evtl. noch Kehrwert nehmen
+}