+/* --- @fgoldi_quosqrt@ --- *
+ *
+ * Arguments: @fgoldi *z@ = where to put the result (may alias @x@ or @y@)
+ * @const fgoldi *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.
+ */
+
+int fgoldi_quosqrt(fgoldi *z, const fgoldi *x, const fgoldi *y)
+{
+ fgoldi t, u, v;
+ octet xb[56], b0[56];
+ int32 rc = -1;
+ mask32 m;
+ unsigned i;
+
+#define SQRN(z, x, n) do { \
+ fgoldi_sqr((z), (x)); \
+ for (i = 1; i < (n); i++) fgoldi_sqr((z), (z)); \
+} while (0)
+
+ /* This is, fortunately, significantly easier than the equivalent problem
+ * in GF(2^255 - 19), since p == 3 (mod 4).
+ *
+ * If x/y is square, then so is v = y^2 x/y = x y, and therefore u has
+ * order r = (p - 1)/2. Let w = v^{(p-3)/4}. Then w^2 = v^{(p-3)/2} =
+ * u^{r-1} = 1/v = 1/x y. Clearly, then, (x w)^2 = x^2/x y = x/y, so x w
+ * is one of the square roots we seek.
+ *
+ * The addition chain, then, is a prefix of the previous one.
+ */
+ fgoldi_mul(&v, x, y);
+
+ fgoldi_sqr(&u, &v); /* 1 | 2 */
+ fgoldi_mul(&t, &u, &v); /* 2 | 3 */
+ SQRN(&u, &t, 2); /* 4 | 12 */
+ fgoldi_mul(&t, &u, &t); /* 5 | 15 */
+ SQRN(&u, &t, 4); /* 9 | 240 */
+ fgoldi_mul(&u, &u, &t); /* 10 | 255 = 2^8 - 1 */
+ SQRN(&u, &u, 4); /* 14 | 2^12 - 16 */
+ fgoldi_mul(&t, &u, &t); /* 15 | 2^12 - 1 */
+ SQRN(&u, &t, 12); /* 27 | 2^24 - 2^12 */
+ fgoldi_mul(&u, &u, &t); /* 28 | 2^24 - 1 */
+ SQRN(&u, &u, 12); /* 40 | 2^36 - 2^12 */
+ fgoldi_mul(&t, &u, &t); /* 41 | 2^36 - 1 */
+ fgoldi_sqr(&t, &t); /* 42 | 2^37 - 2 */
+ fgoldi_mul(&t, &t, &v); /* 43 | 2^37 - 1 */
+ SQRN(&u, &t, 37); /* 80 | 2^74 - 2^37 */
+ fgoldi_mul(&u, &u, &t); /* 81 | 2^74 - 1 */
+ SQRN(&u, &u, 37); /* 118 | 2^111 - 2^37 */
+ fgoldi_mul(&t, &u, &t); /* 119 | 2^111 - 1 */
+ SQRN(&u, &t, 111); /* 230 | 2^222 - 2^111 */
+ fgoldi_mul(&t, &u, &t); /* 231 | 2^222 - 1 */
+ fgoldi_sqr(&u, &t); /* 232 | 2^223 - 2 */
+ fgoldi_mul(&u, &u, &v); /* 233 | 2^223 - 1 */
+ SQRN(&u, &u, 223); /* 456 | 2^446 - 2^223 */
+ fgoldi_mul(&t, &u, &t); /* 457 | 2^446 - 2^222 - 1 */
+
+#undef SQRN
+
+ /* Now we must decide whether the answer was right. We should have z^2 =
+ * x/y, so y z^2 = x.
+ *
+ * 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.
+ */
+ fgoldi_mul(z, x, &t);
+ fgoldi_sqr(&t, z);
+ fgoldi_mul(&t, &t, y);
+ fgoldi_store(xb, x);
+ fgoldi_store(b0, &t);
+ m = -ct_memeq(xb, b0, 56);
+ rc = PICK2(0, rc, m);
+ return (rc);
+}
+