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

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