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