1 /*
2 Copyright (C) 2009 William Hart
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include <gmp.h>
13 #include "flint.h"
14 #include "ulong_extras.h"
15 #include "fmpz.h"
16
17 void
fmpz_submul_ui(fmpz_t f,const fmpz_t g,ulong x)18 fmpz_submul_ui(fmpz_t f, const fmpz_t g, ulong x)
19 {
20 fmpz c1 = *g, r;
21 if ((x == 0) || (c1 == 0))
22 return; /* product is zero */
23
24 r = *f;
25 if (r == 0)
26 {
27 fmpz_mul_ui(f, g, x); /* we are subtracting product from 0 */
28 fmpz_neg(f, f);
29 return;
30 }
31
32 if (!COEFF_IS_MPZ(c1)) /* c1 is small */
33 {
34 mp_limb_t prod[2];
35 ulong uc1 = FLINT_ABS(c1);
36
37 __mpz_struct * mpz_ptr;
38 mpz_t temp;
39
40 umul_ppmm(prod[1], prod[0], uc1, x); /* compute product */
41
42 if (prod[1] == 0) /* product fits in one limb */
43 {
44 if (c1 < WORD(0))
45 fmpz_add_ui(f, f, prod[0]);
46 else
47 fmpz_sub_ui(f, f, prod[0]);
48 return;
49 }
50 else if ((prod[1] == 1) && (!COEFF_IS_MPZ(r)) && ((r ^ c1) >= WORD(0)))
51 {
52 /*
53 only chance at cancellation is if product is one bit past
54 a limb and f is small and same sign as this product
55 */
56 ulong ur = FLINT_ABS(r);
57 if (ur > prod[0]) /* cancellation will occur */
58 {
59 fmpz_set_ui(f, prod[0] - ur);
60 if (r > WORD(0))
61 fmpz_neg(f, f);
62 return;
63 }
64 }
65
66 /*
67 in all remaining cases res is either big already,
68 or will be big in the end
69 */
70 mpz_ptr = _fmpz_promote_val(f);
71 /* set up a temporary, cheap mpz_t to contain prod */
72 temp->_mp_d = prod;
73 temp->_mp_size = (c1 < WORD(0) ? -2 : 2);
74 mpz_sub(mpz_ptr, mpz_ptr, temp);
75 _fmpz_demote_val(f); /* cancellation may have occurred */
76 }
77 else /* c1 is large */
78 {
79 __mpz_struct *mpz_ptr = _fmpz_promote_val(f);
80 flint_mpz_submul_ui(mpz_ptr, COEFF_TO_PTR(c1), x);
81 _fmpz_demote_val(f); /* cancellation may have occurred */
82 }
83 }
84