1 /*
2 Copyright (C) 2009 William Hart
3 Copyright (C) 2011 Fredrik Johansson
4 Copyright (C) 2014 William Hart
5
6 This file is part of FLINT.
7
8 FLINT is free software: you can redistribute it and/or modify it under
9 the terms of the GNU Lesser General Public License (LGPL) as published
10 by the Free Software Foundation; either version 2.1 of the License, or
11 (at your option) any later version. See <http://www.gnu.org/licenses/>.
12 */
13
14 #include <stdlib.h>
15 #include <stdio.h>
16 #include <gmp.h>
17 #include "flint.h"
18 #include "ulong_extras.h"
19 #include "fmpz.h"
20
21 void
_fmpz_CRT(fmpz_t out,const fmpz_t r1,const fmpz_t m1,fmpz_t r2,fmpz_t m2,const fmpz_t m1m2,fmpz_t c,int sign)22 _fmpz_CRT(fmpz_t out, const fmpz_t r1, const fmpz_t m1, fmpz_t r2,
23 fmpz_t m2, const fmpz_t m1m2, fmpz_t c, int sign)
24 {
25 fmpz_t r1normal, tmp, r1mod, s;
26
27 fmpz_init(tmp);
28 fmpz_init(r1mod);
29 fmpz_init(s);
30
31 /* FIXME: assume r1 moved to [0, m1); add tests for this */
32 if (fmpz_sgn(r1) < 0)
33 {
34 fmpz_init(r1normal);
35 fmpz_add(r1normal, r1, m1);
36 }
37 else
38 {
39 *r1normal = *r1;
40 }
41
42 fmpz_mod(r1mod, r1normal, m2);
43 fmpz_sub(s, r2, r1mod);
44 if (fmpz_sgn(s) < 0)
45 fmpz_add(s, s, m2);
46 fmpz_mul(s, s, c);
47 fmpz_mod(s, s, m2);
48 fmpz_mul(tmp, m1, s);
49 fmpz_add(tmp, tmp, r1normal);
50
51 if (fmpz_sgn(r1) < 0)
52 fmpz_clear(r1normal);
53
54 if (sign)
55 {
56 fmpz_sub(out, tmp, m1m2);
57 if (fmpz_cmpabs(tmp, out) <= 0)
58 fmpz_set(out, tmp);
59 }
60 else
61 {
62 fmpz_set(out, tmp);
63 }
64
65 fmpz_clear(tmp);
66 fmpz_clear(r1mod);
67 fmpz_clear(s);
68 }
69
fmpz_CRT(fmpz_t out,const fmpz_t r1,const fmpz_t m1,fmpz_t r2,fmpz_t m2,int sign)70 void fmpz_CRT(fmpz_t out, const fmpz_t r1, const fmpz_t m1,
71 fmpz_t r2, fmpz_t m2, int sign)
72 {
73 fmpz_t m1m2, c;
74
75 fmpz_init(c);
76
77 fmpz_mod(c, m1, m2);
78 fmpz_invmod(c, c, m2);
79
80 if (fmpz_is_zero(c))
81 {
82 flint_printf("Exception (fmpz_CRT). m1 not invertible modulo m2.\n");
83 flint_abort();
84 }
85
86 fmpz_init(m1m2);
87 fmpz_mul(m1m2, m1, m2);
88
89 _fmpz_CRT(out, r1, m1, r2, m2, m1m2, c, sign);
90
91 fmpz_clear(m1m2);
92 fmpz_clear(c);
93 }
94