1 /*
2 Copyright (C) 2010 Sebastian Pancratz
3 Copyright (C) 2010 William Hart
4 Copyright (C) 2011 Fredrik Johansson
5 Copyright (C) 2012 Lina Kulakova
6 Copyright (C) 2013 Martin Lee
7
8 This file is part of FLINT.
9
10 FLINT is free software: you can redistribute it and/or modify it under
11 the terms of the GNU Lesser General Public License (LGPL) as published
12 by the Free Software Foundation; either version 2.1 of the License, or
13 (at your option) any later version. See <http://www.gnu.org/licenses/>.
14 */
15
16 #undef ulong
17 #define ulong ulongxx/* interferes with system includes */
18
19 #include <stdlib.h>
20
21 #undef ulong
22
23 #include <gmp.h>
24
25 #define ulong mp_limb_t
26
27 #include "flint.h"
28 #include "fmpz_vec.h"
29 #include "fmpz_mod_poly.h"
30 #include "ulong_extras.h"
31
32 void
_fmpz_mod_poly_powmod_ui_binexp_preinv(fmpz * res,const fmpz * poly,ulong e,const fmpz * f,slong lenf,const fmpz * finv,slong lenfinv,const fmpz_t p)33 _fmpz_mod_poly_powmod_ui_binexp_preinv(fmpz * res, const fmpz * poly,
34 ulong e, const fmpz * f, slong lenf,
35 const fmpz * finv, slong lenfinv, const fmpz_t p)
36 {
37 fmpz * T, * Q;
38 slong lenT, lenQ;
39 int i;
40
41 if (lenf == 2)
42 {
43 fmpz_powm_ui(res, poly, e, p);
44 return;
45 }
46
47 lenT = 2 * lenf - 3;
48 lenQ = FLINT_MAX(lenT - lenf + 1, 1);
49
50 T = _fmpz_vec_init(lenT + lenQ);
51 Q = T + lenT;
52
53 _fmpz_vec_set(res, poly, lenf - 1);
54
55 for (i = ((int) FLINT_BIT_COUNT(e) - 2); i >= 0; i--)
56 {
57 _fmpz_mod_poly_sqr(T, res, lenf - 1, p);
58 _fmpz_mod_poly_divrem_newton_n_preinv(Q, res, T, 2 * lenf - 3, f, lenf,
59 finv, lenfinv, p);
60
61 if (e & (UWORD (1) << i))
62 {
63 _fmpz_mod_poly_mul(T, res, lenf - 1, poly, lenf - 1, p);
64 _fmpz_mod_poly_divrem_newton_n_preinv(Q, res, T, 2 * lenf - 3, f,
65 lenf, finv, lenfinv, p);
66 }
67 }
68
69 _fmpz_vec_clear(T, lenT + lenQ);
70 }
71
72
73 void
fmpz_mod_poly_powmod_ui_binexp_preinv(fmpz_mod_poly_t res,const fmpz_mod_poly_t poly,ulong e,const fmpz_mod_poly_t f,const fmpz_mod_poly_t finv)74 fmpz_mod_poly_powmod_ui_binexp_preinv(fmpz_mod_poly_t res,
75 const fmpz_mod_poly_t poly, ulong e,
76 const fmpz_mod_poly_t f, const fmpz_mod_poly_t finv)
77 {
78 fmpz * q;
79 slong len = poly->length;
80 slong lenf = f->length;
81 slong trunc = lenf - 1;
82 int qcopy = 0;
83
84 if (lenf == 0)
85 {
86 flint_printf("Exception (fmpz_mod_poly_powmod_ui_binexp_preinv)."
87 "Divide by zero\n");
88 flint_abort();
89 }
90
91 if (lenf == 1)
92 {
93 fmpz_mod_poly_zero(res);
94 return;
95 }
96
97 if (len >= lenf)
98 {
99 fmpz_mod_poly_t t, r;
100 fmpz_mod_poly_init(t, &res->p);
101 fmpz_mod_poly_init(r, &res->p);
102 fmpz_mod_poly_divrem(t, r, poly, f);
103 fmpz_mod_poly_powmod_ui_binexp_preinv(res, r, e, f, finv);
104 fmpz_mod_poly_clear(t);
105 fmpz_mod_poly_clear(r);
106 return;
107 }
108
109 if (e <= 2)
110 {
111 if (e == UWORD (0))
112 {
113 fmpz_mod_poly_fit_length(res, 1);
114 fmpz_one(res->coeffs);
115 _fmpz_mod_poly_set_length(res, 1);
116 }
117 else if (e == UWORD (1))
118 {
119 fmpz_mod_poly_set(res, poly);
120 }
121 else
122 fmpz_mod_poly_mulmod_preinv(res, poly, poly, f, finv);
123 return;
124 }
125
126 if (len == 0)
127 {
128 fmpz_mod_poly_zero(res);
129 return;
130 }
131
132 if (len < trunc)
133 {
134 q = _fmpz_vec_init(trunc);
135 _fmpz_vec_set(q, poly->coeffs, len);
136 _fmpz_vec_zero(q + len, trunc - len);
137 qcopy = 1;
138 } else
139 q = poly->coeffs;
140
141 if ((res == poly && !qcopy) || (res == f) || (res == finv))
142 {
143 fmpz_mod_poly_t t;
144 fmpz_mod_poly_init2(t, &poly->p, 2 * lenf - 3);
145 _fmpz_mod_poly_powmod_ui_binexp_preinv(t->coeffs,
146 q, e, f->coeffs, lenf, finv->coeffs, finv->length, &poly->p);
147 fmpz_mod_poly_swap(res, t);
148 fmpz_mod_poly_clear(t);
149 }
150 else
151 {
152 fmpz_mod_poly_fit_length(res, 2 * lenf - 3);
153 _fmpz_mod_poly_powmod_ui_binexp_preinv(res->coeffs,
154 q, e, f->coeffs, lenf, finv->coeffs, finv->length, &poly->p);
155 }
156
157 if (qcopy)
158 _fmpz_vec_clear(q, trunc);
159
160 _fmpz_mod_poly_set_length(res, trunc);
161 _fmpz_mod_poly_normalise(res);
162 }
163