progs/perftest.c: Use from Glibc syscall numbers.

[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;
  size_t i, s;
  mp *b, *c;
  mp *ainv;
  mp *r, *A, *aa;
  mp *t;
  grand *gr;
  int j;

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

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

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

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

  /* --- Some initial setup --- */

  mpmont_create(&mm, p);
  ainv = mp_modinv(MP_NEW, a, p);	/* %$a^{-1} \bmod p$% */
  ainv = mpmont_mul(&mm, ainv, ainv, mm.r2);
  t = mp_sub(MP_NEW, p, MP_ONE);
  t = mp_odd(t, t, &s);			/* %$2^s t = p - 1$% */
  b = mpmont_mul(&mm, b, b, mm.r2);
  c = mpmont_expr(&mm, b, b, t);	/* %$b^t \bmod p$% */
  t = mp_add(t, t, MP_ONE);
  t = mp_lsr(t, t, 1);			/* %$(t + 1)/2$% */
  a = mpmont_mul(&mm, MP_NEW, a, mm.r2);
  r = mpmont_expr(&mm, a, a, t);	/* %$a^{(t+1)/2} \bmod p$% */

  /* --- Now for the main loop --- *
   *
   * Let %$g = c^{-2}$%; we know that %$g$% is a generator of the order-
   * %$2^{s-1}$% subgroup mod %$p$%.  We also know that %$A = a^t = r^2/a$%
   * is an element of this group.  If we can determine %$m$% such that
   * %$g^m = A$% then %$a^{(t+1)/2}/g^{m/2} = r c^m$% is the square root we
   * seek.
   *
   * Write %$m = m_0 + 2 m'$%.  Then %$A^{2^{s-1}} = g^{m_0 2^{s-1}}$%, which
   * is %$1$% if %$m_0 = 0$% or %$-1$% if %$m_0 = 1$% (modulo %$p$%).  Then
   * %$A/g^{m_0} = (g^2)^{m'}$% and we can proceed inductively.  The end
   * result will me %$A/g^m$%.
   *
   * Note that this loop keeps track of (what will be) %$r c^m$%, since this
   * is the result we want, and computes $A/g^m = r^2/a$% on demand.
   */

  A = mp_sqr(t, r); A = mpmont_reduce(&mm, A, A);
  A = mpmont_mul(&mm, A, A, ainv);	/* %$x^t/g^m$% */

  while (s-- > 1) {
    aa = MP_COPY(A);
    for (i = 1; i < s; i++)
      { aa = mp_sqr(aa, aa); aa = mpmont_reduce(&mm, aa, aa); }
    if (!MP_EQ(aa, mm.r)) {
      r = mpmont_mul(&mm, r, r, c);
      A = mp_sqr(A, r); A = mpmont_reduce(&mm, A, A);
      A = mpmont_mul(&mm, A, A, ainv);	/* %$x^t/g^m$% */
    }
    c = mp_sqr(c, c); c = mpmont_reduce(&mm, c, c);
    MP_DROP(aa);
  }

  /* --- We want the smaller square root --- */

  d = mpmont_reduce(&mm, d, r);
  r = mp_sub(r, p, d);
  if (MP_CMP(r, <, d)) { mp *tt = r; r = d; d = tt; }

  /* --- Clear away all the temporaries --- */

  mp_drop(ainv);
  mp_drop(r); mp_drop(c);
  mp_drop(A);
  mpmont_destroy(&mm);

  /* --- Done --- */

  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 -------------------------------------------------*/