1 /*
2 Copyright (C) 2012 Sebastian Pancratz
3 Copyright (C) 2012 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 "padic.h"
14
15 /*
16 Computes the sum
17 \begin{equation*}
18 (a-1)! x^{1-a} \sum_{i=a}^{b-1} \frac{x^i}{i!}.
19 \end{equation*}
20 in the rational $(T, Q)$.
21
22 Assumes that $1 \leq a < b$.
23
24 If $a + 1 = b$, sets $P = x$, $Q = a$, and $T = x$.
25 If $a + 2 = b$, sets $P = x^2$, $Q = a (a + 1)$, $T = x (a + 1) + x^2$.
26 In general, sets
27 \begin{align*}
28 P & = x^{b-a}, \\
29 Q & = \frac{(b-1)!}{(a-1)!}, \\
30 T & = (b-1)! x^{1-a} \sum_{i=a}^{b-1} \frac{x^i}{i!}.
31 \end{align*}
32 */
33
34 static void
_padic_exp_bsplit_series(fmpz_t P,fmpz_t Q,fmpz_t T,const fmpz_t x,slong a,slong b)35 _padic_exp_bsplit_series(fmpz_t P, fmpz_t Q, fmpz_t T,
36 const fmpz_t x, slong a, slong b)
37 {
38 if (b - a == 1)
39 {
40 fmpz_set(P, x);
41 fmpz_set_ui(Q, a);
42 fmpz_set(T, x);
43 }
44 else if (b - a == 2)
45 {
46 fmpz_mul(P, x, x);
47 fmpz_set_ui(Q, a);
48 fmpz_mul_ui(Q, Q, a + 1);
49 fmpz_mul_ui(T, x, a + 1);
50 fmpz_add(T, T, P);
51 }
52 else
53 {
54 const slong m = (a + b) / 2;
55
56 fmpz_t PR, QR, TR;
57
58 fmpz_init(PR);
59 fmpz_init(QR);
60 fmpz_init(TR);
61
62 _padic_exp_bsplit_series(P, Q, T, x, a, m);
63 _padic_exp_bsplit_series(PR, QR, TR, x, m, b);
64
65 fmpz_mul(T, T, QR);
66 fmpz_addmul(T, P, TR);
67 fmpz_mul(P, P, PR);
68 fmpz_mul(Q, Q, QR);
69
70 fmpz_clear(PR);
71 fmpz_clear(QR);
72 fmpz_clear(TR);
73 }
74 }
75
76 /*
77 Assumes that $x$ is such that $\exp(x)$ converges.
78
79 Assumes that $v = \ord_p(x)$ with $v < N$,
80 which also forces $N$ to positive.
81
82 The result $y$ might not be reduced modulo $p^N$.
83
84 Supports aliasing between $x$ and $y$.
85 */
86
87 static void
_padic_exp_bsplit(fmpz_t y,const fmpz_t x,slong v,const fmpz_t p,slong N)88 _padic_exp_bsplit(fmpz_t y, const fmpz_t x, slong v, const fmpz_t p, slong N)
89 {
90 const slong n = _padic_exp_bound(v, N, p);
91
92 if (n == 1)
93 {
94 fmpz_one(y);
95 }
96 else
97 {
98 fmpz_t P, Q, T;
99
100 fmpz_init(P);
101 fmpz_init(Q);
102 fmpz_init(T);
103
104 _padic_exp_bsplit_series(P, Q, T, x, 1, n);
105
106 fmpz_add(T, T, Q); /* (T,Q) := (T,Q) + 1 */
107
108 /* Note exp(x) is a unit so val(T) == val(Q). */
109 if (fmpz_remove(T, T, p))
110 fmpz_remove(Q, Q, p);
111
112 _padic_inv(Q, Q, p, N);
113 fmpz_mul(y, T, Q);
114
115 fmpz_clear(P);
116 fmpz_clear(Q);
117 fmpz_clear(T);
118 }
119 }
120
_padic_exp_balanced_2(fmpz_t rop,const fmpz_t xu,slong xv,slong N)121 void _padic_exp_balanced_2(fmpz_t rop, const fmpz_t xu, slong xv, slong N)
122 {
123 const fmpz_t p = {WORD(2)};
124
125 fmpz_t r, t;
126 slong w;
127
128 fmpz_init(r);
129 fmpz_init(t);
130
131 w = 1;
132
133 fmpz_mul_2exp(t, xu, xv);
134 fmpz_fdiv_r_2exp(t, t, N);
135
136 fmpz_one(rop);
137
138 while (!fmpz_is_zero(t))
139 {
140 fmpz_fdiv_r_2exp(r, t, 2*w);
141 fmpz_sub(t, t, r);
142
143 if (!fmpz_is_zero(r))
144 {
145 _padic_exp_bsplit(r, r, w, p, N);
146 fmpz_mul(rop, rop, r);
147 fmpz_fdiv_r_2exp(rop, rop, N);
148 }
149
150 w *= 2;
151 }
152
153 fmpz_clear(r);
154 fmpz_clear(t);
155 }
156
_padic_exp_balanced_p(fmpz_t rop,const fmpz_t xu,slong xv,const fmpz_t p,slong N)157 void _padic_exp_balanced_p(fmpz_t rop, const fmpz_t xu, slong xv,
158 const fmpz_t p, slong N)
159 {
160 fmpz_t r, t, pw, pN;
161 slong w;
162
163 fmpz_init(r);
164 fmpz_init(t);
165 fmpz_init(pw);
166 fmpz_init(pN);
167
168 fmpz_set(pw, p);
169 fmpz_pow_ui(pN, p, N);
170 w = 1;
171
172 fmpz_pow_ui(t, p, xv);
173 fmpz_mul(t, t, xu);
174 fmpz_mod(t, t, pN);
175
176 fmpz_one(rop);
177
178 while (!fmpz_is_zero(t))
179 {
180 fmpz_mul(pw, pw, pw);
181
182 fmpz_fdiv_r(r, t, pw);
183 fmpz_sub(t, t, r);
184
185 if (!fmpz_is_zero(r))
186 {
187 _padic_exp_bsplit(r, r, w, p, N);
188 fmpz_mul(rop, rop, r);
189 fmpz_mod(rop, rop, pN);
190 }
191
192 w *= 2;
193 }
194
195 fmpz_clear(r);
196 fmpz_clear(t);
197 fmpz_clear(pw);
198 fmpz_clear(pN);
199 }
200
201 /*
202 Assumes that the exponential series converges at $x \neq 0$,
203 and that $\ord_p(x) < N$.
204
205 Supports aliasing between $x$ and $y$.
206
207 TODO: Take advantage of additional factors of $p$ in $x$.
208 */
209
_padic_exp_balanced(fmpz_t rop,const fmpz_t u,slong v,const fmpz_t p,slong N)210 void _padic_exp_balanced(fmpz_t rop, const fmpz_t u, slong v,
211 const fmpz_t p, slong N)
212 {
213 if (fmpz_equal_ui(p, 2))
214 _padic_exp_balanced_2(rop, u, v, N);
215 else
216 _padic_exp_balanced_p(rop, u, v, p, N);
217 }
218
padic_exp_balanced(padic_t rop,const padic_t op,const padic_ctx_t ctx)219 int padic_exp_balanced(padic_t rop, const padic_t op, const padic_ctx_t ctx)
220 {
221 const slong N = padic_prec(rop);
222 const slong v = padic_val(op);
223 const fmpz *p = ctx->p;
224
225 if (padic_is_zero(op))
226 {
227 padic_one(rop);
228 return 1;
229 }
230
231 if ((fmpz_equal_ui(p, 2) && v <= 1) || (v <= 0))
232 {
233 return 0;
234 }
235 else
236 {
237 if (v < N)
238 {
239 _padic_exp_balanced(padic_unit(rop),
240 padic_unit(op), padic_val(op), p, N);
241 padic_val(rop) = 0;
242 }
243 else
244 {
245 padic_one(rop);
246 }
247 return 1;
248 }
249 }
250