Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /********************************************************************/
16 : /** **/
17 : /** TRANSCENDENTAL FUNCTIONS **/
18 : /** **/
19 : /********************************************************************/
20 : #include "pari.h"
21 : #include "paripriv.h"
22 :
23 : #define DEBUGLEVEL DEBUGLEVEL_trans
24 :
25 : #ifdef LONG_IS_64BIT
26 : static const long SQRTVERYBIGINT = 3037000500L; /* ceil(sqrt(LONG_MAX)) */
27 : #else
28 : static const long SQRTVERYBIGINT = 46341L;
29 : #endif
30 :
31 : static THREAD GEN gcatalan, geuler, glog2, gpi;
32 : void
33 337377 : pari_init_floats(void)
34 : {
35 337377 : gcatalan = geuler = gpi = zetazone = bernzone = glog2 = eulerzone = NULL;
36 337377 : }
37 :
38 : void
39 334920 : pari_close_floats(void)
40 : {
41 334920 : guncloneNULL(gcatalan);
42 333821 : guncloneNULL(geuler);
43 331537 : guncloneNULL(gpi);
44 330995 : guncloneNULL(glog2);
45 330664 : guncloneNULL(zetazone);
46 330360 : guncloneNULL_deep(bernzone);
47 330211 : guncloneNULL_deep(eulerzone);
48 329903 : }
49 :
50 : /********************************************************************/
51 : /** GENERIC BINARY SPLITTING **/
52 : /** (Haible, Papanikolaou) **/
53 : /********************************************************************/
54 : void
55 5579018 : abpq_init(struct abpq *A, long n)
56 : {
57 5579018 : A->a = (GEN*)new_chunk(n+1);
58 5586797 : A->b = (GEN*)new_chunk(n+1);
59 5590048 : A->p = (GEN*)new_chunk(n+1);
60 5590029 : A->q = (GEN*)new_chunk(n+1);
61 5590205 : }
62 : static GEN
63 189991858 : mulii3(GEN a, GEN b, GEN c) { return mulii(mulii(a,b),c); }
64 :
65 : /* T_{n1,n1+1} */
66 : static GEN
67 42950489 : T2(struct abpq *A, long n1)
68 : {
69 42950489 : GEN u = mulii3(A->a[n1], A->b[n1+1], A->q[n1+1]);
70 42990577 : GEN v = mulii3(A->b[n1], A->a[n1+1], A->p[n1+1]);
71 42987566 : return mulii(A->p[n1], addii(u, v));
72 : }
73 :
74 : /* assume n2 > n1. Compute sum_{n1 <= n < n2} a/b(n) p/q(n1)... p/q(n) */
75 : void
76 81050817 : abpq_sum(struct abpq_res *r, long n1, long n2, struct abpq *A)
77 : {
78 : struct abpq_res L, R;
79 : GEN u1, u2;
80 : pari_sp av;
81 : long n;
82 81050817 : switch(n2 - n1)
83 : {
84 : GEN b, q;
85 63 : case 1:
86 63 : r->P = A->p[n1];
87 63 : r->Q = A->q[n1];
88 63 : r->B = A->b[n1];
89 63 : r->T = mulii(A->a[n1], A->p[n1]);
90 43078723 : return;
91 26012203 : case 2:
92 26012203 : r->P = mulii(A->p[n1], A->p[n1+1]);
93 25882765 : r->Q = mulii(A->q[n1], A->q[n1+1]);
94 25870503 : r->B = mulii(A->b[n1], A->b[n1+1]);
95 25872784 : av = avma;
96 25872784 : r->T = gc_INT(av, T2(A, n1));
97 25886465 : return;
98 :
99 17394707 : case 3:
100 17394707 : q = mulii(A->q[n1+1], A->q[n1+2]);
101 17345858 : b = mulii(A->b[n1+1], A->b[n1+2]);
102 17338827 : r->P = mulii3(A->p[n1], A->p[n1+1], A->p[n1+2]);
103 17341387 : r->Q = mulii(A->q[n1], q);
104 17338343 : r->B = mulii(A->b[n1], b);
105 17341785 : av = avma;
106 17341785 : u1 = mulii3(b, q, A->a[n1]);
107 17325656 : u2 = mulii(A->b[n1], T2(A, n1+1));
108 17341286 : r->T = gc_INT(av, mulii(A->p[n1], addii(u1, u2)));
109 17192195 : return;
110 : }
111 :
112 37643844 : av = avma;
113 37643844 : n = (n1 + n2) >> 1;
114 37643844 : abpq_sum(&L, n1, n, A);
115 37696350 : abpq_sum(&R, n, n2, A);
116 :
117 37723899 : r->P = mulii(L.P, R.P);
118 37555665 : r->Q = mulii(L.Q, R.Q);
119 37527051 : r->B = mulii(L.B, R.B);
120 37547341 : u1 = mulii3(R.B, R.Q, L.T);
121 37485852 : u2 = mulii3(L.B, L.P, R.T);
122 37504176 : r->T = addii(u1,u2);
123 37509278 : set_avma(av);
124 37490008 : r->P = icopy(r->P);
125 37732105 : r->Q = icopy(r->Q);
126 37774974 : r->B = icopy(r->B);
127 37776335 : r->T = icopy(r->T);
128 : }
129 :
130 : /********************************************************************/
131 : /** **/
132 : /** PI **/
133 : /** **/
134 : /********************************************************************/
135 : /* replace *old clone by c. Protect against SIGINT */
136 : static void
137 83983 : swap_clone(GEN *old, GEN c)
138 83983 : { GEN tmp = *old; *old = c; guncloneNULL(tmp); }
139 :
140 : /* ----
141 : * 53360 (640320)^(1/2) \ (6n)! (545140134 n + 13591409)
142 : * -------------------- = / ------------------------------
143 : * Pi ---- (n!)^3 (3n)! (-640320)^(3n)
144 : * n>=0
145 : *
146 : * Ramanujan's formula + binary splitting */
147 : static GEN
148 41360 : pi_ramanujan(long prec)
149 : {
150 41360 : const ulong B = 545140134, A = 13591409, C = 640320;
151 41360 : const double alpha2 = 47.11041314; /* 3log(C/12) / log(2) */
152 : long n, nmax, prec2;
153 : struct abpq_res R;
154 : struct abpq S;
155 : GEN D, u;
156 :
157 41360 : nmax = (long)(1 + prec2nbits(prec)/alpha2);
158 : #ifdef LONG_IS_64BIT
159 40833 : D = utoipos(10939058860032000UL); /* C^3/24 */
160 : #else
161 527 : D = uutoi(2546948UL,495419392UL);
162 : #endif
163 41359 : abpq_init(&S, nmax);
164 41362 : S.a[0] = utoipos(A);
165 41361 : S.b[0] = S.p[0] = S.q[0] = gen_1;
166 323916 : for (n = 1; n <= nmax; n++)
167 : {
168 282553 : S.a[n] = addiu(muluu(B, n), A);
169 282463 : S.b[n] = gen_1;
170 282463 : S.p[n] = mulis(muluu(6*n-5, 2*n-1), 1-6*n);
171 282534 : S.q[n] = mulii(sqru(n), muliu(D,n));
172 : }
173 41363 : abpq_sum(&R, 0, nmax, &S); prec2 = prec+EXTRAPREC64;
174 41370 : u = itor(muliu(R.Q,C/12), prec2);
175 41366 : return rtor(mulrr(divri(u, R.T), sqrtr_abs(utor(C,prec2))), prec);
176 : }
177 :
178 : #if 0 /* Much slower than binary splitting at least up to prec = 10^8 */
179 : /* Gauss - Brent-Salamin AGM iteration */
180 : static GEN
181 : pi_brent_salamin(long prec)
182 : {
183 : GEN A, B, C;
184 : pari_sp av2;
185 : long i, G;
186 :
187 : G = - prec2nbits(prec);
188 : incrprec(prec);
189 :
190 : A = real2n(-1, prec);
191 : B = sqrtr_abs(A); /* = 1/sqrt(2) */
192 : setexpo(A, 0);
193 : C = real2n(-2, prec); av2 = avma;
194 : for (i = 0;; i++)
195 : {
196 : GEN y, a, b, B_A = subrr(B, A);
197 : pari_sp av3 = avma;
198 : if (expo(B_A) < G) break;
199 : a = addrr(A,B); shiftr_inplace(a, -1);
200 : b = mulrr(A,B);
201 : affrr(a, A);
202 : affrr(sqrtr_abs(b), B); set_avma(av3);
203 : y = sqrr(B_A); shiftr_inplace(y, i - 2);
204 : affrr(subrr(C, y), C); set_avma(av2);
205 : }
206 : shiftr_inplace(C, 2);
207 : return divrr(sqrr(addrr(A,B)), C);
208 : }
209 : #endif
210 :
211 : GEN
212 46172803 : constpi(long prec)
213 : {
214 : pari_sp av;
215 : GEN tmp;
216 46172803 : if (gpi && realprec(gpi) >= prec) return gpi;
217 :
218 41119 : av = avma;
219 41119 : tmp = gclone(pi_ramanujan(prec));
220 41378 : swap_clone(&gpi,tmp);
221 41374 : return gc_const(av, gpi);
222 : }
223 :
224 : GEN
225 46172583 : mppi(long prec) { return rtor(constpi(prec), prec); }
226 :
227 : /* Pi * 2^n */
228 : GEN
229 31283676 : Pi2n(long n, long prec)
230 : {
231 31283676 : GEN x = mppi(prec); shiftr_inplace(x, n);
232 31282975 : return x;
233 : }
234 :
235 : /* I * Pi * 2^n */
236 : GEN
237 280449 : PiI2n(long n, long prec) { retmkcomplex(gen_0, Pi2n(n, prec)); }
238 :
239 : /* 2I * Pi */
240 : GEN
241 267723 : PiI2(long prec) { return PiI2n(1, prec); }
242 :
243 : /********************************************************************/
244 : /** **/
245 : /** EULER CONSTANT **/
246 : /** **/
247 : /********************************************************************/
248 :
249 : GEN
250 57681 : consteuler(long prec)
251 : {
252 : GEN u,v,a,b,tmpeuler;
253 : long l, n1, n, k, x;
254 : pari_sp av1, av2;
255 :
256 57681 : if (geuler && realprec(geuler) >= prec) return geuler;
257 :
258 508 : av1 = avma; tmpeuler = cgetr_block(prec);
259 :
260 508 : incrprec(prec);
261 :
262 508 : l = prec+EXTRAPREC64; x = (long) (1 + prec2nbits_mul(l, M_LN2/4));
263 508 : a = utor(x,l); u=logr_abs(a); setsigne(u,-1); affrr(u,a);
264 508 : b = real_1(l);
265 508 : v = real_1(l);
266 508 : n = (long)(1+3.591*x); /* z=3.591: z*[ ln(z)-1 ]=1 */
267 508 : n1 = minss(n, SQRTVERYBIGINT);
268 508 : if (x < SQRTVERYBIGINT)
269 : {
270 508 : ulong xx = x*x;
271 508 : av2 = avma;
272 168189 : for (k=1; k<n1; k++)
273 : {
274 167681 : affrr(divru(mulur(xx,b),k*k), b);
275 167689 : affrr(divru(addrr(divru(mulur(xx,a),k),b),k), a);
276 167695 : affrr(addrr(u,a), u);
277 167685 : affrr(addrr(v,b), v); set_avma(av2);
278 : }
279 1016 : for ( ; k<=n; k++)
280 : {
281 508 : affrr(divru(divru(mulur(xx,b),k),k), b);
282 508 : affrr(divru(addrr(divru(mulur(xx,a),k),b),k), a);
283 508 : affrr(addrr(u,a), u);
284 508 : affrr(addrr(v,b), v); set_avma(av2);
285 : }
286 : }
287 : else
288 : {
289 0 : GEN xx = sqru(x);
290 0 : av2 = avma;
291 0 : for (k=1; k<n1; k++)
292 : {
293 0 : affrr(divru(mulir(xx,b),k*k), b);
294 0 : affrr(divru(addrr(divru(mulir(xx,a),k),b),k), a);
295 0 : affrr(addrr(u,a), u);
296 0 : affrr(addrr(v,b), v); set_avma(av2);
297 : }
298 0 : for ( ; k<=n; k++)
299 : {
300 0 : affrr(divru(divru(mulir(xx,b),k),k), b);
301 0 : affrr(divru(addrr(divru(mulir(xx,a),k),b),k), a);
302 0 : affrr(addrr(u,a), u);
303 0 : affrr(addrr(v,b), v); set_avma(av2);
304 : }
305 : }
306 508 : affrr(divrr(u,v),tmpeuler);
307 508 : swap_clone(&geuler,tmpeuler);
308 508 : return gc_const(av1, geuler);
309 : }
310 :
311 : GEN
312 57681 : mpeuler(long prec) { return rtor(consteuler(prec), prec); }
313 :
314 : /********************************************************************/
315 : /** **/
316 : /** CATALAN CONSTANT **/
317 : /** **/
318 : /********************************************************************/
319 : /* inf 256^i (580i^2 - 184i + 15) (2i)!^3 (3i)!^2
320 : * 64 G = SUM ------------------------------------------
321 : * i=1 i^3 (2i-1) (6i)!^2 */
322 : static GEN
323 14 : catalan(long prec)
324 : {
325 14 : long i, nmax = 1 + prec2nbits(prec) / 7.509; /* / log2(729/4) */
326 : struct abpq_res R;
327 : struct abpq A;
328 : GEN u;
329 14 : abpq_init(&A, nmax);
330 14 : A.a[0] = gen_0; A.b[0] = A.p[0] = A.q[0] = gen_1;
331 1750 : for (i = 1; i <= nmax; i++)
332 : {
333 1736 : A.a[i] = addiu(muluu(580*i - 184, i), 15);
334 1736 : A.b[i] = muliu(powuu(i, 3), 2*i - 1);
335 1736 : A.p[i] = mului(64*i-32, powuu(i,3));
336 1736 : A.q[i] = sqri(muluu(6*i - 1, 18*i - 15));
337 : }
338 14 : abpq_sum(&R, 0, nmax, &A);
339 14 : u = rdivii(R.T, mulii(R.B,R.Q),prec);
340 14 : shiftr_inplace(u, -6); return u;
341 : }
342 :
343 : GEN
344 14 : constcatalan(long prec)
345 : {
346 14 : pari_sp av = avma;
347 : GEN tmp;
348 14 : if (gcatalan && realprec(gcatalan) >= prec) return gcatalan;
349 14 : tmp = gclone(catalan(prec));
350 14 : swap_clone(&gcatalan,tmp);
351 14 : return gc_const(av, gcatalan);
352 : }
353 :
354 : GEN
355 14 : mpcatalan(long prec) { return rtor(constcatalan(prec), prec); }
356 :
357 : /********************************************************************/
358 : /** **/
359 : /** TYPE CONVERSION FOR TRANSCENDENTAL FUNCTIONS **/
360 : /** **/
361 : /********************************************************************/
362 : static GEN
363 2006343 : transvec(GEN (*f)(GEN,long), GEN x, long prec)
364 6512824 : { pari_APPLY_same(f(gel(x,i), prec)); }
365 : static GEN
366 329 : transvecgen(void *E, GEN (*f)(void *,GEN,long), GEN x, long prec)
367 735 : { pari_APPLY_same(f(E, gel(x,i), prec)); }
368 :
369 : GEN
370 3927064 : trans_eval(const char *fun, GEN (*f)(GEN,long), GEN x, long prec)
371 : {
372 3927064 : pari_sp av = avma;
373 3927064 : if (prec < LOWDEFAULTPREC) pari_err_BUG("trans_eval [prec < 3]");
374 3927090 : switch(typ(x))
375 : {
376 1663186 : case t_INT: x = f(itor(x,prec),prec); break;
377 257502 : case t_FRAC: x = f(fractor(x, prec),prec); break;
378 7 : case t_QUAD: x = f(quadtofp(x,prec),prec); break;
379 14 : case t_POLMOD: x = transvec(f, polmod_to_embed(x,prec), prec); break;
380 2006332 : case t_VEC:
381 : case t_COL:
382 2006332 : case t_MAT: return transvec(f, x, prec);
383 49 : default: pari_err_TYPE(fun,x);
384 : return NULL;/*LCOV_EXCL_LINE*/
385 : }
386 1920696 : return gc_upto(av, x);
387 : }
388 :
389 : GEN
390 1967 : trans_evalgen(const char *fun, void *E, GEN (*f)(void*,GEN,long),
391 : GEN x, long prec)
392 : {
393 1967 : pari_sp av = avma;
394 1967 : if (prec < LOWDEFAULTPREC) pari_err_BUG("trans_eval [prec < 3]");
395 1967 : switch(typ(x))
396 : {
397 343 : case t_INT: x = f(E, itor(x,prec),prec); break;
398 1260 : case t_FRAC: x = f(E, fractor(x, prec),prec); break;
399 0 : case t_QUAD: x = f(E, quadtofp(x,prec),prec); break;
400 70 : case t_POLMOD: x = transvecgen(E, f, polmod_to_embed(x,prec), prec); break;
401 259 : case t_VEC:
402 : case t_COL:
403 259 : case t_MAT: return transvecgen(E, f, x, prec);
404 35 : default: pari_err_TYPE(fun,x);
405 : return NULL;/*LCOV_EXCL_LINE*/
406 : }
407 1673 : return gc_upto(av, x);
408 : }
409 :
410 : /*******************************************************************/
411 : /* */
412 : /* POWERING */
413 : /* */
414 : /*******************************************************************/
415 : /* x a t_REAL 0, return exp(x) */
416 : static GEN
417 147118 : mpexp0(GEN x)
418 : {
419 147118 : long e = expo(x);
420 147118 : return e >= 0? real_0_bit(e): real_1_bit(-e);
421 : }
422 : static GEN
423 21140 : powr0(GEN x)
424 21140 : { return signe(x)? real_1(realprec(x)): mpexp0(x); }
425 :
426 : /* assume typ(x) = t_VEC */
427 : static int
428 49 : is_ext_qfr(GEN x)
429 35 : { return lg(x) == 3 && typ(gel(x,1)) == t_QFB && !qfb_is_qfi(gel(x,1))
430 84 : && typ(gel(x,2)) == t_REAL; }
431 :
432 : /* x t_POL or t_SER, return scalarpol(Rg_get_1(x)) */
433 : static GEN
434 375055 : scalarpol_get_1(GEN x)
435 : {
436 375055 : GEN y = cgetg(3,t_POL);
437 375055 : y[1] = evalvarn(varn(x)) | evalsigne(1);
438 375055 : gel(y,2) = Rg_get_1(x); return y;
439 : }
440 : /* to be called by the generic function gpowgs(x,s) when s = 0 */
441 : static GEN
442 2725343 : gpowg0(GEN x)
443 : {
444 : long lx, i;
445 : GEN y;
446 :
447 2725343 : switch(typ(x))
448 : {
449 2304587 : case t_INT: case t_REAL: case t_FRAC: case t_PADIC:
450 2304587 : return gen_1;
451 :
452 7 : case t_QUAD: x++; /*fall through*/
453 38107 : case t_COMPLEX: {
454 38107 : pari_sp av = avma;
455 38107 : GEN a = gpowg0(gel(x,1));
456 38107 : GEN b = gpowg0(gel(x,2));
457 38107 : if (a == gen_1) return b;
458 14 : if (b == gen_1) return a;
459 7 : return gc_upto(av, gmul(a,b));
460 : }
461 133 : case t_INTMOD:
462 133 : y = cgetg(3,t_INTMOD);
463 133 : gel(y,1) = icopy(gel(x,1));
464 133 : gel(y,2) = is_pm1(gel(x,1))? gen_0: gen_1;
465 133 : return y;
466 :
467 7287 : case t_FFELT: return FF_1(x);
468 :
469 1550 : case t_POLMOD:
470 1550 : retmkpolmod(scalarpol_get_1(gel(x,1)), gcopy(gel(x,1)));
471 :
472 28 : case t_RFRAC:
473 28 : return scalarpol_get_1(gel(x,2));
474 373477 : case t_POL: case t_SER:
475 373477 : return scalarpol_get_1(x);
476 :
477 84 : case t_MAT:
478 84 : lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
479 77 : if (lx != lgcols(x)) pari_err_DIM("gpow");
480 77 : y = matid(lx-1);
481 252 : for (i=1; i<lx; i++) gcoeff(y,i,i) = gpowg0(gcoeff(x,i,i));
482 77 : return y;
483 21 : case t_VEC: if (!is_ext_qfr(x)) break;
484 : /* fall through handle extended t_QFB */
485 28 : case t_QFB: return qfbpow(x, gen_0);
486 49 : case t_VECSMALL: return identity_perm(lg(x) - 1);
487 : }
488 13 : pari_err_TYPE("gpow",x);
489 : return NULL; /* LCOV_EXCL_LINE */
490 : }
491 :
492 : static GEN
493 6340139 : _sqr(void *data /* ignored */, GEN x) { (void)data; return gsqr(x); }
494 : static GEN
495 4111224 : _mul(void *data /* ignored */, GEN x, GEN y) { (void)data; return gmul(x,y); }
496 : static GEN
497 782753 : _one(void *x) { return gpowg0((GEN) x); }
498 : static GEN
499 81886684 : _sqri(void *data /* ignored */, GEN x) { (void)data; return sqri(x); }
500 : static GEN
501 29990124 : _muli(void *data /* ignored */, GEN x, GEN y) { (void)data; return mulii(x,y); }
502 : static GEN
503 16297186 : _sqrr(void *data /* ignored */, GEN x) { (void)data; return sqrr(x); }
504 : static GEN
505 7128563 : _mulr(void *data /* ignored */, GEN x, GEN y) { (void)data; return mulrr(x,y); }
506 : static GEN
507 14196 : _oner(void *data /* prec */) { return real_1( *(long*) data); }
508 :
509 : /* INTEGER POWERING (a^n for integer a != 0 and integer n > 0)
510 : *
511 : * Use left shift binary algorithm (RS is wasteful: multiplies big numbers,
512 : * with LS one of them is the base, hence small). Sign of result is set
513 : * to s (= 1,-1). Makes life easier for caller, which otherwise might do a
514 : * setsigne(gen_1 / gen_m1) */
515 : static GEN
516 118706850 : powiu_sign(GEN a, ulong N, long s)
517 : {
518 : pari_sp av;
519 : GEN y;
520 :
521 118706850 : if (lgefint(a) == 3)
522 : { /* easy if |a| < 3 */
523 117149401 : ulong q = a[2];
524 117149401 : if (q == 1) return (s>0)? gen_1: gen_m1;
525 100818026 : if (q == 2) { a = int2u(N); setsigne(a,s); return a; }
526 75361852 : q = upowuu(q, N);
527 75364694 : if (q) return s>0? utoipos(q): utoineg(q);
528 : }
529 32704220 : if (N <= 2) {
530 1830978 : if (N == 2) return sqri(a);
531 18658 : a = icopy(a); setsigne(a,s); return a;
532 : }
533 30873242 : av = avma;
534 30873242 : y = gen_powu_i(a, N, NULL, &_sqri, &_muli);
535 30873485 : setsigne(y,s); return gc_INT(av, y);
536 : }
537 : /* a^n */
538 : GEN
539 118555700 : powiu(GEN a, ulong n)
540 : {
541 : long s;
542 118555700 : if (!n) return gen_1;
543 117404028 : s = signe(a);
544 117404028 : if (!s) return gen_0;
545 117322885 : return powiu_sign(a, n, (s < 0 && odd(n))? -1: 1);
546 : }
547 : GEN
548 28734880 : powis(GEN a, long n)
549 : {
550 : long s;
551 : GEN t, y;
552 28734880 : if (n >= 0) return powiu(a, n);
553 630418 : s = signe(a);
554 630418 : if (!s) pari_err_INV("powis",gen_0);
555 630418 : t = (s < 0 && odd(n))? gen_m1: gen_1;
556 630418 : if (is_pm1(a)) return t;
557 : /* n < 0, |a| > 1 */
558 627898 : y = cgetg(3,t_FRAC);
559 627899 : gel(y,1) = t;
560 627899 : gel(y,2) = powiu_sign(a, -n, 1); /* force denominator > 0 */
561 627896 : return y;
562 : }
563 : GEN
564 45322605 : powuu(ulong p, ulong N)
565 : {
566 : pari_sp av;
567 : ulong pN;
568 : GEN y;
569 45322605 : if (!p) return gen_0;
570 45322528 : if (N <= 2)
571 : {
572 39959258 : if (N == 2) return sqru(p);
573 37900746 : if (N == 1) return utoipos(p);
574 4931081 : return gen_1;
575 : }
576 5363270 : pN = upowuu(p, N);
577 5363284 : if (pN) return utoipos(pN);
578 986431 : if (p == 2) return int2u(N);
579 973011 : av = avma;
580 973011 : y = gen_powu_i(utoipos(p), N, NULL, &_sqri, &_muli);
581 973008 : return gc_INT(av, y);
582 : }
583 :
584 : /* return 0 if overflow */
585 : static ulong
586 21932879 : usqru(ulong p) { return p & HIGHMASK? 0: p*p; }
587 : ulong
588 112357452 : upowuu(ulong p, ulong k)
589 : {
590 : #ifdef LONG_IS_64BIT
591 96378117 : const ulong CUTOFF3 = 2642245;
592 96378117 : const ulong CUTOFF4 = 65535;
593 96378117 : const ulong CUTOFF5 = 7131;
594 96378117 : const ulong CUTOFF6 = 1625;
595 96378117 : const ulong CUTOFF7 = 565;
596 96378117 : const ulong CUTOFF8 = 255;
597 96378117 : const ulong CUTOFF9 = 138;
598 96378117 : const ulong CUTOFF10 = 84;
599 96378117 : const ulong CUTOFF11 = 56;
600 96378117 : const ulong CUTOFF12 = 40;
601 96378117 : const ulong CUTOFF13 = 30;
602 96378117 : const ulong CUTOFF14 = 23;
603 96378117 : const ulong CUTOFF15 = 19;
604 96378117 : const ulong CUTOFF16 = 15;
605 96378117 : const ulong CUTOFF17 = 13;
606 96378117 : const ulong CUTOFF18 = 11;
607 96378117 : const ulong CUTOFF19 = 10;
608 96378117 : const ulong CUTOFF20 = 9;
609 : #else
610 15979335 : const ulong CUTOFF3 = 1625;
611 15979335 : const ulong CUTOFF4 = 255;
612 15979335 : const ulong CUTOFF5 = 84;
613 15979335 : const ulong CUTOFF6 = 40;
614 15979335 : const ulong CUTOFF7 = 23;
615 15979335 : const ulong CUTOFF8 = 15;
616 15979335 : const ulong CUTOFF9 = 11;
617 15979335 : const ulong CUTOFF10 = 9;
618 15979335 : const ulong CUTOFF11 = 7;
619 15979335 : const ulong CUTOFF12 = 6;
620 15979335 : const ulong CUTOFF13 = 5;
621 15979335 : const ulong CUTOFF14 = 4;
622 15979335 : const ulong CUTOFF15 = 4;
623 15979335 : const ulong CUTOFF16 = 3;
624 15979335 : const ulong CUTOFF17 = 3;
625 15979335 : const ulong CUTOFF18 = 3;
626 15979335 : const ulong CUTOFF19 = 3;
627 15979335 : const ulong CUTOFF20 = 3;
628 : #endif
629 :
630 112357452 : if (p <= 2)
631 : {
632 9729564 : if (p < 2) return p;
633 9173217 : return k < BITS_IN_LONG? 1UL<<k: 0;
634 : }
635 102627888 : switch(k)
636 : {
637 : ulong p2, p3, p4, p5, p8;
638 8618538 : case 0: return 1;
639 26792571 : case 1: return p;
640 21932875 : case 2: return usqru(p);
641 4160589 : case 3: if (p > CUTOFF3) return 0; return p*p*p;
642 11544354 : case 4: if (p > CUTOFF4) return 0; p2=p*p; return p2*p2;
643 2334703 : case 5: if (p > CUTOFF5) return 0; p2=p*p; return p2*p2*p;
644 6982854 : case 6: if (p > CUTOFF6) return 0; p2=p*p; return p2*p2*p2;
645 614806 : case 7: if (p > CUTOFF7) return 0; p2=p*p; return p2*p2*p2*p;
646 914924 : case 8: if (p > CUTOFF8) return 0; p2=p*p; p4=p2*p2; return p4*p4;
647 650532 : case 9: if (p > CUTOFF9) return 0; p2=p*p; p4=p2*p2; return p4*p4*p;
648 4992183 : case 10: if (p > CUTOFF10)return 0; p2=p*p; p4=p2*p2; return p4*p4*p2;
649 376634 : case 11: if (p > CUTOFF11)return 0; p2=p*p; p4=p2*p2; return p4*p4*p2*p;
650 4795753 : case 12: if (p > CUTOFF12)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4;
651 107190 : case 13: if (p > CUTOFF13)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4*p;
652 4753481 : case 14: if (p > CUTOFF14)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4*p2;
653 165778 : case 15: if (p > CUTOFF15)return 0;
654 105498 : p2=p*p; p3=p2*p; p5=p3*p2; return p5*p5*p5;
655 107303 : case 16: if (p > CUTOFF16)return 0;
656 53639 : p2=p*p; p4=p2*p2; p8=p4*p4; return p8*p8;
657 80833 : case 17: if (p > CUTOFF17)return 0;
658 42143 : p2=p*p; p4=p2*p2; p8=p4*p4; return p*p8*p8;
659 70253 : case 18: if (p > CUTOFF18)return 0;
660 39403 : p2=p*p; p4=p2*p2; p8=p4*p4; return p2*p8*p8;
661 827289 : case 19: if (p > CUTOFF19)return 0;
662 773004 : p2=p*p; p4=p2*p2; p8=p4*p4; return p*p2*p8*p8;
663 81860 : case 20: if (p > CUTOFF20)return 0;
664 39293 : p2=p*p; p4=p2*p2; p8=p4*p4; return p4*p8*p8;
665 : }
666 : #ifdef LONG_IS_64BIT
667 1500296 : switch(p)
668 : {
669 221587 : case 3: if (k > 40) return 0;
670 161687 : break;
671 17034 : case 4: if (k > 31) return 0;
672 780 : return 1UL<<(2*k);
673 637471 : case 5: if (k > 27) return 0;
674 20383 : break;
675 49650 : case 6: if (k > 24) return 0;
676 9180 : break;
677 56696 : case 7: if (k > 22) return 0;
678 3074 : break;
679 517858 : default: return 0;
680 : }
681 : /* no overflow */
682 : {
683 194324 : ulong q = upowuu(p, k >> 1);
684 194324 : q *= q ;
685 194324 : return odd(k)? q*p: q;
686 : }
687 : #else
688 222289 : return 0;
689 : #endif
690 : }
691 :
692 : GEN
693 12017 : upowers(ulong x, long n)
694 : {
695 : long i;
696 12017 : GEN p = cgetg(n + 2, t_VECSMALL);
697 12017 : uel(p,1) = 1; if (n == 0) return p;
698 12017 : uel(p,2) = x;
699 91465 : for (i = 3; i <= n; i++)
700 79448 : uel(p,i) = uel(p,i-1)*x;
701 12017 : return p;
702 : }
703 :
704 : typedef struct {
705 : long prec, a;
706 : GEN (*sqr)(GEN);
707 : GEN (*mulug)(ulong,GEN);
708 : } sr_muldata;
709 :
710 : static GEN
711 1620028 : _rpowuu_sqr(void *data, GEN x)
712 : {
713 1620028 : sr_muldata *D = (sr_muldata *)data;
714 1620028 : if (typ(x) == t_INT && lg2prec(lgefint(x)) >= D->prec)
715 : { /* switch to t_REAL */
716 158021 : D->sqr = &sqrr;
717 158021 : D->mulug = &mulur; x = itor(x, D->prec);
718 : }
719 1620028 : return D->sqr(x);
720 : }
721 :
722 : static GEN
723 628884 : _rpowuu_msqr(void *data, GEN x)
724 : {
725 628884 : GEN x2 = _rpowuu_sqr(data, x);
726 628884 : sr_muldata *D = (sr_muldata *)data;
727 628884 : return D->mulug(D->a, x2);
728 : }
729 :
730 : /* return a^n as a t_REAL of precision prec. Assume a > 0, n > 0 */
731 : GEN
732 445493 : rpowuu(ulong a, ulong n, long prec)
733 : {
734 : pari_sp av;
735 : GEN y, z;
736 : sr_muldata D;
737 :
738 445493 : if (a == 1) return real_1(prec);
739 445493 : if (a == 2) return real2n(n, prec);
740 445493 : if (n == 1) return utor(a, prec);
741 440338 : z = cgetr(prec);
742 440338 : av = avma;
743 440338 : D.sqr = &sqri;
744 440338 : D.mulug = &mului;
745 440338 : D.prec = prec;
746 440338 : D.a = (long)a;
747 440338 : y = gen_powu_fold_i(utoipos(a), n, (void*)&D, &_rpowuu_sqr, &_rpowuu_msqr);
748 440338 : affgr(y, z); return gc_const(av,z);
749 : }
750 :
751 : GEN
752 5095273 : powrs(GEN x, long n)
753 : {
754 5095273 : pari_sp av = avma;
755 : GEN y;
756 5095273 : if (!n) return powr0(x);
757 5095273 : y = gen_powu_i(x, (ulong)labs(n), NULL, &_sqrr, &_mulr);
758 5095616 : if (n < 0) y = invr(y);
759 5095548 : return gc_uptoleaf(av,y);
760 : }
761 : GEN
762 6137160 : powru(GEN x, ulong n)
763 : {
764 6137160 : pari_sp av = avma;
765 : GEN y;
766 6137160 : if (!n) return powr0(x);
767 6116538 : y = gen_powu_i(x, n, NULL, &_sqrr, &_mulr);
768 6116460 : return gc_uptoleaf(av,y);
769 : }
770 :
771 : GEN
772 14196 : powersr(GEN x, long n)
773 : {
774 14196 : long prec = realprec(x);
775 14196 : return gen_powers(x, n, 1, &prec, &_sqrr, &_mulr, &_oner);
776 : }
777 :
778 : /* x^(s/2), assume x t_REAL */
779 : GEN
780 0 : powrshalf(GEN x, long s)
781 : {
782 0 : if (s & 1) return sqrtr(powrs(x, s));
783 0 : return powrs(x, s>>1);
784 : }
785 : /* x^(s/2), assume x t_REAL */
786 : GEN
787 120224 : powruhalf(GEN x, ulong s)
788 : {
789 120224 : if (s & 1) return sqrtr(powru(x, s));
790 8023 : return powru(x, s>>1);
791 : }
792 : /* x^(n/d), assume x t_REAL, return t_REAL */
793 : GEN
794 518 : powrfrac(GEN x, long n, long d)
795 : {
796 : long z;
797 518 : if (!n) return powr0(x);
798 0 : z = cgcd(n, d); if (z > 1) { n /= z; d /= z; }
799 0 : if (d == 1) return powrs(x, n);
800 0 : x = powrs(x, n);
801 0 : if (d == 2) return sqrtr(x);
802 0 : return sqrtnr(x, d);
803 : }
804 :
805 : /* assume x != 0 */
806 : static GEN
807 641269 : pow_monome(GEN x, long n)
808 : {
809 641269 : long i, d, dx = degpol(x);
810 : GEN A, b, y;
811 :
812 641269 : if (n < 0) { n = -n; y = cgetg(3, t_RFRAC); } else y = NULL;
813 :
814 641269 : if (HIGHWORD(dx) || HIGHWORD(n))
815 8 : {
816 : LOCAL_HIREMAINDER;
817 9 : d = (long)mulll((ulong)dx, (ulong)n);
818 9 : if (hiremainder || (d &~ LGBITS)) d = LGBITS; /* overflow */
819 9 : d += 2;
820 : }
821 : else
822 641260 : d = dx*n + 2;
823 641269 : if ((d + 1) & ~LGBITS) pari_err(e_OVERFLOW,"pow_monome [degree]");
824 641262 : A = cgetg(d+1, t_POL); A[1] = x[1];
825 5999995 : for (i=2; i < d; i++) gel(A,i) = gen_0;
826 641262 : b = gpowgs(gel(x,dx+2), n); /* not memory clean if (n < 0) */
827 641262 : if (!y) y = A;
828 : else {
829 20482 : GEN c = denom_i(b);
830 20482 : gel(y,1) = c; if (c != gen_1) b = gmul(b,c);
831 20482 : gel(y,2) = A;
832 : }
833 641262 : gel(A,d) = b; return y;
834 : }
835 :
836 : /* q*p^v */
837 : static GEN
838 1317010 : mulpowu(GEN q, GEN p, ulong v)
839 : {
840 1317010 : pari_sp av = avma;
841 1317010 : if (v == 0) return icopy(q);
842 86758 : return gc_INT(av, mulii(q, powiu(p,v)));
843 : }
844 :
845 : /* x t_PADIC, n != 0 */
846 : static GEN
847 1316945 : powps(GEN x, long n)
848 : {
849 1316945 : long e = valp(x), v;
850 1316945 : GEN p = padic_p(x);
851 :
852 1316945 : if (e)
853 : {
854 658 : pari_sp av = avma;
855 658 : e = itos_or_0(mulss(e, n));
856 658 : if (!e) pari_err_OVERFLOW("valp()");
857 658 : set_avma(av);
858 : }
859 1316945 : if (!signe(padic_u(x))) {
860 84 : if (n < 0) pari_err_INV("powps",x);
861 77 : return zeropadic(p, e);
862 : }
863 1316861 : v = z_pval(n, p);
864 1316863 : if (v && precp(x) == 1 && absequaliu(p, 2)) v++;
865 1316863 : retmkpadic_i(Fp_pows(padic_u(x), n, _pd), icopy(p),
866 : mulpowu(padic_pd(x), p, v), e, precp(x) + v);
867 : }
868 : /* x t_PADIC */
869 : static GEN
870 161 : powp(GEN x, GEN n)
871 : {
872 161 : GEN p = padic_p(x);
873 : long v;
874 :
875 161 : if (valp(x)) pari_err_OVERFLOW("valp()");
876 161 : if (!signe(padic_u(x))) {
877 14 : if (signe(n) < 0) pari_err_INV("powp",x);
878 7 : return zeropadic(p, 0);
879 : }
880 147 : v = Z_pval(n, p);
881 147 : retmkpadic_i(Fp_pow(padic_u(x), n, _pd), icopy(p),
882 : mulpowu(padic_pd(x), p, v), 0, precp(x) + v);
883 : }
884 : static GEN
885 23767 : pow_polmod(GEN x, GEN n)
886 : {
887 23767 : GEN z = cgetg(3, t_POLMOD), a = gel(x,2), T = gel(x,1);
888 23767 : gel(z,1) = gcopy(T);
889 23767 : if (typ(a) != t_POL || varn(a) != varn(T) || lg(a) <= 3)
890 1143 : a = powgi(a, n);
891 : else {
892 22624 : pari_sp av = avma;
893 22624 : GEN p = NULL;
894 22624 : if (RgX_is_FpX(T, &p) && RgX_is_FpX(a, &p) && p)
895 : {
896 8771 : T = RgX_to_FpX(T, p); a = RgX_to_FpX(a, p);
897 8771 : if (lgefint(p) == 3)
898 : {
899 8764 : ulong pp = p[2];
900 8764 : a = Flxq_pow(ZX_to_Flx(a, pp), n, ZX_to_Flx(T, pp), pp);
901 8764 : a = Flx_to_ZX(a);
902 : }
903 : else
904 7 : a = FpXQ_pow(a, n, T, p);
905 8771 : a = FpX_to_mod(a, p);
906 8771 : a = gc_upto(av, a);
907 : }
908 : else
909 : {
910 13853 : set_avma(av);
911 13853 : a = RgXQ_pow(a, n, gel(z,1));
912 : }
913 : }
914 23767 : gel(z,2) = a; return z;
915 : }
916 :
917 : GEN
918 124487898 : gpowgs(GEN x, long n)
919 : {
920 : long m;
921 : pari_sp av;
922 : GEN y;
923 :
924 124487898 : if (n == 0) return gpowg0(x);
925 122621705 : if (n == 1)
926 : {
927 73973837 : long t = typ(x);
928 73973837 : if (is_scalar_t(t)) return gcopy(x);
929 717472 : switch(t)
930 : {
931 664424 : case t_POL: case t_SER: case t_RFRAC: case t_MAT: case t_VECSMALL:
932 664424 : return gcopy(x);
933 21 : case t_VEC: if (!is_ext_qfr(x)) break;
934 : /* fall through handle extended t_QFB */
935 53035 : case t_QFB: return qfbred(x);
936 : }
937 13 : pari_err_TYPE("gpow", x);
938 : }
939 48648012 : if (n ==-1) return ginv(x);
940 40171800 : switch(typ(x))
941 : {
942 28553308 : case t_INT: return powis(x,n);
943 5086377 : case t_REAL: return powrs(x,n);
944 29263 : case t_INTMOD:
945 29263 : y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
946 29263 : gel(y,2) = Fp_pows(gel(x,2), n, gel(x,1));
947 29263 : return y;
948 377840 : case t_FRAC:
949 : {
950 377840 : GEN a = gel(x,1), b = gel(x,2);
951 377840 : long s = (signe(a) < 0 && odd(n))? -1: 1;
952 377840 : if (n < 0) {
953 3045 : n = -n;
954 3045 : if (is_pm1(a)) return powiu_sign(b, n, s); /* +-1/x[2] inverts to t_INT */
955 2821 : swap(a, b);
956 : }
957 377616 : y = cgetg(3, t_FRAC);
958 377616 : gel(y,1) = powiu_sign(a, n, s);
959 377616 : gel(y,2) = powiu_sign(b, n, 1);
960 377615 : return y;
961 : }
962 1316945 : case t_PADIC: return powps(x, n);
963 248997 : case t_RFRAC:
964 : {
965 248997 : av = avma; y = cgetg(3, t_RFRAC); m = labs(n);
966 248997 : gel(y,1) = gpowgs(gel(x,1),m);
967 248997 : gel(y,2) = gpowgs(gel(x,2),m);
968 248997 : if (n < 0) y = ginv(y);
969 248997 : return gc_upto(av,y);
970 : }
971 23760 : case t_POLMOD: {
972 23760 : long N[] = {evaltyp(t_INT) | _evallg(3),0,0};
973 23760 : affsi(n,N); return pow_polmod(x, N);
974 : }
975 7 : case t_VEC: if (!is_ext_qfr(x)) pari_err_TYPE("gpow", x);
976 : /* fall through handle extended t_QFB */
977 1313376 : case t_QFB: return qfbpows(x, n);
978 1346686 : case t_POL:
979 1346686 : if (RgX_is_monomial(x)) return pow_monome(x, n);
980 : default: {
981 2580665 : pari_sp av = avma;
982 2580665 : y = gen_powu_i(x, (ulong)labs(n), NULL, &_sqr, &_mul);
983 2580687 : if (n < 0) y = ginv(y);
984 2580686 : return gc_upto(av,y);
985 : }
986 : }
987 : }
988 :
989 : /* n a t_INT */
990 : GEN
991 104605662 : powgi(GEN x, GEN n)
992 : {
993 : GEN y;
994 :
995 104605662 : if (!is_bigint(n)) return gpowgs(x, itos(n));
996 : /* probable overflow for nonmodular types (typical exception: (X^0)^N) */
997 25583 : switch(typ(x))
998 : {
999 25261 : case t_INTMOD:
1000 25261 : y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
1001 25265 : gel(y,2) = Fp_pow(gel(x,2), n, gel(x,1));
1002 25274 : return y;
1003 101 : case t_FFELT: return FF_pow(x,n);
1004 161 : case t_PADIC: return powp(x, n);
1005 :
1006 35 : case t_INT:
1007 35 : if (is_pm1(x)) return (signe(x) < 0 && mpodd(n))? gen_m1: gen_1;
1008 14 : if (signe(x)) pari_err_OVERFLOW("lg()");
1009 7 : if (signe(n) < 0) pari_err_INV("powgi",gen_0);
1010 7 : return gen_0;
1011 7 : case t_FRAC:
1012 7 : pari_err_OVERFLOW("lg()");
1013 :
1014 0 : case t_VEC: if (!is_ext_qfr(x)) pari_err_TYPE("gpow",x);
1015 : /* fall through handle extended t_QFB */
1016 18 : case t_QFB: return qfbpow(x, n);
1017 7 : case t_POLMOD: return pow_polmod(x, n);
1018 7 : default: {
1019 7 : pari_sp av = avma;
1020 7 : y = gen_pow_i(x, n, NULL, &_sqr, &_mul);
1021 7 : if (signe(n) < 0) return gc_upto(av, ginv(y));
1022 7 : return gc_GEN(av,y);
1023 : }
1024 : }
1025 : }
1026 :
1027 : /* Assume x = 1 + O(t), n a scalar. Return x^n */
1028 : static GEN
1029 8022 : ser_pow_1(GEN x, GEN n)
1030 : {
1031 : long lx, mi, i, j, d;
1032 8022 : GEN y = cgetg_copy(x, &lx), X = x+2, Y = y + 2;
1033 8022 : y[1] = evalsigne(1) | _evalvalser(0) | evalvarn(varn(x));
1034 74347 : d = mi = lx-3; while (mi>=1 && isrationalzero(gel(X,mi))) mi--;
1035 8022 : gel(Y,0) = gen_1;
1036 111741 : for (i=1; i<=d; i++)
1037 : {
1038 103719 : pari_sp av = avma;
1039 103719 : GEN s = gen_0;
1040 494319 : for (j=1; j<=minss(i,mi); j++)
1041 : {
1042 390600 : GEN t = gsubgs(gmulgu(n,j),i-j);
1043 390600 : s = gadd(s, gmul(gmul(t, gel(X,j)), gel(Y,i-j)));
1044 : }
1045 103719 : gel(Y,i) = gc_upto(av, gdivgu(s,i));
1046 : }
1047 8022 : return y;
1048 : }
1049 :
1050 : /* we suppose n != 0, valser(x) = 0 and leading-term(x) != 0. Not stack clean */
1051 : static GEN
1052 8127 : ser_pow(GEN x, GEN n, long prec)
1053 : {
1054 : GEN y, c, lead;
1055 8127 : if (varncmp(gvar(n), varn(x)) <= 0) return gexp(gmul(n, glog(x,prec)), prec);
1056 8022 : lead = gel(x,2);
1057 8022 : if (gequal1(lead)) return ser_pow_1(x, n);
1058 7567 : x = ser_normalize(x);
1059 7567 : if (typ(n) == t_FRAC && !isinexact(lead) && ispower(lead, gel(n,2), &c))
1060 210 : c = powgi(c, gel(n,1));
1061 : else
1062 7357 : c = gpow(lead,n, prec);
1063 7567 : y = gmul(c, ser_pow_1(x, n));
1064 : /* gpow(t_POLMOD,n) can be a t_COL [conjvec] */
1065 7567 : if (typ(y) != t_SER) pari_err_TYPE("gpow", y);
1066 7567 : return y;
1067 : }
1068 :
1069 : static long
1070 8036 : val_from_i(GEN E)
1071 : {
1072 8036 : if (is_bigint(E)) pari_err_OVERFLOW("sqrtn [valuation]");
1073 8029 : return itos(E);
1074 : }
1075 :
1076 : /* return x^q, assume typ(x) = t_SER, typ(q) = t_INT/t_FRAC and q != 0 */
1077 : static GEN
1078 8043 : ser_powfrac(GEN x, GEN q, long prec)
1079 : {
1080 8043 : GEN y, E = gmulsg(valser(x), q);
1081 : long e;
1082 :
1083 8043 : if (!signe(x))
1084 : {
1085 21 : if (gsigne(q) < 0) pari_err_INV("gpow", x);
1086 21 : return zeroser(varn(x), val_from_i(gfloor(E)));
1087 : }
1088 8022 : if (typ(E) != t_INT)
1089 7 : pari_err_DOMAIN("sqrtn", "valuation", "!=", mkintmod(gen_0, gel(q,2)), x);
1090 8015 : e = val_from_i(E);
1091 8015 : y = leafcopy(x); setvalser(y, 0);
1092 8015 : y = ser_pow(y, q, prec);
1093 8015 : setvalser(y, e); return y;
1094 : }
1095 :
1096 : static GEN
1097 126 : gpow0(GEN z, GEN x, long prec)
1098 : {
1099 126 : pari_sp av = avma;
1100 126 : switch(typ(x))
1101 : {
1102 84 : case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD:
1103 84 : break;
1104 35 : case t_VEC: case t_COL: case t_MAT:
1105 105 : pari_APPLY_same(gpow0(z,gel(x,i),prec));
1106 7 : default: pari_err_TYPE("gpow(0,x)", x);
1107 : }
1108 84 : x = real_i(x);
1109 84 : if (gsigne(x) <= 0) pari_err_DOMAIN("gpow(0,x)", "x", "<=", gen_0, x);
1110 77 : if (!precision(z)) return gcopy(z);
1111 :
1112 14 : z = ground(gmulsg(gexpo(z),x));
1113 14 : if (is_bigint(z) || uel(z,2) >= HIGHEXPOBIT)
1114 7 : pari_err_OVERFLOW("gpow");
1115 7 : set_avma(av); return real_0_bit(itos(z));
1116 : }
1117 :
1118 : /* centermod(x, log(2)), set *sh to the quotient */
1119 : static GEN
1120 19024044 : modlog2(GEN x, long *sh)
1121 : {
1122 19024044 : double d = rtodbl(x), qd = (fabs(d) + M_LN2/2)/M_LN2;
1123 : long q;
1124 19024036 : if (dblexpo(qd) >= BITS_IN_LONG-1) pari_err_OVERFLOW("expo()");
1125 19023937 : q = d < 0 ? - (long) qd: (long) qd;
1126 19023937 : *sh = q;
1127 19023937 : if (q) {
1128 15809630 : long l = realprec(x) + EXTRAPRECWORD;
1129 15809630 : x = subrr(rtor(x,l), mulsr(q, mplog2(l)));
1130 15809596 : if (!signe(x)) return NULL;
1131 : }
1132 19023903 : return x;
1133 : }
1134 :
1135 : /* x^n, n a t_FRAC */
1136 : static GEN
1137 10529413 : powfrac(GEN x, GEN n, long prec)
1138 : {
1139 10529413 : GEN a = gel(n,1), d = gel(n,2);
1140 10529413 : long D = itos_or_0(d);
1141 10528901 : if (D == 2)
1142 : {
1143 8930562 : GEN y = gsqrt(x,prec);
1144 8933582 : if (!equali1(a)) y = gmul(y, powgi(x, shifti(subiu(a,1), -1)));
1145 8930992 : return y;
1146 : }
1147 1598339 : if (D && is_real_t(typ(x)) && gsigne(x) > 0)
1148 : { /* x^n = x^q * x^(r/D) */
1149 1594178 : GEN z, r, q = truedvmdis(a, D, &r);
1150 1594185 : if (typ(x) == t_REAL)
1151 : {
1152 171840 : z = sqrtnr(x, D);
1153 171840 : if (!equali1(r)) z = powgi(z, r);
1154 171840 : if (signe(q)) z = gmul(z, powgi(x, q));
1155 : }
1156 : else
1157 : {
1158 1422345 : GEN X = x;
1159 1422345 : x = gtofp(x, prec + nbits2extraprec(expi(r)));
1160 1422345 : z = sqrtnr(x, D);
1161 1422345 : if (!equali1(r)) z = powgi(z, r);
1162 1422345 : if (signe(q))
1163 : {
1164 17068 : long e = typ(X)==t_INT? expi(X): maxuu(expi(gel(X,1)), expi(gel(X,2)));
1165 17068 : z = gmul(z, powgi(cmpiu(muliu(q,e), realprec(x)) > 0? x: X, q));
1166 : }
1167 : }
1168 1594185 : return z;
1169 : }
1170 4161 : return NULL;
1171 : }
1172 :
1173 : /* n = a+ib, x > 0 real, ex ~ |log2(x)|; return precision at which
1174 : * log(x) must be computed to evaluate x^n */
1175 : long
1176 192854 : powcx_prec(long ex, GEN n, long prec)
1177 : {
1178 192854 : GEN a = gel(n,1), b = gel(n,2);
1179 192854 : long e = (ex < 2)? 0: expu(ex);
1180 192854 : e += gexpo_safe(is_rational_t(typ(a))? b: n);
1181 192854 : return e > 2? prec + nbits2extraprec(e): prec;
1182 : }
1183 : GEN
1184 5520522 : powcx(GEN x, GEN logx, GEN n, long prec)
1185 : {
1186 5520522 : GEN sxb, cxb, xa, a = gel(n,1), xb = gmul(gel(n,2), logx);
1187 5522003 : long sh, p = realprec(logx);
1188 5522003 : switch(typ(a))
1189 : {
1190 49882 : case t_INT: xa = powgi(x, a); break;
1191 5378399 : case t_FRAC: xa = powfrac(x, a, prec);
1192 5377811 : if (xa) break;
1193 : default:
1194 93773 : xa = modlog2(gmul(gel(n,1), logx), &sh);
1195 93785 : if (!xa) xa = real2n(sh, prec);
1196 : else
1197 : {
1198 93785 : if (signe(xa) && realprec(xa) > prec) setprec(xa, prec);
1199 93785 : xa = mpexp(xa); shiftr_inplace(xa, sh);
1200 : }
1201 : }
1202 5521383 : if (typ(xb) != t_REAL) return xa;
1203 5521383 : if (gexpo(xb) > 30)
1204 : {
1205 5184155 : GEN q, P = Pi2n(-2, p), z = addrr(xb,P); /* = x + Pi/4 */
1206 5182916 : shiftr_inplace(P, 1);
1207 5181371 : q = floorr(divrr(z, P)); /* round ( x / (Pi/2) ) */
1208 5165103 : xb = subrr(xb, mulir(q, P)); /* x mod Pi/2 */
1209 5182191 : sh = Mod4(q);
1210 : }
1211 : else
1212 : {
1213 336913 : long q = floor(rtodbl(xb) / (M_PI/2) + 0.5);
1214 336927 : if (q) xb = subrr(xb, mulsr(q, Pi2n(-1,p))); /* x mod Pi/2 */
1215 336922 : sh = q & 3;
1216 : }
1217 5518040 : if (signe(xb) && realprec(xb) > prec) setprec(xb, prec);
1218 5518040 : mpsincos(xb, &sxb, &cxb);
1219 5524801 : return gmul(xa, mulcxpowIs(mkcomplex(cxb, sxb), sh));
1220 : }
1221 :
1222 : GEN
1223 21699587 : gpow(GEN x, GEN n, long prec)
1224 : {
1225 21699587 : long i, prec0, tx, tn = typ(n);
1226 : pari_sp av;
1227 : GEN y;
1228 :
1229 21699587 : if (tn == t_INT) return powgi(x,n);
1230 6162195 : tx = typ(x);
1231 6162279 : if (is_matvec_t(tx)) pari_APPLY_same(gpow(gel(x,i),n,prec));
1232 6162215 : av = avma;
1233 6162215 : switch (tx)
1234 : {
1235 28 : case t_POL: case t_RFRAC: x = toser_i(x); /* fall through */
1236 7560 : case t_SER:
1237 7560 : if (tn == t_FRAC) return gc_upto(av, ser_powfrac(x, n, prec));
1238 140 : if (valser(x))
1239 21 : pari_err_DOMAIN("gpow [irrational exponent]",
1240 : "valuation", "!=", gen_0, x);
1241 119 : if (lg(x) == 2) return gc_GEN(av, x); /* O(1) */
1242 112 : return gc_upto(av, ser_pow(x, n, prec));
1243 : }
1244 6154649 : if (gequal0(x)) return gpow0(x, n, prec);
1245 6154595 : if (tn == t_FRAC)
1246 : {
1247 5154672 : GEN p, z, a = gel(n,1), d = gel(n,2);
1248 5154672 : switch (tx)
1249 : {
1250 1481235 : case t_INT:
1251 1481235 : if (signe(x) < 0)
1252 : {
1253 42 : if (equaliu(d, 2) && Z_issquareall(negi(x), &z))
1254 : {
1255 21 : z = powgi(z, a);
1256 21 : if (Mod4(a) == 3) z = gneg(z);
1257 5150423 : return gc_GEN(av, mkcomplex(gen_0, z));
1258 : }
1259 21 : break;
1260 : }
1261 1481193 : if (ispower(x, d, &z)) return powgi(z, a);
1262 1479233 : break;
1263 70030 : case t_FRAC:
1264 70030 : if (signe(gel(x,1)) < 0)
1265 : {
1266 28 : if (equaliu(d, 2) && ispower(absfrac(x), d, &z))
1267 7 : return gc_GEN(av, mkcomplex(gen_0, powgi(z, a)));
1268 21 : break;
1269 : }
1270 70002 : if (ispower(x, d, &z)) return powgi(z, a);
1271 68630 : break;
1272 :
1273 21 : case t_INTMOD:
1274 21 : p = gel(x,1);
1275 21 : if (!BPSW_psp(p)) pari_err_PRIME("gpow",p);
1276 14 : y = cgetg(3,t_INTMOD); gel(y,1) = icopy(p);
1277 14 : av = avma;
1278 14 : z = Fp_sqrtn(gel(x,2), d, p, NULL);
1279 14 : if (!z) pari_err_SQRTN("gpow",x);
1280 7 : gel(y,2) = gc_INT(av, Fp_pow(z, a, p));
1281 7 : return y;
1282 :
1283 14 : case t_PADIC:
1284 14 : z = Qp_sqrtn(x, d, NULL); if (!z) pari_err_SQRTN("gpow",x);
1285 7 : return gc_upto(av, powgi(z, a));
1286 :
1287 21 : case t_FFELT:
1288 21 : return gc_upto(av,FF_pow(FF_sqrtn(x,d,NULL),a));
1289 : }
1290 5151256 : z = powfrac(x, n, prec);
1291 5151392 : if (z) return gc_upto(av, z);
1292 : }
1293 1004031 : if (tn == t_COMPLEX && is_real_t(typ(x)) && gsigne(x) > 0)
1294 : {
1295 180997 : long p = powcx_prec(fabs(dbllog2(x)), n, prec);
1296 180996 : return gc_upto(av, powcx(x, glog(x, p), n, prec));
1297 : }
1298 823033 : if (tn == t_PADIC) x = gcvtop(x, padic_p(n), precp(n));
1299 823033 : i = precision(n);
1300 823034 : if (i) prec = i;
1301 823034 : prec0 = prec;
1302 823034 : if (!gprecision(x))
1303 : {
1304 91615 : long e = gexpo_safe(n); /* avoided if n = 0 or gexpo not defined */
1305 91615 : if (e > 2) prec += nbits2extraprec(e);
1306 : }
1307 823034 : y = gmul(n, glog(x,prec));
1308 823006 : y = gexp(y,prec);
1309 823006 : if (prec0 == prec) return gc_upto(av, y);
1310 79730 : return gc_GEN(av, gprec_wtrunc(y,prec0));
1311 : }
1312 : GEN
1313 11683 : powPis(GEN s, long prec)
1314 : {
1315 11683 : pari_sp av = avma;
1316 : GEN x;
1317 11683 : if (typ(s) != t_COMPLEX) return gpow(mppi(prec), s, prec);
1318 490 : x = mppi(powcx_prec(1, s, prec));
1319 490 : return gc_upto(av, powcx(x, logr_abs(x), s, prec));
1320 : }
1321 : GEN
1322 12208 : pow2Pis(GEN s, long prec)
1323 : {
1324 12208 : pari_sp av = avma;
1325 : GEN x;
1326 12208 : if (typ(s) != t_COMPLEX) return gpow(Pi2n(1,prec), s, prec);
1327 1876 : x = Pi2n(1, powcx_prec(2, s, prec));
1328 1876 : return gc_upto(av, powcx(x, logr_abs(x), s, prec));
1329 : }
1330 :
1331 : GEN
1332 208228 : gpowers0(GEN x, long n, GEN x0)
1333 : {
1334 : long i, l;
1335 : GEN V;
1336 208228 : if (!x0) return gpowers(x,n);
1337 193739 : if (n < 0) return cgetg(1,t_VEC);
1338 193739 : l = n+2; V = cgetg(l, t_VEC); gel(V,1) = gcopy(x0);
1339 7605256 : for (i = 2; i < l; i++) gel(V,i) = gmul(gel(V,i-1),x);
1340 193774 : return V;
1341 : }
1342 :
1343 : GEN
1344 782763 : gpowers(GEN x, long n)
1345 : {
1346 782763 : if (n < 0) return cgetg(1,t_VEC);
1347 782756 : return gen_powers(x, n, 0, (void*)x, &_sqr, &_mul, &_one);
1348 : }
1349 :
1350 : /* return [q^1,q^4,...,q^{n^2}] */
1351 : GEN
1352 39711 : gsqrpowers(GEN q, long n)
1353 : {
1354 39711 : pari_sp av = avma;
1355 39711 : GEN L = gpowers0(gsqr(q), n, q); /* L[i] = q^(2i - 1), i <= n+1 */
1356 39711 : GEN v = cgetg(n+1, t_VEC);
1357 : long i;
1358 39711 : gel(v, 1) = gcopy(q);
1359 6737729 : for (i = 2; i <= n ; ++i) gel(v, i) = q = gmul(q, gel(L,i)); /* q^(i^2) */
1360 39711 : return gc_upto(av, v);
1361 : }
1362 :
1363 : /* 4 | N. returns a vector RU which contains exp(2*i*k*Pi/N), k=0..N-1 */
1364 : static GEN
1365 1005544 : grootsof1_4(long N, long prec)
1366 : {
1367 1005544 : GEN z, RU = cgetg(N+1,t_COL), *v = ((GEN*)RU) + 1;
1368 1005544 : long i, N2 = (N>>1), N4 = (N>>2), N8 = (N>>3);
1369 : /* z^N2 = -1, z^N4 = I; if z^k = a+I*b, then z^(N4-k) = I*conj(z) = b+a*I */
1370 :
1371 1005544 : v[0] = gen_1; v[1] = z = rootsof1u_cx(N, prec);
1372 1005543 : if (odd(N4)) N8++;
1373 1114828 : for (i=1; i<N8; i++)
1374 : {
1375 109285 : GEN t = v[i];
1376 109285 : v[i+1] = gmul(z, t);
1377 109285 : v[N4-i] = mkcomplex(gel(t,2), gel(t,1));
1378 : }
1379 2541131 : for (i=0; i<N4; i++) v[i+N4] = mulcxI(v[i]);
1380 4076712 : for (i=0; i<N2; i++) v[i+N2] = gneg(v[i]);
1381 1005540 : return RU;
1382 : }
1383 :
1384 : /* as above, N arbitrary */
1385 : GEN
1386 1171466 : grootsof1(long N, long prec)
1387 : {
1388 : GEN z, RU, *v;
1389 : long i, k;
1390 :
1391 1171466 : if (N <= 0) pari_err_DOMAIN("rootsof1", "N", "<=", gen_0, stoi(N));
1392 1171453 : if ((N & 3) == 0) return grootsof1_4(N, prec);
1393 165910 : if (N <= 2) return N == 1? mkcol(gen_1): mkcol2(gen_1, gen_m1);
1394 45613 : k = (N+1)>>1;
1395 45613 : RU = cgetg(N+1,t_COL);
1396 45613 : v = ((GEN*)RU) + 1;
1397 45613 : v[0] = gen_1; v[1] = z = rootsof1u_cx(N, prec);
1398 108658 : for (i=2; i<k; i++) v[i] = gmul(z, v[i-1]);
1399 45613 : if (!odd(N)) v[i++] = gen_m1; /*avoid loss of accuracy*/
1400 154271 : for ( ; i<N; i++) v[i] = gconj(v[N-i]);
1401 45613 : return RU;
1402 : }
1403 :
1404 : /********************************************************************/
1405 : /** **/
1406 : /** RACINE CARREE **/
1407 : /** **/
1408 : /********************************************************************/
1409 : /* assume x unit, e = precp(x) */
1410 : GEN
1411 144690 : Z2_sqrt(GEN x, long e)
1412 : {
1413 144690 : ulong r = signe(x)>=0?mod16(x):16-mod16(x);
1414 : GEN z;
1415 : long ez;
1416 : pari_sp av;
1417 :
1418 144690 : switch(e)
1419 : {
1420 7 : case 1: return gen_1;
1421 161 : case 2: return (r & 3UL) == 1? gen_1: NULL;
1422 28 : case 3: return (r & 7UL) == 1? gen_1: NULL;
1423 71064 : case 4: if (r == 1) return gen_1;
1424 35133 : else return (r == 9)? utoipos(3): NULL;
1425 73430 : default: if ((r&7UL) != 1) return NULL;
1426 : }
1427 73430 : av = avma; z = (r==1)? gen_1: utoipos(3);
1428 73430 : ez = 3; /* number of correct bits in z (compared to sqrt(x)) */
1429 : for(;;)
1430 47978 : {
1431 : GEN mod;
1432 121408 : ez = (ez<<1) - 1;
1433 121408 : if (ez > e) ez = e;
1434 121408 : mod = int2n(ez);
1435 121408 : z = addii(z, remi2n(mulii(x, Fp_inv(z,mod)), ez));
1436 121408 : z = shifti(z, -1); /* (z + x/z) / 2 */
1437 121408 : if (e == ez) return gc_INT(av, z);
1438 47978 : if (ez < e) ez--;
1439 47978 : if (gc_needed(av,2))
1440 : {
1441 0 : if (DEBUGMEM > 1) pari_warn(warnmem,"Qp_sqrt");
1442 0 : z = gc_INT(av,z);
1443 : }
1444 : }
1445 : }
1446 :
1447 : static GEN
1448 1834 : Up_sqrt(GEN u, GEN p, GEN pd, long d, long D)
1449 1834 : { retmkpadic(Zp_sqrt(u, p, d), icopy(p), icopy(pd), 0, D); }
1450 :
1451 : /* x unit defined modulo p^e, e > 0 */
1452 : GEN
1453 1862 : Qp_sqrt(GEN x)
1454 : {
1455 1862 : long D, d, e = valp(x);
1456 1862 : GEN y, mod, p = padic_p(x);
1457 :
1458 1862 : if (gequal0(x)) return zeropadic(p, (e+1) >> 1);
1459 1848 : if (e & 1) return NULL;
1460 :
1461 1834 : mod = padic_pd(x); D = d = precp(x); e >>= 1;
1462 1834 : if (absequaliu(p,2)) { D = (d <= 3) ? 1 : d-1; mod = int2n(D); }
1463 1834 : y = Up_sqrt(padic_u(x), p, mod, d, D);
1464 1834 : if (!padic_u(y)) return NULL;
1465 1778 : setvalp(y, e); return y;
1466 : }
1467 :
1468 : GEN
1469 420 : Zn_sqrt(GEN d, GEN fn)
1470 : {
1471 420 : pari_sp ltop = avma, btop;
1472 420 : GEN b = gen_0, m = gen_1;
1473 : long j, np;
1474 420 : if (typ(d) != t_INT) pari_err_TYPE("Zn_sqrt",d);
1475 420 : if (typ(fn) == t_INT)
1476 0 : fn = absZ_factor(fn);
1477 420 : else if (!is_Z_factorpos(fn))
1478 0 : pari_err_TYPE("Zn_sqrt",fn);
1479 420 : np = nbrows(fn);
1480 420 : btop = avma;
1481 1680 : for (j = 1; j <= np; ++j)
1482 : {
1483 : GEN bp, mp, pr, r;
1484 1260 : GEN p = gcoeff(fn, j, 1);
1485 1260 : long e = itos(gcoeff(fn, j, 2));
1486 1260 : long v = Z_pvalrem(d,p,&r);
1487 1260 : if (v >= e) bp =gen_0;
1488 : else
1489 : {
1490 1134 : if (odd(v)) return NULL;
1491 1134 : bp = Zp_sqrt(r, p, e-v);
1492 1134 : if (!bp) return NULL;
1493 1134 : if (v) bp = mulii(bp, powiu(p, v>>1L));
1494 : }
1495 1260 : mp = powiu(p, e);
1496 1260 : pr = mulii(m, mp);
1497 1260 : b = Z_chinese_coprime(b, bp, m, mp, pr);
1498 1260 : m = pr;
1499 1260 : if (gc_needed(btop, 1))
1500 0 : (void)gc_all(btop, 2, &b, &m);
1501 : }
1502 420 : return gc_upto(ltop, b);
1503 : }
1504 :
1505 : static GEN
1506 18739 : sqrt_ser(GEN b, long prec)
1507 : {
1508 18739 : long e = valser(b), vx = varn(b), lx, lold, j;
1509 : ulong mask;
1510 : GEN a, x, lta, ltx;
1511 :
1512 18739 : if (!signe(b)) return zeroser(vx, e>>1);
1513 18739 : a = leafcopy(b);
1514 18739 : x = cgetg_copy(b, &lx);
1515 18739 : if (e & 1)
1516 14 : pari_err_DOMAIN("sqrtn", "valuation", "!=", mkintmod(gen_0, gen_2), b);
1517 18725 : a[1] = x[1] = evalsigne(1) | evalvarn(0) | _evalvalser(0);
1518 18725 : lta = gel(a,2);
1519 18725 : if (gequal1(lta)) ltx = lta;
1520 14833 : else if (!issquareall(lta,<x)) ltx = gsqrt(lta,prec);
1521 18718 : gel(x,2) = ltx;
1522 316771 : for (j = 3; j < lx; j++) gel(x,j) = gen_0;
1523 18718 : setlg(x,3);
1524 18718 : mask = quadratic_prec_mask(lx - 2);
1525 18718 : lold = 1;
1526 96715 : while (mask > 1)
1527 : {
1528 77997 : GEN y, x2 = gmul2n(x,1);
1529 77997 : long l = lold << 1, lx;
1530 :
1531 77997 : if (mask & 1) l--;
1532 77997 : mask >>= 1;
1533 77997 : setlg(a, l + 2);
1534 77997 : setlg(x, l + 2);
1535 77997 : y = sqr_ser_part(x, lold, l-1) - lold;
1536 376050 : for (j = lold+2; j < l+2; j++) gel(y,j) = gsub(gel(y,j), gel(a,j));
1537 77997 : y += lold; setvalser(y, lold);
1538 77997 : y = normalizeser(y);
1539 77997 : y = gsub(x, gdiv(y, x2)); /* = gmul2n(gsub(x, gdiv(a,x)), -1); */
1540 77997 : lx = minss(l+2, lg(y));
1541 376043 : for (j = lold+2; j < lx; j++) gel(x,j) = gel(y,j);
1542 77997 : lold = l;
1543 : }
1544 18718 : x[1] = evalsigne(1) | evalvarn(vx) | _evalvalser(e >> 1);
1545 18718 : return x;
1546 : }
1547 :
1548 : GEN
1549 62734369 : gsqrt(GEN x, long prec)
1550 : {
1551 : pari_sp av;
1552 : GEN y;
1553 :
1554 62734369 : switch(typ(x))
1555 : {
1556 5535134 : case t_INT:
1557 5535134 : if (!signe(x)) return real_0(prec); /* no loss of accuracy */
1558 5535064 : x = itor(x,prec); /* fall through */
1559 55986031 : case t_REAL: return sqrtr(x);
1560 :
1561 35 : case t_INTMOD:
1562 : {
1563 35 : GEN p = gel(x,1), a;
1564 35 : y = cgetg(3,t_INTMOD); gel(y,1) = icopy(p);
1565 35 : a = Fp_sqrt(gel(x,2),p);
1566 21 : if (!a)
1567 : {
1568 7 : if (!BPSW_psp(p)) pari_err_PRIME("sqrt [modulus]",p);
1569 7 : pari_err_SQRTN("gsqrt",x);
1570 : }
1571 14 : gel(y,2) = a; return y;
1572 : }
1573 :
1574 6464387 : case t_COMPLEX:
1575 : { /* (u+iv)^2 = a+ib <=> u^2+v^2 = sqrt(a^2+b^2), u^2-v^2=a, 2uv=b */
1576 6464387 : GEN a = gel(x,1), b = gel(x,2), r, u, v;
1577 6464387 : if (isrationalzero(b)) return gsqrt(a, prec);
1578 6464387 : y = cgetg(3,t_COMPLEX); av = avma;
1579 :
1580 6464388 : r = cxnorm(x);
1581 6464384 : if (typ(r) == t_INTMOD || typ(r) == t_PADIC)
1582 0 : pari_err_IMPL("sqrt(complex of t_INTMODs)");
1583 6464384 : r = gsqrt(r, prec); /* t_REAL, |a+Ib| */
1584 6464389 : if (!signe(r))
1585 67 : u = v = gc_uptoleaf(av, sqrtr(r));
1586 6464322 : else if (gsigne(a) < 0)
1587 : {
1588 : /* v > 0 since r > 0, a < 0, rounding errors can't make the sum of two
1589 : * positive numbers = 0 */
1590 192020 : v = sqrtr( gmul2n(gsub(r,a), -1) );
1591 192020 : if (gsigne(b) < 0) togglesign(v);
1592 192020 : v = gc_uptoleaf(av, v); av = avma;
1593 : /* v = 0 is impossible */
1594 192020 : u = gc_uptoleaf(av, gdiv(b, shiftr(v,1)));
1595 : } else {
1596 6272302 : u = sqrtr( gmul2n(gadd(r,a), -1) );
1597 6272302 : u = gc_uptoleaf(av, u); av = avma;
1598 6272302 : if (!signe(u)) /* possible if a = 0.0, e.g. sqrt(0.e-10+1e-10*I) */
1599 7 : v = u;
1600 : else
1601 6272295 : v = gc_uptoleaf(av, gdiv(b, shiftr(u,1)));
1602 : }
1603 6464384 : gel(y,1) = u;
1604 6464384 : gel(y,2) = v; return y;
1605 : }
1606 :
1607 63 : case t_PADIC:
1608 63 : y = Qp_sqrt(x);
1609 63 : if (!y) pari_err_SQRTN("Qp_sqrt",x);
1610 42 : return y;
1611 :
1612 7161 : case t_FFELT: return FF_sqrt(x);
1613 :
1614 274498 : default:
1615 274498 : av = avma; if (!(y = toser_i(x))) break;
1616 18739 : return gc_GEN(av, sqrt_ser(y, prec));
1617 : }
1618 255759 : return trans_eval("sqrt",gsqrt,x,prec);
1619 : }
1620 : /********************************************************************/
1621 : /** **/
1622 : /** N-th ROOT **/
1623 : /** **/
1624 : /********************************************************************/
1625 :
1626 : static GEN
1627 304145 : Z_to_padic(GEN a, GEN p, long e)
1628 : {
1629 304145 : if (!signe(a))
1630 1316 : return zeropadic(p, e);
1631 : else
1632 : {
1633 302829 : long v = Z_pvalrem(a, p, &a), d = e - v;
1634 302827 : retmkpadic(icopy(a), icopy(p), powiu(p, d), v, d);
1635 : }
1636 : }
1637 :
1638 : GEN
1639 196210 : Qp_log(GEN x)
1640 : {
1641 196210 : pari_sp av = avma;
1642 196210 : GEN y, p = padic_p(x), a = padic_u(x);
1643 196210 : long e = precp(x);
1644 :
1645 196210 : if (!signe(a)) pari_err_DOMAIN("Qp_log", "argument", "=", gen_0, x);
1646 196189 : if (absequaliu(p,2) || equali1(modii(a, p)))
1647 75508 : y = Zp_log(a, p, e);
1648 : else
1649 : { /* compute log(x^(p-1)) / (p-1) */
1650 120682 : GEN q = padic_pd(x), t = subiu(p, 1);
1651 120682 : a = Fp_pow(a, t, q);
1652 120682 : y = Fp_mul(Zp_log(a, p, e), diviiexact(subsi(1, q), t), q);
1653 : }
1654 196190 : return gc_upto(av, Z_to_padic(y, p, e));
1655 : }
1656 :
1657 : static GEN Qp_exp_safe(GEN x);
1658 :
1659 : /*compute the p^e th root of x p-adic, ignoring valuation; assume x != 0 */
1660 : static GEN
1661 784 : Up_sqrtn_ram(GEN x, long e)
1662 : {
1663 784 : GEN a, p = padic_p(x), n = powiu(p,e);
1664 : long va;
1665 : /*If p = 2, -1 is a root of 1 in U1: need extra check*/
1666 784 : if (absequaliu(p, 2) && mod8(padic_u(x)) != 1) return NULL;
1667 749 : a = Qp_log(x); va = valp(a) - e;
1668 749 : if (va <= 0)
1669 : {
1670 287 : if (signe(padic_u(a))) return NULL;
1671 : /* all accuracy lost */
1672 119 : a = cvtop(remii(padic_u(x),p), p, 1);
1673 : }
1674 : else
1675 : {
1676 462 : setvalp(a, va); /* divide by p^e */
1677 462 : a = Qp_exp_safe(a);
1678 462 : if (!a) return NULL;
1679 : /* n=p^e and a^n=z*x where z is a (p-1)th-root of 1.
1680 : * Since z^n=z, we have (a/z)^n = x. */
1681 462 : a = gdiv(x, powgi(a,subiu(n,1))); /* = a/z = x/a^(n-1)*/
1682 : }
1683 581 : return a;
1684 : }
1685 :
1686 : /* set q = s/n if n | s */
1687 : static int
1688 2296 : dvdsi_quot(long s, GEN n, long *q)
1689 : {
1690 : long r;
1691 2296 : *q = sdivsi_rem(s, n, &r);
1692 2296 : return !r;
1693 : }
1694 :
1695 : GEN
1696 2534 : Qp_sqrtn(GEN x, GEN n, GEN *zetan)
1697 : {
1698 2534 : pari_sp av = avma;
1699 2534 : GEN a, u, q, p = padic_p(x);
1700 : long e, v, prec;
1701 2534 : if (absequaliu(n, 2))
1702 : {
1703 35 : if (zetan) *zetan = gen_m1;
1704 35 : if (signe(n) < 0) x = ginv(x);
1705 28 : return Qp_sqrt(x);
1706 : }
1707 2499 : if (!signe(padic_u(x)))
1708 : {
1709 203 : if (signe(n) < 0) pari_err_INV("Qp_sqrtn", x);
1710 203 : q = divii(addis(n, valp(x)-1), n);
1711 203 : if (zetan) *zetan = gen_1;
1712 203 : set_avma(av); return zeropadic(p, itos(q));
1713 : }
1714 2296 : if (!dvdsi_quot(valp(x), n, &v)) return NULL;
1715 : /* treat the ramified part using logarithms */
1716 2219 : e = Z_pvalrem(n, p, &q);
1717 2219 : if (e) { x = Up_sqrtn_ram(x,e); if (!x) return NULL; }
1718 2016 : if (is_pm1(q))
1719 : { /* finished */
1720 21 : if (signe(q) < 0) x = ginv(x);
1721 21 : if (v) setvalp(x, v);
1722 21 : x = gc_upto(av, x);
1723 49 : if (zetan) *zetan = (e && absequaliu(p, 2))? gen_m1 /*-1 in Q_2*/
1724 28 : : gen_1;
1725 21 : return x;
1726 : }
1727 : /* Hensel lift for unramified case */
1728 1995 : u = padic_u(x); prec = precp(x);
1729 1995 : a = Fp_sqrtn(u, q, p, zetan); if (!a) return NULL;
1730 1995 : u = Zp_sqrtnlift(u, q, a, p, prec);
1731 1995 : x = mkpadic(u, p, padic_pd(x), v, prec);
1732 1995 : if (!zetan) return gc_GEN(av, x);
1733 :
1734 14 : u = Zp_sqrtnlift(gen_1, q, *zetan, p, prec);
1735 14 : *zetan = mkpadic(u, p, padic_pd(x), 0, prec);
1736 14 : if (e && absequaliu(p, 2)) *zetan = gneg(*zetan); /*-1 in Q_2*/
1737 14 : return gc_all(av, 2, &x, zetan);
1738 : }
1739 :
1740 : GEN
1741 26070 : sqrtnint(GEN a, long n)
1742 : {
1743 26070 : pari_sp av = avma;
1744 : GEN x, b, q;
1745 : long s, k, e;
1746 26070 : const ulong nm1 = n - 1;
1747 26070 : if (n == 2) return sqrtint(a);
1748 22171 : if (typ(a) != t_INT)
1749 : {
1750 35 : if (typ(a) == t_REAL)
1751 : {
1752 : long e;
1753 14 : switch(signe(a))
1754 : {
1755 0 : case 0: return gen_0;
1756 7 : case -1: pari_err_DOMAIN("sqrtnint", "argument", "<", gen_0,a);
1757 : }
1758 7 : e = expo(a); if (e < 0) return gen_0;
1759 7 : if (nbits2lg(e+1) > lg(a))
1760 0 : a = floorr(sqrtnr(a,n)); /* try to avoid precision loss in truncation */
1761 : else
1762 7 : a = sqrtnint(truncr(a),n);
1763 : }
1764 : else
1765 : {
1766 21 : GEN b = gfloor(a);
1767 21 : if (typ(b) != t_INT) pari_err_TYPE("sqrtint",a);
1768 14 : if (signe(b) < 0) pari_err_DOMAIN("sqrtnint", "argument", "<", gen_0,b);
1769 7 : a = sqrtnint(b, n);
1770 : }
1771 14 : return gc_INT(av, a);
1772 : }
1773 22136 : if (n <= 0) pari_err_DOMAIN("sqrtnint", "n", "<=", gen_0, stoi(n));
1774 22129 : if (n == 1) return icopy(a);
1775 19973 : s = signe(a);
1776 19973 : if (s < 0) pari_err_DOMAIN("sqrtnint", "x", "<", gen_0, a);
1777 19973 : if (!s) return gen_0;
1778 19896 : if (lgefint(a) == 3) return utoi(usqrtn(itou(a), n));
1779 13665 : e = expi(a); k = e/(2*n);
1780 13665 : if (k == 0)
1781 : {
1782 : long flag;
1783 291 : if (n > e) return gc_const(av, gen_1);
1784 291 : flag = cmpii(a, powuu(3, n)); set_avma(av);
1785 291 : return (flag < 0) ? gen_2: stoi(3);
1786 : }
1787 13374 : if (e < n*BITS_IN_LONG - 1)
1788 : {
1789 : ulong xs, qs;
1790 6785 : b = itor(a, (2*e < n*BITS_IN_LONG)? DEFAULTPREC: MEDDEFAULTPREC);
1791 6785 : x = mpexp(divru(logr_abs(b), n));
1792 6785 : xs = itou(floorr(x)) + 1; /* >= a^(1/n) */
1793 : for(;;) {
1794 13392 : q = divii(a, powuu(xs, nm1));
1795 13392 : if (lgefint(q) > 3) break;
1796 13385 : qs = itou(q); if (qs >= xs) break;
1797 6607 : xs -= (xs - qs + nm1)/n;
1798 : }
1799 6785 : return utoi(xs);
1800 : }
1801 6589 : b = addui(1, shifti(a, -n*k));
1802 6589 : x = shifti(addui(1, sqrtnint(b, n)), k);
1803 6589 : q = divii(a, powiu(x, nm1));
1804 14775 : while (cmpii(q, x) < 0) /* a priori one iteration, no GC necessary */
1805 : {
1806 8186 : x = subii(x, divis(addui(nm1, subii(x, q)), n));
1807 8186 : q = divii(a, powiu(x, nm1));
1808 : }
1809 6589 : return gc_uptoleaf(av, x);
1810 : }
1811 :
1812 : ulong
1813 8114 : usqrtn(ulong a, ulong n)
1814 : {
1815 : ulong x, s, q;
1816 8114 : const ulong nm1 = n - 1;
1817 8114 : if (!n) pari_err_DOMAIN("sqrtnint", "n", "=", gen_0, utoi(n));
1818 8114 : if (n == 1 || a == 0) return a;
1819 8114 : s = 1 + expu(a)/n; x = 1UL << s;
1820 8114 : q = (nm1*s >= BITS_IN_LONG)? 0: a >> (nm1*s);
1821 21073 : while (q < x) {
1822 : ulong X;
1823 12959 : x -= (x - q + nm1)/n;
1824 12959 : X = upowuu(x, nm1);
1825 12959 : q = X? a/X: 0;
1826 : }
1827 8114 : return x;
1828 : }
1829 :
1830 : static ulong
1831 1735977 : cubic_prec_mask(long n)
1832 : {
1833 1735977 : long a = n, i;
1834 1735977 : ulong mask = 0;
1835 1735977 : for(i = 1;; i++, mask *= 3)
1836 8254173 : {
1837 9990150 : long c = a%3;
1838 9990150 : if (c) mask += 3 - c;
1839 9990150 : a = (a+2)/3;
1840 9990150 : if (a==1) return mask + upowuu(3, i);
1841 : }
1842 : }
1843 :
1844 : /* cubic Newton iteration, |a|^(1/n), assuming a != 0 */
1845 : GEN
1846 2799246 : sqrtnr_abs(GEN a, long n)
1847 : {
1848 : pari_sp av;
1849 : GEN x, b;
1850 : long eextra, eold, n1, n2, prec, B, v;
1851 : ulong mask;
1852 2799246 : double K = n, X;
1853 :
1854 2799246 : if (n == 1) return mpabs(a);
1855 2798551 : if (n == 2) return sqrtr_abs(a);
1856 :
1857 2451671 : prec = realprec(a); v = expo(a) / n; av = avma;
1858 2451671 : if (v) a = shiftr(a, -n*v);
1859 2451676 : b = rtor(a, DEFAULTPREC);
1860 2451698 : x = mpexp(divru(logr_abs(b), n));
1861 2451706 : if (prec == DEFAULTPREC)
1862 : {
1863 754144 : if (v) shiftr_inplace(x, v);
1864 754145 : return gc_uptoleaf(av, x);
1865 : }
1866 1697562 : X = rtodbl(x);
1867 1697562 : K = (K*K-1) / (12*X*X); /* |x_{n+1} - x| < K |x_n - x|^3 */
1868 1697562 : eextra = dblexpo(K);
1869 1697562 : n1 = n+1;
1870 1697562 : n2 = 2*n;
1871 1697562 : B = prec2nbits(prec);
1872 1697562 : mask = cubic_prec_mask(B + 63);
1873 1697562 : eold = 1;
1874 : for(;;)
1875 6768337 : { /* reach 64 */
1876 8465899 : long enew = eold * 3;
1877 8465899 : enew -= mask % 3;
1878 8465899 : if (enew > 64) break; /* back up one step */
1879 6768337 : mask /= 3;
1880 6768337 : eold = enew;
1881 : }
1882 : for(;;)
1883 1318731 : {
1884 3016293 : long pr, enew = eold * 3;
1885 : GEN y, z;
1886 3016293 : enew -= mask % 3;
1887 3016293 : mask /= 3;
1888 3016293 : pr = nbits2prec(enew + eextra);
1889 3016293 : b = rtor(a, pr); setsigne(b,1);
1890 3016293 : x = rtor(x, pr);
1891 3016293 : y = subrr(powru(x, n), b);
1892 3016293 : z = divrr(y, addrr(mulur(n1, y), mulur(n2, b)));
1893 3016293 : shiftr_inplace(z,1);
1894 3016293 : x = subrr(x, mulrr(x,z));
1895 3016293 : if (mask == 1)
1896 : {
1897 1697562 : if (v) shiftr_inplace(x, v);
1898 1697562 : return gc_uptoleaf(av, gprec_wtrunc(x,prec));
1899 : }
1900 1318731 : eold = enew;
1901 : }
1902 : }
1903 :
1904 : static void
1905 55428 : shiftc_inplace(GEN z, long d)
1906 : {
1907 55428 : shiftr_inplace(gel(z,1), d);
1908 55428 : shiftr_inplace(gel(z,2), d);
1909 55428 : }
1910 :
1911 : /* exp(2*Pi*I/n), same iteration as sqrtnr_abs, different initial point */
1912 : static GEN
1913 554096 : sqrtnof1(ulong n, long prec)
1914 : {
1915 : pari_sp av;
1916 : GEN x;
1917 : long eold, n1, n2, B;
1918 : ulong mask;
1919 :
1920 554096 : B = prec2nbits(prec);
1921 554096 : n1 = n+1;
1922 554096 : n2 = 2*n; av = avma;
1923 :
1924 554096 : x = expIr(divru(Pi2n(1, LOWDEFAULTPREC), n));
1925 554096 : if (prec == LOWDEFAULTPREC) return gc_upto(av, x);
1926 38415 : mask = cubic_prec_mask(B + BITS_IN_LONG-1);
1927 38415 : eold = 1;
1928 : for(;;)
1929 150092 : { /* reach BITS_IN_LONG */
1930 188507 : long enew = eold * 3;
1931 188507 : enew -= mask % 3;
1932 188507 : if (enew > BITS_IN_LONG) break; /* back up one step */
1933 150092 : mask /= 3;
1934 150092 : eold = enew;
1935 : }
1936 : for(;;)
1937 17013 : {
1938 55428 : long pr, enew = eold * 3;
1939 : GEN y, z;
1940 55428 : enew -= mask % 3;
1941 55428 : mask /= 3;
1942 55428 : pr = nbits2prec(enew);
1943 55428 : x = cxtofp(x, pr);
1944 55428 : y = gsub(gpowgs(x, n), gen_1);
1945 55428 : z = gdiv(y, gaddgs(gmulsg(n1, y), n2));
1946 55428 : shiftc_inplace(z,1);
1947 55428 : x = gmul(x, gsubsg(1, z));
1948 55428 : if (mask == 1) return gc_GEN(av, gprec_w(x,prec));
1949 17013 : eold = enew;
1950 : }
1951 : }
1952 :
1953 : /* exp(2iPi/d) */
1954 : GEN
1955 2161017 : rootsof1u_cx(ulong n, long prec)
1956 : {
1957 2161017 : switch(n)
1958 : {
1959 15421 : case 1: return gen_1;
1960 4081 : case 2: return gen_m1;
1961 695881 : case 4: return gen_I();
1962 42260 : case 3: case 6: case 12:
1963 : {
1964 42260 : pari_sp av = avma;
1965 42260 : GEN a = (n == 3)? mkfrac(gen_m1,gen_2): ghalf;
1966 42260 : GEN sq3 = sqrtr_abs(utor(3, prec));
1967 42260 : shiftr_inplace(sq3, -1);
1968 42260 : a = (n == 12)? mkcomplex(sq3, a): mkcomplex(a, sq3);
1969 42260 : return gc_GEN(av, a);
1970 : }
1971 849289 : case 8:
1972 : {
1973 849289 : pari_sp av = avma;
1974 849289 : GEN sq2 = sqrtr_abs(utor(2, prec));
1975 849269 : shiftr_inplace(sq2,-1);
1976 849278 : return gc_GEN(av, mkcomplex(sq2, sq2));
1977 : }
1978 : }
1979 554085 : return sqrtnof1(n, prec);
1980 : }
1981 : /* e(a/b) */
1982 : GEN
1983 14616 : rootsof1q_cx(long a, long b, long prec)
1984 : {
1985 14616 : long g = cgcd(a,b);
1986 : GEN z;
1987 14616 : if (g != 1) { a /= g; b /= g; }
1988 14616 : if (b < 0) { b = -b; a = -a; }
1989 14616 : z = rootsof1u_cx(b, prec);
1990 14616 : if (a < 0) { z = conj_i(z); a = -a; }
1991 14616 : return gpowgs(z, a);
1992 : }
1993 :
1994 : /* initializes powers of e(a/b) */
1995 : GEN
1996 15603 : rootsof1powinit(long a, long b, long prec)
1997 : {
1998 15603 : long g = cgcd(a,b);
1999 15603 : if (g != 1) { a /= g; b /= g; }
2000 15603 : if (b < 0) { b = -b; a = -a; }
2001 15603 : a %= b; if (a < 0) a += b;
2002 15603 : return mkvec2(grootsof1(b,prec), mkvecsmall2(a,b));
2003 : }
2004 : /* T = rootsof1powinit(a,b); return e(a/b)^c */
2005 : GEN
2006 12939507 : rootsof1pow(GEN T, long c)
2007 : {
2008 12939507 : GEN vz = gel(T,1), ab = gel(T,2);
2009 12939507 : long a = ab[1], b = ab[2]; /* a >= 0, b > 0 */
2010 12939507 : c %= b; if (c < 0) c += b;
2011 12939507 : a = Fl_mul(a, c, b);
2012 12939507 : return gel(vz, a + 1);
2013 : }
2014 :
2015 : /* exp(2iPi/d), assume d a t_INT */
2016 : GEN
2017 4536 : rootsof1_cx(GEN d, long prec)
2018 : {
2019 4536 : if (lgefint(d) == 3) return rootsof1u_cx((ulong)d[2], prec);
2020 0 : return expIr(divri(Pi2n(1,prec), d));
2021 : }
2022 :
2023 : GEN
2024 42986 : gsqrtn(GEN x, GEN n, GEN *zetan, long prec)
2025 : {
2026 : long i, tx;
2027 : pari_sp av;
2028 : GEN y, z;
2029 42986 : if (typ(n)!=t_INT) pari_err_TYPE("sqrtn",n);
2030 42986 : if (!signe(n)) pari_err_DOMAIN("sqrtn", "n", "=", gen_0, n);
2031 42986 : if (is_pm1(n))
2032 : {
2033 70 : if (zetan) *zetan = gen_1;
2034 70 : return (signe(n) > 0)? gcopy(x): ginv(x);
2035 : }
2036 42916 : if (zetan) *zetan = gen_0;
2037 42916 : tx = typ(x);
2038 42930 : if (is_matvec_t(tx)) pari_APPLY_same(gsqrtn(gel(x,i),n,NULL,prec));
2039 42909 : av = avma;
2040 42909 : switch(tx)
2041 : {
2042 182 : case t_INTMOD:
2043 : {
2044 182 : GEN p = gel(x,1), s;
2045 182 : z = gen_0;
2046 182 : y = cgetg(3,t_INTMOD); gel(y,1) = icopy(p);
2047 182 : if (zetan) { z = cgetg(3,t_INTMOD); gel(z,1) = gel(y,1); }
2048 182 : s = Fp_sqrtn(gel(x,2),n,p,zetan);
2049 161 : if (!s) {
2050 35 : if (zetan) return gc_const(av,gen_0);
2051 28 : if (!BPSW_psp(p)) pari_err_PRIME("sqrtn [modulus]",p);
2052 14 : pari_err_SQRTN("gsqrtn",x);
2053 : }
2054 126 : gel(y,2) = s;
2055 126 : if (zetan) { gel(z,2) = *zetan; *zetan = z; }
2056 126 : return y;
2057 : }
2058 :
2059 56 : case t_PADIC:
2060 56 : y = Qp_sqrtn(x,n,zetan);
2061 49 : if (!y) {
2062 7 : if (zetan) return gen_0;
2063 7 : pari_err_SQRTN("gsqrtn",x);
2064 : }
2065 42 : return y;
2066 :
2067 616 : case t_FFELT: return FF_sqrtn(x,n,zetan);
2068 :
2069 41425 : case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
2070 41425 : i = precision(x); if (i) prec = i;
2071 41425 : if (isint1(x))
2072 7 : y = real_1(prec);
2073 41418 : else if (gequal0(x))
2074 : {
2075 : long b;
2076 21 : if (signe(n) < 0) pari_err_INV("gsqrtn",x);
2077 21 : if (isinexactreal(x))
2078 14 : b = sdivsi(gexpo(x), n);
2079 : else
2080 7 : b = -prec2nbits(prec);
2081 21 : if (typ(x) == t_COMPLEX)
2082 : {
2083 7 : y = cgetg(3,t_COMPLEX);
2084 7 : gel(y,1) = gel(y,2) = real_0_bit(b);
2085 : }
2086 : else
2087 14 : y = real_0_bit(b);
2088 : }
2089 : else
2090 : {
2091 41397 : long nn = itos_or_0(n);
2092 41397 : if (tx == t_INT) { x = itor(x,prec); tx = t_REAL; }
2093 41397 : if (nn > 0 && tx == t_REAL && signe(x) > 0)
2094 31160 : y = sqrtnr(x, nn);
2095 : else
2096 10237 : y = gexp(gdiv(glog(x,prec), n), prec);
2097 41397 : y = gc_upto(av, y);
2098 : }
2099 41425 : if (zetan) *zetan = rootsof1_cx(n, prec);
2100 41425 : return y;
2101 :
2102 7 : case t_QUAD:
2103 7 : return gsqrtn(quadtofp(x, prec), n, zetan, prec);
2104 :
2105 623 : default:
2106 623 : av = avma; if (!(y = toser_i(x))) break;
2107 623 : return gc_upto(av, ser_powfrac(y, ginv(n), prec));
2108 : }
2109 0 : pari_err_TYPE("sqrtn",x);
2110 : return NULL;/* LCOV_EXCL_LINE */
2111 : }
2112 :
2113 : /********************************************************************/
2114 : /** **/
2115 : /** EXP(X) - 1 **/
2116 : /** **/
2117 : /********************************************************************/
2118 : /* exp(|x|) - 1, assume x != 0.
2119 : * For efficiency, x should be reduced mod log(2): if so, we have a < 0 */
2120 : GEN
2121 19006855 : exp1r_abs(GEN x)
2122 : {
2123 19006855 : long l = realprec(x), a = expo(x), b = prec2nbits(l), L, i, n, m, B;
2124 : GEN y, p2, X;
2125 : pari_sp av;
2126 : double d;
2127 :
2128 19006696 : if (b + a <= 0) return mpabs(x);
2129 :
2130 18990480 : y = cgetr(l); av = avma;
2131 18989954 : B = b/3 + BITS_IN_LONG + (BITS_IN_LONG*BITS_IN_LONG)/ b;
2132 18989954 : d = a/2.; m = (long)(d + sqrt(d*d + B)); /* >= 0 */
2133 18989954 : if (m < (-a) * 0.1) m = 0; /* not worth it */
2134 : /* Multiplication is quadratic in this range (l is small, otherwise we
2135 : * use logAGM + Newton). Set Y = 2^(-e-a) x, compute truncated series
2136 : * sum_{k <= n} Y^k/k!: this costs roughly
2137 : * m b^2 + sum_{k <= n} (k e + BITS_IN_LONG)^2
2138 : * bit operations with n ~ b/e, |x| < 2^(1+a), |Y| < 2^(1-e) , m = e+a and
2139 : * b bits of accuracy needed, so
2140 : * B := (b / 3 + BITS_IN_LONG + BITS_IN_LONG^2 / b) ~ m(m-a)
2141 : * we want b ~ 3 m (m-a) or m~b+a hence
2142 : * m = min( a/2 + sqrt(a^2/4 + B), b + a )
2143 : * NB: e ~ (b/3)^(1/2) as b -> oo
2144 : *
2145 : * Truncate the sum at k = n (>= 1), the remainder is
2146 : * sum_{k >= n+1} Y^k / k! < Y^(n+1) / (n+1)! (1-Y) < Y^(n+1) / n!
2147 : * We want Y^(n+1) / n! <= Y 2^-b, hence -n log_2 |Y| + log_2 n! >= b
2148 : * log n! ~ (n + 1/2) log(n+1) - (n+1) + log(2Pi)/2,
2149 : * error bounded by 1/6(n+1) <= 1/12. Finally, we want
2150 : * n (-1/log(2) -log_2 |Y| + log_2(n+1)) >= b */
2151 18989954 : d = m-dbllog2(x)-1/M_LN2; /* ~ -log_2 Y - 1/log(2) */
2152 18990694 : while (d <= 0) { d++; m++; } /* d < 0 can occur from expm1 */
2153 18990688 : L = l + nbits2extraprec(m);
2154 18990660 : b += m;
2155 18990660 : n = (long)(b / d); /* > 0 */
2156 18990660 : if (n == 1)
2157 748112 : n = (long)(b / (d + log2((double)n+1))); /* log ~ const in small ranges */
2158 20359220 : while (n*(d+log2((double)n+1)) < b) n++; /* expect few corrections */
2159 :
2160 18990660 : X = rtor(x,L); shiftr_inplace(X, -m); setsigne(X, 1);
2161 18991867 : if (n == 1) p2 = X;
2162 : else
2163 : {
2164 18991867 : long s = 0, l1 = nbits2prec((long)(d + n + 16));
2165 18991757 : GEN unr = real_1(L);
2166 : pari_sp av2;
2167 :
2168 18991114 : p2 = cgetr(L); av2 = avma;
2169 360164902 : for (i=n; i>=2; i--, set_avma(av2))
2170 : { /* compute X^(n-1)/n! + ... + X/2 + 1 */
2171 : GEN p1, p3;
2172 341246580 : setprec(X,l1); p3 = divru(X,i);
2173 341616858 : l1 += nbits2extraprec(dvmdsBIL(s - expo(p3), &s)<<TWOPOTBITS_IN_LONG);
2174 341467893 : if (l1>L) l1=L;
2175 341467893 : setprec(unr,l1); p1 = addrr_sign(unr,1, i == n? p3: mulrr(p3,p2),1);
2176 341079111 : setprec(p2,l1); affrr(p1,p2); /* p2 <- 1 + (X/i)*p2 */
2177 : }
2178 18989560 : setprec(X,L); p2 = mulrr(X,p2);
2179 : }
2180 :
2181 18991458 : B = prec2nbits(L);
2182 205625409 : for (i = 1; i <= m; i++)
2183 : {
2184 186632832 : if (realprec(p2) > L) setprec(p2,L);
2185 186632832 : if (expo(p2) < -B)
2186 0 : shiftr_inplace(p2, 1); /* 2 + p2 ~ 2 and may blow up accuracy */
2187 : else
2188 186632832 : p2 = mulrr(p2, addsr(2,p2));
2189 : }
2190 18992577 : affrr_fixlg(p2,y); return gc_const(av,y);
2191 : }
2192 :
2193 : GEN
2194 24683 : mpexpm1(GEN x)
2195 : {
2196 24683 : const long s = 6;
2197 24683 : long B, l, sx = signe(x);
2198 : GEN y, z;
2199 : pari_sp av;
2200 24683 : if (!sx) return real_0_bit(expo(x));
2201 24676 : l = realprec(x);
2202 24676 : if (l > maxss(EXPNEWTON_LIMIT, BITS_IN_LONG<<s))
2203 : {
2204 6 : long e = expo(x);
2205 6 : if (e < 0) x = rtor(x, l + nbits2extraprec(-e));
2206 6 : return subrs(mpexp(x), 1);
2207 : }
2208 24670 : if (sx > 0) return exp1r_abs(x);
2209 10298 : B = prec2nbits(l);
2210 10298 : if (cmpsr(-B, x) > 0) return real_m1(l);
2211 : /* compute exp(x) * (1 - exp(-x)) */
2212 10291 : av = avma; y = exp1r_abs(x); /* > 0 */
2213 10291 : if (expo(y) >= -B) { z = addsr(1, y); y = divrr(y, z); }
2214 10291 : setsigne(y, -1);
2215 10291 : return gc_uptoleaf(av, y);
2216 : }
2217 :
2218 : static GEN serexp(GEN x, long prec);
2219 : GEN
2220 26506 : gexpm1(GEN x, long prec)
2221 : {
2222 26506 : switch(typ(x))
2223 : {
2224 4220 : case t_REAL: return mpexpm1(x);
2225 20172 : case t_COMPLEX: return cxexpm1(x,prec);
2226 14 : case t_PADIC: return gsubgs(Qp_exp(x), 1);
2227 2100 : default:
2228 : {
2229 2100 : pari_sp av = avma;
2230 : long ey;
2231 : GEN y;
2232 2100 : if (!(y = toser_i(x))) break;
2233 2079 : ey = valser(y);
2234 2079 : if (ey < 0) pari_err_DOMAIN("expm1","valuation", "<", gen_0, x);
2235 2079 : if (gequal0(y)) return gcopy(y);
2236 2072 : if (ey)
2237 511 : return gc_upto(av, gsubgs(serexp(y,prec), 1));
2238 : else
2239 : {
2240 1561 : GEN e1 = gexpm1(gel(y,2), prec), e = gaddgs(e1,1);
2241 1561 : y = gmul(e, serexp(serchop0(y),prec));
2242 1561 : gel(y,2) = e1;
2243 1561 : return gc_GEN(av, y);
2244 : }
2245 : }
2246 : }
2247 21 : return trans_eval("expm1",gexpm1,x,prec);
2248 : }
2249 : /********************************************************************/
2250 : /** **/
2251 : /** EXP(X) **/
2252 : /** **/
2253 : /********************************************************************/
2254 : static GEN
2255 18929980 : mpexp_basecase(GEN x)
2256 : {
2257 18929980 : pari_sp av = avma;
2258 18929980 : long sh, l = realprec(x);
2259 : GEN y, z;
2260 :
2261 18929980 : y = modlog2(x, &sh);
2262 18930052 : if (!y) { set_avma(av); return real2n(sh, l); }
2263 18930052 : z = addsr(1, exp1r_abs(y));
2264 18929632 : if (signe(y) < 0) z = invr(z);
2265 18930006 : if (sh) {
2266 15716705 : shiftr_inplace(z, sh);
2267 15716602 : if (realprec(z) > l) z = rtor(z, l); /* spurious precision increase */
2268 : }
2269 : #ifdef DEBUG
2270 : {
2271 : GEN t = mplog(z), u = divrr(subrr(x, t),x);
2272 : if (signe(u) && expo(u) > 5-prec2nbits(minss(l,realprec(t))))
2273 : pari_err_BUG("exp");
2274 : }
2275 : #endif
2276 18929839 : return gc_uptoleaf(av, z); /* NOT affrr, precision often increases */
2277 : }
2278 :
2279 : GEN
2280 19076975 : mpexp(GEN x)
2281 : {
2282 19076975 : const long s = 6; /*Initial steps using basecase*/
2283 19076975 : long i, p, l = realprec(x), sh;
2284 : GEN a, t, z;
2285 : ulong mask;
2286 :
2287 19076975 : if (l <= maxss(EXPNEWTON_LIMIT, (BITS_IN_LONG<<s) + 2))
2288 : {
2289 19077045 : if (!signe(x)) return mpexp0(x);
2290 18929927 : return mpexp_basecase(x);
2291 : }
2292 11 : z = cgetr(l); /* room for result */
2293 13 : x = modlog2(x, &sh);
2294 13 : if (!x) { set_avma((pari_sp)(z+lg(z))); return real2n(sh, l); }
2295 13 : constpi(l); /* precompute for later logr_abs() */
2296 13 : mask = quadratic_prec_mask(prec2nbits(l)+BITS_IN_LONG);
2297 168 : for(i=0, p=1; i<s+TWOPOTBITS_IN_LONG; i++) { p <<= 1; if (mask & 1) p-=1; mask >>= 1; }
2298 13 : a = mpexp_basecase(rtor(x, nbits2prec(p)));
2299 13 : x = addrs(x,1);
2300 13 : if (realprec(x) < l+EXTRAPREC64) x = rtor(x, l+EXTRAPREC64);
2301 13 : a = rtor(a, l+EXTRAPREC64); /*append 0s */
2302 13 : t = NULL;
2303 : for(;;)
2304 : {
2305 14 : p <<= 1; if (mask & 1) p--;
2306 14 : mask >>= 1;
2307 14 : setprec(x, nbits2prec(p));
2308 14 : setprec(a, nbits2prec(p));
2309 14 : t = mulrr(a, subrr(x, logr_abs(a))); /* a (x - log(a)) */
2310 14 : if (mask == 1) break;
2311 1 : affrr(t, a); set_avma((pari_sp)a);
2312 : }
2313 13 : affrr(t,z);
2314 13 : if (sh) shiftr_inplace(z, sh);
2315 13 : return gc_const((pari_sp)z, z);
2316 : }
2317 :
2318 : /* x != 0; k = ceil(tn / (te-1)), t = p-1 */
2319 : long
2320 98 : Qp_exp_prec(GEN x)
2321 : {
2322 98 : long e = valp(x), n = precp(x);
2323 : ulong a, b, q, r, p, t;
2324 :
2325 98 : if (e < 1) return -1;
2326 77 : if (e > n) return 1;
2327 77 : p = itos_or_0(padic_p(x));
2328 77 : if (!p) return n / e + 1;
2329 77 : if (p == 2) return e < 2? -1: ceildivuu(n, e - 1);
2330 : /* n >= e > 0, n = qe + r */
2331 : /* tn = q (te-1) + rt + q = (q+1)(te-1) - t(e-r) + q + 1 */
2332 63 : t = p - 1;
2333 63 : if (e == 1) return n + ceildivuu(n, t - 1);
2334 0 : q = n / e;
2335 0 : r = n % e; /* k = q + 1 if rt + q < te */
2336 0 : a = umuluu_or_0(e - r, t); if (!a || a > q) return q + 1;
2337 0 : b = umuluu_or_0(e, t); if (!b) return q + 2;
2338 0 : return q + 1 + ceildivuu(q + 1 - a, b - 1);
2339 : }
2340 :
2341 : static GEN
2342 109557 : Qp_exp_safe(GEN x)
2343 : {
2344 109557 : pari_sp av = avma;
2345 109557 : GEN p = padic_p(x), a = padic_u(x), z;
2346 109557 : long d = precp(x), v = valp(x), e = d+v;
2347 109557 : if (gequal0(x)) return gaddgs(x,1);
2348 107961 : if (v < (equaliu(p,2)? 2:1)) return NULL;
2349 107954 : z = Zp_exp(mulii(a,powiu(p,v)), p, e);
2350 107955 : return gc_upto(av, Z_to_padic(z, p, e));
2351 : }
2352 :
2353 : GEN
2354 109095 : Qp_exp(GEN x)
2355 : {
2356 109095 : GEN y = Qp_exp_safe(x);
2357 109096 : if (!y) pari_err_DOMAIN("gexp(t_PADIC)","argument","",gen_0,x);
2358 109089 : return y;
2359 : }
2360 :
2361 : static GEN
2362 49 : cos_p(GEN x)
2363 : {
2364 : long k;
2365 : pari_sp av;
2366 : GEN x2, y;
2367 :
2368 49 : if (gequal0(x)) return gaddgs(x,1);
2369 28 : k = Qp_exp_prec(x);
2370 28 : if (k < 0) return NULL;
2371 21 : av = avma; x2 = gsqr(x);
2372 21 : if (k & 1) k--;
2373 105 : for (y=gen_1; k; k-=2)
2374 : {
2375 84 : GEN t = gdiv(gmul(y,x2), muluu(k, k-1));
2376 84 : y = gsubsg(1, t);
2377 : }
2378 21 : return gc_upto(av, y);
2379 : }
2380 : static GEN
2381 63 : sin_p(GEN x)
2382 : {
2383 : long k;
2384 : pari_sp av;
2385 : GEN x2, y;
2386 :
2387 63 : if (gequal0(x)) return gcopy(x);
2388 42 : k = Qp_exp_prec(x);
2389 42 : if (k < 0) return NULL;
2390 28 : av = avma; x2 = gsqr(x);
2391 28 : if (k & 1) k--;
2392 133 : for (y=gen_1; k; k-=2)
2393 : {
2394 105 : GEN t = gdiv(gmul(y,x2), muluu(k, k+1));
2395 105 : y = gsubsg(1, t);
2396 : }
2397 28 : return gc_upto(av, gmul(y, x));
2398 : }
2399 :
2400 : static GEN
2401 4688205 : cxexp(GEN x, long prec)
2402 : {
2403 4688205 : GEN r, p1, p2, y = cgetg(3,t_COMPLEX);
2404 4688152 : pari_sp av = avma, tetpil;
2405 : long l;
2406 4688152 : l = precision(x); if (l > prec) prec = l;
2407 4688208 : if (gequal0(gel(x,1)))
2408 : {
2409 346662 : gsincos(gel(x,2),&gel(y,2),&gel(y,1),prec);
2410 346666 : return y;
2411 : }
2412 4341529 : r = gexp(gel(x,1),prec);
2413 4341676 : gsincos(gel(x,2),&p2,&p1,prec);
2414 4341904 : tetpil = avma;
2415 4341904 : gel(y,1) = gmul(r,p1);
2416 4341791 : gel(y,2) = gmul(r,p2);
2417 4341786 : gc_slice_unsafe(av,tetpil,y+1,2);
2418 4341858 : return y;
2419 : }
2420 :
2421 : /* given a t_SER x^v s(x), with s(0) != 0, return x^v(s - s(0)), shallow */
2422 : GEN
2423 37576 : serchop0(GEN s)
2424 : {
2425 37576 : long i, l = lg(s);
2426 : GEN y;
2427 37576 : if (l == 2) return s;
2428 37576 : if (l == 3 && isexactzero(gel(s,2))) return s;
2429 37576 : y = cgetg(l, t_SER); y[1] = s[1];
2430 164990 : gel(y,2) = gen_0; for (i=3; i <l; i++) gel(y,i) = gel(s,i);
2431 37576 : return normalizeser(y);
2432 : }
2433 :
2434 : GEN
2435 42 : serchop_i(GEN s, long n)
2436 : {
2437 42 : long i, m, l = lg(s);
2438 : GEN y;
2439 42 : if (l == 2 || (l == 3 && isexactzero(gel(s,2))))
2440 : {
2441 14 : if (valser(s) < n) { s = shallowcopy(s); setvalser(s,n); }
2442 14 : return s;
2443 : }
2444 28 : m = n - valser(s); if (m < 0) return s;
2445 21 : if (l-m <= 2) return zeroser(varn(s), n);
2446 14 : y = cgetg(l-m, t_SER); y[1] = s[1]; setvalser(y, valser(y)+m);
2447 42 : for (i=m+2; i < l; i++) gel(y,i-m) = gel(s,i);
2448 14 : return normalizeser(y);
2449 : }
2450 : GEN
2451 42 : serchop(GEN s, long n)
2452 : {
2453 42 : pari_sp av = avma;
2454 42 : if (typ(s) != t_SER) pari_err_TYPE("serchop",s);
2455 42 : return gc_GEN(av, serchop_i(s,n));
2456 : }
2457 :
2458 : static GEN
2459 83433 : serexp(GEN x, long prec)
2460 : {
2461 83433 : long i, j, lx, ly, mi, e = valser(x);
2462 : GEN y, xd, yd;
2463 : pari_sp av;
2464 :
2465 83433 : if (e < 0) pari_err_DOMAIN("exp","valuation", "<", gen_0, x);
2466 83426 : if (gequal0(x)) return gaddsg(1,x);
2467 70497 : lx = lg(x);
2468 70497 : if (e)
2469 : {
2470 : GEN X;
2471 55699 : ly = lx+e; y = cgetg(ly,t_SER);
2472 566888 : mi = lx-1; while (mi>=3 && isrationalzero(gel(x,mi))) mi--;
2473 55699 : mi += e-2;
2474 55699 : y[1] = evalsigne(1) | _evalvalser(0) | evalvarn(varn(x));
2475 : /* zd[i] = coefficient of X^i in z */
2476 55699 : xd = x+2-e; yd = y+2; ly -= 2;
2477 55699 : X = gel(xd,e); if (e != 1) X = gmulgu(X, e); /* left on stack */
2478 55699 : X = isint1(X)? NULL: X;
2479 55699 : gel(yd,0) = gen_1;
2480 56070 : for (i = 1; i < e; i++) gel(yd,i) = gen_0;
2481 664615 : for ( ; i < ly; i++)
2482 : {
2483 608916 : GEN t = gel(yd,i-e);
2484 608916 : long J = minss(i, mi);
2485 608916 : av = avma; if (X) t = gmul(t, X);
2486 2578765 : for (j = e + 1; j <= J; j++)
2487 1969849 : t = gadd(t, gmulgu(gmul(gel(xd,j),gel(yd,i-j)), j));
2488 608916 : gel(yd,i) = gc_upto(av, gdivgu(t, i));
2489 : }
2490 55699 : return y;
2491 : }
2492 14798 : av = avma;
2493 14798 : return gc_upto(av, gmul(gexp(gel(x,2),prec), serexp(serchop0(x),prec)));
2494 : }
2495 :
2496 : static GEN
2497 1468573 : expQ(GEN x, long prec)
2498 : {
2499 1468573 : GEN p, q, z, z0 = NULL;
2500 : pari_sp av;
2501 1468573 : long n, nmax, s, e, b = prec2nbits(prec);
2502 : double ex;
2503 : struct abpq_res R;
2504 : struct abpq S;
2505 :
2506 1468573 : if (typ(x) == t_INT)
2507 : {
2508 24683 : if (!signe(x)) return real_1(prec);
2509 24612 : p = x; q = gen_1;
2510 24612 : e = expi(p);
2511 24612 : if (e > b) return mpexp(itor(x, prec));
2512 : }
2513 : else
2514 : {
2515 1443890 : long ep, eq, B = usqrt(b) / 2;
2516 1443890 : p = gel(x,1); ep = expi(p);
2517 1443890 : q = gel(x,2); eq = expi(q);
2518 1443890 : if (ep > B || eq > B) return mpexp(fractor(x, prec));
2519 14637 : e = ep - eq;
2520 14637 : if (e < -3) prec += nbits2extraprec(-e); /* see addrr 'extend' rule */
2521 : }
2522 39249 : if (e > 2) { z0 = cgetr(prec); prec += EXTRAPREC64; b += BITS_IN_LONG; }
2523 39249 : z = cgetr(prec); av = avma;
2524 39248 : if (e > 0)
2525 : { /* simplify x/2^e = p / (q * 2^e) */
2526 2478 : long v = minss(e, vali(p));
2527 2478 : if (v) p = shifti(p, -v);
2528 2478 : if (e - v) q = shifti(q, e - v);
2529 : }
2530 39248 : s = signe(p);
2531 39248 : if (s < 0) p = negi(p);
2532 39251 : ex = exp2(dbllog2(x) - e) * 2.718281828; /* exp(1) * x / 2^e, x / 2^e < 2 */
2533 39254 : nmax = (long)(1 + exp(dbllambertW0(M_LN2 * b / ex)) * ex);
2534 39252 : abpq_init(&S, nmax);
2535 39276 : S.a[0] = S.b[0] = S.p[0] = S.q[0] = gen_1;
2536 3368343 : for (n = 1; n <= nmax; n++)
2537 : {
2538 3329100 : S.a[n] = gen_1;
2539 3329100 : S.b[n] = gen_1;
2540 3329100 : S.p[n] = p;
2541 3329100 : S.q[n] = muliu(q, n);
2542 : }
2543 39243 : abpq_sum(&R, 0, nmax, &S);
2544 39256 : if (s > 0) rdiviiz(R.T, R.Q, z); else rdiviiz(R.Q, R.T, z);
2545 39257 : if (e > 0)
2546 : {
2547 17136 : q = z; while (e--) q = sqrr(q);
2548 2478 : if (z0) { affrr(q, z0); z = z0; } else affrr(q,z);
2549 : }
2550 39257 : return gc_const(av,z);
2551 : }
2552 :
2553 : GEN
2554 18635163 : gexp(GEN x, long prec)
2555 : {
2556 18635163 : switch(typ(x))
2557 : {
2558 1468573 : case t_INT: case t_FRAC: return expQ(x, prec);
2559 11128375 : case t_REAL: return mpexp(x);
2560 4688168 : case t_COMPLEX: return cxexp(x,prec);
2561 70 : case t_PADIC: return Qp_exp(x);
2562 1349977 : default:
2563 : {
2564 1349977 : pari_sp av = avma;
2565 : GEN y;
2566 1349977 : if (!(y = toser_i(x))) break;
2567 66563 : return gc_upto(av, serexp(y,prec));
2568 : }
2569 : }
2570 1283987 : return trans_eval("exp",gexp,x,prec);
2571 : }
2572 :
2573 : /********************************************************************/
2574 : /** **/
2575 : /** AGM(X, Y) **/
2576 : /** **/
2577 : /********************************************************************/
2578 : static int
2579 16365937 : agmr_gap(GEN a, GEN b, long L)
2580 : {
2581 16365937 : GEN d = subrr(b, a);
2582 16365787 : return (signe(d) && expo(d) - expo(b) >= L);
2583 : }
2584 : /* assume x > 0 */
2585 : static GEN
2586 1104958 : agm1r_abs(GEN x)
2587 : {
2588 1104958 : long l = realprec(x), L = 5-prec2nbits(l);
2589 1104958 : GEN a1, b1, y = cgetr(l);
2590 1104957 : pari_sp av = avma;
2591 :
2592 1104957 : a1 = addrr(real_1(l), x); shiftr_inplace(a1, -1);
2593 1104956 : b1 = sqrtr_abs(x);
2594 16365958 : while (agmr_gap(a1,b1,L))
2595 : {
2596 15260863 : GEN a = a1;
2597 15260863 : a1 = addrr(a,b1); shiftr_inplace(a1, -1);
2598 15260958 : b1 = sqrtr_abs(mulrr(a,b1));
2599 : }
2600 1104927 : affrr_fixlg(a1,y); return gc_const(av,y);
2601 : }
2602 :
2603 : struct agmcx_gap_t { long L, ex, cnt; };
2604 :
2605 : static void
2606 366075 : agmcx_init(GEN x, long *prec, struct agmcx_gap_t *S)
2607 : {
2608 366075 : long l = precision(x);
2609 366075 : if (l) *prec = l;
2610 366075 : S->L = 1-prec2nbits(*prec);
2611 366075 : S->cnt = 0;
2612 366075 : S->ex = LONG_MAX;
2613 366075 : }
2614 :
2615 : static long
2616 366075 : agmcx_a_b(GEN x, GEN *a1, GEN *b1, long prec)
2617 : {
2618 366075 : long rotate = 0;
2619 366075 : if (gsigne(real_i(x))<0)
2620 : { /* Rotate by +/-Pi/2, so that the choice of the principal square
2621 : * root gives the optimal AGM. So a1 = +/-I*a1, b1=sqrt(-x). */
2622 11655 : if (gsigne(imag_i(x))<0) { *a1=mulcxI(*a1); rotate=-1; }
2623 11137 : else { *a1=mulcxmI(*a1); rotate=1; }
2624 11655 : x = gneg(x);
2625 : }
2626 366075 : *b1 = gsqrt(x, prec);
2627 366075 : return rotate;
2628 : }
2629 : /* return 0 if we must stop the AGM loop (a=b or a ~ b), 1 otherwise */
2630 : static int
2631 5540946 : agmcx_gap(GEN a, GEN b, struct agmcx_gap_t *S)
2632 : {
2633 5540946 : GEN d = gsub(b, a);
2634 5540946 : long ex = S->ex;
2635 5540946 : S->ex = gexpo(d);
2636 5540946 : if (gequal0(d) || S->ex - gexpo(b) < S->L) return 0;
2637 : /* if (S->ex >= ex) we're no longer making progress; twice in a row */
2638 5279147 : if (S->ex < ex) S->cnt = 0;
2639 : else
2640 209085 : if (S->cnt++) return 0;
2641 5174871 : return 1;
2642 : }
2643 : static GEN
2644 337326 : agm1cx(GEN x, long prec)
2645 : {
2646 : struct agmcx_gap_t S;
2647 : GEN a1, b1;
2648 337326 : pari_sp av = avma;
2649 : long rotate;
2650 337326 : agmcx_init(x, &prec, &S);
2651 337326 : a1 = gtofp(gmul2n(gadd(real_1(prec), x), -1), prec);
2652 337326 : rotate = agmcx_a_b(x, &a1, &b1, prec);
2653 5358931 : while (agmcx_gap(a1,b1,&S))
2654 : {
2655 5021605 : GEN a = a1;
2656 5021605 : a1 = gmul2n(gadd(a,b1),-1);
2657 5021605 : b1 = gsqrt(gmul(a,b1), prec);
2658 : }
2659 337326 : if (rotate) a1 = rotate>0 ? mulcxI(a1):mulcxmI(a1);
2660 337326 : return gc_GEN(av,a1);
2661 : }
2662 :
2663 : GEN
2664 28749 : zellagmcx(GEN a0, GEN b0, GEN r, GEN t, long prec)
2665 : {
2666 : struct agmcx_gap_t S;
2667 28749 : pari_sp av = avma;
2668 28749 : GEN x = gdiv(a0, b0), a1, b1;
2669 : long rotate;
2670 28749 : agmcx_init(x, &prec, &S);
2671 28749 : a1 = gtofp(gmul2n(gadd(real_1(prec), x), -1), prec);
2672 28749 : r = gsqrt(gdiv(gmul(a1,gaddgs(r, 1)),gadd(r, x)), prec);
2673 28749 : t = gmul(r, t);
2674 28749 : rotate = agmcx_a_b(x, &a1, &b1, prec);
2675 182015 : while (agmcx_gap(a1,b1,&S))
2676 : {
2677 153266 : GEN a = a1, b = b1;
2678 153266 : a1 = gmul2n(gadd(a,b),-1);
2679 153266 : b1 = gsqrt(gmul(a,b), prec);
2680 153266 : r = gsqrt(gdiv(gmul(a1,gaddgs(r, 1)),gadd(gmul(b, r), a )), prec);
2681 153266 : t = gmul(r, t);
2682 : }
2683 28749 : if (rotate) a1 = rotate>0 ? mulcxI(a1):mulcxmI(a1);
2684 28749 : a1 = gmul(a1, b0);
2685 28749 : t = gatan(gdiv(a1,t), prec);
2686 : /* send t to the fundamental domain if necessary */
2687 28749 : if (gsigne(real_i(t))<0) t = gadd(t, mppi(prec));
2688 28749 : return gc_upto(av,gdiv(t,a1));
2689 : }
2690 :
2691 : static long
2692 49 : ser_cmp_expo(GEN A, GEN B)
2693 : {
2694 49 : long e = -(long)HIGHEXPOBIT, d = valser(B) - valser(A);
2695 49 : long i, la = lg(A), v = varn(B);
2696 9849 : for (i = 2; i < la; i++)
2697 : {
2698 9800 : GEN a = gel(A,i), b;
2699 : long ei;
2700 9800 : if (isexactzero(a)) continue;
2701 9800 : b = polcoef_i(B, i-2 + d, v);
2702 9800 : ei = gexpo(a);
2703 9800 : if (!isexactzero(b)) ei -= gexpo(b);
2704 9800 : e = maxss(e, ei);
2705 : }
2706 49 : return e;
2707 : }
2708 :
2709 : static GEN
2710 21 : ser_agm1(GEN y, long prec)
2711 : {
2712 21 : GEN a1 = y, b1 = gen_1;
2713 21 : long l = lg(y)-2, l2 = 6-prec2nbits(prec), eold = LONG_MAX;
2714 : for(;;)
2715 84 : {
2716 105 : GEN a = a1, p1;
2717 105 : a1 = gmul2n(gadd(a,b1),-1);
2718 105 : b1 = gsqrt(gmul(a,b1), prec);
2719 105 : p1 = gsub(b1,a1);
2720 105 : if (isinexactreal(p1))
2721 : {
2722 49 : long e = ser_cmp_expo(p1, b1);
2723 49 : if (e < l2 || e >= eold) break;
2724 42 : eold = e;
2725 : }
2726 56 : else if (valser(p1)-valser(b1) >= l || gequal0(p1)) break;
2727 : }
2728 21 : return a1;
2729 : }
2730 :
2731 : /* agm(1,x) */
2732 : static GEN
2733 112280 : agm1(GEN x, long prec)
2734 : {
2735 : GEN y;
2736 : pari_sp av;
2737 :
2738 112280 : if (gequal0(x)) return gcopy(x);
2739 112280 : switch(typ(x))
2740 : {
2741 28 : case t_INT:
2742 28 : if (!is_pm1(x)) break;
2743 21 : return (signe(x) > 0)? real_1(prec): real_0(prec);
2744 :
2745 74732 : case t_REAL: return signe(x) > 0? agm1r_abs(x): agm1cx(x, prec);
2746 :
2747 37380 : case t_COMPLEX:
2748 37380 : if (gequal0(gel(x,2))) return agm1(gel(x,1), prec);
2749 37345 : return agm1cx(x, prec);
2750 :
2751 14 : case t_PADIC:
2752 : {
2753 14 : GEN a1 = x, b1 = gen_1;
2754 14 : long l = precp(x);
2755 14 : av = avma;
2756 : for(;;)
2757 14 : {
2758 28 : GEN a = a1, p1;
2759 : long ep;
2760 28 : a1 = gmul2n(gadd(a,b1),-1);
2761 28 : a = gmul(a,b1);
2762 28 : b1 = Qp_sqrt(a); if (!b1) pari_err_SQRTN("Qp_sqrt",a);
2763 21 : p1 = gsub(b1,a1); ep = valp(p1)-valp(b1);
2764 21 : if (ep<=0) { b1 = gneg_i(b1); p1 = gsub(b1,a1); ep=valp(p1)-valp(b1); }
2765 21 : if (ep >= l || gequal0(p1)) return gc_GEN(av,a1);
2766 : }
2767 : }
2768 :
2769 126 : default:
2770 126 : av = avma; if (!(y = toser_i(x))) break;
2771 21 : return gc_GEN(av, ser_agm1(y, prec));
2772 : }
2773 112 : return trans_eval("agm",agm1,x,prec);
2774 : }
2775 :
2776 : GEN
2777 111986 : agm(GEN x, GEN y, long prec)
2778 : {
2779 : pari_sp av;
2780 111986 : if (is_matvec_t(typ(y)))
2781 : {
2782 14 : if (is_matvec_t(typ(x))) pari_err_TYPE2("agm",x,y);
2783 7 : swap(x, y);
2784 : }
2785 111979 : if (gequal0(y)) return gcopy(y);
2786 111979 : av = avma;
2787 111979 : return gc_upto(av, gmul(y, agm1(gdiv(x,y), prec)));
2788 : }
2789 :
2790 : /* b2 != 0 */
2791 : static GEN
2792 147 : ellK_i(GEN b2, long prec)
2793 147 : { return gdiv(Pi2n(-1, prec), agm1(gsqrt(b2, prec), prec)); }
2794 : GEN
2795 140 : ellK(GEN k, long prec)
2796 : {
2797 140 : pari_sp av = avma;
2798 140 : GEN k2 = gsqr(k), b2 = gsubsg(1, k2);
2799 140 : if (gequal0(b2)) pari_err_DOMAIN("ellK", "k^2", "=", gen_1, k2);
2800 133 : return gc_upto(av, ellK_i(b2, prec));
2801 : }
2802 :
2803 : static int
2804 84 : magm_gap(GEN a, GEN b, long L)
2805 : {
2806 84 : GEN d = gsub(b, a);
2807 84 : return !gequal0(d) && gexpo(d) - gexpo(b) >= L;
2808 : }
2809 :
2810 : /* http://www.ams.org/notices/201208/rtx120801094p.pdf
2811 : * An Eloquent Formula for the Perimeter of an Ellipse
2812 : * Semjon Adlaj, Notices of the AMS */
2813 : static GEN
2814 14 : magm(GEN a, GEN b, long prec)
2815 : {
2816 14 : long L = -prec2nbits(prec) + 16;
2817 14 : GEN c = gen_0;
2818 84 : while (magm_gap(a, b, L))
2819 : {
2820 70 : GEN u = gsqrt(gmul(gsub(a, c), gsub(b, c)), prec);
2821 70 : a = gmul2n(gadd(a, b), -1);
2822 70 : b = gadd(c, u); c = gsub(c, u);
2823 : }
2824 14 : return gmul2n(gadd(a, b), -1);
2825 : }
2826 :
2827 : GEN
2828 21 : ellE(GEN k, long prec)
2829 : {
2830 21 : pari_sp av = avma;
2831 21 : GEN b2 = gsubsg(1, gsqr(k));
2832 21 : if (gequal0(b2)) { set_avma(av); return real_1(prec); }
2833 14 : return gc_upto(av, gmul(ellK_i(b2, prec), magm(gen_1, b2, prec)));
2834 : }
2835 :
2836 : /********************************************************************/
2837 : /** **/
2838 : /** LOG(X) **/
2839 : /** **/
2840 : /********************************************************************/
2841 : /* log(2) = 18*atanh(1/26)-2*atanh(1/4801)+8*atanh(1/8749)
2842 : * faster than 10*atanh(1/17)+4*atanh(13/499) for all precisions,
2843 : * and than Pi/2M(1,4/2^n) ~ n log(2) for bitprec at least up to 10^8 */
2844 : static GEN
2845 42098 : log2_split(long prec)
2846 : {
2847 42098 : GEN u = atanhuu(1, 26, prec);
2848 42101 : GEN v = atanhuu(1, 4801, prec);
2849 42099 : GEN w = atanhuu(1, 8749, prec);
2850 42094 : shiftr_inplace(v, 1); setsigne(v, -1);
2851 42096 : shiftr_inplace(w, 3);
2852 42095 : return addrr(mulur(18, u), addrr(v, w));
2853 : }
2854 : GEN
2855 29049765 : constlog2(long prec)
2856 : {
2857 : pari_sp av;
2858 : GEN tmp;
2859 29049765 : if (glog2 && realprec(glog2) >= prec) return glog2;
2860 :
2861 41959 : tmp = cgetr_block(prec);
2862 42098 : av = avma;
2863 42098 : affrr(log2_split(prec+EXTRAPREC64), tmp);
2864 42091 : swap_clone(&glog2,tmp);
2865 42099 : return gc_const(av,glog2);
2866 : }
2867 :
2868 : GEN
2869 29049727 : mplog2(long prec) { return rtor(constlog2(prec), prec); }
2870 :
2871 : /* dont check that q != 2^expo(q), done in logr_abs */
2872 : static GEN
2873 1030264 : logagmr_abs(GEN q)
2874 : {
2875 1030264 : long prec = realprec(q), e = expo(q), lim;
2876 1030264 : GEN z = cgetr(prec), y, Q, _4ovQ;
2877 1030260 : pari_sp av = avma;
2878 :
2879 1030260 : incrprec(prec);
2880 1030260 : lim = prec2nbits(prec) >> 1;
2881 1030260 : Q = rtor(q,prec);
2882 1030267 : shiftr_inplace(Q,lim-e); setsigne(Q,1);
2883 :
2884 1030265 : _4ovQ = invr(Q); shiftr_inplace(_4ovQ, 2); /* 4/Q */
2885 : /* Pi / 2agm(1, 4/Q) ~ log(Q), q = Q * 2^(e-lim) */
2886 1030268 : y = divrr(Pi2n(-1, prec), agm1r_abs(_4ovQ));
2887 1030268 : y = addrr(y, mulsr(e - lim, mplog2(prec)));
2888 1030266 : affrr_fixlg(y, z); return gc_const(av,z);
2889 : }
2890 :
2891 : /* sum_{k >= 0} y^(2k+1) / (2k+1), y close to 0 */
2892 : static GEN
2893 11923756 : logr_aux(GEN y)
2894 : {
2895 11923756 : long k, L = realprec(y); /* should be ~ l+1 - (k-2) */
2896 : /* log(x) = log(1+y) - log(1-y) = 2 sum_{k odd} y^k / k
2897 : * Truncate the sum at k = 2n+1, the remainder is
2898 : * 2 sum_{k >= 2n+3} y^k / k < 2y^(2n+3) / (2n+3)(1-y) < y^(2n+3)
2899 : * We want y^(2n+3) < y 2^(-prec2nbits(L)), hence
2900 : * n+1 > -prec2nbits(L) /-log_2(y^2) */
2901 11923756 : double d = -2*dbllog2r(y); /* ~ -log_2(y^2) */
2902 11923717 : k = (long)(2*(prec2nbits(L) / d));
2903 11923678 : k |= 1;
2904 11923678 : if (k >= 3)
2905 : {
2906 11884389 : GEN T, S = cgetr(L), y2 = sqrr(y), unr = real_1(L);
2907 11884379 : pari_sp av = avma;
2908 11884379 : long s = 0, incs = (long)d, l1 = nbits2prec((long)d);
2909 11884400 : setprec(S, l1);
2910 11884353 : setprec(unr,l1); affrr(divru(unr,k), S);
2911 214410971 : for (k -= 2;; k -= 2) /* k = 2n+1, ..., 1 */
2912 : { /* S = y^(2n+1-k)/(2n+1) + ... + 1 / k */
2913 214410971 : setprec(y2, l1); T = mulrr(S,y2);
2914 214511368 : if (k == 1) break;
2915 :
2916 202626912 : l1 += nbits2extraprec(dvmdsBIL(s + incs, &s)<<TWOPOTBITS_IN_LONG);
2917 202608947 : if (l1>L) l1=L;
2918 202608947 : setprec(S, l1);
2919 202599419 : setprec(unr,l1);
2920 202575194 : affrr(addrr(divru(unr, k), T), S); set_avma(av);
2921 : }
2922 : /* k = 1 special-cased for eficiency */
2923 11884456 : y = mulrr(y, addsr(1,T)); /* = log(X)/2 */
2924 : }
2925 11923763 : return y;
2926 : }
2927 : /*return log(|x|), assuming x != 0 */
2928 : GEN
2929 13795058 : logr_abs(GEN X)
2930 : {
2931 13795058 : long EX, L, m, k, a, b, l = lg(X), p = realprec(X);
2932 : GEN z, x, y;
2933 : ulong u;
2934 : double d;
2935 :
2936 : /* Assuming 1 < x < 2, we want delta = x-1, 1-x/2, 1-1/x, or 2/x-1 small.
2937 : * We have 2/x-1 > 1-x/2, 1-1/x < x-1. So one should be choosing between
2938 : * 1-1/x and 1-x/2 ( crossover sqrt(2), worse ~ 0.29 ). To avoid an inverse,
2939 : * we choose between x-1 and 1-x/2 ( crossover 4/3, worse ~ 0.33 ) */
2940 13795058 : EX = expo(X);
2941 13795058 : u = uel(X,2);
2942 13795058 : k = 2;
2943 13795058 : if (u > (~0UL / 3) * 2) { /* choose 1-x/2 */
2944 7794100 : EX++; u = ~u;
2945 7907279 : while (!u && ++k < l) { u = uel(X,k); u = ~u; }
2946 : } else { /* choose x - 1 */
2947 6000958 : u &= ~HIGHBIT; /* u - HIGHBIT, assuming HIGHBIT set */
2948 7754445 : while (!u && ++k < l) u = uel(X,k);
2949 : }
2950 13795058 : if (k == l) return EX? mulsr(EX, mplog2(p)): real_0(p);
2951 12953847 : a = bit_accuracy(k) + bfffo(u); /* ~ -log2 |1-x| */
2952 12953910 : L = p+EXTRAPRECWORD;
2953 12953910 : b = prec2nbits(L - (bit_accuracy(k))); /* take loss of accuracy into account */
2954 12953877 : if (b > 24*a*log2(prec2lg(L)) && p > LOGAGM_LIMIT) return logagmr_abs(X);
2955 :
2956 11923654 : z = cgetr(EX? p: p - bit_accuracy(k));
2957 :
2958 : /* Multiplication is quadratic in this range (l is small, otherwise we
2959 : * use AGM). Set Y = x^(1/2^m), y = (Y - 1) / (Y + 1) and compute truncated
2960 : * series sum y^(2k+1)/(2k+1): the costs is less than
2961 : * m b^2 + sum_{k <= n} ((2k+1) e + BITS_IN_LONG)^2
2962 : * bit operations with |x-1| < 2^(1-a), |Y| < 2^(1-e) , m = e-a and b bits of
2963 : * accuracy needed (+ BITS_IN_LONG since bit accuracies increase by
2964 : * increments of BITS_IN_LONG), so
2965 : * 4n^3/3 e^2 + n^2 2e BITS_IN_LONG+ n BITS_IN_LONG ~ m b^2, with n ~ b/2e
2966 : * or b/6e + BITS_IN_LONG/2e + BITS_IN_LONG/2be ~ m
2967 : * B := (b / 6 + BITS_IN_LONG/2 + BITS_IN_LONG^2 / 2b) ~ m(m+a)
2968 : * m = min( -a/2 + sqrt(a^2/4 + B), b - a )
2969 : * NB: e ~ (b/6)^(1/2) as b -> oo
2970 : * Instead of the above pessimistic estimate for the cost of the sum, use
2971 : * optimistic estimate (BITS_IN_LONG -> 0) */
2972 11923644 : d = -a/2.; m = (long)(d + sqrt(d*d + b/6)); /* >= 0 */
2973 :
2974 11923644 : if (m > b-a) m = b-a;
2975 11923644 : if (m < 0.2*a) m = 0; else L += nbits2extraprec(m);
2976 11923633 : x = rtor(X,L);
2977 11923716 : setsigne(x,1); shiftr_inplace(x,-EX);
2978 : /* 2/3 < x < 4/3 */
2979 70318734 : for (k=1; k<=m; k++) x = sqrtr_abs(x);
2980 :
2981 11923855 : y = divrr(subrs(x,1), addrs(x,1)); /* = (x-1) / (x+1), close to 0 */
2982 11923721 : y = logr_aux(y); /* log(1+y) - log(1-y) = log(x) */
2983 11923694 : shiftr_inplace(y, m + 1);
2984 11923617 : if (EX) y = addrr(y, mulsr(EX, mplog2(p+EXTRAPRECWORD)));
2985 11923455 : affrr_fixlg(y, z); return gc_const((pari_sp)z, z);
2986 : }
2987 :
2988 : /* assume Im(q) != 0 and precision(q) >= prec. Compute log(q) with accuracy
2989 : * prec [disregard input accuracy] */
2990 : GEN
2991 299939 : logagmcx(GEN q, long prec)
2992 : {
2993 299939 : GEN z = cgetc(prec), y, Q, a, b;
2994 : long lim, e, ea, eb;
2995 299939 : pari_sp av = avma;
2996 299939 : int neg = 0;
2997 :
2998 299939 : incrprec(prec);
2999 299939 : if (gsigne(gel(q,1)) < 0) { q = gneg(q); neg = 1; }
3000 299939 : lim = prec2nbits(prec) >> 1;
3001 299939 : Q = gtofp(q, prec);
3002 299939 : a = gel(Q,1);
3003 299939 : b = gel(Q,2);
3004 299939 : if (gequal0(a)) {
3005 0 : affrr_fixlg(logr_abs(b), gel(z,1));
3006 0 : y = Pi2n(-1, prec);
3007 0 : if (signe(b) < 0) setsigne(y, -1);
3008 0 : affrr_fixlg(y, gel(z,2)); return gc_const(av,z);
3009 : }
3010 299939 : ea = expo(a);
3011 299939 : eb = expo(b);
3012 299939 : e = ea <= eb ? lim - eb : lim - ea;
3013 299939 : shiftr_inplace(a, e);
3014 299939 : shiftr_inplace(b, e);
3015 :
3016 : /* Pi / 2agm(1, 4/Q) ~ log(Q), q = Q * 2^e */
3017 299939 : y = gdiv(Pi2n(-1, prec), agm1cx( gdivsg(4, Q), prec ));
3018 299939 : a = gel(y,1);
3019 299939 : b = gel(y,2);
3020 299939 : a = addrr(a, mulsr(-e, mplog2(prec)));
3021 299939 : if (realprec(a) <= LOWDEFAULTPREC) a = real_0_bit(expo(a));
3022 418458 : if (neg) b = gsigne(b) <= 0? gadd(b, mppi(prec))
3023 118519 : : gsub(b, mppi(prec));
3024 299939 : affrr_fixlg(a, gel(z,1));
3025 299939 : affrr_fixlg(b, gel(z,2)); return gc_const(av,z);
3026 : }
3027 :
3028 : GEN
3029 203666 : mplog(GEN x)
3030 : {
3031 203666 : if (signe(x)<=0) pari_err_DOMAIN("mplog", "argument", "<=", gen_0, x);
3032 203666 : return logr_abs(x);
3033 : }
3034 :
3035 : /* pe = p^e, p prime, 0 < x < pe a t_INT coprime to p. Return the (p-1)-th
3036 : * root of 1 in (Z/pe)^* congruent to x mod p, resp x mod 4 if p = 2.
3037 : * Simplified form of Zp_sqrtnlift: 1/(p-1) is trivial to compute */
3038 : GEN
3039 10815 : Zp_teichmuller(GEN x, GEN p, long e, GEN pe)
3040 : {
3041 : GEN q, z, p1;
3042 : pari_sp av;
3043 : ulong mask;
3044 10815 : if (absequaliu(p,2)) return (mod4(x) & 2)? subiu(pe,1): gen_1;
3045 10136 : if (e == 1) return icopy(x);
3046 10136 : av = avma;
3047 10136 : p1 = subiu(p, 1);
3048 10135 : mask = quadratic_prec_mask(e);
3049 10135 : q = p; z = remii(x, p);
3050 35504 : while (mask > 1)
3051 : { /* Newton iteration solving z^{1 - p} = 1, z = x (mod p) */
3052 25368 : GEN w, t, qold = q;
3053 25368 : if (mask <= 3) /* last iteration */
3054 10136 : q = pe;
3055 : else
3056 : {
3057 15232 : q = sqri(q);
3058 15231 : if (mask & 1) q = diviiexact(q, p);
3059 : }
3060 25367 : mask >>= 1;
3061 : /* q <= qold^2 */
3062 25367 : if (lgefint(q) == 3)
3063 : {
3064 24381 : ulong Z = uel(z,2), Q = uel(q,2), P1 = uel(p1,2);
3065 24381 : ulong W = (Q-1) / P1; /* -1/(p-1) + O(qold) */
3066 24381 : ulong T = Fl_mul(W, Fl_powu(Z,P1,Q) - 1, Q);
3067 24380 : Z = Fl_mul(Z, 1 + T, Q);
3068 24379 : z = utoi(Z);
3069 : }
3070 : else
3071 : {
3072 986 : w = diviiexact(subiu(qold,1),p1); /* -1/(p-1) + O(qold) */
3073 986 : t = Fp_mul(w, subiu(Fp_pow(z,p1,q), 1), q);
3074 986 : z = Fp_mul(z, addui(1,t), q);
3075 : }
3076 : }
3077 10136 : return gc_INT(av, z);
3078 : }
3079 :
3080 : GEN
3081 1225 : teichmullerinit(long p, long n)
3082 : {
3083 : GEN t, pn, g, v;
3084 : ulong gp, tp;
3085 : long a, m;
3086 :
3087 1225 : if (p == 2) return mkvec(gen_1);
3088 1225 : if (!uisprime(p)) pari_err_PRIME("teichmullerinit",utoipos(p));
3089 :
3090 1225 : m = p >> 1; /* (p-1)/2 */
3091 1225 : tp= gp= pgener_Fl(p); /* order (p-1), gp^m = -1 */
3092 1225 : pn = powuu(p, n);
3093 1225 : v = cgetg(p, t_VEC);
3094 1225 : t = g = Zp_teichmuller(utoipos(gp), utoipos(p), n, pn);
3095 1225 : gel(v, 1) = gen_1;
3096 1225 : gel(v, p-1) = subiu(pn,1);
3097 3031 : for (a = 1; a < m; a++)
3098 : {
3099 1806 : gel(v, tp) = t;
3100 1806 : gel(v, p - tp) = Fp_neg(t, pn); /* g^(m+a) = -g^a */
3101 1806 : if (a < m-1)
3102 : {
3103 1029 : t = Fp_mul(t, g, pn); /* g^(a+1) */
3104 1029 : tp = Fl_mul(tp, gp, p); /* t mod p */
3105 : }
3106 : }
3107 1225 : return v;
3108 : }
3109 :
3110 : /* tab from teichmullerinit or NULL */
3111 : GEN
3112 5803 : teichmuller(GEN x, GEN tab)
3113 : {
3114 : GEN p, q, z;
3115 5803 : long n, tx = typ(x);
3116 :
3117 5803 : if (!tab)
3118 : {
3119 5691 : if (tx == t_VEC && lg(x) == 3)
3120 : {
3121 7 : p = gel(x,1);
3122 7 : q = gel(x,2);
3123 7 : if (typ(p) == t_INT && typ(q) == t_INT)
3124 7 : return teichmullerinit(itos(p), itos(q));
3125 : }
3126 : }
3127 112 : else if (typ(tab) != t_VEC) pari_err_TYPE("teichmuller",tab);
3128 5796 : if (tx!=t_PADIC) pari_err_TYPE("teichmuller",x);
3129 5796 : z = padic_u(x);
3130 5796 : if (!signe(z)) return gcopy(x);
3131 5796 : p = padic_p(x);
3132 5796 : q = padic_pd(x);
3133 5796 : n = precp(x);
3134 5796 : if (tab)
3135 : {
3136 112 : ulong pp = itou_or_0(p);
3137 112 : if (lg(tab) != (long)pp) pari_err_TYPE("teichmuller",tab);
3138 112 : z = gel(tab, umodiu(z, pp));
3139 112 : if (typ(z) != t_INT) pari_err_TYPE("teichmuller",tab);
3140 : }
3141 5796 : retmkpadic(tab? remii(z, q): Zp_teichmuller(z, p, n, q),
3142 : icopy(p), icopy(q), 0, n);
3143 : }
3144 : GEN
3145 5565 : teich(GEN x) { return teichmuller(x, NULL); }
3146 :
3147 : GEN
3148 17997168 : glog(GEN x, long prec)
3149 : {
3150 : pari_sp av;
3151 : GEN y, p1;
3152 : long l;
3153 :
3154 17997168 : switch(typ(x))
3155 : {
3156 10364007 : case t_REAL:
3157 10364007 : if (signe(x) >= 0)
3158 : {
3159 8695953 : if (!signe(x)) pari_err_DOMAIN("log", "argument", "=", gen_0, x);
3160 8695946 : return logr_abs(x);
3161 : }
3162 1668054 : retmkcomplex(logr_abs(x), mppi(realprec(x)));
3163 :
3164 517117 : case t_FRAC:
3165 : {
3166 : GEN a, b;
3167 : long e1, e2;
3168 517117 : av = avma;
3169 517117 : a = gel(x,1);
3170 517117 : b = gel(x,2);
3171 517117 : e1 = expi(subii(a,b)); e2 = expi(b);
3172 517114 : if (e2 > e1) prec += nbits2extraprec(e2 - e1);
3173 517114 : x = fractor(x, prec);
3174 517116 : return gc_upto(av, glog(x, prec));
3175 : }
3176 4731731 : case t_COMPLEX:
3177 4731731 : if (ismpzero(gel(x,2))) return glog(gel(x,1), prec);
3178 4721075 : l = precision(x); if (l > prec) prec = l;
3179 4721078 : if (ismpzero(gel(x,1)))
3180 : {
3181 70314 : GEN a = gel(x,2), b;
3182 70314 : av = avma; b = Pi2n(-1,prec);
3183 70313 : if (gsigne(a) < 0) { setsigne(b, -1); a = gabs(a,prec); }
3184 70313 : a = isint1(a) ? gen_0: glog(a,prec);
3185 70313 : return gc_GEN(av, mkcomplex(a, b));
3186 : }
3187 4650771 : if (prec >= LOGAGMCX_LIMIT) return logagmcx(x, prec);
3188 4351016 : y = cgetg(3,t_COMPLEX);
3189 4351021 : gel(y,2) = garg(x,prec); av = avma;
3190 4351018 : gel(y,1) = gc_upto(av, gmul2n(glog(cxnorm(x),prec),-1)); return y;
3191 :
3192 322 : case t_PADIC: return Qp_log(x);
3193 2383991 : default:
3194 2383991 : av = avma; if (!(y = toser_i(x))) break;
3195 140 : if (!signe(y)) pari_err_DOMAIN("log", "argument", "=", gen_0, x);
3196 140 : if (valser(y)) pari_err_DOMAIN("log", "series valuation", "!=", gen_0, x);
3197 133 : p1 = integser(gdiv(derivser(y), y)); /* log(y)' = y'/y */
3198 133 : if (!gequal1(gel(y,2))) p1 = gadd(p1, glog(gel(y,2),prec));
3199 133 : return gc_upto(av, p1);
3200 : }
3201 2384151 : return trans_eval("log",glog,x,prec);
3202 : }
3203 :
3204 : static GEN
3205 63 : mplog1p(GEN x)
3206 : {
3207 : long ex, a, b, l, L;
3208 63 : if (!signe(x)) return rcopy(x);
3209 63 : ex = expo(x); if (ex >= -3) return glog(addrs(x,1), 0);
3210 42 : a = -ex;
3211 42 : b = realprec(x); L = b+1;
3212 42 : if (b > a*log2(L) && b > LOGAGM_LIMIT)
3213 : {
3214 0 : x = addrs(x,1); l = b + nbits2extraprec(a);
3215 0 : if (realprec(x) < l) x = rtor(x,l);
3216 0 : return logagmr_abs(x);
3217 : }
3218 42 : x = rtor(x, L);
3219 42 : x = logr_aux(divrr(x, addrs(x,2)));
3220 42 : if (realprec(x) > b) fixlg(x, b);
3221 42 : shiftr_inplace(x,1); return x;
3222 : }
3223 :
3224 : static GEN log1p_i(GEN x, long prec);
3225 : static GEN
3226 14 : cxlog1p(GEN x, long prec)
3227 : {
3228 : pari_sp av;
3229 14 : GEN z, a, b = gel(x,2);
3230 : long l;
3231 14 : if (ismpzero(b)) return log1p_i(gel(x,1), prec);
3232 14 : l = precision(x); if (l > prec) prec = l;
3233 14 : if (prec >= LOGAGMCX_LIMIT) return logagmcx(gaddgs(x,1), prec);
3234 14 : a = gel(x,1);
3235 14 : z = cgetg(3,t_COMPLEX); av = avma;
3236 14 : a = gadd(gadd(gmul2n(a,1), gsqr(a)), gsqr(b));
3237 14 : a = log1p_i(a, prec); shiftr_inplace(a,-1);
3238 14 : gel(z,1) = gc_upto(av, a);
3239 14 : gel(z,2) = garg(gaddgs(x,1),prec); return z;
3240 : }
3241 : static GEN
3242 133 : log1p_i(GEN x, long prec)
3243 : {
3244 133 : switch(typ(x))
3245 : {
3246 63 : case t_REAL: return mplog1p(x);
3247 14 : case t_COMPLEX: return cxlog1p(x, prec);
3248 7 : case t_PADIC: return Qp_log(gaddgs(x,1));
3249 49 : default:
3250 : {
3251 : long ey;
3252 : GEN y;
3253 49 : if (!(y = toser_i(x))) break;
3254 21 : ey = valser(y);
3255 21 : if (ey < 0) pari_err_DOMAIN("log1p","valuation", "<", gen_0, x);
3256 21 : if (gequal0(y)) return gcopy(y);
3257 14 : if (ey)
3258 7 : return glog(gaddgs(y,1),prec);
3259 : else
3260 : {
3261 7 : GEN a = gel(y,2), a1 = gaddgs(a,1);
3262 7 : y = gdiv(y, a1); gel(y,2) = gen_1;
3263 7 : return gadd(glog1p(a,prec), glog(y, prec));
3264 : }
3265 : }
3266 : }
3267 28 : return trans_eval("log1p",glog1p,x,prec);
3268 : }
3269 : GEN
3270 119 : glog1p(GEN x, long prec)
3271 : {
3272 119 : pari_sp av = avma;
3273 119 : return gc_upto(av, log1p_i(x, prec));
3274 : }
3275 : /********************************************************************/
3276 : /** **/
3277 : /** SINE, COSINE **/
3278 : /** **/
3279 : /********************************************************************/
3280 :
3281 : /* Reduce x0 mod Pi/2 to x in [-Pi/4, Pi/4]. Return cos(x)-1 */
3282 : static GEN
3283 17424766 : mpcosm1(GEN x, long *ptmod8)
3284 : {
3285 17424766 : long a = expo(x), l = realprec(x), b, L, i, n, m, B;
3286 : GEN y, u, x2;
3287 : double d;
3288 :
3289 17424766 : n = 0;
3290 17424766 : if (a >= 0)
3291 : {
3292 : long p;
3293 : GEN q;
3294 10194393 : if (a > 30)
3295 : {
3296 684652 : GEN z, P = Pi2n(-2, nbits2prec(a + 32));
3297 684652 : z = addrr(x,P); /* = x + Pi/4 */
3298 684652 : if (expo(z) >= bit_prec(z) + 3) pari_err_PREC("mpcosm1");
3299 684652 : shiftr_inplace(P, 1);
3300 684652 : q = floorr(divrr(z, P)); /* round ( x / (Pi/2) ) */
3301 684652 : p = l+EXTRAPREC64; x = rtor(x,p);
3302 : } else {
3303 9509741 : q = stoi((long)floor(rtodbl(x) / (M_PI/2) + 0.5));
3304 9509757 : p = l;
3305 : }
3306 10198612 : if (signe(q))
3307 : {
3308 10194409 : GEN y = subrr(x, mulir(q, Pi2n(-1,p))); /* x mod Pi/2 */
3309 10194311 : long b = expo(y);
3310 10194311 : if (a - b < 7) x = y;
3311 : else
3312 : {
3313 6114138 : p += nbits2extraprec(a-b); x = rtor(x, p);
3314 6114144 : x = subrr(x, mulir(q, Pi2n(-1,p)));
3315 : }
3316 10194291 : a = b;
3317 10194291 : if (!signe(x) && a >= 0) pari_err_PREC("mpcosm1");
3318 10194291 : n = Mod4(q);
3319 : }
3320 : }
3321 : /* a < 0 */
3322 17428835 : b = signe(x); *ptmod8 = (b < 0)? 4 + n: n;
3323 17428835 : if (!b) return real_0_bit(expo(x)*2 - 1);
3324 :
3325 17428835 : b = prec2nbits(l);
3326 17423916 : if (b + 2*a <= 0) {
3327 1382951 : y = sqrr(x); shiftr_inplace(y, -1); setsigne(y, -1);
3328 1382952 : return y;
3329 : }
3330 :
3331 16040965 : y = cgetr(l);
3332 16042811 : B = b/6 + BITS_IN_LONG/2 + (BITS_IN_LONG*BITS_IN_LONG/2)/ b;
3333 16042811 : d = a/2.; m = (long)(d + sqrt(d*d + B)); /* >= 0 ,*/
3334 16042811 : if (m < (-a) * 0.1) m = 0; /* not worth it */
3335 16042811 : L = l + nbits2extraprec(m);
3336 :
3337 16042503 : b += m;
3338 16042503 : d = 2.0 * (m-dbllog2r(x)-1/M_LN2); /* ~ 2( - log_2 Y - 1/log(2) ) */
3339 16042340 : n = (long)(b / d);
3340 16042340 : if (n > 1)
3341 15984363 : n = (long)(b / (d + log2((double)n+1))); /* log~constant in small ranges */
3342 34561022 : while (n*(d+log2((double)n+1)) < b) n++; /* expect few corrections */
3343 :
3344 : /* Multiplication is quadratic in this range (l is small, otherwise we
3345 : * use logAGM + Newton). Set Y = 2^(-e-a) x, compute truncated series
3346 : * sum Y^2k/(2k)!: this costs roughly
3347 : * m b^2 + sum_{k <= n} (2k e + BITS_IN_LONG)^2
3348 : * ~ (b/2e) b^2 / 3 + m b^2
3349 : * bit operations with n ~ b/2e, |x| < 2^(1+a), |Y| < 2^(1-e) , m = e+a and
3350 : * b bits of accuracy needed, so
3351 : * B := (b / 6 + BITS_IN_LONG/2 + BITS_IN_LONG^2 / 2b) ~ m(m-a)
3352 : * we want b ~ 6 m (m-a) or m~b+a hence
3353 : * m = min( a/2 + sqrt(a^2/4 + b/6), b/2 + a )
3354 : * NB: e ~ (b/6)^(1/2) or b/2.
3355 : *
3356 : * Truncate the sum at k = n (>= 1), the remainder is
3357 : * < sum_{k >= n+1} Y^2k / 2k! < Y^(2n+2) / (2n+2)!(1-Y^2) < Y^(2n+2)/(2n+1)!
3358 : * We want ... <= Y^2 2^-b, hence -2n log_2 |Y| + log_2 (2n+1)! >= b
3359 : * log n! ~ (n + 1/2) log(n+1) - (n+1) + log(2Pi)/2,
3360 : * error bounded by 1/6(n+1) <= 1/12. Finally, we want
3361 : * 2n (-1/log(2) - log_2 |Y| + log_2(2n+2)) >= b */
3362 16042340 : x = rtor(x, L); shiftr_inplace(x, -m); setsigne(x, 1);
3363 16043578 : x2 = sqrr(x);
3364 16045329 : if (n == 1) { u = x2; shiftr_inplace(u, -1); setsigne(u, -1); } /*-Y^2/2*/
3365 : else
3366 : {
3367 16045329 : GEN un = real_1(L);
3368 : pari_sp av;
3369 16041793 : long s = 0, l1 = nbits2prec((long)(d + n + 16));
3370 :
3371 16041675 : u = cgetr(L); av = avma;
3372 240286624 : for (i = n; i >= 2; i--)
3373 : {
3374 : GEN t;
3375 224246372 : setprec(x2,l1); t = divrunextu(x2, 2*i-1);
3376 224821781 : l1 += nbits2extraprec(dvmdsBIL(s - expo(t), &s)<<TWOPOTBITS_IN_LONG);
3377 224595673 : if (l1 > L) l1 = L;
3378 224595673 : if (i != n) t = mulrr(t,u);
3379 224964288 : setprec(un,l1); t = addrr_sign(un,1, t,-signe(t));
3380 224424805 : setprec(u,l1); affrr(t,u); set_avma(av);
3381 : }
3382 16040252 : shiftr_inplace(u, -1); togglesign(u); /* u := -u/2 */
3383 16042085 : setprec(x2,L); u = mulrr(x2,u);
3384 : }
3385 : /* Now u = sum {1<= i <=n} (-1)^i x^(2i) / (2i)! ~ cos(x) - 1 */
3386 143127685 : for (i = 1; i <= m; i++)
3387 : { /* u = cos(x)-1 <- cos(2x)-1 = 2cos(x)^2 - 2 = 4u + 2u^2*/
3388 127119198 : GEN q = sqrr(u);
3389 127323794 : shiftr_inplace(u, 1); u = addrr(u, q);
3390 127175724 : shiftr_inplace(u, 1);
3391 127078602 : if ((i & 31) == 0) u = gc_uptoleaf((pari_sp)y, u);
3392 : }
3393 16008487 : affrr_fixlg(u, y); return y;
3394 : }
3395 :
3396 : /* sqrt (|1 - (1+x)^2|) = sqrt(|x*(x+2)|). Sends cos(x)-1 to |sin(x)| */
3397 : static GEN
3398 15722973 : mpaut(GEN x)
3399 : {
3400 15722973 : GEN t = mulrr(x, addsr(2,x)); /* != 0 */
3401 15728708 : if (!signe(t)) return real_0_bit(expo(t) >> 1);
3402 15728708 : return sqrtr_abs(t);
3403 : }
3404 :
3405 : /********************************************************************/
3406 : /** COSINE **/
3407 : /********************************************************************/
3408 :
3409 : GEN
3410 2745068 : mpcos(GEN x)
3411 : {
3412 : long mod8;
3413 : pari_sp av;
3414 : GEN y, z;
3415 :
3416 2745068 : if (!signe(x)) {
3417 75 : long l = nbits2prec(-expo(x));
3418 75 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
3419 75 : return real_1(l);
3420 : }
3421 2744993 : av = avma; z = mpcosm1(x,&mod8);
3422 2744986 : switch(mod8)
3423 : {
3424 760059 : case 0: case 4: y = addsr(1,z); break;
3425 688607 : case 1: case 7: y = mpaut(z); togglesign(y); break;
3426 682287 : case 2: case 6: y = subsr(-1,z); break;
3427 614033 : default: y = mpaut(z); break; /* case 3: case 5: */
3428 : }
3429 2745017 : return gc_uptoleaf(av, y);
3430 : }
3431 :
3432 : /* convert INT or FRAC to REAL, which is later reduced mod 2Pi : avoid
3433 : * cancellation */
3434 : static GEN
3435 13252 : tofp_safe(GEN x, long prec)
3436 : {
3437 13252 : return (typ(x) == t_INT || gexpo(x) > 0)? gadd(x, real_0(prec))
3438 26489 : : fractor(x, prec);
3439 : }
3440 :
3441 : GEN
3442 154859 : gcos(GEN x, long prec)
3443 : {
3444 : pari_sp av;
3445 : GEN a, b, u, v, y, u1, v1;
3446 : long i;
3447 :
3448 154859 : switch(typ(x))
3449 : {
3450 153577 : case t_REAL: return mpcos(x);
3451 42 : case t_COMPLEX:
3452 42 : a = gel(x,1);
3453 42 : b = gel(x,2);
3454 42 : if (isintzero(a)) return gcosh(b, prec);
3455 28 : i = precision(x); if (i) prec = i;
3456 28 : y = cgetc(prec); av = avma;
3457 28 : if (typ(b) != t_REAL) b = gtofp(b, prec);
3458 28 : mpsinhcosh(b, &u1, &v1); u1 = mpneg(u1);
3459 28 : if (typ(a) != t_REAL) a = gtofp(a, prec);
3460 28 : mpsincos(a, &u, &v);
3461 28 : affrr_fixlg(gmul(v1,v), gel(y,1));
3462 28 : affrr_fixlg(gmul(u1,u), gel(y,2)); return gc_const(av,y);
3463 :
3464 1156 : case t_INT: case t_FRAC:
3465 1156 : y = cgetr(prec); av = avma;
3466 1156 : affrr_fixlg(mpcos(tofp_safe(x,prec)), y); return gc_const(av,y);
3467 :
3468 49 : case t_PADIC: y = cos_p(x);
3469 49 : if (!y) pari_err_DOMAIN("gcos(t_PADIC)","argument","",gen_0,x);
3470 42 : return y;
3471 :
3472 35 : default:
3473 35 : av = avma; if (!(y = toser_i(x))) break;
3474 28 : if (gequal0(y)) return gc_upto(av, gaddsg(1,y));
3475 28 : if (valser(y) < 0)
3476 7 : pari_err_DOMAIN("cos","valuation", "<", gen_0, x);
3477 21 : gsincos(y,&u,&v,prec);
3478 21 : return gc_GEN(av,v);
3479 : }
3480 7 : return trans_eval("cos",gcos,x,prec);
3481 : }
3482 : /********************************************************************/
3483 : /** SINE **/
3484 : /********************************************************************/
3485 :
3486 : GEN
3487 834962 : mpsin(GEN x)
3488 : {
3489 : long mod8;
3490 : pari_sp av;
3491 : GEN y, z;
3492 :
3493 834962 : if (!signe(x)) return real_0_bit(expo(x));
3494 834753 : av = avma; z = mpcosm1(x,&mod8);
3495 834749 : switch(mod8)
3496 : {
3497 311654 : case 0: case 6: y = mpaut(z); break;
3498 131344 : case 1: case 5: y = addsr(1,z); break;
3499 264277 : case 2: case 4: y = mpaut(z); togglesign(y); break;
3500 127474 : default: y = subsr(-1,z); break; /* case 3: case 7: */
3501 : }
3502 834773 : return gc_uptoleaf(av, y);
3503 : }
3504 :
3505 : GEN
3506 1252137 : gsin(GEN x, long prec)
3507 : {
3508 : pari_sp av;
3509 : GEN a, b, u, v, y, v1, u1;
3510 : long i;
3511 :
3512 1252137 : switch(typ(x))
3513 : {
3514 829780 : case t_REAL: return mpsin(x);
3515 416959 : case t_COMPLEX:
3516 416959 : a = gel(x,1);
3517 416959 : b = gel(x,2);
3518 416959 : if (isintzero(a)) retmkcomplex(gen_0,gsinh(b,prec));
3519 410974 : i = precision(x); if (i) prec = i;
3520 410974 : y = cgetc(prec); av = avma;
3521 410974 : if (typ(b) != t_REAL) b = gtofp(b, prec);
3522 410974 : mpsinhcosh(b, &u1, &v1);
3523 410974 : if (typ(a) != t_REAL) a = gtofp(a, prec);
3524 410974 : mpsincos(a, &u, &v);
3525 410974 : affrr_fixlg(gmul(v1,u), gel(y,1));
3526 410974 : affrr_fixlg(gmul(u1,v), gel(y,2)); return gc_const(av,y);
3527 :
3528 5118 : case t_INT: case t_FRAC:
3529 5118 : y = cgetr(prec); av = avma;
3530 5118 : affrr_fixlg(mpsin(tofp_safe(x,prec)), y); return gc_const(av,y);
3531 :
3532 49 : case t_PADIC: y = sin_p(x);
3533 49 : if (!y) pari_err_DOMAIN("gsin(t_PADIC)","argument","",gen_0,x);
3534 42 : return y;
3535 :
3536 231 : default:
3537 231 : av = avma; if (!(y = toser_i(x))) break;
3538 224 : if (gequal0(y)) return gc_GEN(av, y);
3539 224 : if (valser(y) < 0)
3540 7 : pari_err_DOMAIN("sin","valuation", "<", gen_0, x);
3541 217 : gsincos(y,&u,&v,prec);
3542 217 : return gc_GEN(av,u);
3543 : }
3544 7 : return trans_eval("sin",gsin,x,prec);
3545 : }
3546 : /********************************************************************/
3547 : /** SINE, COSINE together **/
3548 : /********************************************************************/
3549 :
3550 : void
3551 13846028 : mpsincos(GEN x, GEN *s, GEN *c)
3552 : {
3553 : long mod8;
3554 : pari_sp av, tetpil;
3555 : GEN z;
3556 :
3557 13846028 : if (!signe(x))
3558 : {
3559 20952 : long e = expo(x);
3560 20952 : *s = real_0_bit(e);
3561 20952 : *c = e >= 0? real_0_bit(e): real_1_bit(-e);
3562 20952 : return;
3563 : }
3564 :
3565 13825076 : av = avma; z = mpcosm1(x, &mod8); tetpil = avma;
3566 13828130 : switch(mod8)
3567 : {
3568 4457414 : case 0: *c = addsr( 1,z); *s = mpaut(z); break;
3569 688171 : case 1: *s = addsr( 1,z); *c = mpaut(z); togglesign(*c); break;
3570 1041342 : case 2: *c = subsr(-1,z); *s = mpaut(z); togglesign(*s); break;
3571 644902 : case 3: *s = subsr(-1,z); *c = mpaut(z); break;
3572 3750713 : case 4: *c = addsr( 1,z); *s = mpaut(z); togglesign(*s); break;
3573 677032 : case 5: *s = addsr( 1,z); *c = mpaut(z); break;
3574 1893522 : case 6: *c = subsr(-1,z); *s = mpaut(z); break;
3575 677862 : case 7: *s = subsr(-1,z); *c = mpaut(z); togglesign(*c); break;
3576 : }
3577 13829511 : (void)gc_all_unsafe(av,tetpil,2,s,c);
3578 : }
3579 :
3580 : /* SINE and COSINE - 1 */
3581 : void
3582 20130 : mpsincosm1(GEN x, GEN *s, GEN *c)
3583 : {
3584 : long mod8;
3585 : pari_sp av, tetpil;
3586 : GEN z;
3587 :
3588 20130 : if (!signe(x))
3589 : {
3590 0 : long e = expo(x);
3591 0 : *s = real_0_bit(e);
3592 0 : *c = real_0_bit(2*e-1);
3593 0 : return;
3594 : }
3595 20130 : av = avma; z = mpcosm1(x,&mod8); tetpil = avma;
3596 20130 : switch(mod8)
3597 : {
3598 6912 : case 0: *c = rcopy(z); *s = mpaut(z); break;
3599 1844 : case 1: *s = addsr(1,z); *c = addrs(mpaut(z),1); togglesign(*c); break;
3600 1591 : case 2: *c = subsr(-2,z); *s = mpaut(z); togglesign(*s); break;
3601 1881 : case 3: *s = subsr(-1,z); *c = subrs(mpaut(z),1); break;
3602 2295 : case 4: *c = rcopy(z); *s = mpaut(z); togglesign(*s); break;
3603 2062 : case 5: *s = addsr( 1,z); *c = subrs(mpaut(z),1); break;
3604 1650 : case 6: *c = subsr(-2,z); *s = mpaut(z); break;
3605 1895 : case 7: *s = subsr(-1,z); *c = subsr(-1,mpaut(z)); break;
3606 : }
3607 20130 : (void)gc_all_unsafe(av,tetpil,2,s,c);
3608 : }
3609 :
3610 : /* return exp(ix), x a t_REAL */
3611 : GEN
3612 899259 : expIr(GEN x)
3613 : {
3614 899259 : pari_sp av = avma;
3615 899259 : GEN v = cgetg(3,t_COMPLEX);
3616 899259 : mpsincos(x, (GEN*)(v+2), (GEN*)(v+1));
3617 899264 : if (!signe(gel(v,2))) return gc_GEN(av, gel(v,1));
3618 879638 : return v;
3619 : }
3620 :
3621 : /* return exp(ix)-1, x a t_REAL */
3622 : static GEN
3623 20130 : expm1_Ir(GEN x)
3624 : {
3625 20130 : pari_sp av = avma;
3626 20130 : GEN v = cgetg(3,t_COMPLEX);
3627 20130 : mpsincosm1(x, (GEN*)(v+2), (GEN*)(v+1));
3628 20130 : if (!signe(gel(v,2))) return gc_GEN(av, gel(v,1));
3629 20130 : return v;
3630 : }
3631 :
3632 : /* return exp(z)-1, z complex */
3633 : GEN
3634 20186 : cxexpm1(GEN z, long prec)
3635 : {
3636 20186 : pari_sp av = avma;
3637 20186 : GEN X, Y, x = real_i(z), y = imag_i(z);
3638 20186 : long l = precision(z);
3639 20186 : if (l) prec = l;
3640 20186 : if (typ(x) != t_REAL) x = gtofp(x, prec);
3641 20186 : if (typ(y) != t_REAL) y = gtofp(y, prec);
3642 20186 : if (gequal0(y)) return mpexpm1(x);
3643 20130 : if (gequal0(x)) return expm1_Ir(y);
3644 19997 : X = mpexpm1(x); /* t_REAL */
3645 19997 : Y = expm1_Ir(y);
3646 : /* exp(x+iy) - 1 = (exp(x)-1)(exp(iy)-1) + exp(x)-1 + exp(iy)-1 */
3647 19997 : return gc_upto(av, gadd(gadd(X,Y), gmul(X,Y)));
3648 : }
3649 :
3650 : void
3651 4697537 : gsincos(GEN x, GEN *s, GEN *c, long prec)
3652 : {
3653 : long i, j, ex, ex2, lx, ly, mi;
3654 : pari_sp av, tetpil;
3655 : GEN y, r, u, v, u1, v1, p1, p2, p3, p4, ps, pc;
3656 :
3657 4697537 : switch(typ(x))
3658 : {
3659 6948 : case t_INT: case t_FRAC:
3660 6948 : *s = cgetr(prec);
3661 6945 : *c = cgetr(prec); av = avma;
3662 6943 : mpsincos(tofp_safe(x, prec), &ps, &pc);
3663 6952 : affrr_fixlg(ps,*s);
3664 4697742 : affrr_fixlg(pc,*c); set_avma(av); return;
3665 :
3666 4685960 : case t_REAL:
3667 4685960 : mpsincos(x,s,c); return;
3668 :
3669 4123 : case t_COMPLEX:
3670 4123 : i = precision(x); if (i) prec = i;
3671 4123 : ps = cgetc(prec); *s = ps;
3672 4123 : pc = cgetc(prec); *c = pc; av = avma;
3673 4123 : r = gexp(gel(x,2),prec);
3674 4123 : v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
3675 4123 : u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
3676 4123 : gsincos(gel(x,1), &u,&v, prec);
3677 4123 : affrr_fixlg(mulrr(v1,u), gel(ps,1));
3678 4123 : affrr_fixlg(mulrr(u1,v), gel(ps,2));
3679 4123 : affrr_fixlg(mulrr(v1,v), gel(pc,1));
3680 4123 : affrr_fixlg(mulrr(u1,u), gel(pc,2)); togglesign(gel(pc,2));
3681 4123 : set_avma(av); return;
3682 :
3683 0 : case t_QUAD:
3684 0 : av = avma; gsincos(quadtofp(x, prec), s, c, prec);
3685 0 : (void)gc_all(av, 2, s, c); return;
3686 :
3687 506 : default:
3688 506 : av = avma; if (!(y = toser_i(x))) break;
3689 518 : if (gequal0(y)) { *s = gc_GEN(av,y); *c = gaddsg(1,*s); return; }
3690 :
3691 518 : ex = valser(y); lx = lg(y); ex2 = 2*ex+2;
3692 518 : if (ex < 0) pari_err_DOMAIN("gsincos","valuation", "<", gen_0, x);
3693 518 : if (ex2 > lx)
3694 : {
3695 98 : *s = x == y? gcopy(y): gc_GEN(av, y); av = avma;
3696 98 : *c = gc_upto(av, gsubsg(1, gdivgu(gsqr(y),2)));
3697 98 : return;
3698 : }
3699 420 : if (!ex)
3700 : {
3701 105 : gsincos(serchop0(y),&u,&v,prec);
3702 105 : gsincos(gel(y,2),&u1,&v1,prec);
3703 105 : p1 = gmul(v1,v);
3704 105 : p2 = gmul(u1,u);
3705 105 : p3 = gmul(v1,u);
3706 105 : p4 = gmul(u1,v); tetpil = avma;
3707 105 : *c = gsub(p1,p2);
3708 105 : *s = gadd(p3,p4);
3709 105 : gc_all_unsafe(av,tetpil,2,s,c); return;
3710 : }
3711 :
3712 315 : ly = lx+2*ex;
3713 3066 : mi = lx-1; while (mi>=3 && isrationalzero(gel(y,mi))) mi--;
3714 315 : mi += ex-2;
3715 315 : pc = cgetg(ly,t_SER); *c = pc;
3716 315 : ps = cgetg(lx,t_SER); *s = ps;
3717 315 : pc[1] = evalsigne(1) | _evalvalser(0) | evalvarn(varn(y));
3718 315 : gel(pc,2) = gen_1; ps[1] = y[1];
3719 637 : for (i=2; i<ex+2; i++) gel(ps,i) = gcopy(gel(y,i));
3720 644 : for (i=3; i< ex2; i++) gel(pc,i) = gen_0;
3721 3577 : for (i=ex2; i<ly; i++)
3722 : {
3723 3262 : long ii = i-ex;
3724 3262 : av = avma; p1 = gen_0;
3725 7476 : for (j=ex; j<=minss(ii-2,mi); j++)
3726 4214 : p1 = gadd(p1, gmulgu(gmul(gel(y,j-ex+2),gel(ps,ii-j)),j));
3727 3262 : gel(pc,i) = gc_upto(av, gdivgs(p1,2-i));
3728 3262 : if (ii < lx)
3729 : {
3730 2940 : av = avma; p1 = gen_0;
3731 6202 : for (j=ex; j<=minss(i-ex2,mi); j++)
3732 3262 : p1 = gadd(p1,gmulgu(gmul(gel(y,j-ex+2),gel(pc,i-j)),j));
3733 2940 : p1 = gdivgu(p1,i-2);
3734 2940 : gel(ps,ii) = gc_upto(av, gadd(p1,gel(y,ii)));
3735 : }
3736 : }
3737 315 : return;
3738 : }
3739 0 : pari_err_TYPE("gsincos",x);
3740 : }
3741 :
3742 : /********************************************************************/
3743 : /** **/
3744 : /** SINC **/
3745 : /** **/
3746 : /********************************************************************/
3747 : static GEN
3748 2319450 : mpsinc(GEN x)
3749 : {
3750 2319450 : pari_sp av = avma;
3751 : GEN s, c;
3752 :
3753 2319450 : if (!signe(x)) {
3754 0 : long l = nbits2prec(-expo(x));
3755 0 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
3756 0 : return real_1(l);
3757 : }
3758 2319450 : mpsincos(x,&s,&c);
3759 2319450 : return gc_uptoleaf(av, divrr(s,x));
3760 : }
3761 :
3762 : GEN
3763 2319562 : gsinc(GEN x, long prec)
3764 : {
3765 : pari_sp av;
3766 : GEN r, u, v, y, u1, v1;
3767 : long i;
3768 :
3769 2319562 : switch(typ(x))
3770 : {
3771 2319429 : case t_REAL: return mpsinc(x);
3772 49 : case t_COMPLEX:
3773 49 : if (isintzero(gel(x,1)))
3774 : {
3775 28 : av = avma; x = gel(x,2);
3776 28 : if (gequal0(x)) return gcosh(x,prec);
3777 14 : return gc_uptoleaf(av,gdiv(gsinh(x,prec),x));
3778 : }
3779 21 : i = precision(x); if (i) prec = i;
3780 21 : y = cgetc(prec); av = avma;
3781 21 : r = gexp(gel(x,2),prec);
3782 21 : v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
3783 21 : u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
3784 21 : gsincos(gel(x,1),&u,&v,prec);
3785 21 : affc_fixlg(gdiv(mkcomplex(gmul(v1,u), gmul(u1,v)), x), y);
3786 21 : return gc_const(av,y);
3787 :
3788 14 : case t_INT:
3789 14 : if (!signe(x)) return real_1(prec); /*fall through*/
3790 : case t_FRAC:
3791 21 : y = cgetr(prec); av = avma;
3792 21 : affrr_fixlg(mpsinc(tofp_safe(x,prec)), y); return gc_const(av,y);
3793 :
3794 21 : case t_PADIC:
3795 21 : if (gequal0(x)) return cvtop(gen_1, padic_p(x), valp(x));
3796 14 : av = avma; y = sin_p(x);
3797 14 : if (!y) pari_err_DOMAIN("gsinc(t_PADIC)","argument","",gen_0,x);
3798 7 : return gc_upto(av, gdiv(y,x));
3799 :
3800 35 : default:
3801 : {
3802 : long ex;
3803 35 : av = avma; if (!(y = toser_i(x))) break;
3804 35 : if (gequal0(y)) return gc_upto(av, gaddsg(1,y));
3805 35 : ex = valser(y);
3806 35 : if (ex < 0) pari_err_DOMAIN("sinc","valuation", "<", gen_0, x);
3807 28 : if (ex)
3808 : {
3809 28 : gsincos(y,&u,&v,prec);
3810 28 : y = gc_upto(av, gdiv(u,y));
3811 28 : if (lg(y) > 2) gel(y,2) = gen_1;
3812 28 : return y;
3813 : }
3814 : else
3815 : {
3816 0 : GEN z0, y0 = gel(y,2), y1 = serchop0(y), y10 = y1;
3817 0 : if (!gequal1(y0)) y10 = gdiv(y10, y0);
3818 0 : gsincos(y1,&u,&v,prec);
3819 0 : z0 = gdiv(gcos(y0,prec), y0);
3820 0 : y = gaddsg(1, y10);
3821 0 : u = gadd(gmul(gsinc(y0, prec),v), gmul(z0, u));
3822 0 : return gc_upto(av,gdiv(u,y));
3823 : }
3824 : }
3825 : }
3826 0 : return trans_eval("sinc",gsinc,x,prec);
3827 : }
3828 :
3829 : /********************************************************************/
3830 : /** **/
3831 : /** TANGENT and COTANGENT **/
3832 : /** **/
3833 : /********************************************************************/
3834 : static GEN
3835 133 : mptan(GEN x)
3836 : {
3837 133 : pari_sp av = avma;
3838 : GEN s, c;
3839 :
3840 133 : mpsincos(x,&s,&c);
3841 133 : if (!signe(c))
3842 0 : pari_err_DOMAIN("tan", "argument", "=", strtoGENstr("Pi/2 + kPi"),x);
3843 133 : return gc_uptoleaf(av, divrr(s,c));
3844 : }
3845 :
3846 : /* If exp(-|im(x)|) << 1, avoid overflow in sincos(x) */
3847 : static int
3848 4018 : tan_huge_im(GEN ix, long prec)
3849 : {
3850 4018 : long b, p = precision(ix);
3851 4018 : if (!p) p = prec;
3852 4018 : b = prec2nbits(p);
3853 4018 : return (gexpo(ix) > b || fabs(gtodouble(ix)) > (M_LN2 / 2) * b);
3854 : }
3855 : /* \pm I */
3856 : static GEN
3857 35 : real_I(long s, long prec)
3858 : {
3859 35 : GEN z = cgetg(3, t_COMPLEX);
3860 35 : gel(z,1) = real_0(prec);
3861 35 : gel(z,2) = s > 0? real_1(prec): real_m1(prec); return z;
3862 : }
3863 :
3864 : GEN
3865 224 : gtan(GEN x, long prec)
3866 : {
3867 : pari_sp av;
3868 : GEN y, s, c;
3869 :
3870 224 : switch(typ(x))
3871 : {
3872 126 : case t_REAL: return mptan(x);
3873 :
3874 42 : case t_COMPLEX: {
3875 42 : if (isintzero(gel(x,1))) retmkcomplex(gen_0,gtanh(gel(x,2),prec));
3876 28 : if (tan_huge_im(gel(x,2), prec)) return real_I(gsigne(gel(x,2)), prec);
3877 14 : av = avma; y = mulcxmI(gtanh(mulcxI(x), prec)); /* tan x = -I th(I x) */
3878 14 : gel(y,1) = gcopy(gel(y,1)); return gc_upto(av, y);
3879 : }
3880 7 : case t_INT: case t_FRAC:
3881 7 : y = cgetr(prec); av = avma;
3882 7 : affrr_fixlg(mptan(tofp_safe(x,prec)), y); return gc_const(av,y);
3883 :
3884 14 : case t_PADIC:
3885 14 : av = avma;
3886 14 : return gc_upto(av, gdiv(gsin(x,prec), gcos(x,prec)));
3887 :
3888 35 : default:
3889 35 : av = avma; if (!(y = toser_i(x))) break;
3890 28 : if (gequal0(y)) return gc_GEN(av, y);
3891 28 : if (valser(y) < 0)
3892 7 : pari_err_DOMAIN("tan","valuation", "<", gen_0, x);
3893 21 : gsincos(y,&s,&c,prec);
3894 21 : return gc_upto(av, gdiv(s,c));
3895 : }
3896 7 : return trans_eval("tan",gtan,x,prec);
3897 : }
3898 :
3899 : static GEN
3900 70 : mpcotan(GEN x)
3901 : {
3902 70 : pari_sp av = avma;
3903 : GEN s,c;
3904 70 : mpsincos(x,&s,&c); return gc_uptoleaf(av, divrr(c,s));
3905 : }
3906 :
3907 : GEN
3908 4214 : gcotan(GEN x, long prec)
3909 : {
3910 : pari_sp av;
3911 : GEN y, s, c;
3912 :
3913 4214 : switch(typ(x))
3914 : {
3915 63 : case t_REAL:
3916 63 : return mpcotan(x);
3917 :
3918 4011 : case t_COMPLEX:
3919 4011 : if (isintzero(gel(x,1))) {
3920 21 : GEN z = cgetg(3, t_COMPLEX);
3921 21 : gel(z,1) = gen_0; av = avma;
3922 21 : gel(z,2) = gc_upto(av, gneg(ginv(gtanh(gel(x,2),prec))));
3923 21 : return z;
3924 : }
3925 3990 : if (tan_huge_im(gel(x,2), prec)) return real_I(-gsigne(gel(x,2)), prec);
3926 3969 : av = avma; gsincos(x,&s,&c,prec);
3927 3969 : return gc_upto(av, gdiv(c,s));
3928 :
3929 7 : case t_INT: case t_FRAC:
3930 7 : y = cgetr(prec); av = avma;
3931 7 : affrr_fixlg(mpcotan(tofp_safe(x,prec)), y); return gc_const(av,y);
3932 :
3933 14 : case t_PADIC:
3934 14 : av = avma;
3935 14 : return gc_upto(av, gdiv(gcos(x,prec), gsin(x,prec)));
3936 :
3937 119 : default:
3938 119 : av = avma; if (!(y = toser_i(x))) break;
3939 112 : if (gequal0(y)) pari_err_DOMAIN("cotan", "argument", "=", gen_0, y);
3940 112 : if (valser(y) < 0) pari_err_DOMAIN("cotan","valuation", "<", gen_0, x);
3941 105 : gsincos(y,&s,&c,prec);
3942 105 : return gc_upto(av, gdiv(c,s));
3943 : }
3944 7 : return trans_eval("cotan",gcotan,x,prec);
3945 : }
|
