|
Packit |
5c3484 |
/* hgcd_step.c.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY
|
|
Packit |
5c3484 |
SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
|
|
Packit |
5c3484 |
GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Copyright 2003-2005, 2008, 2011, 2012 Free Software Foundation, Inc.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
This file is part of the GNU MP Library.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The GNU MP Library is free software; you can redistribute it and/or modify
|
|
Packit |
5c3484 |
it under the terms of either:
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
* the GNU Lesser General Public License as published by the Free
|
|
Packit |
5c3484 |
Software Foundation; either version 3 of the License, or (at your
|
|
Packit |
5c3484 |
option) any later version.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
or
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
* the GNU General Public License as published by the Free Software
|
|
Packit |
5c3484 |
Foundation; either version 2 of the License, or (at your option) any
|
|
Packit |
5c3484 |
later version.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
or both in parallel, as here.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The GNU MP Library is distributed in the hope that it will be useful, but
|
|
Packit |
5c3484 |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
|
Packit |
5c3484 |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
Packit |
5c3484 |
for more details.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
You should have received copies of the GNU General Public License and the
|
|
Packit |
5c3484 |
GNU Lesser General Public License along with the GNU MP Library. If not,
|
|
Packit |
5c3484 |
see https://www.gnu.org/licenses/. */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#include "gmp.h"
|
|
Packit |
5c3484 |
#include "gmp-impl.h"
|
|
Packit |
5c3484 |
#include "longlong.h"
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
static void
|
|
Packit |
5c3484 |
hgcd_hook (void *p, mp_srcptr gp, mp_size_t gn,
|
|
Packit |
5c3484 |
mp_srcptr qp, mp_size_t qn, int d)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
ASSERT (!gp);
|
|
Packit |
5c3484 |
ASSERT (d >= 0);
|
|
Packit |
5c3484 |
ASSERT (d <= 1);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
MPN_NORMALIZE (qp, qn);
|
|
Packit |
5c3484 |
if (qn > 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
struct hgcd_matrix *M = (struct hgcd_matrix *) p;
|
|
Packit |
5c3484 |
/* NOTES: This is a bit ugly. A tp area is passed to
|
|
Packit |
5c3484 |
gcd_subdiv_step, which stores q at the start of that area. We
|
|
Packit |
5c3484 |
now use the rest. */
|
|
Packit |
5c3484 |
mp_ptr tp = (mp_ptr) qp + qn;
|
|
Packit |
5c3484 |
mpn_hgcd_matrix_update_q (M, qp, qn, d, tp);
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Perform a few steps, using some of mpn_hgcd2, subtraction and
|
|
Packit |
5c3484 |
division. Reduces the size by almost one limb or more, but never
|
|
Packit |
5c3484 |
below the given size s. Return new size for a and b, or 0 if no
|
|
Packit |
5c3484 |
more steps are possible.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
If hgcd2 succeeds, needs temporary space for hgcd_matrix_mul_1, M->n
|
|
Packit |
5c3484 |
limbs, and hgcd_mul_matrix1_inverse_vector, n limbs. If hgcd2
|
|
Packit |
5c3484 |
fails, needs space for the quotient, qn <= n - s limbs, for and
|
|
Packit |
5c3484 |
hgcd_matrix_update_q, qn + (size of the appropriate column of M) <=
|
|
Packit |
5c3484 |
(resulting size of M) + 1.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
If N is the input size to the calling hgcd, then s = floor(N/2) +
|
|
Packit |
5c3484 |
1, M->n < N, qn + product size <= n - s + n - s + 1 = 2 (n - s) + 1
|
|
Packit |
5c3484 |
<= N.
|
|
Packit |
5c3484 |
*/
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
mp_size_t
|
|
Packit |
5c3484 |
mpn_hgcd_step (mp_size_t n, mp_ptr ap, mp_ptr bp, mp_size_t s,
|
|
Packit |
5c3484 |
struct hgcd_matrix *M, mp_ptr tp)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
struct hgcd_matrix1 M1;
|
|
Packit |
5c3484 |
mp_limb_t mask;
|
|
Packit |
5c3484 |
mp_limb_t ah, al, bh, bl;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
ASSERT (n > s);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
mask = ap[n-1] | bp[n-1];
|
|
Packit |
5c3484 |
ASSERT (mask > 0);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (n == s + 1)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
if (mask < 4)
|
|
Packit |
5c3484 |
goto subtract;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
ah = ap[n-1]; al = ap[n-2];
|
|
Packit |
5c3484 |
bh = bp[n-1]; bl = bp[n-2];
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else if (mask & GMP_NUMB_HIGHBIT)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
ah = ap[n-1]; al = ap[n-2];
|
|
Packit |
5c3484 |
bh = bp[n-1]; bl = bp[n-2];
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
int shift;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
count_leading_zeros (shift, mask);
|
|
Packit |
5c3484 |
ah = MPN_EXTRACT_NUMB (shift, ap[n-1], ap[n-2]);
|
|
Packit |
5c3484 |
al = MPN_EXTRACT_NUMB (shift, ap[n-2], ap[n-3]);
|
|
Packit |
5c3484 |
bh = MPN_EXTRACT_NUMB (shift, bp[n-1], bp[n-2]);
|
|
Packit |
5c3484 |
bl = MPN_EXTRACT_NUMB (shift, bp[n-2], bp[n-3]);
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Try an mpn_hgcd2 step */
|
|
Packit |
5c3484 |
if (mpn_hgcd2 (ah, al, bh, bl, &M1))
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* Multiply M <- M * M1 */
|
|
Packit |
5c3484 |
mpn_hgcd_matrix_mul_1 (M, &M1, tp);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Can't swap inputs, so we need to copy. */
|
|
Packit |
5c3484 |
MPN_COPY (tp, ap, n);
|
|
Packit |
5c3484 |
/* Multiply M1^{-1} (a;b) */
|
|
Packit |
5c3484 |
return mpn_matrix22_mul1_inverse_vector (&M1, ap, tp, bp, n);
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
subtract:
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
return mpn_gcd_subdiv_step (ap, bp, n, s, hgcd_hook, M, tp);
|
|
Packit |
5c3484 |
}
|