chiark - git - mdw - catacomb/blob - mp-modsqrt.c

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