math/mp.h, math/mp-jacobi.c: Whitespace fix.

[u/mdw/catacomb] / math / mp-modsqrt.c

/* -*-c-*-
 *
 * Compute square roots modulo a prime
 *
 * (c) 2000 Straylight/Edgeware
 */

/*----- Licensing notice --------------------------------------------------*
 *
 * This file is part of Catacomb.
 *
 * Catacomb is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Library General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 *
 * Catacomb is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with Catacomb; if not, write to the Free
 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 * MA 02111-1307, USA.
 */

/*----- Header files ------------------------------------------------------*/

#include "fibrand.h"
#include "grand.h"
#include "mp.h"
#include "mpmont.h"
#include "mprand.h"

/*----- Main code ---------------------------------------------------------*/

/* --- @mp_modsqrt@ --- *
 *
 * Arguments:	@mp *d@ = destination integer
 *		@mp *a@ = source integer
 *		@mp *p@ = modulus (must be prime)
 *
 * Returns:	If %$a$% is a quadratic residue, a square root of %$a$%; else
 *		a null pointer.
 *
 * Use:		Returns an integer %$x$% such that %$x^2 \equiv a \pmod{p}$%,
 *		if one exists; else a null pointer.  This function will not
 *		work if %$p$% is composite: you must factor the modulus, take
 *		a square root mod each factor, and recombine the results
 *		using the Chinese Remainder Theorem.
 *
 *		We guarantee that the square root returned is the smallest
 *		one (i.e., the `positive' square root).
 */

mp *mp_modsqrt(mp *d, mp *a, mp *p)
{
  mpmont mm;
  mp *t;
  size_t s;
  mp *b;
  mp *ainv;
  mp *c, *r;
  size_t i, j;
  mp *dd, *mone;

  /* --- Cope if %$a \not\in Q_p$% --- */

  if (mp_jacobi(a, p) != 1) {
    mp_drop(d);
    return (0);
  }

  /* --- Choose some quadratic non-residue --- */

  {
    grand *g = fibrand_create(0);

    b = MP_NEW;
    do
      b = mprand_range(b, p, g, 0);
    while (mp_jacobi(b, p) != -1);
    g->ops->destroy(g);
  }

  /* --- Find the inverse of %$a$% --- */

  ainv = mp_modinv(MP_NEW, a, p);

  /* --- Split %$p - 1$% into a power of two and an odd number --- */

  t = mp_sub(MP_NEW, p, MP_ONE);
  t = mp_odd(t, t, &s);

  /* --- Now to really get going --- */

  mpmont_create(&mm, p);
  b = mpmont_mul(&mm, b, b, mm.r2);
  c = mpmont_expr(&mm, b, b, t);
  t = mp_add(t, t, MP_ONE);
  t = mp_lsr(t, t, 1);
  dd = mpmont_mul(&mm, MP_NEW, a, mm.r2);
  r = mpmont_expr(&mm, t, dd, t);
  mp_drop(dd);
  ainv = mpmont_mul(&mm, ainv, ainv, mm.r2);

  mone = mp_sub(MP_NEW, p, mm.r);

  dd = MP_NEW;

  for (i = 1; i < s; i++) {

    /* --- Compute %$d_0 = r^2a^{-1}$% --- */

    dd = mp_sqr(dd, r);
    dd = mpmont_reduce(&mm, dd, dd);
    dd = mpmont_mul(&mm, dd, dd, ainv);

    /* --- Now %$d = d_0^{2^{s - i - 1}}$% --- */

    for (j = i; j < s - 1; j++) {
      dd = mp_sqr(dd, dd);
      dd = mpmont_reduce(&mm, dd, dd);
    }

    /* --- Fiddle at the end --- */

    if (MP_EQ(dd, mone))
      r = mpmont_mul(&mm, r, r, c);
    c = mp_sqr(c, c);
    c = mpmont_reduce(&mm, c, c);
  }

  /* --- Done, so tidy up --- *
   *
   * Canonify the answer.
   */

  d = mpmont_reduce(&mm, d, r);
  r = mp_sub(r, p, d);
  if (MP_CMP(r, <, d)) { mp *tt = r; r = d; d = tt; }
  mp_drop(ainv);
  mp_drop(r); mp_drop(c);
  mp_drop(dd);
  mp_drop(mone);
  mpmont_destroy(&mm);

  return (d);
}

/*----- Test rig ----------------------------------------------------------*/

#ifdef TEST_RIG

#include <mLib/testrig.h>

static int verify(dstr *v)
{
  mp *a = *(mp **)v[0].buf;
  mp *p = *(mp **)v[1].buf;
  mp *rr = *(mp **)v[2].buf;
  mp *r = mp_modsqrt(MP_NEW, a, p);
  int ok = 0;

  if (!r)
    ok = 0;
  else if (MP_EQ(r, rr))
    ok = 1;

  if (!ok) {
    fputs("\n*** fail\n", stderr);
    fputs("a  = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr);
    fputs("p  = ", stderr); mp_writefile(p, stderr, 10); fputc('\n', stderr);
    if (r) {
      fputs("r	= ", stderr);
      mp_writefile(r, stderr, 10);
      fputc('\n', stderr);
    } else
      fputs("r	= <undef>\n", stderr);
    fputs("rr = ", stderr); mp_writefile(rr, stderr, 10); fputc('\n', stderr);
    ok = 0;
  }

  mp_drop(a);
  mp_drop(p);
  mp_drop(r);
  mp_drop(rr);
  assert(mparena_count(MPARENA_GLOBAL) == 0);
  return (ok);
}

static test_chunk tests[] = {
  { "modsqrt", verify, { &type_mp, &type_mp, &type_mp, 0 } },
  { 0, 0, { 0 } }
};

int main(int argc, char *argv[])
{
  sub_init();
  test_run(argc, argv, tests, SRCDIR "/t/mp");
  return (0);
}

#endif

/*----- That's all, folks -------------------------------------------------*/