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   mp *t;
  61   size_t s;
  62   mp *b;
  63   mp *ainv;
  64   mp *c, *r;
  65   size_t i, j;
  66   mp *dd, *mone;
  67
  68   /* --- Cope if %$a \not\in Q_p$% --- */
  69
  70   if (mp_jacobi(a, p) != 1) {
  71     mp_drop(d);
  72     return (0);
  73   }
  74
  75   /* --- Choose some quadratic non-residue --- */
  76
  77   {
  78     grand *g = fibrand_create(0);
  79
  80     b = MP_NEW;
  81     do
  82       b = mprand_range(b, p, g, 0);
  83     while (mp_jacobi(b, p) != -1);
  84     g->ops->destroy(g);
  85   }
  86
  87   /* --- Find the inverse of %$a$% --- */
  88
  89   ainv = mp_modinv(MP_NEW, a, p);
  90
  91   /* --- Split %$p - 1$% into a power of two and an odd number --- */
  92
  93   t = mp_sub(MP_NEW, p, MP_ONE);
  94   t = mp_odd(t, t, &s);
  95
  96   /* --- Now to really get going --- */
  97
  98   mpmont_create(&mm, p);
  99   b = mpmont_mul(&mm, b, b, mm.r2);
 100   c = mpmont_expr(&mm, b, b, t);
 101   t = mp_add(t, t, MP_ONE);
 102   t = mp_lsr(t, t, 1);
 103   dd = mpmont_mul(&mm, MP_NEW, a, mm.r2);
 104   r = mpmont_expr(&mm, t, dd, t);
 105   mp_drop(dd);
 106   ainv = mpmont_mul(&mm, ainv, ainv, mm.r2);
 107
 108   mone = mp_sub(MP_NEW, p, mm.r);
 109
 110   dd = MP_NEW;
 111
 112   for (i = 1; i < s; i++) {
 113
 114     /* --- Compute %$d_0 = r^2a^{-1}$% --- */
 115
 116     dd = mp_sqr(dd, r);
 117     dd = mpmont_reduce(&mm, dd, dd);
 118     dd = mpmont_mul(&mm, dd, dd, ainv);
 119
 120     /* --- Now %$d = d_0^{2^{s - i - 1}}$% --- */
 121
 122     for (j = i; j < s - 1; j++) {
 123       dd = mp_sqr(dd, dd);
 124       dd = mpmont_reduce(&mm, dd, dd);
 125     }
 126
 127     /* --- Fiddle at the end --- */
 128
 129     if (MP_EQ(dd, mone))
 130       r = mpmont_mul(&mm, r, r, c);
 131     c = mp_sqr(c, c);
 132     c = mpmont_reduce(&mm, c, c);
 133   }
 134
 135   /* --- Done, so tidy up --- *
 136    *
 137    * Canonify the answer.
 138    */
 139
 140   d = mpmont_reduce(&mm, d, r);
 141   r = mp_sub(r, p, d);
 142   if (MP_CMP(r, <, d)) { mp *tt = r; r = d; d = tt; }
 143   mp_drop(ainv);
 144   mp_drop(r); mp_drop(c);
 145   mp_drop(dd);
 146   mp_drop(mone);
 147   mpmont_destroy(&mm);
 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 -------------------------------------------------*/