zpow.c - libzahl - big integer library
 (HTM) git clone git://git.suckless.org/libzahl
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
       zpow.c (1157B)
       ---
            1 /* See LICENSE file for copyright and license details. */
            2 #include "internals.h"
            3 
            4 #define tb  libzahl_tmp_pow_b
            5 #define tc  libzahl_tmp_pow_c
            6 
            7 
            8 void
            9 zpow(z_t a, z_t b, z_t c)
           10 {
           11         /*
           12          * Exponentiation by squaring.
           13          * 
           14          * 7↑19 = 7↑10011₂ = 7↑2⁰ ⋅ 7↑2¹ ⋅ 7↑2⁴ where a↑2↑(n + 1) = (a↑2↑n)².
           15          */
           16 
           17         /* TODO use zpowu when possible */
           18 
           19         size_t i, j, n, bits;
           20         zahl_char_t x;
           21         int neg;
           22 
           23         if (unlikely(zsignum(c) <= 0)) {
           24                 if (zzero(c)) {
           25                         if (check(zzero(b)))
           26                                 libzahl_failure(-ZERROR_0_POW_0);
           27                         zsetu(a, 1);
           28                 } else if (check(zzero(b))) {
           29                         libzahl_failure(-ZERROR_DIV_0);
           30                 } else {
           31                         SET_SIGNUM(a, 0);
           32                 }
           33                 return;
           34         } else if (unlikely(zzero(b))) {
           35                 SET_SIGNUM(a, 0);
           36                 return;
           37         }
           38 
           39         bits = zbits(c);
           40         n = FLOOR_BITS_TO_CHARS(bits);
           41 
           42         neg = znegative(b) && zodd(c);
           43         zabs(tb, b);
           44         zset(tc, c);
           45         zsetu(a, 1);
           46 
           47         for (i = 0; i < n; i++) { /* Remember, n is floored. */
           48                 x = tc->chars[i];
           49                 for (j = BITS_PER_CHAR; j--; x >>= 1) {
           50                         if (x & 1)
           51                                 zmul_ll(a, a, tb);
           52                         zsqr_ll(tb, tb);
           53                 }
           54         }
           55         x = tc->chars[i];
           56         for (; x; x >>= 1) {
           57                 if (x & 1)
           58                         zmul_ll(a, a, tb);
           59                 zsqr_ll(tb, tb);
           60         }
           61 
           62         if (neg)
           63                 zneg(a, a);
           64 }