git.distorted.org.uk Git - u/mdw/catacomb/blob - mp-modsqrt.c

   1 /* -*-c-*-
   2  *
   3  * $Id: mp-modsqrt.c,v 1.4 2001/06/16 12:56:38 mdw Exp $
   4  *
   5  * Compute square roots modulo a prime
   6  *
   7  * (c) 2000 Straylight/Edgeware
   8  */
   9
  10 /*----- Licensing notice --------------------------------------------------*
  11  *
  12  * This file is part of Catacomb.
  13  *
  14  * Catacomb is free software; you can redistribute it and/or modify
  15  * it under the terms of the GNU Library General Public License as
  16  * published by the Free Software Foundation; either version 2 of the
  17  * License, or (at your option) any later version.
  18  *
  19  * Catacomb is distributed in the hope that it will be useful,
  20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22  * GNU Library General Public License for more details.
  23  *
  24  * You should have received a copy of the GNU Library General Public
  25  * License along with Catacomb; if not, write to the Free
  26  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  27  * MA 02111-1307, USA.
  28  */
  29
  30 /*----- Revision history --------------------------------------------------*
  31  *
  32  * $Log: mp-modsqrt.c,v $
  33  * Revision 1.4  2001/06/16 12:56:38  mdw
  34  * Fixes for interface change to @mpmont_expr@ and @mpmont_mexpr@.
  35  *
  36  * Revision 1.3  2001/02/03 12:00:29  mdw
  37  * Now @mp_drop@ checks its argument is non-NULL before attempting to free
  38  * it.  Note that the macro version @MP_DROP@ doesn't do this.
  39  *
  40  * Revision 1.2  2000/10/08 12:02:21  mdw
  41  * Use @MP_EQ@ instead of @MP_CMP@.
  42  *
  43  * Revision 1.1  2000/06/22 19:01:31  mdw
  44  * Compute square roots in a prime field.
  45  *
  46  */
  47
  48 /*----- Header files ------------------------------------------------------*/
  49
  50 #include "fibrand.h"
  51 #include "grand.h"
  52 #include "mp.h"
  53 #include "mpmont.h"
  54 #include "mprand.h"
  55
  56 /*----- Main code ---------------------------------------------------------*/
  57
  58 /* --- @mp_modsqrt@ --- *
  59  *
  60  * Arguments:   @mp *d@ = destination integer
  61  *              @mp *a@ = source integer
  62  *              @mp *p@ = modulus (must be prime)
  63  *
  64  * Returns:     If %$a$% is a quadratic residue, a square root of %$a$%; else
  65  *              a null pointer.
  66  *
  67  * Use:         Returns an integer %$x$% such that %$x^2 \equiv a \pmod{p}$%,
  68  *              if one exists; else a null pointer.  This function will not
  69  *              work if %$p$% is composite: you must factor the modulus, take
  70  *              a square root mod each factor, and recombine the results
  71  *              using the Chinese Remainder Theorem.
  72  */
  73
  74 mp *mp_modsqrt(mp *d, mp *a, mp *p)
  75 {
  76   mpmont mm;
  77   mp *t;
  78   size_t s;
  79   mp *b;
  80   mp *ainv;
  81   mp *c, *r;
  82   size_t i, j;
  83   mp *dd, *mone;
  84
  85   /* --- Cope if %$a \not\in Q_p$% --- */
  86
  87   if (mp_jacobi(a, p) != 1) {
  88     mp_drop(d);
  89     return (0);
  90   }
  91
  92   /* --- Choose some quadratic non-residue --- */
  93
  94   {
  95     grand *g = fibrand_create(0);
  96
  97     b = MP_NEW;
  98     do
  99       b = mprand_range(b, p, g, 0);
 100     while (mp_jacobi(b, p) != -1);
 101     g->ops->destroy(g);
 102   }
 103
 104   /* --- Find the inverse of %$a$% --- */
 105
 106   ainv = MP_NEW;
 107   mp_gcd(0, &ainv, 0, a, p);
 108
 109   /* --- Split %$p - 1$% into a power of two and an odd number --- */
 110
 111   t = mp_sub(MP_NEW, p, MP_ONE);
 112   t = mp_odd(t, t, &s);
 113
 114   /* --- Now to really get going --- */
 115
 116   mpmont_create(&mm, p);
 117   b = mpmont_mul(&mm, b, b, mm.r2);
 118   c = mpmont_expr(&mm, b, b, t);
 119   t = mp_add(t, t, MP_ONE);
 120   t = mp_lsr(t, t, 1);
 121   dd = mpmont_mul(&mm, MP_NEW, a, mm.r2);
 122   r = mpmont_expr(&mm, t, dd, t);
 123   mp_drop(dd);
 124   ainv = mpmont_mul(&mm, ainv, ainv, mm.r2);
 125
 126   mone = mp_sub(MP_NEW, p, mm.r);
 127
 128   dd = MP_NEW;
 129
 130   for (i = 1; i < s; i++) {
 131
 132     /* --- Compute %$d_0 = r^2a^{-1}$% --- */
 133
 134     dd = mp_sqr(dd, r);
 135     dd = mpmont_reduce(&mm, dd, dd);
 136     dd = mpmont_mul(&mm, dd, dd, ainv);
 137
 138     /* --- Now %$d = d_0^{s - i - 1}$% --- */
 139
 140     for (j = i; j < s - 1; j++) {
 141       dd = mp_sqr(dd, dd);
 142       dd = mpmont_reduce(&mm, dd, dd);
 143     }
 144
 145     /* --- Fiddle at the end --- */
 146
 147     if (MP_EQ(dd, mone))
 148       r = mpmont_mul(&mm, r, r, c);
 149     c = mp_sqr(c, c);
 150     c = mpmont_reduce(&mm, c, c);
 151   }
 152
 153   /* --- Done, so tidy up --- */
 154
 155   d = mpmont_reduce(&mm, d, r);
 156   mp_drop(ainv);
 157   mp_drop(r); mp_drop(c);
 158   mp_drop(dd);
 159   mp_drop(mone);
 160   mpmont_destroy(&mm);
 161
 162   return (d);
 163 }
 164
 165 /*----- Test rig ----------------------------------------------------------*/
 166
 167 #ifdef TEST_RIG
 168
 169 #include <mLib/testrig.h>
 170
 171 static int verify(dstr *v)
 172 {
 173   mp *a = *(mp **)v[0].buf;
 174   mp *p = *(mp **)v[1].buf;
 175   mp *rr = *(mp **)v[2].buf;
 176   mp *r = mp_modsqrt(MP_NEW, a, p);
 177   int ok = 0;
 178
 179   if (!r)
 180     ok = 0;
 181   else if (MP_EQ(r, rr))
 182     ok = 1;
 183   else {
 184     r = mp_sub(r, p, r);
 185     if (MP_EQ(r, rr))
 186       ok = 1;
 187   }
 188
 189   if (!ok) {
 190     fputs("\n*** fail\n", stderr);
 191     fputs("a  = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr);
 192     fputs("p  = ", stderr); mp_writefile(p, stderr, 10); fputc('\n', stderr);
 193     if (r) {
 194       fputs("r  = ", stderr);
 195       mp_writefile(r, stderr, 10);
 196       fputc('\n', stderr);
 197     } else
 198       fputs("r  = <undef>\n", stderr);
 199     fputs("rr = ", stderr); mp_writefile(rr, stderr, 10); fputc('\n', stderr);
 200     ok = 0;
 201   }
 202
 203   mp_drop(a);
 204   mp_drop(p);
 205   mp_drop(r);
 206   mp_drop(rr);
 207   assert(mparena_count(MPARENA_GLOBAL) == 0);
 208   return (ok);
 209 }
 210
 211 static test_chunk tests[] = {
 212   { "modsqrt", verify, { &type_mp, &type_mp, &type_mp, 0 } },
 213   { 0, 0, { 0 } }
 214 };
 215
 216 int main(int argc, char *argv[])
 217 {
 218   sub_init();
 219   test_run(argc, argv, tests, SRCDIR "/tests/mp");
 220   return (0);
 221 }
 222
 223 #endif
 224
 225 /*----- That's all, folks -------------------------------------------------*/