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
  67   /* --- Cope if %$a \not\in Q_p$% --- */
  68
  69   if (mp_jacobi(a, p) != 1) {
  70     mp_drop(d);
  71     return (0);
  72   }
  73
  74   /* --- Choose some quadratic non-residue --- */
  75
  76   gr = fibrand_create(0);
  77   b = MP_NEW;
  78   do b = mprand_range(b, p, gr, 0); while (mp_jacobi(b, p) != -1);
  79   gr->ops->destroy(gr);
  80
  81   /* --- Some initial setup --- */
  82
  83   mpmont_create(&mm, p);
  84   ainv = mp_modinv(MP_NEW, a, p);       /* %$a^{-1} \bmod p$% */
  85   ainv = mpmont_mul(&mm, ainv, ainv, mm.r2);
  86   t = mp_sub(MP_NEW, p, MP_ONE);
  87   t = mp_odd(t, t, &s);                 /* %$2^s t = p - 1$% */
  88   b = mpmont_mul(&mm, b, b, mm.r2);
  89   c = mpmont_expr(&mm, b, b, t);        /* %$b^t \bmod p$% */
  90   t = mp_add(t, t, MP_ONE);
  91   t = mp_lsr(t, t, 1);                  /* %$(t + 1)/2$% */
  92   a = mpmont_mul(&mm, MP_NEW, a, mm.r2);
  93   r = mpmont_expr(&mm, a, a, t);        /* %$a^{(t+1)/2} \bmod p$% */
  94
  95   /* --- Now for the main loop --- *
  96    *
  97    * Let %$g = c^{-2}$%; we know that %$g$% is a generator of the order-
  98    * %$2^{s-1}$% subgroup mod %$p$%.  We also know that %$A = a^t = r^2/a$%
  99    * is an element of this group.  If we can determine %$m$% such that
 100    * %$g^m = A$% then %$a^{(t+1)/2}/g^{m/2} = r c^m$% is the square root we
 101    * seek.
 102    *
 103    * Write %$m = m_0 + 2 m'$%.  Then %$A^{2^{s-1}} = g^{m_0 2^{s-1}}$%, which
 104    * is %$1$% if %$m_0 = 0$% or %$-1$% if %$m_0 = 1$% (modulo %$p$%).  Then
 105    * %$A/g^{m_0} = (g^2)^{m'}$% and we can proceed inductively.  The end
 106    * result will me %$A/g^m$%.
 107    *
 108    * Note that this loop keeps track of (what will be) %$r c^m$%, since this
 109    * is the result we want, and computes $A/g^m = r^2/a$% on demand.
 110    */
 111
 112   A = mp_sqr(t, r); A = mpmont_reduce(&mm, A, A);
 113   A = mpmont_mul(&mm, A, A, ainv);      /* %$x^t/g^m$% */
 114
 115   while (s-- > 1) {
 116     aa = MP_COPY(A);
 117     for (i = 1; i < s; i++)
 118       { aa = mp_sqr(aa, aa); aa = mpmont_reduce(&mm, aa, aa); }
 119     if (!MP_EQ(aa, mm.r)) {
 120       r = mpmont_mul(&mm, r, r, c);
 121       A = mp_sqr(A, r); A = mpmont_reduce(&mm, A, A);
 122       A = mpmont_mul(&mm, A, A, ainv);  /* %$x^t/g^m$% */
 123     }
 124     c = mp_sqr(c, c); c = mpmont_reduce(&mm, c, c);
 125     MP_DROP(aa);
 126   }
 127
 128   /* --- We want the smaller square root --- */
 129
 130   d = mpmont_reduce(&mm, d, r);
 131   r = mp_sub(r, p, d);
 132   if (MP_CMP(r, <, d)) { mp *tt = r; r = d; d = tt; }
 133
 134   /* --- Clear away all the temporaries --- */
 135
 136   mp_drop(ainv);
 137   mp_drop(r); mp_drop(c);
 138   mp_drop(A);
 139   mpmont_destroy(&mm);
 140
 141   /* --- Done --- */
 142
 143   return (d);
 144 }
 145
 146 /*----- Test rig ----------------------------------------------------------*/
 147
 148 #ifdef TEST_RIG
 149
 150 #include <mLib/testrig.h>
 151
 152 static int verify(dstr *v)
 153 {
 154   mp *a = *(mp **)v[0].buf;
 155   mp *p = *(mp **)v[1].buf;
 156   mp *rr = *(mp **)v[2].buf;
 157   mp *r = mp_modsqrt(MP_NEW, a, p);
 158   int ok = 0;
 159
 160   if (!r)
 161     ok = 0;
 162   else if (MP_EQ(r, rr))
 163     ok = 1;
 164
 165   if (!ok) {
 166     fputs("\n*** fail\n", stderr);
 167     fputs("a  = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr);
 168     fputs("p  = ", stderr); mp_writefile(p, stderr, 10); fputc('\n', stderr);
 169     if (r) {
 170       fputs("r  = ", stderr);
 171       mp_writefile(r, stderr, 10);
 172       fputc('\n', stderr);
 173     } else
 174       fputs("r  = <undef>\n", stderr);
 175     fputs("rr = ", stderr); mp_writefile(rr, stderr, 10); fputc('\n', stderr);
 176     ok = 0;
 177   }
 178
 179   mp_drop(a);
 180   mp_drop(p);
 181   mp_drop(r);
 182   mp_drop(rr);
 183   assert(mparena_count(MPARENA_GLOBAL) == 0);
 184   return (ok);
 185 }
 186
 187 static test_chunk tests[] = {
 188   { "modsqrt", verify, { &type_mp, &type_mp, &type_mp, 0 } },
 189   { 0, 0, { 0 } }
 190 };
 191
 192 int main(int argc, char *argv[])
 193 {
 194   sub_init();
 195   test_run(argc, argv, tests, SRCDIR "/t/mp");
 196   return (0);
 197 }
 198
 199 #endif
 200
 201 /*----- That's all, folks -------------------------------------------------*/