libsail: optimise real_power

The Arm spec uses the value 2.0^1000000 to represent infinity
so it is worth making real_power take logarithmic time.

1 files changed, 19 insertions, 12 deletions

diff --git a/lib/sail.c b/lib/sail.c
index 1fa3fa62..a5ddffe0 100644
--- a/lib/sail.c
+++ b/lib/sail.c
@@ -932,20 +932,27 @@ void real_power(real *rop, const real base, const sail_int exp)
 {
   int64_t exp_si = mpz_get_si(exp);
 
-  if (exp_si == 0) {
-    mpz_set_ui(mpq_numref(*rop), 1);
-    mpz_set_ui(mpq_denref(*rop), 1);
-  } else {
-    mpq_set(*rop, base);
-
-    for (int i = 1; i < abs(exp_si); i++) {
-      mpq_mul(*rop, *rop, base);
-    }
-
-    if (exp_si < 0) {
-      mpq_inv(*rop, *rop);
+  mpz_set_ui(mpq_numref(*rop), 1);
+  mpz_set_ui(mpq_denref(*rop), 1);
+
+  real b;
+  mpq_init(b);
+  mpq_set(b, base);
+  int64_t pexp = llabs(exp_si);
+  while (pexp != 0) {
+    // invariant: rop * b^pexp == base^abs(exp)
+    if (pexp & 1) { // b^(e+1) = b * b^e
+      mpq_mul(*rop, *rop, b);
+      pexp -= 1;
+    } else { // b^(2e) = (b*b)^e
+      mpq_mul(b, b, b);
+      pexp >>= 1;
     }
   }
+  if (exp_si < 0) {
+    mpq_inv(*rop, *rop);
+  }
+  mpq_clear(b);
 }
 
 void CREATE_OF(real, sail_string)(real *rop, const sail_string op)