1 /*
2 Copyright (C) 2010, 2011 Sebastian Pancratz
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 <stdlib.h>
13 #include <gmp.h>
14 #include "flint.h"
15 #include "fmpz.h"
16 #include "fmpz_vec.h"
17 #include "fmpz_poly.h"
18 #include "fmpz_mod_poly.h"
19
_fmpz_mod_poly_pow(fmpz * res,const fmpz * poly,slong len,ulong e,const fmpz_t p)20 void _fmpz_mod_poly_pow(fmpz *res, const fmpz *poly, slong len, ulong e,
21 const fmpz_t p)
22 {
23 ulong bit = ~((~UWORD(0)) >> 1);
24 slong rlen;
25 slong alloc = (slong) e * (len - 1) + 1;
26 fmpz *v = _fmpz_vec_init(alloc);
27 fmpz *R, *S, *T;
28
29 /*
30 Set bits to the bitmask with a 1 one place lower than the msb of e
31 */
32
33 while ((bit & e) == UWORD(0))
34 bit >>= 1;
35
36 bit >>= 1;
37
38 /*
39 Trial run without any polynomial arithmetic to determine the parity
40 of the number of swaps; then set R and S accordingly
41 */
42
43 {
44 unsigned int swaps = 0U;
45 ulong bit2 = bit;
46 if ((bit2 & e))
47 swaps = ~swaps;
48 while (bit2 >>= 1)
49 if ((bit2 & e) == UWORD(0))
50 swaps = ~swaps;
51
52 if (swaps == 0U)
53 {
54 R = res;
55 S = v;
56 }
57 else
58 {
59 R = v;
60 S = res;
61 }
62 }
63
64 /*
65 We unroll the first step of the loop, referring to {poly, len}
66 */
67
68 _fmpz_mod_poly_sqr(R, poly, len, p);
69 rlen = 2 * len - 1;
70 if ((bit & e))
71 {
72 _fmpz_mod_poly_mul(S, R, rlen, poly, len, p);
73 rlen += len - 1;
74 T = R;
75 R = S;
76 S = T;
77 }
78
79 while ((bit >>= 1))
80 {
81 if ((bit & e))
82 {
83 _fmpz_mod_poly_sqr(S, R, rlen, p);
84 rlen += rlen - 1;
85 _fmpz_mod_poly_mul(R, S, rlen, poly, len, p);
86 rlen += len - 1;
87 }
88 else
89 {
90 _fmpz_mod_poly_sqr(S, R, rlen, p);
91 rlen += rlen - 1;
92 T = R;
93 R = S;
94 S = T;
95 }
96 }
97
98 _fmpz_vec_clear(v, alloc);
99 }
100
fmpz_mod_poly_pow(fmpz_mod_poly_t rop,const fmpz_mod_poly_t op,ulong e)101 void fmpz_mod_poly_pow(fmpz_mod_poly_t rop, const fmpz_mod_poly_t op, ulong e)
102 {
103 const slong len = op->length;
104 slong rlen;
105
106 if ((len < 2) || (e < UWORD(3)))
107 {
108 if (e == UWORD(0))
109 fmpz_mod_poly_set_ui(rop, 1);
110 else if (len == 0)
111 fmpz_mod_poly_zero(rop);
112 else if (len == 1)
113 {
114 fmpz_mod_poly_fit_length(rop, 1);
115 fmpz_powm_ui(rop->coeffs, op->coeffs, e, &(rop->p));
116 _fmpz_mod_poly_set_length(rop, 1);
117 _fmpz_mod_poly_normalise(rop);
118 }
119 else if (e == UWORD(1))
120 fmpz_mod_poly_set(rop, op);
121 else /* e == UWORD(2) */
122 fmpz_mod_poly_sqr(rop, op);
123 return;
124 }
125
126 rlen = (slong) e * (len - 1) + 1;
127
128 if (rop != op)
129 {
130 fmpz_mod_poly_fit_length(rop, rlen);
131 _fmpz_mod_poly_pow(rop->coeffs, op->coeffs, len, e, &(rop->p));
132 _fmpz_mod_poly_set_length(rop, rlen);
133 }
134 else
135 {
136 fmpz *t = _fmpz_vec_init(rlen);
137
138 _fmpz_mod_poly_pow(t, op->coeffs, len, e, &(rop->p));
139
140 _fmpz_vec_clear(rop->coeffs, rop->alloc);
141 rop->coeffs = t;
142 rop->alloc = rlen;
143 rop->length = rlen;
144 }
145
146 _fmpz_mod_poly_normalise(rop);
147 }
148
149