matrix22_mul1_inverse_vector.c revision 1.1.1.1 
1/* matrix22_mul1_inverse_vector.c
2
3   THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
4   SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
5   GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
6
7Copyright 2008, 2010 Free Software Foundation, Inc.
8
9This file is part of the GNU MP Library.
10
11The GNU MP Library is free software; you can redistribute it and/or modify
12it under the terms of the GNU Lesser General Public License as published by
13the Free Software Foundation; either version 3 of the License, or (at your
14option) any later version.
15
16The GNU MP Library is distributed in the hope that it will be useful, but
17WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
19License for more details.
20
21You should have received a copy of the GNU Lesser General Public License
22along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
23
24#include "gmp.h"
25#include "gmp-impl.h"
26#include "longlong.h"
27
28/* Sets (r;b) = M^{-1}(a;b), with M^{-1} = (u11, -u01; -u10, u00) from
29   the left. Uses three buffers, to avoid a copy. */
30mp_size_t
31mpn_matrix22_mul1_inverse_vector (const struct hgcd_matrix1 *M,
32				  mp_ptr rp, mp_srcptr ap, mp_ptr bp, mp_size_t n)
33{
34  mp_limb_t h0, h1;
35
36  /* Compute (r;b) <-- (u11 a - u01 b; -u10 a + u00 b) as
37
38     r  = u11 * a
39     r -= u01 * b
40     b *= u00
41     b -= u10 * a
42  */
43
44  h0 =    mpn_mul_1 (rp, ap, n, M->u[1][1]);
45  h1 = mpn_submul_1 (rp, bp, n, M->u[0][1]);
46  ASSERT (h0 == h1);
47
48  h0 =    mpn_mul_1 (bp, bp, n, M->u[0][0]);
49  h1 = mpn_submul_1 (bp, ap, n, M->u[1][0]);
50  ASSERT (h0 == h1);
51
52  n -= (rp[n-1] | bp[n-1]) == 0;
53  return n;
54}
55