1 /* $OpenBSD: bn_isqrt.c,v 1.2 2022/07/13 11:20:00 tb Exp $ */ 2 /* 3 * Copyright (c) 2022 Theo Buehler <tb@openbsd.org> 4 * 5 * Permission to use, copy, modify, and distribute this software for any 6 * purpose with or without fee is hereby granted, provided that the above 7 * copyright notice and this permission notice appear in all copies. 8 * 9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 11 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 12 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 13 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 14 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 15 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 16 */ 17 18 #include <stddef.h> 19 #include <stdint.h> 20 21 #include <openssl/bn.h> 22 #include <openssl/err.h> 23 24 #include "bn_lcl.h" 25 26 #define CTASSERT(x) extern char _ctassert[(x) ? 1 : -1 ] \ 27 __attribute__((__unused__)) 28 29 /* 30 * Calculate integer square root of |n| using a variant of Newton's method. 31 * 32 * Returns the integer square root of |n| in the caller-provided |out_sqrt|; 33 * |*out_perfect| is set to 1 if and only if |n| is a perfect square. 34 * One of |out_sqrt| and |out_perfect| can be NULL; |in_ctx| can be NULL. 35 * 36 * Returns 0 on error, 1 on success. 37 * 38 * Adapted from pure Python describing cpython's math.isqrt(), without bothering 39 * with any of the optimizations in the C code. A correctness proof is here: 40 * https://github.com/mdickinson/snippets/blob/master/proofs/isqrt/src/isqrt.lean 41 * The comments in the Python code also give a rather detailed proof. 42 */ 43 44 int 45 bn_isqrt(BIGNUM *out_sqrt, int *out_perfect, const BIGNUM *n, BN_CTX *in_ctx) 46 { 47 BN_CTX *ctx = NULL; 48 BIGNUM *a, *b; 49 int c, d, e, s; 50 int cmp, perfect; 51 int ret = 0; 52 53 if (out_perfect == NULL && out_sqrt == NULL) { 54 BNerror(ERR_R_PASSED_NULL_PARAMETER); 55 goto err; 56 } 57 58 if (BN_is_negative(n)) { 59 BNerror(BN_R_INVALID_RANGE); 60 goto err; 61 } 62 63 if ((ctx = in_ctx) == NULL) 64 ctx = BN_CTX_new(); 65 if (ctx == NULL) 66 goto err; 67 68 BN_CTX_start(ctx); 69 70 if ((a = BN_CTX_get(ctx)) == NULL) 71 goto err; 72 if ((b = BN_CTX_get(ctx)) == NULL) 73 goto err; 74 75 if (BN_is_zero(n)) { 76 perfect = 1; 77 if (!BN_zero(a)) 78 goto err; 79 goto done; 80 } 81 82 if (!BN_one(a)) 83 goto err; 84 85 c = (BN_num_bits(n) - 1) / 2; 86 d = 0; 87 88 /* Calculate s = floor(log(c)). */ 89 if (!BN_set_word(b, c)) 90 goto err; 91 s = BN_num_bits(b) - 1; 92 93 /* 94 * By definition, the loop below is run <= floor(log(log(n))) times. 95 * Comments in the cpython code establish the loop invariant that 96 * 97 * (a - 1)^2 < n / 4^(c - d) < (a + 1)^2 98 * 99 * holds true in every iteration. Once this is proved via induction, 100 * correctness of the algorithm is easy. 101 * 102 * Roughly speaking, A = (a << (d - e)) is used for one Newton step 103 * "a = (A >> 1) + (m >> 1) / A" approximating m = (n >> 2 * (c - d)). 104 */ 105 106 for (; s >= 0; s--) { 107 e = d; 108 d = c >> s; 109 110 if (!BN_rshift(b, n, 2 * c - d - e + 1)) 111 goto err; 112 113 if (!BN_div_ct(b, NULL, b, a, ctx)) 114 goto err; 115 116 if (!BN_lshift(a, a, d - e - 1)) 117 goto err; 118 119 if (!BN_add(a, a, b)) 120 goto err; 121 } 122 123 /* 124 * The loop invariant implies that either a or a - 1 is isqrt(n). 125 * Figure out which one it is. The invariant also implies that for 126 * a perfect square n, a must be the square root. 127 */ 128 129 if (!BN_sqr(b, a, ctx)) 130 goto err; 131 132 /* If a^2 > n, we must have isqrt(n) == a - 1. */ 133 if ((cmp = BN_cmp(b, n)) > 0) { 134 if (!BN_sub_word(a, 1)) 135 goto err; 136 } 137 138 perfect = cmp == 0; 139 140 done: 141 if (out_perfect != NULL) 142 *out_perfect = perfect; 143 144 if (out_sqrt != NULL) { 145 if (!BN_copy(out_sqrt, a)) 146 goto err; 147 } 148 149 ret = 1; 150 151 err: 152 BN_CTX_end(ctx); 153 154 if (ctx != in_ctx) 155 BN_CTX_free(ctx); 156 157 return ret; 158 } 159 160 /* 161 * is_square_mod_N[r % N] indicates whether r % N has a square root modulo N. 162 * The tables are generated in regress/lib/libcrypto/bn/bn_isqrt.c. 163 */ 164 165 const uint8_t is_square_mod_11[] = { 166 1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 167 }; 168 CTASSERT(sizeof(is_square_mod_11) == 11); 169 170 const uint8_t is_square_mod_63[] = { 171 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 172 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 173 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 174 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 175 }; 176 CTASSERT(sizeof(is_square_mod_63) == 63); 177 178 const uint8_t is_square_mod_64[] = { 179 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 180 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 181 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 182 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 183 }; 184 CTASSERT(sizeof(is_square_mod_64) == 64); 185 186 const uint8_t is_square_mod_65[] = { 187 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 188 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 189 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 190 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 191 1, 192 }; 193 CTASSERT(sizeof(is_square_mod_65) == 65); 194 195 /* 196 * Determine whether n is a perfect square or not. 197 * 198 * Returns 1 on success and 0 on error. In case of success, |*out_perfect| is 199 * set to 1 if and only if |n| is a perfect square. 200 */ 201 202 int 203 bn_is_perfect_square(int *out_perfect, const BIGNUM *n, BN_CTX *ctx) 204 { 205 BN_ULONG r; 206 207 *out_perfect = 0; 208 209 if (BN_is_negative(n)) 210 return 1; 211 212 /* 213 * Before performing an expensive bn_isqrt() operation, weed out many 214 * obvious non-squares. See H. Cohen, "A course in computational 215 * algebraic number theory", Algorithm 1.7.3. 216 * 217 * The idea is that a square remains a square when reduced modulo any 218 * number. The moduli are chosen in such a way that a non-square has 219 * probability < 1% of passing the four table lookups. 220 */ 221 222 /* n % 64 */ 223 r = BN_lsw(n) & 0x3f; 224 225 if (!is_square_mod_64[r % 64]) 226 return 1; 227 228 if ((r = BN_mod_word(n, 11 * 63 * 65)) == (BN_ULONG)-1) 229 return 0; 230 231 if (!is_square_mod_63[r % 63] || 232 !is_square_mod_65[r % 65] || 233 !is_square_mod_11[r % 11]) 234 return 1; 235 236 return bn_isqrt(NULL, out_perfect, n, ctx); 237 } 238