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

   1 /* -*-c-*-
   2  *
   3  * $Id: mp-modsqrt.c,v 1.2 2000/10/08 12:02:21 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.2  2000/10/08 12:02:21  mdw
  34  * Use @MP_EQ@ instead of @MP_CMP@.
  35  *
  36  * Revision 1.1  2000/06/22 19:01:31  mdw
  37  * Compute square roots in a prime field.
  38  *
  39  */
  40
  41 /*----- Header files ------------------------------------------------------*/
  42
  43 #include "fibrand.h"
  44 #include "grand.h"
  45 #include "mp.h"
  46 #include "mpmont.h"
  47 #include "mprand.h"
  48
  49 /*----- Main code ---------------------------------------------------------*/
  50
  51 /* --- @mp_modsqrt@ --- *
  52  *
  53  * Arguments:   @mp *d@ = destination integer
  54  *              @mp *a@ = source integer
  55  *              @mp *p@ = modulus (must be prime)
  56  *
  57  * Returns:     If %$a$% is a quadratic residue, a square root of %$a$%; else
  58  *              a null pointer.
  59  *
  60  * Use:         Returns an integer %$x$% such that %$x^2 \equiv a \pmod{p}$%,
  61  *              if one exists; else a null pointer.  This function will not
  62  *              work if %$p$% is composite: you must factor the modulus, take
  63  *              a square root mod each factor, and recombine the results
  64  *              using the Chinese Remainder Theorem.
  65  */
  66
  67 mp *mp_modsqrt(mp *d, mp *a, mp *p)
  68 {
  69   mpmont mm;
  70   mp *t;
  71   size_t s;
  72   mp *b;
  73   mp *ainv;
  74   mp *c, *r;
  75   size_t i, j;
  76   mp *dd, *mone;
  77
  78   /* --- Cope if %$a \not\in Q_p$% --- */
  79
  80   if (mp_jacobi(a, p) != 1) {
  81     if (d)
  82       mp_drop(d);
  83     return (0);
  84   }
  85
  86   /* --- Choose some quadratic non-residue --- */
  87
  88   {
  89     grand *g = fibrand_create(0);
  90
  91     b = MP_NEW;
  92     do
  93       b = mprand_range(b, p, g, 0);
  94     while (mp_jacobi(b, p) != -1);
  95     g->ops->destroy(g);
  96   }
  97
  98   /* --- Find the inverse of %$a$% --- */
  99
 100   ainv = MP_NEW;
 101   mp_gcd(0, &ainv, 0, a, p);
 102
 103   /* --- Split %$p - 1$% into a power of two and an odd number --- */
 104
 105   t = mp_sub(MP_NEW, p, MP_ONE);
 106   t = mp_odd(t, t, &s);
 107
 108   /* --- Now to really get going --- */
 109
 110   mpmont_create(&mm, p);
 111   c = mpmont_expr(&mm, b, b, t);
 112   t = mp_add(t, t, MP_ONE);
 113   t = mp_lsr(t, t, 1);
 114   r = mpmont_expr(&mm, t, a, t);
 115   ainv = mpmont_mul(&mm, ainv, ainv, mm.r2);
 116
 117   mone = mp_sub(MP_NEW, p, mm.r);
 118
 119   dd = MP_NEW;
 120
 121   for (i = 1; i < s; i++) {
 122
 123     /* --- Compute %$d_0 = r^2a^{-1}$% --- */
 124
 125     dd = mp_sqr(dd, r);
 126     dd = mpmont_reduce(&mm, dd, dd);
 127     dd = mpmont_mul(&mm, dd, dd, ainv);
 128
 129     /* --- Now %$d = d_0^{s - i - 1}$% --- */
 130
 131     for (j = i; j < s - 1; j++) {
 132       dd = mp_sqr(dd, dd);
 133       dd = mpmont_reduce(&mm, dd, dd);
 134     }
 135
 136     /* --- Fiddle at the end --- */
 137
 138     if (MP_EQ(dd, mone))
 139       r = mpmont_mul(&mm, r, r, c);
 140     c = mp_sqr(c, c);
 141     c = mpmont_reduce(&mm, c, c);
 142   }
 143
 144   /* --- Done, so tidy up --- */
 145
 146   d = mpmont_reduce(&mm, d, r);
 147   mp_drop(ainv);
 148   mp_drop(r); mp_drop(c);
 149   if (dd)
 150     mp_drop(dd);
 151   mp_drop(mone);
 152   mpmont_destroy(&mm);
 153
 154   return (d);
 155 }
 156
 157 /*----- Test rig ----------------------------------------------------------*/
 158
 159 #ifdef TEST_RIG
 160
 161 #include <mLib/testrig.h>
 162
 163 static int verify(dstr *v)
 164 {
 165   mp *a = *(mp **)v[0].buf;
 166   mp *p = *(mp **)v[1].buf;
 167   mp *rr = *(mp **)v[2].buf;
 168   mp *r = mp_modsqrt(MP_NEW, a, p);
 169   int ok = 0;
 170
 171   if (!r)
 172     ok = 0;
 173   else if (MP_EQ(r, rr))
 174     ok = 1;
 175   else {
 176     r = mp_sub(r, p, r);
 177     if (MP_EQ(r, rr))
 178       ok = 1;
 179   }
 180
 181   if (!ok) {
 182     fputs("\n*** fail\n", stderr);
 183     fputs("a  = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr);
 184     fputs("p  = ", stderr); mp_writefile(p, stderr, 10); fputc('\n', stderr);
 185     if (r) {
 186       fputs("r  = ", stderr);
 187       mp_writefile(r, stderr, 10);
 188       fputc('\n', stderr);
 189     } else
 190       fputs("r  = <undef>\n", stderr);
 191     fputs("rr = ", stderr); mp_writefile(rr, stderr, 10); fputc('\n', stderr);
 192     ok = 0;
 193   }
 194
 195   mp_drop(a);
 196   mp_drop(p);
 197   if (r)
 198     mp_drop(r);
 199   mp_drop(rr);
 200   assert(mparena_count(MPARENA_GLOBAL) == 0);
 201   return (ok);
 202 }
 203
 204 static test_chunk tests[] = {
 205   { "modsqrt", verify, { &type_mp, &type_mp, &type_mp, 0 } },
 206   { 0, 0, { 0 } }
 207 };
 208
 209 int main(int argc, char *argv[])
 210 {
 211   sub_init();
 212   test_run(argc, argv, tests, SRCDIR "/tests/mp");
 213   return (0);
 214 }
 215
 216 #endif
 217
 218 /*----- That's all, folks -------------------------------------------------*/