chiark - git - mdw - catacomb/blob - math/mp-modsqrt.c

   1 /* -*-c-*-
   2  *
   3  * Compute square roots modulo a prime
   4  *
   5  * (c) 2000 Straylight/Edgeware
   6  */
   7
   8 /*----- Licensing notice --------------------------------------------------*
   9  *
  10  * This file is part of Catacomb.
  11  *
  12  * Catacomb is free software; you can redistribute it and/or modify
  13  * it under the terms of the GNU Library General Public License as
  14  * published by the Free Software Foundation; either version 2 of the
  15  * License, or (at your option) any later version.
  16  *
  17  * Catacomb is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20  * GNU Library General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Library General Public
  23  * License along with Catacomb; if not, write to the Free
  24  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  25  * MA 02111-1307, USA.
  26  */
  27
  28 /*----- Header files ------------------------------------------------------*/
  29
  30 #include "fibrand.h"
  31 #include "grand.h"
  32 #include "mp.h"
  33 #include "mpmont.h"
  34 #include "mprand.h"
  35
  36 /*----- Main code ---------------------------------------------------------*/
  37
  38 /* --- @mp_modsqrt@ --- *
  39  *
  40  * Arguments:   @mp *d@ = destination integer
  41  *              @mp *a@ = source integer
  42  *              @mp *p@ = modulus (must be prime)
  43  *
  44  * Returns:     If %$a$% is a quadratic residue, a square root of %$a$%; else
  45  *              a null pointer.
  46  *
  47  * Use:         Returns an integer %$x$% such that %$x^2 \equiv a \pmod{p}$%,
  48  *              if one exists; else a null pointer.  This function will not
  49  *              work if %$p$% is composite: you must factor the modulus, take
  50  *              a square root mod each factor, and recombine the results
  51  *              using the Chinese Remainder Theorem.
  52  *
  53  *              We guarantee that the square root returned is the smallest
  54  *              one (i.e., the `positive' square root).
  55  */
  56
  57 mp *mp_modsqrt(mp *d, mp *a, mp *p)
  58 {
  59   mpmont mm;
  60   size_t i, s;
  61   mp *b, *c;
  62   mp *ainv;
  63   mp *r, *A, *aa;
  64   mp *t;
  65   grand *gr;
  66   int j;
  67
  68   /* --- Cope if %$a \not\in Q_p$% --- */
  69
  70   j = mp_jacobi(a, p);
  71   if (j == -1) {
  72     mp_drop(d);
  73     return (0);
  74   } else if (j == 0) {
  75     if (d != a) mp_drop(d);
  76     d = MP_COPY(a);
  77     return (d);
  78   }
  79
  80   /* --- Choose some quadratic non-residue --- */
  81
  82   gr = fibrand_create(0);
  83   b = MP_NEW;
  84   do b = mprand_range(b, p, gr, 0); while (mp_jacobi(b, p) != -1);
  85   gr->ops->destroy(gr);
  86
  87   /* --- Some initial setup --- */
  88
  89   mpmont_create(&mm, p);
  90   ainv = mp_modinv(MP_NEW, a, p);       /* %$a^{-1} \bmod p$% */
  91   ainv = mpmont_mul(&mm, ainv, ainv, mm.r2);
  92   t = mp_sub(MP_NEW, p, MP_ONE);
  93   t = mp_odd(t, t, &s);                 /* %$2^s t = p - 1$% */
  94   b = mpmont_mul(&mm, b, b, mm.r2);
  95   c = mpmont_expr(&mm, b, b, t);        /* %$b^t \bmod p$% */
  96   t = mp_add(t, t, MP_ONE);
  97   t = mp_lsr(t, t, 1);                  /* %$(t + 1)/2$% */
  98   a = mpmont_mul(&mm, MP_NEW, a, mm.r2);
  99   r = mpmont_expr(&mm, a, a, t);        /* %$a^{(t+1)/2} \bmod p$% */
 100
 101   /* --- Now for the main loop --- *
 102    *
 103    * Let %$g = c^{-2}$%; we know that %$g$% is a generator of the order-
 104    * %$2^{s-1}$% subgroup mod %$p$%.  We also know that %$A = a^t = r^2/a$%
 105    * is an element of this group.  If we can determine %$m$% such that
 106    * %$g^m = A$% then %$a^{(t+1)/2}/g^{m/2} = r c^m$% is the square root we
 107    * seek.
 108    *
 109    * Write %$m = m_0 + 2 m'$%.  Then %$A^{2^{s-1}} = g^{m_0 2^{s-1}}$%, which
 110    * is %$1$% if %$m_0 = 0$% or %$-1$% if %$m_0 = 1$% (modulo %$p$%).  Then
 111    * %$A/g^{m_0} = (g^2)^{m'}$% and we can proceed inductively.  The end
 112    * result will me %$A/g^m$%.
 113    *
 114    * Note that this loop keeps track of (what will be) %$r c^m$%, since this
 115    * is the result we want, and computes $A/g^m = r^2/a$% on demand.
 116    */
 117
 118   A = mp_sqr(t, r); A = mpmont_reduce(&mm, A, A);
 119   A = mpmont_mul(&mm, A, A, ainv);      /* %$x^t/g^m$% */
 120
 121   while (s-- > 1) {
 122     aa = MP_COPY(A);
 123     for (i = 1; i < s; i++)
 124       { aa = mp_sqr(aa, aa); aa = mpmont_reduce(&mm, aa, aa); }
 125     if (!MP_EQ(aa, mm.r)) {
 126       r = mpmont_mul(&mm, r, r, c);
 127       A = mp_sqr(A, r); A = mpmont_reduce(&mm, A, A);
 128       A = mpmont_mul(&mm, A, A, ainv);  /* %$x^t/g^m$% */
 129     }
 130     c = mp_sqr(c, c); c = mpmont_reduce(&mm, c, c);
 131     MP_DROP(aa);
 132   }
 133
 134   /* --- We want the smaller square root --- */
 135
 136   d = mpmont_reduce(&mm, d, r);
 137   r = mp_sub(r, p, d);
 138   if (MP_CMP(r, <, d)) { mp *tt = r; r = d; d = tt; }
 139
 140   /* --- Clear away all the temporaries --- */
 141
 142   mp_drop(ainv);
 143   mp_drop(r); mp_drop(c);
 144   mp_drop(A);
 145   mpmont_destroy(&mm);
 146
 147   /* --- Done --- */
 148
 149   return (d);
 150 }
 151
 152 /*----- Test rig ----------------------------------------------------------*/
 153
 154 #ifdef TEST_RIG
 155
 156 #include <mLib/testrig.h>
 157
 158 static int verify(dstr *v)
 159 {
 160   mp *a = *(mp **)v[0].buf;
 161   mp *p = *(mp **)v[1].buf;
 162   mp *rr = *(mp **)v[2].buf;
 163   mp *r = mp_modsqrt(MP_NEW, a, p);
 164   int ok = 0;
 165
 166   if (!r)
 167     ok = 0;
 168   else if (MP_EQ(r, rr))
 169     ok = 1;
 170
 171   if (!ok) {
 172     fputs("\n*** fail\n", stderr);
 173     fputs("a  = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr);
 174     fputs("p  = ", stderr); mp_writefile(p, stderr, 10); fputc('\n', stderr);
 175     if (r) {
 176       fputs("r  = ", stderr);
 177       mp_writefile(r, stderr, 10);
 178       fputc('\n', stderr);
 179     } else
 180       fputs("r  = <undef>\n", stderr);
 181     fputs("rr = ", stderr); mp_writefile(rr, stderr, 10); fputc('\n', stderr);
 182     ok = 0;
 183   }
 184
 185   mp_drop(a);
 186   mp_drop(p);
 187   mp_drop(r);
 188   mp_drop(rr);
 189   assert(mparena_count(MPARENA_GLOBAL) == 0);
 190   return (ok);
 191 }
 192
 193 static test_chunk tests[] = {
 194   { "modsqrt", verify, { &type_mp, &type_mp, &type_mp, 0 } },
 195   { 0, 0, { 0 } }
 196 };
 197
 198 int main(int argc, char *argv[])
 199 {
 200   sub_init();
 201   test_run(argc, argv, tests, SRCDIR "/t/mp");
 202   return (0);
 203 }
 204
 205 #endif
 206
 207 /*----- That's all, folks -------------------------------------------------*/