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 }