+/* --- @f25519_quosqrt@ --- *
+ *
+ * Arguments: @f25519 *z@ = where to put the result (may alias @x@ or @y@)
+ * @const f25519 *x, *y@ = two operands
+ *
+ * Returns: Zero if successful, @-1@ if %$x/y$% is not a square.
+ *
+ * Use: Stores in @z@ the one of the square roots %$\pm\sqrt{x/y}$%.
+ * If %$x = y = 0% then the result is zero; if %$y = 0$% but %$x
+ * \ne 0$% then the operation fails. If you wanted a specific
+ * square root then you'll have to pick it yourself.
+ */
+
+static const piece sqrtm1_pieces[NPIECE] = {
+#if F25519_IMPL == 26
+ -32595792, -7943725, 9377950, 3500415, 12389472,
+ -272473, -25146209, -2005654, 326686, 11406482
+#elif F25519_IMPL == 10
+ 176, -88, 161, 157, -485, -196, -231, -220, -416,
+ -169, -255, 50, 189, -89, -266, -32, 202, -511,
+ 423, 357, 248, -249, 80, 288, 50, 174
+#endif
+};
+#define SQRTM1 ((const f25519 *)sqrtm1_pieces)
+
+int f25519_quosqrt(f25519 *z, const f25519 *x, const f25519 *y)
+{
+ f25519 t, u, v, w, t15;
+ octet xb[32], b0[32], b1[32];
+ int32 rc = -1;
+ mask32 m;
+ unsigned i;
+
+#define SQRN(z, x, n) do { \
+ f25519_sqr((z), (x)); \
+ for (i = 1; i < (n); i++) f25519_sqr((z), (z)); \
+} while (0)
+
+ /* This is a bit tricky; the algorithm is loosely based on Bernstein, Duif,
+ * Lange, Schwabe, and Yang, `High-speed high-security signatures',
+ * 2011-09-26, https://ed25519.cr.yp.to/ed25519-20110926.pdf.
+ */
+ f25519_mul(&v, x, y);
+
+ /* Now for an addition chain. */ /* step | value */
+ f25519_sqr(&u, &v); /* 1 | 2 */
+ f25519_mul(&t, &u, &v); /* 2 | 3 */
+ SQRN(&u, &t, 2); /* 4 | 12 */
+ f25519_mul(&t15, &u, &t); /* 5 | 15 */
+ f25519_sqr(&u, &t15); /* 6 | 30 */
+ f25519_mul(&t, &u, &v); /* 7 | 31 = 2^5 - 1 */
+ SQRN(&u, &t, 5); /* 12 | 2^10 - 2^5 */
+ f25519_mul(&t, &u, &t); /* 13 | 2^10 - 1 */
+ SQRN(&u, &t, 10); /* 23 | 2^20 - 2^10 */
+ f25519_mul(&u, &u, &t); /* 24 | 2^20 - 1 */
+ SQRN(&u, &u, 10); /* 34 | 2^30 - 2^10 */
+ f25519_mul(&t, &u, &t); /* 35 | 2^30 - 1 */
+ f25519_sqr(&u, &t); /* 36 | 2^31 - 2 */
+ f25519_mul(&t, &u, &v); /* 37 | 2^31 - 1 */
+ SQRN(&u, &t, 31); /* 68 | 2^62 - 2^31 */
+ f25519_mul(&t, &u, &t); /* 69 | 2^62 - 1 */
+ SQRN(&u, &t, 62); /* 131 | 2^124 - 2^62 */
+ f25519_mul(&t, &u, &t); /* 132 | 2^124 - 1 */
+ SQRN(&u, &t, 124); /* 256 | 2^248 - 2^124 */
+ f25519_mul(&t, &u, &t); /* 257 | 2^248 - 1 */
+ f25519_sqr(&u, &t); /* 258 | 2^249 - 2 */
+ f25519_mul(&t, &u, &v); /* 259 | 2^249 - 1 */
+ SQRN(&t, &t, 3); /* 262 | 2^252 - 8 */
+ f25519_sqr(&u, &t); /* 263 | 2^253 - 16 */
+ f25519_mul(&t, &u, &t); /* 264 | 3*2^252 - 24 */
+ f25519_mul(&t, &t, &t15); /* 265 | 3*2^252 - 9 */
+ f25519_mul(&w, &t, &v); /* 266 | 3*2^252 - 8 */
+
+ /* Awesome. Now let me explain. Let v be a square in GF(p), and let w =
+ * v^(3*2^252 - 8). In particular, let's consider
+ *
+ * v^2 w^4 = v^2 v^{3*2^254 - 32} = (v^{2^254 - 10})^3
+ *
+ * But 2^254 - 10 = ((2^255 - 19) - 1)/2 = (p - 1)/2. Since v is a square,
+ * it has order dividing (p - 1)/2, and therefore v^2 w^4 = 1 and
+ *
+ * w^4 = 1/v^2
+ *
+ * That in turn implies that w^2 = ±1/v. Now, recall that v = x y, and let
+ * w' = w x. Then w'^2 = ±x^2/v = ±x/y. If y w'^2 = x then we set
+ * z = w', since we have z^2 = x/y; otherwise let z = i w', where i^2 = -1,
+ * so z^2 = -w^2 = x/y, and we're done.
+ *
+ * The easiest way to compare is to encode. This isn't as wasteful as it
+ * sounds: the hard part is normalizing the representations, which we have
+ * to do anyway.
+ */
+ f25519_mul(&w, &w, x);
+ f25519_sqr(&t, &w);
+ f25519_mul(&t, &t, y);
+ f25519_neg(&u, &t);
+ f25519_store(xb, x);
+ f25519_store(b0, &t);
+ f25519_store(b1, &u);
+ f25519_mul(&u, &w, SQRTM1);
+
+ m = -ct_memeq(b0, xb, 32);
+ rc = PICK2(0, rc, m);
+ f25519_pick2(z, &w, &u, m);
+ m = -ct_memeq(b1, xb, 32);
+ rc = PICK2(0, rc, m);
+
+ /* And we're done. */
+ return (rc);
+}
+