1 /*
2 Copyright (C) 2008, 2009, William Hart
3 Copyright (C) 2010 Fredrik Johansson
4
5 This file is part of FLINT.
6
7 FLINT is free software: you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License (LGPL) as published
9 by the Free Software Foundation; either version 2.1 of the License, or
10 (at your option) any later version. See <http://www.gnu.org/licenses/>.
11 */
12
13 #include <gmp.h>
14 #include "flint.h"
15 #include "ulong_extras.h"
16 #include "fmpz.h"
17 #include "nmod_vec.h"
18
19
20 void
__fmpz_multi_CRT_ui_sign(fmpz_t output,const fmpz_t input,const fmpz_comb_t comb,fmpz_t temp)21 __fmpz_multi_CRT_ui_sign(fmpz_t output, const fmpz_t input,
22 const fmpz_comb_t comb, fmpz_t temp)
23 {
24 slong n = comb->n;
25 slong p;
26
27 if (n == WORD(0))
28 {
29 if (fmpz_is_zero(input))
30 {
31 fmpz_zero(output);
32 return;
33 }
34
35 /* XXX: overflow possible? */
36 p = comb->primes[0];
37 if ((p - (*input)) < (*input))
38 fmpz_set_si(output, (slong) ((*input) - p));
39 else
40 fmpz_set_ui(output, (*input));
41 return;
42 }
43
44 fmpz_sub(temp, input, comb->comb[comb->n - 1]);
45
46 if (fmpz_cmpabs(temp, input) <= 0)
47 fmpz_set(output, temp);
48 else
49 fmpz_set(output, input);
50
51 return;
52 }
53
fmpz_multi_CRT_ui(fmpz_t output,mp_srcptr residues,const fmpz_comb_t comb,fmpz_comb_temp_t ctemp,int sign)54 void fmpz_multi_CRT_ui(fmpz_t output, mp_srcptr residues,
55 const fmpz_comb_t comb, fmpz_comb_temp_t ctemp, int sign)
56 {
57 slong i, j;
58 slong n = comb->n;
59 slong num;
60 slong log_res;
61 slong num_primes = comb->num_primes;
62
63 fmpz ** comb_temp = ctemp->comb_temp;
64 fmpz * temp = ctemp->temp;
65 fmpz * temp2 = ctemp->temp2;
66
67 /* The output is less than a single prime, so just output the result */
68 if (num_primes == 1)
69 {
70 if (sign)
71 {
72 mp_limb_t p = comb->primes[0];
73
74 if ((p - residues[0]) < residues[0])
75 fmpz_set_si(output, residues[0] - p);
76 else
77 fmpz_set_ui(output, residues[0]);
78 }
79 else
80 {
81 fmpz_set_ui(output, residues[0]);
82 }
83 return;
84 }
85
86 /* First layer of reconstruction */
87 num = (WORD(1) << n);
88
89 for (i = 0, j = 0; i + 2 <= num_primes; i += 2, j++)
90 {
91 fmpz_set_ui(temp, residues[i]);
92 fmpz_mod_ui(temp2, temp, comb->primes[i+1]);
93 fmpz_sub_ui(temp2, temp2, residues[i + 1]);
94 fmpz_neg(temp2, temp2);
95 fmpz_mul(temp, temp2, comb->res[0] + j);
96 fmpz_mod_ui(temp2, temp, comb->primes[i+1]);
97 fmpz_mul_ui(temp, temp2, comb->primes[i]);
98 fmpz_add_ui(comb_temp[0] + j, temp, residues[i]);
99 }
100
101 if (i < num_primes)
102 fmpz_set_ui(comb_temp[0] + j, residues[i]);
103
104 /* Compute other layers of reconstruction */
105 num /= 2;
106 log_res = 1;
107
108 while (log_res < n)
109 {
110 for (i = 0, j = 0; i < num; i += 2, j++)
111 {
112 if (fmpz_is_one(comb->comb[log_res-1] + i + 1))
113 {
114 if (!fmpz_is_one(comb->comb[log_res-1] + i))
115 fmpz_set(comb_temp[log_res] + j, comb_temp[log_res-1] + i);
116 }
117 else
118 {
119 fmpz_mod(temp2, comb_temp[log_res-1] + i,
120 comb->comb[log_res-1] + i + 1);
121 fmpz_sub(temp, comb_temp[log_res-1] + i + 1, temp2);
122 fmpz_mul(temp2, temp, comb->res[log_res] + j);
123 fmpz_mod(temp, temp2, comb->comb[log_res-1] + i + 1);
124 fmpz_mul(temp2, temp, comb->comb[log_res-1] + i);
125 fmpz_add(comb_temp[log_res] + j, temp2,
126 comb_temp[log_res-1] + i);
127 }
128 }
129 log_res++;
130 num /= 2;
131 }
132
133 /* Write out the output */
134 if (sign)
135 __fmpz_multi_CRT_ui_sign(output, comb_temp[log_res - 1], comb, temp);
136 else
137 fmpz_set(output, comb_temp[log_res - 1]);
138 }
139