Fix chinese_remainder() so modular GCD actually works.
[ginac.git] / ginac / polynomial / poly_cra.h
1 #ifndef GINAC_POLY_CRA_H
2 #define GINAC_POLY_CRA_H
3 #include "ex.h"
4 #include <cln/integer.h>
5 #include "smod_helpers.h"
6
7 namespace GiNaC
8 {
9
10 /**
11  * @brief Chinese reamainder algorithm for polynomials.
12  *
13  * Given two polynomials \f$e_1 \in Z_{q_1}[x_1, \ldots, x_n]\f$ and
14  * \f$e_2 \in Z_{q_2}[x_1, \ldots, x_n]\f$, compute the polynomial
15  * \f$r \in Z_{q_1 q_2}[x_1, \ldots, x_n]\f$ such that \f$ r mod q_1 = e_1\f$
16  * and \f$ r mod q_2 = e_2 \f$ 
17  */
18 ex chinese_remainder(const ex& e1, const cln::cl_I& q1,
19                      const ex& e2, const long q2)
20 {
21         // res = v_1 + v_2 q_1
22         // v_1 = e_1 mod q_1
23         // v_2 = (e_2 - v_1)/q_1 mod q_2
24         const numeric q2n(q2);
25         const numeric q1n(q1);
26         ex v1 = e1.smod(q1n);
27         ex u = v1.smod(q2n);
28         ex v2 = (e2.smod(q2n) - v1.smod(q2n)).expand().smod(q2n);
29         const numeric q1_1(recip(q1, q2)); // 1/q_1 mod q_2
30         v2 = (v2*q1_1).smod(q2n);
31         ex ret = (v1 + v2*q1n).expand();
32         return ret;
33 }
34
35 } // namespace GiNaC
36
37 #endif /* GINAC_POLY_CRA_H */
38