chiark - git - 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  *              We guarantee that the square root returned is the smallest
  56  *              one (i.e., the `positive' square root).
  57  */
  58
  59 mp *mp_modsqrt(mp *d, mp *a, mp *p)
  60 {
  61   mpmont mm;
  62   mp *t;
  63   size_t s;
  64   mp *b;
  65   mp *ainv;
  66   mp *c, *r;
  67   size_t i, j;
  68   mp *dd, *mone;
  69
  70   /* --- Cope if %$a \not\in Q_p$% --- */
  71
  72   if (mp_jacobi(a, p) != 1) {
  73     mp_drop(d);
  74     return (0);
  75   }
  76
  77   /* --- Choose some quadratic non-residue --- */
  78
  79   {
  80     grand *g = fibrand_create(0);
  81
  82     b = MP_NEW;
  83     do
  84       b = mprand_range(b, p, g, 0);
  85     while (mp_jacobi(b, p) != -1);
  86     g->ops->destroy(g);
  87   }
  88
  89   /* --- Find the inverse of %$a$% --- */
  90
  91   ainv = mp_modinv(MP_NEW, a, p);
  92
  93   /* --- Split %$p - 1$% into a power of two and an odd number --- */
  94
  95   t = mp_sub(MP_NEW, p, MP_ONE);
  96   t = mp_odd(t, t, &s);
  97
  98   /* --- Now to really get going --- */
  99
 100   mpmont_create(&mm, p);
 101   b = mpmont_mul(&mm, b, b, mm.r2);
 102   c = mpmont_expr(&mm, b, b, t);
 103   t = mp_add(t, t, MP_ONE);
 104   t = mp_lsr(t, t, 1);
 105   dd = mpmont_mul(&mm, MP_NEW, a, mm.r2);
 106   r = mpmont_expr(&mm, t, dd, t);
 107   mp_drop(dd);
 108   ainv = mpmont_mul(&mm, ainv, ainv, mm.r2);
 109
 110   mone = mp_sub(MP_NEW, p, mm.r);
 111
 112   dd = MP_NEW;
 113
 114   for (i = 1; i < s; i++) {
 115
 116     /* --- Compute %$d_0 = r^2a^{-1}$% --- */
 117
 118     dd = mp_sqr(dd, r);
 119     dd = mpmont_reduce(&mm, dd, dd);
 120     dd = mpmont_mul(&mm, dd, dd, ainv);
 121
 122     /* --- Now %$d = d_0^{s - i - 1}$% --- */
 123
 124     for (j = i; j < s - 1; j++) {
 125       dd = mp_sqr(dd, dd);
 126       dd = mpmont_reduce(&mm, dd, dd);
 127     }
 128
 129     /* --- Fiddle at the end --- */
 130
 131     if (MP_EQ(dd, mone))
 132       r = mpmont_mul(&mm, r, r, c);
 133     c = mp_sqr(c, c);
 134     c = mpmont_reduce(&mm, c, c);
 135   }
 136
 137   /* --- Done, so tidy up --- *
 138    *
 139    * Canonify the answer.
 140    */
 141
 142   d = mpmont_reduce(&mm, d, r);
 143   r = mp_sub(r, p, d);
 144   if (MP_CMP(r, <, d)) { mp *tt = r; r = d; d = tt; }
 145   mp_drop(ainv);
 146   mp_drop(r); mp_drop(c);
 147   mp_drop(dd);
 148   mp_drop(mone);
 149   mpmont_destroy(&mm);
 150
 151   return (d);
 152 }
 153
 154 /*----- Test rig ----------------------------------------------------------*/
 155
 156 #ifdef TEST_RIG
 157
 158 #include <mLib/testrig.h>
 159
 160 static int verify(dstr *v)
 161 {
 162   mp *a = *(mp **)v[0].buf;
 163   mp *p = *(mp **)v[1].buf;
 164   mp *rr = *(mp **)v[2].buf;
 165   mp *r = mp_modsqrt(MP_NEW, a, p);
 166   int ok = 0;
 167
 168   if (!r)
 169     ok = 0;
 170   else if (MP_EQ(r, rr))
 171     ok = 1;
 172
 173   if (!ok) {
 174     fputs("\n*** fail\n", stderr);
 175     fputs("a  = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr);
 176     fputs("p  = ", stderr); mp_writefile(p, stderr, 10); fputc('\n', stderr);
 177     if (r) {
 178       fputs("r  = ", stderr);
 179       mp_writefile(r, stderr, 10);
 180       fputc('\n', stderr);
 181     } else
 182       fputs("r  = <undef>\n", stderr);
 183     fputs("rr = ", stderr); mp_writefile(rr, stderr, 10); fputc('\n', stderr);
 184     ok = 0;
 185   }
 186
 187   mp_drop(a);
 188   mp_drop(p);
 189   mp_drop(r);
 190   mp_drop(rr);
 191   assert(mparena_count(MPARENA_GLOBAL) == 0);
 192   return (ok);
 193 }
 194
 195 static test_chunk tests[] = {
 196   { "modsqrt", verify, { &type_mp, &type_mp, &type_mp, 0 } },
 197   { 0, 0, { 0 } }
 198 };
 199
 200 int main(int argc, char *argv[])
 201 {
 202   sub_init();
 203   test_run(argc, argv, tests, SRCDIR "/tests/mp");
 204   return (0);
 205 }
 206
 207 #endif
 208
 209 /*----- That's all, folks -------------------------------------------------*/