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 227702 : pari_init_floats(void)
34 : {
35 227702 : gcatalan = geuler = gpi = zetazone = bernzone = glog2 = eulerzone = NULL;
36 227702 : }
37 :
38 : void
39 228823 : pari_close_floats(void)
40 : {
41 228823 : guncloneNULL(gcatalan);
42 228177 : guncloneNULL(geuler);
43 227356 : guncloneNULL(gpi);
44 227176 : guncloneNULL(glog2);
45 226940 : guncloneNULL(zetazone);
46 226778 : guncloneNULL_deep(bernzone);
47 226548 : guncloneNULL_deep(eulerzone);
48 226369 : }
49 :
50 : /********************************************************************/
51 : /** GENERIC BINARY SPLITTING **/
52 : /** (Haible, Papanikolaou) **/
53 : /********************************************************************/
54 : void
55 265416 : abpq_init(struct abpq *A, long n)
56 : {
57 265416 : A->a = (GEN*)new_chunk(n+1);
58 265590 : A->b = (GEN*)new_chunk(n+1);
59 265782 : A->p = (GEN*)new_chunk(n+1);
60 266065 : A->q = (GEN*)new_chunk(n+1);
61 266064 : }
62 : static GEN
63 19441431 : mulii3(GEN a, GEN b, GEN c) { return mulii(mulii(a,b),c); }
64 :
65 : /* T_{n1,n1+1} */
66 : static GEN
67 4193937 : T2(struct abpq *A, long n1)
68 : {
69 4193937 : GEN u = mulii3(A->a[n1], A->b[n1+1], A->q[n1+1]);
70 4193788 : GEN v = mulii3(A->b[n1], A->a[n1+1], A->p[n1+1]);
71 4194258 : 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 8202049 : 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 8202049 : switch(n2 - n1)
83 : {
84 : GEN b, q;
85 56 : case 1:
86 56 : r->P = A->p[n1];
87 56 : r->Q = A->q[n1];
88 56 : r->B = A->b[n1];
89 56 : r->T = mulii(A->a[n1], A->p[n1]);
90 4214781 : return;
91 2370666 : case 2:
92 2370666 : r->P = mulii(A->p[n1], A->p[n1+1]);
93 2354011 : r->Q = mulii(A->q[n1], A->q[n1+1]);
94 2352903 : r->B = mulii(A->b[n1], A->b[n1+1]);
95 2353054 : av = avma;
96 2353054 : r->T = gerepileuptoint(av, T2(A, n1));
97 2361618 : return;
98 :
99 1871196 : case 3:
100 1871196 : q = mulii(A->q[n1+1], A->q[n1+2]);
101 1865297 : b = mulii(A->b[n1+1], A->b[n1+2]);
102 1864824 : r->P = mulii3(A->p[n1], A->p[n1+1], A->p[n1+2]);
103 1864797 : r->Q = mulii(A->q[n1], q);
104 1864907 : r->B = mulii(A->b[n1], b);
105 1865099 : av = avma;
106 1865099 : u1 = mulii3(b, q, A->a[n1]);
107 1864570 : u2 = mulii(A->b[n1], T2(A, n1+1));
108 1864847 : r->T = gerepileuptoint(av, mulii(A->p[n1], addii(u1, u2)));
109 1853107 : return;
110 : }
111 :
112 3960131 : av = avma;
113 3960131 : n = (n1 + n2) >> 1;
114 3960131 : abpq_sum(&L, n1, n, A);
115 3966387 : abpq_sum(&R, n, n2, A);
116 :
117 3967668 : r->P = mulii(L.P, R.P);
118 3941467 : r->Q = mulii(L.Q, R.Q);
119 3943612 : r->B = mulii(L.B, R.B);
120 3941281 : u1 = mulii3(R.B, R.Q, L.T);
121 3941146 : u2 = mulii3(L.B, L.P, R.T);
122 3938560 : r->T = addii(u1,u2);
123 3941073 : set_avma(av);
124 3943339 : r->P = icopy(r->P);
125 3965573 : r->Q = icopy(r->Q);
126 3970285 : r->B = icopy(r->B);
127 3970066 : 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 76833 : swap_clone(GEN *old, GEN c)
138 76833 : { 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 37876 : pi_ramanujan(long prec)
149 : {
150 37876 : const ulong B = 545140134, A = 13591409, C = 640320;
151 37876 : 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 37876 : nmax = (long)(1 + prec2nbits(prec)/alpha2);
158 : #ifdef LONG_IS_64BIT
159 37401 : D = utoipos(10939058860032000UL); /* C^3/24 */
160 : #else
161 478 : D = uutoi(2546948UL,495419392UL);
162 : #endif
163 37882 : abpq_init(&S, nmax);
164 37883 : S.a[0] = utoipos(A);
165 37882 : S.b[0] = S.p[0] = S.q[0] = gen_1;
166 309150 : for (n = 1; n <= nmax; n++)
167 : {
168 271307 : S.a[n] = addiu(muluu(B, n), A);
169 271126 : S.b[n] = gen_1;
170 271126 : S.p[n] = mulis(muluu(6*n-5, 2*n-1), 1-6*n);
171 271134 : S.q[n] = mulii(sqru(n), muliu(D,n));
172 : }
173 37843 : abpq_sum(&R, 0, nmax, &S); prec2 = prec+EXTRAPRECWORD;
174 37900 : u = itor(muliu(R.Q,C/12), prec2);
175 37894 : 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 29434332 : constpi(long prec)
213 : {
214 : pari_sp av;
215 : GEN tmp;
216 29434332 : if (gpi && realprec(gpi) >= prec) return gpi;
217 :
218 37635 : av = avma;
219 37635 : tmp = gclone(pi_ramanujan(prec));
220 37902 : swap_clone(&gpi,tmp);
221 37903 : return gc_const(av, gpi);
222 : }
223 :
224 : GEN
225 29434358 : mppi(long prec) { return rtor(constpi(prec), prec); }
226 :
227 : /* Pi * 2^n */
228 : GEN
229 19739794 : Pi2n(long n, long prec)
230 : {
231 19739794 : GEN x = mppi(prec); shiftr_inplace(x, n);
232 19739827 : return x;
233 : }
234 :
235 : /* I * Pi * 2^n */
236 : GEN
237 262235 : PiI2n(long n, long prec) { retmkcomplex(gen_0, Pi2n(n, prec)); }
238 :
239 : /* 2I * Pi */
240 : GEN
241 261346 : PiI2(long prec) { return PiI2n(1, prec); }
242 :
243 : /********************************************************************/
244 : /** **/
245 : /** EULER CONSTANT **/
246 : /** **/
247 : /********************************************************************/
248 :
249 : GEN
250 49239 : 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 49239 : if (geuler && realprec(geuler) >= prec) return geuler;
257 :
258 505 : av1 = avma; tmpeuler = cgetr_block(prec);
259 :
260 505 : incrprec(prec);
261 :
262 505 : l = prec+EXTRAPRECWORD; x = (long) (1 + prec2nbits_mul(l, M_LN2/4));
263 505 : a = utor(x,l); u=logr_abs(a); setsigne(u,-1); affrr(u,a);
264 505 : b = real_1(l);
265 505 : v = real_1(l);
266 505 : n = (long)(1+3.591*x); /* z=3.591: z*[ ln(z)-1 ]=1 */
267 505 : n1 = minss(n, SQRTVERYBIGINT);
268 505 : if (x < SQRTVERYBIGINT)
269 : {
270 505 : ulong xx = x*x;
271 505 : av2 = avma;
272 161013 : for (k=1; k<n1; k++)
273 : {
274 160508 : affrr(divru(mulur(xx,b),k*k), b);
275 160467 : affrr(divru(addrr(divru(mulur(xx,a),k),b),k), a);
276 160445 : affrr(addrr(u,a), u);
277 160331 : affrr(addrr(v,b), v); set_avma(av2);
278 : }
279 1010 : for ( ; k<=n; k++)
280 : {
281 505 : affrr(divru(divru(mulur(xx,b),k),k), b);
282 505 : affrr(divru(addrr(divru(mulur(xx,a),k),b),k), a);
283 505 : affrr(addrr(u,a), u);
284 505 : 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 505 : divrrz(u,v,tmpeuler);
307 505 : swap_clone(&geuler,tmpeuler);
308 505 : return gc_const(av1, geuler);
309 : }
310 :
311 : GEN
312 49239 : 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 1394573 : transvec(GEN (*f)(GEN,long), GEN x, long prec)
364 4693735 : { 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 2504665 : trans_eval(const char *fun, GEN (*f)(GEN,long), GEN x, long prec)
371 : {
372 2504665 : pari_sp av = avma;
373 2504665 : if (prec < 3) pari_err_BUG("trans_eval [prec < 3]");
374 2504682 : switch(typ(x))
375 : {
376 862276 : case t_INT: x = f(itor(x,prec),prec); break;
377 247773 : 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 1394563 : case t_VEC:
381 : case t_COL:
382 1394563 : case t_MAT: return transvec(f, x, prec);
383 49 : default: pari_err_TYPE(fun,x);
384 : return NULL;/*LCOV_EXCL_LINE*/
385 : }
386 1110026 : return gerepileupto(av, x);
387 : }
388 :
389 : GEN
390 1876 : trans_evalgen(const char *fun, void *E, GEN (*f)(void*,GEN,long),
391 : GEN x, long prec)
392 : {
393 1876 : pari_sp av = avma;
394 1876 : if (prec < 3) pari_err_BUG("trans_eval [prec < 3]");
395 1876 : switch(typ(x))
396 : {
397 266 : case t_INT: x = f(E, itor(x,prec),prec); break;
398 1246 : 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 1582 : return gerepileupto(av, x);
408 : }
409 :
410 : /*******************************************************************/
411 : /* */
412 : /* POWERING */
413 : /* */
414 : /*******************************************************************/
415 : /* x a t_REAL 0, return exp(x) */
416 : static GEN
417 68841 : mpexp0(GEN x)
418 : {
419 68841 : long e = expo(x);
420 68841 : return e >= 0? real_0_bit(e): real_1_bit(-e);
421 : }
422 : static GEN
423 20972 : powr0(GEN x)
424 20972 : { return signe(x)? real_1(realprec(x)): mpexp0(x); }
425 :
426 : /* x t_POL or t_SER, return scalarpol(Rg_get_1(x)) */
427 : static GEN
428 366459 : scalarpol_get_1(GEN x)
429 : {
430 366459 : GEN y = cgetg(3,t_POL);
431 366459 : y[1] = evalvarn(varn(x)) | evalsigne(1);
432 366459 : gel(y,2) = Rg_get_1(x); return y;
433 : }
434 : /* to be called by the generic function gpowgs(x,s) when s = 0 */
435 : static GEN
436 1572957 : gpowg0(GEN x)
437 : {
438 : long lx, i;
439 : GEN y;
440 :
441 1572957 : switch(typ(x))
442 : {
443 1163831 : case t_INT: case t_REAL: case t_FRAC: case t_PADIC:
444 1163831 : return gen_1;
445 :
446 7 : case t_QUAD: x++; /*fall through*/
447 36605 : case t_COMPLEX: {
448 36605 : pari_sp av = avma;
449 36605 : GEN a = gpowg0(gel(x,1));
450 36614 : GEN b = gpowg0(gel(x,2));
451 36614 : if (a == gen_1) return b;
452 14 : if (b == gen_1) return a;
453 7 : return gerepileupto(av, gmul(a,b));
454 : }
455 133 : case t_INTMOD:
456 133 : y = cgetg(3,t_INTMOD);
457 133 : gel(y,1) = icopy(gel(x,1));
458 133 : gel(y,2) = is_pm1(gel(x,1))? gen_0: gen_1;
459 133 : return y;
460 :
461 5768 : case t_FFELT: return FF_1(x);
462 :
463 973 : case t_POLMOD:
464 973 : retmkpolmod(scalarpol_get_1(gel(x,1)), gcopy(gel(x,1)));
465 :
466 7 : case t_RFRAC:
467 7 : return scalarpol_get_1(gel(x,2));
468 365479 : case t_POL: case t_SER:
469 365479 : return scalarpol_get_1(x);
470 :
471 84 : case t_MAT:
472 84 : lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
473 77 : if (lx != lgcols(x)) pari_err_DIM("gpow");
474 77 : y = matid(lx-1);
475 252 : for (i=1; i<lx; i++) gcoeff(y,i,i) = gpowg0(gcoeff(x,i,i));
476 77 : return y;
477 21 : case t_VEC: /* handle extended t_QFB */
478 21 : case t_QFB: return qfbpow(x, gen_0);
479 49 : case t_VECSMALL: return identity_perm(lg(x) - 1);
480 : }
481 7 : pari_err_TYPE("gpow",x);
482 : return NULL; /* LCOV_EXCL_LINE */
483 : }
484 :
485 : static GEN
486 5864178 : _sqr(void *data /* ignored */, GEN x) { (void)data; return gsqr(x); }
487 : static GEN
488 2292139 : _mul(void *data /* ignored */, GEN x, GEN y) { (void)data; return gmul(x,y); }
489 : static GEN
490 318490 : _one(void *x) { return gpowg0((GEN) x); }
491 : static GEN
492 14501121 : _sqri(void *data /* ignored */, GEN x) { (void)data; return sqri(x); }
493 : static GEN
494 9047276 : _muli(void *data /* ignored */, GEN x, GEN y) { (void)data; return mulii(x,y); }
495 : static GEN
496 18018884 : _sqrr(void *data /* ignored */, GEN x) { (void)data; return sqrr(x); }
497 : static GEN
498 4836514 : _mulr(void *data /* ignored */, GEN x, GEN y) { (void)data; return mulrr(x,y); }
499 : static GEN
500 13763 : _oner(void *data /* prec */) { return real_1( *(long*) data); }
501 :
502 : /* INTEGER POWERING (a^n for integer a != 0 and integer n > 0)
503 : *
504 : * Use left shift binary algorithm (RS is wasteful: multiplies big numbers,
505 : * with LS one of them is the base, hence small). Sign of result is set
506 : * to s (= 1,-1). Makes life easier for caller, which otherwise might do a
507 : * setsigne(gen_1 / gen_m1) */
508 : static GEN
509 64142859 : powiu_sign(GEN a, ulong N, long s)
510 : {
511 : pari_sp av;
512 : GEN y;
513 :
514 64142859 : if (lgefint(a) == 3)
515 : { /* easy if |a| < 3 */
516 60205442 : ulong q = a[2];
517 60205442 : if (q == 1) return (s>0)? gen_1: gen_m1;
518 51002143 : if (q == 2) { a = int2u(N); setsigne(a,s); return a; }
519 27030685 : q = upowuu(q, N);
520 27033125 : if (q) return s>0? utoipos(q): utoineg(q);
521 : }
522 6357359 : if (N <= 2) {
523 2978466 : if (N == 2) return sqri(a);
524 9013 : a = icopy(a); setsigne(a,s); return a;
525 : }
526 3378893 : av = avma;
527 3378893 : y = gen_powu_i(a, N, NULL, &_sqri, &_muli);
528 3379462 : setsigne(y,s); return gerepileuptoint(av, y);
529 : }
530 : /* a^n */
531 : GEN
532 64113092 : powiu(GEN a, ulong n)
533 : {
534 : long s;
535 64113092 : if (!n) return gen_1;
536 63040868 : s = signe(a);
537 63040868 : if (!s) return gen_0;
538 62967137 : return powiu_sign(a, n, (s < 0 && odd(n))? -1: 1);
539 : }
540 : GEN
541 23507471 : powis(GEN a, long n)
542 : {
543 : long s;
544 : GEN t, y;
545 23507471 : if (n >= 0) return powiu(a, n);
546 577820 : s = signe(a);
547 577820 : if (!s) pari_err_INV("powis",gen_0);
548 577820 : t = (s < 0 && odd(n))? gen_m1: gen_1;
549 577820 : if (is_pm1(a)) return t;
550 : /* n < 0, |a| > 1 */
551 575419 : y = cgetg(3,t_FRAC);
552 575418 : gel(y,1) = t;
553 575418 : gel(y,2) = powiu_sign(a, -n, 1); /* force denominator > 0 */
554 575421 : return y;
555 : }
556 : GEN
557 48406538 : powuu(ulong p, ulong N)
558 : {
559 : pari_sp av;
560 : ulong pN;
561 : GEN y;
562 48406538 : if (!p) return gen_0;
563 48406461 : if (N <= 2)
564 : {
565 41610201 : if (N == 2) return sqru(p);
566 38126994 : if (N == 1) return utoipos(p);
567 5184071 : return gen_1;
568 : }
569 6796260 : pN = upowuu(p, N);
570 6796402 : if (pN) return utoipos(pN);
571 1125780 : if (p == 2) return int2u(N);
572 1112335 : av = avma;
573 1112335 : y = gen_powu_i(utoipos(p), N, NULL, &_sqri, &_muli);
574 1112330 : return gerepileuptoint(av, y);
575 : }
576 :
577 : /* return 0 if overflow */
578 : static ulong
579 9332127 : usqru(ulong p) { return p & HIGHMASK? 0: p*p; }
580 : ulong
581 59338049 : upowuu(ulong p, ulong k)
582 : {
583 : #ifdef LONG_IS_64BIT
584 50252319 : const ulong CUTOFF3 = 2642245;
585 50252319 : const ulong CUTOFF4 = 65535;
586 50252319 : const ulong CUTOFF5 = 7131;
587 50252319 : const ulong CUTOFF6 = 1625;
588 50252319 : const ulong CUTOFF7 = 565;
589 50252319 : const ulong CUTOFF8 = 255;
590 50252319 : const ulong CUTOFF9 = 138;
591 50252319 : const ulong CUTOFF10 = 84;
592 50252319 : const ulong CUTOFF11 = 56;
593 50252319 : const ulong CUTOFF12 = 40;
594 50252319 : const ulong CUTOFF13 = 30;
595 50252319 : const ulong CUTOFF14 = 23;
596 50252319 : const ulong CUTOFF15 = 19;
597 50252319 : const ulong CUTOFF16 = 15;
598 50252319 : const ulong CUTOFF17 = 13;
599 50252319 : const ulong CUTOFF18 = 11;
600 50252319 : const ulong CUTOFF19 = 10;
601 50252319 : const ulong CUTOFF20 = 9;
602 : #else
603 9085730 : const ulong CUTOFF3 = 1625;
604 9085730 : const ulong CUTOFF4 = 255;
605 9085730 : const ulong CUTOFF5 = 84;
606 9085730 : const ulong CUTOFF6 = 40;
607 9085730 : const ulong CUTOFF7 = 23;
608 9085730 : const ulong CUTOFF8 = 15;
609 9085730 : const ulong CUTOFF9 = 11;
610 9085730 : const ulong CUTOFF10 = 9;
611 9085730 : const ulong CUTOFF11 = 7;
612 9085730 : const ulong CUTOFF12 = 6;
613 9085730 : const ulong CUTOFF13 = 5;
614 9085730 : const ulong CUTOFF14 = 4;
615 9085730 : const ulong CUTOFF15 = 4;
616 9085730 : const ulong CUTOFF16 = 3;
617 9085730 : const ulong CUTOFF17 = 3;
618 9085730 : const ulong CUTOFF18 = 3;
619 9085730 : const ulong CUTOFF19 = 3;
620 9085730 : const ulong CUTOFF20 = 3;
621 : #endif
622 :
623 59338049 : if (p <= 2)
624 : {
625 8384394 : if (p < 2) return p;
626 7857391 : return k < BITS_IN_LONG? 1UL<<k: 0;
627 : }
628 50953655 : switch(k)
629 : {
630 : ulong p2, p3, p4, p5, p8;
631 7863376 : case 0: return 1;
632 17908003 : case 1: return p;
633 9332126 : case 2: return usqru(p);
634 3549800 : case 3: if (p > CUTOFF3) return 0; return p*p*p;
635 2321394 : case 4: if (p > CUTOFF4) return 0; p2=p*p; return p2*p2;
636 1950814 : case 5: if (p > CUTOFF5) return 0; p2=p*p; return p2*p2*p;
637 1858810 : case 6: if (p > CUTOFF6) return 0; p2=p*p; return p2*p2*p2;
638 332437 : case 7: if (p > CUTOFF7) return 0; p2=p*p; return p2*p2*p2*p;
639 349820 : case 8: if (p > CUTOFF8) return 0; p2=p*p; p4=p2*p2; return p4*p4;
640 1175133 : case 9: if (p > CUTOFF9) return 0; p2=p*p; p4=p2*p2; return p4*p4*p;
641 257809 : case 10: if (p > CUTOFF10)return 0; p2=p*p; p4=p2*p2; return p4*p4*p2;
642 166046 : case 11: if (p > CUTOFF11)return 0; p2=p*p; p4=p2*p2; return p4*p4*p2*p;
643 136525 : case 12: if (p > CUTOFF12)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4;
644 135163 : case 13: if (p > CUTOFF13)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4*p;
645 129628 : case 14: if (p > CUTOFF14)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4*p2;
646 163103 : case 15: if (p > CUTOFF15)return 0;
647 87587 : p2=p*p; p3=p2*p; p5=p3*p2; return p5*p5*p5;
648 111475 : case 16: if (p > CUTOFF16)return 0;
649 50774 : p2=p*p; p4=p2*p2; p8=p4*p4; return p8*p8;
650 82526 : case 17: if (p > CUTOFF17)return 0;
651 41542 : p2=p*p; p4=p2*p2; p8=p4*p4; return p*p8*p8;
652 71967 : case 18: if (p > CUTOFF18)return 0;
653 39237 : p2=p*p; p4=p2*p2; p8=p4*p4; return p2*p8*p8;
654 790042 : case 19: if (p > CUTOFF19)return 0;
655 762889 : p2=p*p; p4=p2*p2; p8=p4*p4; return p*p2*p8*p8;
656 45499 : case 20: if (p > CUTOFF20)return 0;
657 21800 : p2=p*p; p4=p2*p2; p8=p4*p4; return p4*p8*p8;
658 : }
659 : #ifdef LONG_IS_64BIT
660 1912426 : switch(p)
661 : {
662 227190 : case 3: if (k > 40) return 0;
663 144539 : break;
664 17028 : case 4: if (k > 31) return 0;
665 774 : return 1UL<<(2*k);
666 1036341 : case 5: if (k > 27) return 0;
667 20260 : break;
668 49686 : case 6: if (k > 24) return 0;
669 9216 : break;
670 55807 : case 7: if (k > 22) return 0;
671 2682 : break;
672 526374 : default: return 0;
673 : }
674 : /* no overflow */
675 : {
676 176697 : ulong q = upowuu(p, k >> 1);
677 176697 : q *= q ;
678 176697 : return odd(k)? q*p: q;
679 : }
680 : #else
681 309733 : return 0;
682 : #endif
683 : }
684 :
685 : GEN
686 12017 : upowers(ulong x, long n)
687 : {
688 : long i;
689 12017 : GEN p = cgetg(n + 2, t_VECSMALL);
690 12017 : uel(p,1) = 1; if (n == 0) return p;
691 12017 : uel(p,2) = x;
692 91465 : for (i = 3; i <= n; i++)
693 79448 : uel(p,i) = uel(p,i-1)*x;
694 12017 : return p;
695 : }
696 :
697 : typedef struct {
698 : long prec, a;
699 : GEN (*sqr)(GEN);
700 : GEN (*mulug)(ulong,GEN);
701 : } sr_muldata;
702 :
703 : static GEN
704 1588551 : _rpowuu_sqr(void *data, GEN x)
705 : {
706 1588551 : sr_muldata *D = (sr_muldata *)data;
707 1588551 : if (typ(x) == t_INT && lgefint(x) >= D->prec)
708 : { /* switch to t_REAL */
709 157685 : D->sqr = &sqrr;
710 157685 : D->mulug = &mulur; x = itor(x, D->prec);
711 : }
712 1588551 : return D->sqr(x);
713 : }
714 :
715 : static GEN
716 621085 : _rpowuu_msqr(void *data, GEN x)
717 : {
718 621085 : GEN x2 = _rpowuu_sqr(data, x);
719 621085 : sr_muldata *D = (sr_muldata *)data;
720 621085 : return D->mulug(D->a, x2);
721 : }
722 :
723 : /* return a^n as a t_REAL of precision prec. Assume a > 0, n > 0 */
724 : GEN
725 426939 : rpowuu(ulong a, ulong n, long prec)
726 : {
727 : pari_sp av;
728 : GEN y, z;
729 : sr_muldata D;
730 :
731 426939 : if (a == 1) return real_1(prec);
732 426939 : if (a == 2) return real2n(n, prec);
733 426939 : if (n == 1) return utor(a, prec);
734 422260 : z = cgetr(prec);
735 422260 : av = avma;
736 422260 : D.sqr = &sqri;
737 422260 : D.mulug = &mului;
738 422260 : D.prec = prec;
739 422260 : D.a = (long)a;
740 422260 : y = gen_powu_fold_i(utoipos(a), n, (void*)&D, &_rpowuu_sqr, &_rpowuu_msqr);
741 422260 : mpaff(y, z); return gc_const(av,z);
742 : }
743 :
744 : GEN
745 9278089 : powrs(GEN x, long n)
746 : {
747 9278089 : pari_sp av = avma;
748 : GEN y;
749 9278089 : if (!n) return powr0(x);
750 9278089 : y = gen_powu_i(x, (ulong)labs(n), NULL, &_sqrr, &_mulr);
751 9278787 : if (n < 0) y = invr(y);
752 9278405 : return gerepileuptoleaf(av,y);
753 : }
754 : GEN
755 3421214 : powru(GEN x, ulong n)
756 : {
757 3421214 : pari_sp av = avma;
758 : GEN y;
759 3421214 : if (!n) return powr0(x);
760 3400753 : y = gen_powu_i(x, n, NULL, &_sqrr, &_mulr);
761 3400671 : return gerepileuptoleaf(av,y);
762 : }
763 :
764 : GEN
765 13763 : powersr(GEN x, long n)
766 : {
767 13763 : long prec = realprec(x);
768 13763 : return gen_powers(x, n, 1, &prec, &_sqrr, &_mulr, &_oner);
769 : }
770 :
771 : /* x^(s/2), assume x t_REAL */
772 : GEN
773 0 : powrshalf(GEN x, long s)
774 : {
775 0 : if (s & 1) return sqrtr(powrs(x, s));
776 0 : return powrs(x, s>>1);
777 : }
778 : /* x^(s/2), assume x t_REAL */
779 : GEN
780 115139 : powruhalf(GEN x, ulong s)
781 : {
782 115139 : if (s & 1) return sqrtr(powru(x, s));
783 6759 : return powru(x, s>>1);
784 : }
785 : /* x^(n/d), assume x t_REAL, return t_REAL */
786 : GEN
787 511 : powrfrac(GEN x, long n, long d)
788 : {
789 : long z;
790 511 : if (!n) return powr0(x);
791 0 : z = cgcd(n, d); if (z > 1) { n /= z; d /= z; }
792 0 : if (d == 1) return powrs(x, n);
793 0 : x = powrs(x, n);
794 0 : if (d == 2) return sqrtr(x);
795 0 : return sqrtnr(x, d);
796 : }
797 :
798 : /* assume x != 0 */
799 : static GEN
800 669296 : pow_monome(GEN x, long n)
801 : {
802 669296 : long i, d, dx = degpol(x);
803 : GEN A, b, y;
804 :
805 669296 : if (n < 0) { n = -n; y = cgetg(3, t_RFRAC); } else y = NULL;
806 :
807 669296 : if (HIGHWORD(dx) || HIGHWORD(n))
808 8 : {
809 : LOCAL_HIREMAINDER;
810 9 : d = (long)mulll((ulong)dx, (ulong)n);
811 9 : if (hiremainder || (d &~ LGBITS)) d = LGBITS; /* overflow */
812 9 : d += 2;
813 : }
814 : else
815 669287 : d = dx*n + 2;
816 669296 : if ((d + 1) & ~LGBITS) pari_err(e_OVERFLOW,"pow_monome [degree]");
817 669289 : A = cgetg(d+1, t_POL); A[1] = x[1];
818 6089376 : for (i=2; i < d; i++) gel(A,i) = gen_0;
819 669289 : b = gpowgs(gel(x,dx+2), n); /* not memory clean if (n < 0) */
820 669290 : if (!y) y = A;
821 : else {
822 20468 : GEN c = denom_i(b);
823 20468 : gel(y,1) = c; if (c != gen_1) b = gmul(b,c);
824 20468 : gel(y,2) = A;
825 : }
826 669290 : gel(A,d) = b; return y;
827 : }
828 :
829 : /* x t_PADIC */
830 : static GEN
831 1308829 : powps(GEN x, long n)
832 : {
833 1308829 : long e = n*valp(x), v;
834 1308829 : GEN t, y, mod, p = gel(x,2);
835 : pari_sp av;
836 :
837 1308829 : if (!signe(gel(x,4))) {
838 84 : if (n < 0) pari_err_INV("powps",x);
839 77 : return zeropadic(p, e);
840 : }
841 1308745 : v = z_pval(n, p);
842 :
843 1308744 : y = cgetg(5,t_PADIC);
844 1308742 : mod = gel(x,3);
845 1308742 : if (v == 0) mod = icopy(mod);
846 : else
847 : {
848 90293 : if (precp(x) == 1 && absequaliu(p, 2)) v++;
849 90293 : mod = mulii(mod, powiu(p,v));
850 90293 : mod = gerepileuptoint((pari_sp)y, mod);
851 : }
852 1308742 : y[1] = evalprecp(precp(x) + v) | evalvalp(e);
853 1308742 : gel(y,2) = icopy(p);
854 1308741 : gel(y,3) = mod;
855 :
856 1308741 : av = avma; t = gel(x,4);
857 1308741 : if (n < 0) { t = Fp_inv(t, mod); n = -n; }
858 1308741 : t = Fp_powu(t, n, mod);
859 1308745 : gel(y,4) = gerepileuptoint(av, t);
860 1308744 : return y;
861 : }
862 : /* x t_PADIC */
863 : static GEN
864 161 : powp(GEN x, GEN n)
865 : {
866 : long v;
867 161 : GEN y, mod, p = gel(x,2);
868 :
869 161 : if (valp(x)) pari_err_OVERFLOW("valp()");
870 :
871 161 : if (!signe(gel(x,4))) {
872 14 : if (signe(n) < 0) pari_err_INV("powp",x);
873 7 : return zeropadic(p, 0);
874 : }
875 147 : v = Z_pval(n, p);
876 :
877 147 : y = cgetg(5,t_PADIC);
878 147 : mod = gel(x,3);
879 147 : if (v == 0) mod = icopy(mod);
880 : else
881 : {
882 70 : mod = mulii(mod, powiu(p,v));
883 70 : mod = gerepileuptoint((pari_sp)y, mod);
884 : }
885 147 : y[1] = evalprecp(precp(x) + v) | _evalvalp(0);
886 147 : gel(y,2) = icopy(p);
887 147 : gel(y,3) = mod;
888 147 : gel(y,4) = Fp_pow(gel(x,4), n, mod);
889 147 : return y;
890 : }
891 : static GEN
892 31192 : pow_polmod(GEN x, GEN n)
893 : {
894 31192 : GEN z = cgetg(3, t_POLMOD), a = gel(x,2), T = gel(x,1);
895 31192 : gel(z,1) = gcopy(T);
896 31192 : if (typ(a) != t_POL || varn(a) != varn(T) || lg(a) <= 3)
897 1554 : a = powgi(a, n);
898 : else {
899 29638 : pari_sp av = avma;
900 29638 : GEN p = NULL;
901 29638 : if (RgX_is_FpX(T, &p) && RgX_is_FpX(a, &p) && p)
902 : {
903 7602 : T = RgX_to_FpX(T, p); a = RgX_to_FpX(a, p);
904 7602 : if (lgefint(p) == 3)
905 : {
906 7595 : ulong pp = p[2];
907 7595 : a = Flxq_pow(ZX_to_Flx(a, pp), n, ZX_to_Flx(T, pp), pp);
908 7595 : a = Flx_to_ZX(a);
909 : }
910 : else
911 7 : a = FpXQ_pow(a, n, T, p);
912 7602 : a = FpX_to_mod(a, p);
913 7602 : a = gerepileupto(av, a);
914 : }
915 : else
916 : {
917 22036 : set_avma(av);
918 22036 : a = RgXQ_pow(a, n, gel(z,1));
919 : }
920 : }
921 31192 : gel(z,2) = a; return z;
922 : }
923 :
924 : GEN
925 115172568 : gpowgs(GEN x, long n)
926 : {
927 : long m;
928 : pari_sp av;
929 : GEN y;
930 :
931 115172568 : if (n == 0) return gpowg0(x);
932 113991497 : if (n == 1)
933 72498173 : switch (typ(x)) {
934 52180 : case t_VEC: /* handle extended t_QFB */
935 52180 : case t_QFB: return qfbred(x);
936 72445993 : default: return gcopy(x);
937 : }
938 41493324 : if (n ==-1) return ginv(x);
939 38381999 : switch(typ(x))
940 : {
941 23350414 : case t_INT: return powis(x,n);
942 9269718 : case t_REAL: return powrs(x,n);
943 22365 : case t_INTMOD:
944 22365 : y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
945 22365 : gel(y,2) = Fp_pows(gel(x,2), n, gel(x,1));
946 22365 : return y;
947 299274 : case t_FRAC:
948 : {
949 299274 : GEN a = gel(x,1), b = gel(x,2);
950 299274 : long s = (signe(a) < 0 && odd(n))? -1: 1;
951 299274 : if (n < 0) {
952 700 : n = -n;
953 700 : if (is_pm1(a)) return powiu_sign(b, n, s); /* +-1/x[2] inverts to t_INT */
954 490 : swap(a, b);
955 : }
956 299064 : y = cgetg(3, t_FRAC);
957 299064 : gel(y,1) = powiu_sign(a, n, s);
958 299064 : gel(y,2) = powiu_sign(b, n, 1);
959 299064 : return y;
960 : }
961 1308829 : case t_PADIC: return powps(x, n);
962 249165 : case t_RFRAC:
963 : {
964 249165 : av = avma; y = cgetg(3, t_RFRAC); m = labs(n);
965 249165 : gel(y,1) = gpowgs(gel(x,1),m);
966 249165 : gel(y,2) = gpowgs(gel(x,2),m);
967 249165 : if (n < 0) y = ginv(y);
968 249165 : return gerepileupto(av,y);
969 : }
970 31185 : case t_POLMOD: {
971 31185 : long N[] = {evaltyp(t_INT) | _evallg(3),0,0};
972 31185 : affsi(n,N); return pow_polmod(x, N);
973 : }
974 1151174 : case t_VEC: /* handle extended t_QFB */
975 1151174 : case t_QFB: return qfbpows(x, n);
976 732155 : case t_POL:
977 732155 : if (RgX_is_monomial(x)) return pow_monome(x, n);
978 : default: {
979 2030579 : pari_sp av = avma;
980 2030579 : y = gen_powu_i(x, (ulong)labs(n), NULL, &_sqr, &_mul);
981 2030726 : if (n < 0) y = ginv(y);
982 2030732 : return gerepileupto(av,y);
983 : }
984 : }
985 : }
986 :
987 : /* n a t_INT */
988 : GEN
989 97509434 : powgi(GEN x, GEN n)
990 : {
991 : GEN y;
992 :
993 97509434 : if (!is_bigint(n)) return gpowgs(x, itos(n));
994 : /* probable overflow for nonmodular types (typical exception: (X^0)^N) */
995 25578 : switch(typ(x))
996 : {
997 25297 : case t_INTMOD:
998 25297 : y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
999 25298 : gel(y,2) = Fp_pow(gel(x,2), n, gel(x,1));
1000 25302 : return y;
1001 59 : case t_FFELT: return FF_pow(x,n);
1002 161 : case t_PADIC: return powp(x, n);
1003 :
1004 35 : case t_INT:
1005 35 : if (is_pm1(x)) return (signe(x) < 0 && mpodd(n))? gen_m1: gen_1;
1006 14 : if (signe(x)) pari_err_OVERFLOW("lg()");
1007 7 : if (signe(n) < 0) pari_err_INV("powgi",gen_0);
1008 7 : return gen_0;
1009 7 : case t_FRAC:
1010 7 : pari_err_OVERFLOW("lg()");
1011 :
1012 12 : case t_VEC: /* handle extended t_QFB */
1013 12 : case t_QFB: return qfbpow(x, n);
1014 7 : case t_POLMOD: return pow_polmod(x, n);
1015 7 : default: {
1016 7 : pari_sp av = avma;
1017 7 : y = gen_pow_i(x, n, NULL, &_sqr, &_mul);
1018 7 : if (signe(n) < 0) return gerepileupto(av, ginv(y));
1019 7 : return gerepilecopy(av,y);
1020 : }
1021 : }
1022 : }
1023 :
1024 : /* Assume x = 1 + O(t), n a scalar. Return x^n */
1025 : static GEN
1026 7805 : ser_pow_1(GEN x, GEN n)
1027 : {
1028 : long lx, mi, i, j, d;
1029 7805 : GEN y = cgetg_copy(x, &lx), X = x+2, Y = y + 2;
1030 7805 : y[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(x));
1031 74130 : d = mi = lx-3; while (mi>=1 && isrationalzero(gel(X,mi))) mi--;
1032 7805 : gel(Y,0) = gen_1;
1033 109970 : for (i=1; i<=d; i++)
1034 : {
1035 102165 : pari_sp av = avma;
1036 102165 : GEN s = gen_0;
1037 485758 : for (j=1; j<=minss(i,mi); j++)
1038 : {
1039 383593 : GEN t = gsubgs(gmulgu(n,j),i-j);
1040 383593 : s = gadd(s, gmul(gmul(t, gel(X,j)), gel(Y,i-j)));
1041 : }
1042 102165 : gel(Y,i) = gerepileupto(av, gdivgu(s,i));
1043 : }
1044 7805 : return y;
1045 : }
1046 :
1047 : /* we suppose n != 0, valp(x) = 0 and leading-term(x) != 0. Not stack clean */
1048 : static GEN
1049 7910 : ser_pow(GEN x, GEN n, long prec)
1050 : {
1051 : GEN y, c, lead;
1052 7910 : if (varncmp(gvar(n), varn(x)) <= 0) return gexp(gmul(n, glog(x,prec)), prec);
1053 7805 : lead = gel(x,2);
1054 7805 : if (gequal1(lead)) return ser_pow_1(x, n);
1055 7462 : x = ser_normalize(x);
1056 7462 : if (typ(n) == t_FRAC && !isinexact(lead) && ispower(lead, gel(n,2), &c))
1057 105 : c = powgi(c, gel(n,1));
1058 : else
1059 7357 : c = gpow(lead,n, prec);
1060 7462 : y = gmul(c, ser_pow_1(x, n));
1061 : /* gpow(t_POLMOD,n) can be a t_COL [conjvec] */
1062 7462 : if (typ(y) != t_SER) pari_err_TYPE("gpow", y);
1063 7462 : return y;
1064 : }
1065 :
1066 : static long
1067 7819 : val_from_i(GEN E)
1068 : {
1069 7819 : if (is_bigint(E)) pari_err_OVERFLOW("sqrtn [valuation]");
1070 7812 : return itos(E);
1071 : }
1072 :
1073 : /* return x^q, assume typ(x) = t_SER, typ(q) = t_INT/t_FRAC and q != 0 */
1074 : static GEN
1075 7826 : ser_powfrac(GEN x, GEN q, long prec)
1076 : {
1077 7826 : GEN y, E = gmulsg(valp(x), q);
1078 : long e;
1079 :
1080 7826 : if (!signe(x))
1081 : {
1082 21 : if (gsigne(q) < 0) pari_err_INV("gpow", x);
1083 21 : return zeroser(varn(x), val_from_i(gfloor(E)));
1084 : }
1085 7805 : if (typ(E) != t_INT)
1086 7 : pari_err_DOMAIN("sqrtn", "valuation", "!=", mkintmod(gen_0, gel(q,2)), x);
1087 7798 : e = val_from_i(E);
1088 7798 : y = leafcopy(x); setvalp(y, 0);
1089 7798 : y = ser_pow(y, q, prec);
1090 7798 : setvalp(y, e); return y;
1091 : }
1092 :
1093 : static GEN
1094 126 : gpow0(GEN x, GEN n, long prec)
1095 : {
1096 126 : pari_sp av = avma;
1097 : long i, lx;
1098 : GEN y;
1099 126 : switch(typ(n))
1100 : {
1101 84 : case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD:
1102 84 : break;
1103 35 : case t_VEC: case t_COL: case t_MAT:
1104 35 : y = cgetg_copy(n, &lx);
1105 105 : for (i=1; i<lx; i++) gel(y,i) = gpow0(x,gel(n,i),prec);
1106 35 : return y;
1107 7 : default: pari_err_TYPE("gpow(0,n)", n);
1108 : }
1109 84 : n = real_i(n);
1110 84 : if (gsigne(n) <= 0) pari_err_DOMAIN("gpow(0,n)", "n", "<=", gen_0, n);
1111 77 : if (!precision(x)) return gcopy(x);
1112 :
1113 14 : x = ground(gmulsg(gexpo(x),n));
1114 14 : if (is_bigint(x) || uel(x,2) >= HIGHEXPOBIT)
1115 7 : pari_err_OVERFLOW("gpow");
1116 7 : set_avma(av); return real_0_bit(itos(x));
1117 : }
1118 :
1119 : /* centermod(x, log(2)), set *sh to the quotient */
1120 : static GEN
1121 12998282 : modlog2(GEN x, long *sh)
1122 : {
1123 12998282 : double d = rtodbl(x), qd = (fabs(d) + M_LN2/2)/M_LN2;
1124 12998539 : long q = (long)qd;
1125 12998539 : if (dblexpo(qd) >= BITS_IN_LONG-1) pari_err_OVERFLOW("expo()");
1126 12998689 : if (d < 0) q = -q;
1127 12998689 : *sh = q;
1128 12998689 : if (q) {
1129 10615489 : long l = realprec(x) + 1;
1130 10615489 : x = subrr(rtor(x,l), mulsr(q, mplog2(l)));
1131 10614270 : if (!signe(x)) return NULL;
1132 : }
1133 12997470 : return x;
1134 : }
1135 :
1136 : /* x^n, n a t_FRAC */
1137 : static GEN
1138 4001804 : powfrac(GEN x, GEN n, long prec)
1139 : {
1140 4001804 : GEN a = gel(n,1), d = gel(n,2);
1141 4001804 : long D = itos_or_0(d);
1142 4001775 : if (D == 2)
1143 : {
1144 3223333 : GEN y = gsqrt(x,prec);
1145 3224724 : if (!equali1(a)) y = gmul(y, powgi(x, shifti(subiu(a,1), -1)));
1146 3223527 : return y;
1147 : }
1148 778442 : if (D && (is_real_t(typ(x)) && gsigne(x) > 0))
1149 : {
1150 : GEN z;
1151 776682 : prec += nbits2extraprec(expi(a));
1152 776686 : if (typ(x) != t_REAL) x = gtofp(x, prec);
1153 776686 : z = sqrtnr(x, D);
1154 776686 : if (!equali1(a)) z = powgi(z, a);
1155 776686 : return z;
1156 : }
1157 1760 : return NULL;
1158 : }
1159 :
1160 : /* n = a+ib, x > 0 real, ex ~ |log2(x)|; return precision at which
1161 : * log(x) must be computed to evaluate x^n */
1162 : long
1163 178676 : powcx_prec(long ex, GEN n, long prec)
1164 : {
1165 178676 : GEN a = gel(n,1), b = gel(n,2);
1166 178676 : long e = (ex < 2)? 0: expu(ex);
1167 178676 : e += gexpo_safe(is_rational_t(typ(a))? b: n);
1168 178676 : return e > 2? prec + nbits2extraprec(e): prec;
1169 : }
1170 : GEN
1171 257449 : powcx(GEN x, GEN logx, GEN n, long prec)
1172 : {
1173 257449 : GEN sxb, cxb, xa, a = gel(n,1), xb = gmul(gel(n,2), logx);
1174 257449 : long sh, p = realprec(logx);
1175 257449 : switch(typ(a))
1176 : {
1177 44563 : case t_INT: xa = powgi(x, a); break;
1178 123690 : case t_FRAC: xa = powfrac(x, a, prec);
1179 123690 : if (xa) break;
1180 : default:
1181 89203 : xa = modlog2(gmul(gel(n,1), logx), &sh);
1182 89203 : if (!xa) xa = real2n(sh, prec);
1183 : else
1184 : {
1185 89203 : if (signe(xa) && realprec(xa) > prec) setlg(xa, prec);
1186 89203 : xa = mpexp(xa); shiftr_inplace(xa, sh);
1187 : }
1188 : }
1189 257449 : if (typ(xb) != t_REAL) return xa;
1190 257449 : if (gexpo(xb) > 30)
1191 : {
1192 0 : GEN q, P = Pi2n(-2, p), z = addrr(xb,P); /* = x + Pi/4 */
1193 0 : shiftr_inplace(P, 1);
1194 0 : q = floorr(divrr(z, P)); /* round ( x / (Pi/2) ) */
1195 0 : xb = subrr(xb, mulir(q, P)); /* x mod Pi/2 */
1196 0 : sh = Mod4(q);
1197 : }
1198 : else
1199 : {
1200 257449 : long q = floor(rtodbl(xb) / (M_PI/2) + 0.5);
1201 257449 : if (q) xb = subrr(xb, mulsr(q, Pi2n(-1,p))); /* x mod Pi/2 */
1202 257449 : sh = q & 3;
1203 : }
1204 257449 : if (signe(xb) && realprec(xb) > prec) setlg(xb, prec);
1205 257449 : mpsincos(xb, &sxb, &cxb);
1206 257449 : return gmul(xa, mulcxpowIs(mkcomplex(cxb, sxb), sh));
1207 : }
1208 :
1209 : GEN
1210 19254470 : gpow(GEN x, GEN n, long prec)
1211 : {
1212 19254470 : long prec0, i, lx, tx, tn = typ(n);
1213 : pari_sp av;
1214 : GEN y;
1215 :
1216 19254470 : if (tn == t_INT) return powgi(x,n);
1217 4214393 : tx = typ(x);
1218 4214393 : if (is_matvec_t(tx))
1219 : {
1220 49 : y = cgetg_copy(x, &lx);
1221 133 : for (i=1; i<lx; i++) gel(y,i) = gpow(gel(x,i),n,prec);
1222 49 : return y;
1223 : }
1224 4214407 : av = avma;
1225 4214407 : switch (tx)
1226 : {
1227 28 : case t_POL: case t_RFRAC: x = toser_i(x); /* fall through */
1228 7560 : case t_SER:
1229 7560 : if (tn == t_FRAC) return gerepileupto(av, ser_powfrac(x, n, prec));
1230 140 : if (valp(x))
1231 21 : pari_err_DOMAIN("gpow [irrational exponent]",
1232 : "valuation", "!=", gen_0, x);
1233 119 : if (lg(x) == 2) return gerepilecopy(av, x); /* O(1) */
1234 112 : return gerepileupto(av, ser_pow(x, n, prec));
1235 : }
1236 4206841 : if (gequal0(x)) return gpow0(x, n, prec);
1237 4206778 : if (tn == t_FRAC)
1238 : {
1239 3881142 : GEN p, z, a = gel(n,1), d = gel(n,2);
1240 3881142 : switch (tx)
1241 : {
1242 674226 : case t_INT:
1243 674226 : if (signe(x) < 0)
1244 : {
1245 42 : if (equaliu(d, 2) && Z_issquareall(negi(x), &z))
1246 : {
1247 21 : z = powgi(z, a);
1248 21 : if (Mod4(a) == 3) z = gneg(z);
1249 3879269 : return gerepilecopy(av, mkcomplex(gen_0, z));
1250 : }
1251 21 : break;
1252 : }
1253 674184 : if (ispower(x, d, &z)) return powgi(z, a);
1254 672627 : break;
1255 69820 : case t_FRAC:
1256 69820 : if (signe(gel(x,1)) < 0)
1257 : {
1258 28 : if (equaliu(d, 2) && ispower(absfrac(x), d, &z))
1259 7 : return gerepilecopy(av, mkcomplex(gen_0, powgi(z, a)));
1260 21 : break;
1261 : }
1262 69792 : if (ispower(x, d, &z)) return powgi(z, a);
1263 68420 : break;
1264 :
1265 21 : case t_INTMOD:
1266 21 : p = gel(x,1);
1267 21 : if (!BPSW_psp(p)) pari_err_PRIME("gpow",p);
1268 14 : y = cgetg(3,t_INTMOD); gel(y,1) = icopy(p);
1269 14 : av = avma;
1270 14 : z = Fp_sqrtn(gel(x,2), d, p, NULL);
1271 14 : if (!z) pari_err_SQRTN("gpow",x);
1272 7 : gel(y,2) = gerepileuptoint(av, Fp_pow(z, a, p));
1273 7 : return y;
1274 :
1275 14 : case t_PADIC:
1276 14 : z = Qp_sqrtn(x, d, NULL); if (!z) pari_err_SQRTN("gpow",x);
1277 7 : return gerepileupto(av, powgi(z, a));
1278 :
1279 21 : case t_FFELT:
1280 21 : return gerepileupto(av,FF_pow(FF_sqrtn(x,d,NULL),a));
1281 : }
1282 3878129 : z = powfrac(x, n, prec);
1283 3878259 : if (z) return gerepileupto(av, z);
1284 : }
1285 327389 : if (tn == t_COMPLEX && is_real_t(typ(x)) && gsigne(x) > 0)
1286 : {
1287 168057 : long p = powcx_prec(fabs(dbllog2(x)), n, prec);
1288 168057 : return gerepileupto(av, powcx(x, glog(x, p), n, prec));
1289 : }
1290 159332 : if (tn == t_PADIC) x = gcvtop(x, gel(n,2), precp(n));
1291 159332 : i = precision(n);
1292 159333 : if (i) prec = i;
1293 159333 : prec0 = prec;
1294 159333 : if (!gprecision(x))
1295 : {
1296 38206 : long e = gexpo_safe(n); /* avoided if n = 0 or gexpo not defined */
1297 38206 : if (e > 2) prec += nbits2extraprec(e);
1298 : }
1299 159333 : y = gmul(n, glog(x,prec));
1300 159305 : y = gexp(y,prec);
1301 159305 : if (prec0 == prec) return gerepileupto(av, y);
1302 29246 : return gerepilecopy(av, gprec_wtrunc(y,prec0));
1303 : }
1304 : GEN
1305 9786 : powPis(GEN s, long prec)
1306 : {
1307 9786 : pari_sp av = avma;
1308 : GEN x;
1309 9786 : if (typ(s) != t_COMPLEX) return gpow(mppi(prec), s, prec);
1310 427 : x = mppi(powcx_prec(1, s, prec));
1311 427 : return gerepileupto(av, powcx(x, logr_abs(x), s, prec));
1312 : }
1313 : GEN
1314 7861 : pow2Pis(GEN s, long prec)
1315 : {
1316 7861 : pari_sp av = avma;
1317 : GEN x;
1318 7861 : if (typ(s) != t_COMPLEX) return gpow(Pi2n(1,prec), s, prec);
1319 1876 : x = Pi2n(1, powcx_prec(2, s, prec));
1320 1876 : return gerepileupto(av, powcx(x, logr_abs(x), s, prec));
1321 : }
1322 :
1323 : GEN
1324 199744 : gpowers0(GEN x, long n, GEN x0)
1325 : {
1326 : long i, l;
1327 : GEN V;
1328 199744 : if (!x0) return gpowers(x,n);
1329 185258 : if (n < 0) return cgetg(1,t_VEC);
1330 185258 : l = n+2; V = cgetg(l, t_VEC); gel(V,1) = gcopy(x0);
1331 7609238 : for (i = 2; i < l; i++) gel(V,i) = gmul(gel(V,i-1),x);
1332 185290 : return V;
1333 : }
1334 :
1335 : GEN
1336 318496 : gpowers(GEN x, long n)
1337 : {
1338 318496 : if (n < 0) return cgetg(1,t_VEC);
1339 318489 : return gen_powers(x, n, 1, (void*)x, &_sqr, &_mul, &_one);
1340 : }
1341 :
1342 : /* return [q^1,q^4,...,q^{n^2}] */
1343 : GEN
1344 37429 : gsqrpowers(GEN q, long n)
1345 : {
1346 37429 : pari_sp av = avma;
1347 37429 : GEN L = gpowers0(gsqr(q), n, q); /* L[i] = q^(2i - 1), i <= n+1 */
1348 37429 : GEN v = cgetg(n+1, t_VEC);
1349 : long i;
1350 37429 : gel(v, 1) = gcopy(q);
1351 6951917 : for (i = 2; i <= n ; ++i) gel(v, i) = q = gmul(q, gel(L,i)); /* q^(i^2) */
1352 37429 : return gerepileupto(av, v);
1353 : }
1354 :
1355 : /* 4 | N. returns a vector RU which contains exp(2*i*k*Pi/N), k=0..N-1 */
1356 : static GEN
1357 568434 : grootsof1_4(long N, long prec)
1358 : {
1359 568434 : GEN z, RU = cgetg(N+1,t_COL), *v = ((GEN*)RU) + 1;
1360 568432 : long i, N2 = (N>>1), N4 = (N>>2), N8 = (N>>3);
1361 : /* z^N2 = -1, z^N4 = I; if z^k = a+I*b, then z^(N4-k) = I*conj(z) = b+a*I */
1362 :
1363 568432 : v[0] = gen_1; v[1] = z = rootsof1u_cx(N, prec);
1364 568434 : if (odd(N4)) N8++;
1365 678301 : for (i=1; i<N8; i++)
1366 : {
1367 109873 : GEN t = v[i];
1368 109873 : v[i+1] = gmul(z, t);
1369 109872 : v[N4-i] = mkcomplex(gel(t,2), gel(t,1));
1370 : }
1371 1669078 : for (i=0; i<N4; i++) v[i+N4] = mulcxI(v[i]);
1372 2769725 : for (i=0; i<N2; i++) v[i+N2] = gneg(v[i]);
1373 568431 : return RU;
1374 : }
1375 :
1376 : /* as above, N arbitrary */
1377 : GEN
1378 700817 : grootsof1(long N, long prec)
1379 : {
1380 : GEN z, RU, *v;
1381 : long i, k;
1382 :
1383 700817 : if (N <= 0) pari_err_DOMAIN("rootsof1", "N", "<=", gen_0, stoi(N));
1384 700804 : if ((N & 3) == 0) return grootsof1_4(N, prec);
1385 132370 : if (N <= 2) return N == 1? mkcol(gen_1): mkcol2(gen_1, gen_m1);
1386 14211 : k = (N+1)>>1;
1387 14211 : RU = cgetg(N+1,t_COL);
1388 14211 : v = ((GEN*)RU) + 1;
1389 14211 : v[0] = gen_1; v[1] = z = rootsof1u_cx(N, prec);
1390 74375 : for (i=2; i<k; i++) v[i] = gmul(z, v[i-1]);
1391 14211 : if (!odd(N)) v[i++] = gen_m1; /*avoid loss of accuracy*/
1392 88586 : for ( ; i<N; i++) v[i] = gconj(v[N-i]);
1393 14211 : return RU;
1394 : }
1395 :
1396 : /********************************************************************/
1397 : /** **/
1398 : /** RACINE CARREE **/
1399 : /** **/
1400 : /********************************************************************/
1401 : /* assume x unit, e = precp(x) */
1402 : GEN
1403 144690 : Z2_sqrt(GEN x, long e)
1404 : {
1405 144690 : ulong r = signe(x)>=0?mod16(x):16-mod16(x);
1406 : GEN z;
1407 : long ez;
1408 : pari_sp av;
1409 :
1410 144690 : switch(e)
1411 : {
1412 7 : case 1: return gen_1;
1413 161 : case 2: return (r & 3UL) == 1? gen_1: NULL;
1414 28 : case 3: return (r & 7UL) == 1? gen_1: NULL;
1415 71064 : case 4: if (r == 1) return gen_1;
1416 35133 : else return (r == 9)? utoipos(3): NULL;
1417 73430 : default: if ((r&7UL) != 1) return NULL;
1418 : }
1419 73430 : av = avma; z = (r==1)? gen_1: utoipos(3);
1420 73430 : ez = 3; /* number of correct bits in z (compared to sqrt(x)) */
1421 : for(;;)
1422 47978 : {
1423 : GEN mod;
1424 121408 : ez = (ez<<1) - 1;
1425 121408 : if (ez > e) ez = e;
1426 121408 : mod = int2n(ez);
1427 121408 : z = addii(z, remi2n(mulii(x, Fp_inv(z,mod)), ez));
1428 121408 : z = shifti(z, -1); /* (z + x/z) / 2 */
1429 121408 : if (e == ez) return gerepileuptoint(av, z);
1430 47978 : if (ez < e) ez--;
1431 47978 : if (gc_needed(av,2))
1432 : {
1433 0 : if (DEBUGMEM > 1) pari_warn(warnmem,"Qp_sqrt");
1434 0 : z = gerepileuptoint(av,z);
1435 : }
1436 : }
1437 : }
1438 :
1439 : /* x unit defined modulo p^e, e > 0 */
1440 : GEN
1441 1897 : Qp_sqrt(GEN x)
1442 : {
1443 1897 : long pp, e = valp(x);
1444 1897 : GEN z,y,mod, p = gel(x,2);
1445 :
1446 1897 : if (gequal0(x)) return zeropadic(p, (e+1) >> 1);
1447 1883 : if (e & 1) return NULL;
1448 :
1449 1869 : y = cgetg(5,t_PADIC);
1450 1869 : pp = precp(x);
1451 1869 : mod = gel(x,3);
1452 1869 : z = gel(x,4); /* lift to t_INT */
1453 1869 : e >>= 1;
1454 1869 : z = Zp_sqrt(z, p, pp);
1455 1869 : if (!z) return NULL;
1456 1806 : if (absequaliu(p,2))
1457 : {
1458 805 : pp = (pp <= 3) ? 1 : pp-1;
1459 805 : mod = int2n(pp);
1460 : }
1461 1001 : else mod = icopy(mod);
1462 1806 : y[1] = evalprecp(pp) | evalvalp(e);
1463 1806 : gel(y,2) = icopy(p);
1464 1806 : gel(y,3) = mod;
1465 1806 : gel(y,4) = z; 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 : gerepileall(btop, 2, &b, &m);
1501 : }
1502 420 : return gerepileupto(ltop, b);
1503 : }
1504 :
1505 : static GEN
1506 22351 : sqrt_ser(GEN b, long prec)
1507 : {
1508 22351 : long e = valp(b), vx = varn(b), lx, lold, j;
1509 : ulong mask;
1510 : GEN a, x, lta, ltx;
1511 :
1512 22351 : if (!signe(b)) return zeroser(vx, e>>1);
1513 22351 : a = leafcopy(b);
1514 22351 : x = cgetg_copy(b, &lx);
1515 22351 : if (e & 1)
1516 14 : pari_err_DOMAIN("sqrtn", "valuation", "!=", mkintmod(gen_0, gen_2), b);
1517 22337 : a[1] = x[1] = evalsigne(1) | evalvarn(0) | _evalvalp(0);
1518 22337 : lta = gel(a,2);
1519 22337 : if (gequal1(lta)) ltx = lta;
1520 14833 : else if (!issquareall(lta,<x)) ltx = gsqrt(lta,prec);
1521 22330 : gel(x,2) = ltx;
1522 335993 : for (j = 3; j < lx; j++) gel(x,j) = gen_0;
1523 22330 : setlg(x,3);
1524 22330 : mask = quadratic_prec_mask(lx - 2);
1525 22330 : lold = 1;
1526 107887 : while (mask > 1)
1527 : {
1528 85557 : GEN y, x2 = gmul2n(x,1);
1529 85557 : long l = lold << 1, lx;
1530 :
1531 85557 : if (mask & 1) l--;
1532 85557 : mask >>= 1;
1533 85557 : setlg(a, l + 2);
1534 85557 : setlg(x, l + 2);
1535 85557 : y = sqr_ser_part(x, lold, l-1) - lold;
1536 399220 : for (j = lold+2; j < l+2; j++) gel(y,j) = gsub(gel(y,j), gel(a,j));
1537 85557 : y += lold; setvalp(y, lold);
1538 85557 : y = normalizeser(y);
1539 85557 : y = gsub(x, gdiv(y, x2)); /* = gmul2n(gsub(x, gdiv(a,x)), -1); */
1540 85557 : lx = minss(l+2, lg(y));
1541 399213 : for (j = lold+2; j < lx; j++) gel(x,j) = gel(y,j);
1542 85557 : lold = l;
1543 : }
1544 22330 : x[1] = evalsigne(1) | evalvarn(vx) | _evalvalp(e >> 1);
1545 22330 : return x;
1546 : }
1547 :
1548 : GEN
1549 42116731 : gsqrt(GEN x, long prec)
1550 : {
1551 : pari_sp av;
1552 : GEN y;
1553 :
1554 42116731 : switch(typ(x))
1555 : {
1556 380278 : case t_INT:
1557 380278 : if (!signe(x)) return real_0(prec); /* no loss of accuracy */
1558 380208 : x = itor(x,prec); /* fall through */
1559 32510445 : 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 9337621 : 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 9337621 : GEN a = gel(x,1), b = gel(x,2), r, u, v;
1577 9337621 : if (isrationalzero(b)) return gsqrt(a, prec);
1578 9337623 : y = cgetg(3,t_COMPLEX); av = avma;
1579 :
1580 9337625 : r = cxnorm(x);
1581 9337617 : if (typ(r) == t_INTMOD || typ(r) == t_PADIC)
1582 0 : pari_err_IMPL("sqrt(complex of t_INTMODs)");
1583 9337617 : r = gsqrt(r, prec); /* t_REAL, |a+Ib| */
1584 9337621 : if (!signe(r))
1585 67 : u = v = gerepileuptoleaf(av, sqrtr(r));
1586 9337554 : 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 133689 : v = sqrtr( gmul2n(gsub(r,a), -1) );
1591 133689 : if (gsigne(b) < 0) togglesign(v);
1592 133689 : v = gerepileuptoleaf(av, v); av = avma;
1593 : /* v = 0 is impossible */
1594 133689 : u = gerepileuptoleaf(av, gdiv(b, shiftr(v,1)));
1595 : } else {
1596 9203865 : u = sqrtr( gmul2n(gadd(r,a), -1) );
1597 9203868 : u = gerepileuptoleaf(av, u); av = avma;
1598 9203870 : if (!signe(u)) /* possible if a = 0.0, e.g. sqrt(0.e-10+1e-10*I) */
1599 7 : v = u;
1600 : else
1601 9203863 : v = gerepileuptoleaf(av, gdiv(b, shiftr(u,1)));
1602 : }
1603 9337626 : gel(y,1) = u;
1604 9337626 : 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 70 : case t_FFELT: return FF_sqrt(x);
1613 :
1614 268423 : default:
1615 268423 : av = avma; if (!(y = toser_i(x))) break;
1616 22351 : return gerepilecopy(av, sqrt_ser(y, prec));
1617 : }
1618 246072 : return trans_eval("sqrt",gsqrt,x,prec);
1619 : }
1620 : /********************************************************************/
1621 : /** **/
1622 : /** N-th ROOT **/
1623 : /** **/
1624 : /********************************************************************/
1625 :
1626 : static GEN
1627 424937 : Z_to_padic(GEN a, GEN p, long e)
1628 : {
1629 424937 : if (signe(a)==0)
1630 1029 : return zeropadic(p, e);
1631 : else
1632 : {
1633 423908 : GEN z = cgetg(5, t_PADIC);
1634 423908 : long v = Z_pvalrem(a, p, &a), d = e - v;
1635 423908 : z[1] = evalprecp(d) | evalvalp(v);
1636 423908 : gel(z,2) = icopy(p);
1637 423908 : gel(z,3) = powiu(p, d);
1638 423908 : gel(z,4) = a;
1639 423908 : return z;
1640 : }
1641 : }
1642 :
1643 : GEN
1644 317729 : Qp_log(GEN x)
1645 : {
1646 317729 : pari_sp av = avma;
1647 317729 : GEN y, p = gel(x,2), a = gel(x,4);
1648 317729 : long e = precp(x);
1649 :
1650 317729 : if (!signe(a)) pari_err_DOMAIN("Qp_log", "argument", "=", gen_0, x);
1651 317708 : if (absequaliu(p,2) || equali1(modii(a, p)))
1652 108703 : y = Zp_log(a, p, e);
1653 : else
1654 : { /* compute log(x^(p-1)) / (p-1) */
1655 209005 : GEN q = gel(x,3), t = subiu(p, 1);
1656 209005 : a = Fp_pow(a, t, q);
1657 209005 : y = Fp_mul(Zp_log(a, p, e), diviiexact(subsi(1, q), t), q);
1658 : }
1659 317708 : return gerepileupto(av, Z_to_padic(y, p, e));
1660 : }
1661 :
1662 : static GEN Qp_exp_safe(GEN x);
1663 :
1664 : /*compute the p^e th root of x p-adic, assume x != 0 */
1665 : static GEN
1666 854 : Qp_sqrtn_ram(GEN x, long e)
1667 : {
1668 854 : pari_sp ltop=avma;
1669 854 : GEN a, p = gel(x,2), n = powiu(p,e);
1670 854 : long v = valp(x), va;
1671 854 : if (v)
1672 : {
1673 : long z;
1674 161 : v = sdivsi_rem(v, n, &z);
1675 161 : if (z) return NULL;
1676 91 : x = leafcopy(x);
1677 91 : setvalp(x,0);
1678 : }
1679 : /*If p = 2, -1 is a root of 1 in U1: need extra check*/
1680 784 : if (absequaliu(p, 2) && mod8(gel(x,4)) != 1) return NULL;
1681 749 : a = Qp_log(x);
1682 749 : va = valp(a) - e;
1683 749 : if (va <= 0)
1684 : {
1685 287 : if (signe(gel(a,4))) return NULL;
1686 : /* all accuracy lost */
1687 119 : a = cvtop(remii(gel(x,4),p), p, 1);
1688 : }
1689 : else
1690 : {
1691 462 : setvalp(a, va); /* divide by p^e */
1692 462 : a = Qp_exp_safe(a);
1693 462 : if (!a) return NULL;
1694 : /* n=p^e and a^n=z*x where z is a (p-1)th-root of 1.
1695 : * Since z^n=z, we have (a/z)^n = x. */
1696 462 : a = gdiv(x, powgi(a,subiu(n,1))); /* = a/z = x/a^(n-1)*/
1697 462 : if (v) setvalp(a,v);
1698 : }
1699 581 : return gerepileupto(ltop,a);
1700 : }
1701 :
1702 : /*compute the nth root of x p-adic p prime with n*/
1703 : static GEN
1704 2037 : Qp_sqrtn_unram(GEN x, GEN n, GEN *zetan)
1705 : {
1706 : pari_sp av;
1707 2037 : GEN Z, a, r, p = gel(x,2);
1708 2037 : long v = valp(x);
1709 2037 : if (v)
1710 : {
1711 : long z;
1712 84 : v = sdivsi_rem(v,n,&z);
1713 84 : if (z) return NULL;
1714 : }
1715 2030 : r = cgetp(x); setvalp(r,v);
1716 2030 : Z = NULL; /* -Wall */
1717 2030 : if (zetan) Z = cgetp(x);
1718 2030 : av = avma; a = Fp_sqrtn(gel(x,4), n, p, zetan);
1719 2030 : if (!a) return NULL;
1720 2016 : affii(Zp_sqrtnlift(gel(x,4), n, a, p, precp(x)), gel(r,4));
1721 2016 : if (zetan)
1722 : {
1723 14 : affii(Zp_sqrtnlift(gen_1, n, *zetan, p, precp(x)), gel(Z,4));
1724 14 : *zetan = Z;
1725 : }
1726 2016 : return gc_const(av,r);
1727 : }
1728 :
1729 : GEN
1730 2604 : Qp_sqrtn(GEN x, GEN n, GEN *zetan)
1731 : {
1732 : pari_sp av, tetpil;
1733 : GEN q, p;
1734 : long e;
1735 2604 : if (absequaliu(n, 2))
1736 : {
1737 70 : if (zetan) *zetan = gen_m1;
1738 70 : if (signe(n) < 0) x = ginv(x);
1739 63 : return Qp_sqrt(x);
1740 : }
1741 2534 : av = avma; p = gel(x,2);
1742 2534 : if (!signe(gel(x,4)))
1743 : {
1744 203 : if (signe(n) < 0) pari_err_INV("Qp_sqrtn", x);
1745 203 : q = divii(addis(n, valp(x)-1), n);
1746 203 : if (zetan) *zetan = gen_1;
1747 203 : set_avma(av); return zeropadic(p, itos(q));
1748 : }
1749 : /* treat the ramified part using logarithms */
1750 2331 : e = Z_pvalrem(n, p, &q);
1751 2331 : if (e) { x = Qp_sqrtn_ram(x,e); if (!x) return NULL; }
1752 2058 : if (is_pm1(q))
1753 : { /* finished */
1754 21 : if (signe(q) < 0) x = ginv(x);
1755 21 : x = gerepileupto(av, x);
1756 21 : if (zetan)
1757 28 : *zetan = (e && absequaliu(p, 2))? gen_m1 /*-1 in Q_2*/
1758 28 : : gen_1;
1759 21 : return x;
1760 : }
1761 2037 : tetpil = avma;
1762 : /* use hensel lift for unramified case */
1763 2037 : x = Qp_sqrtn_unram(x, q, zetan);
1764 2037 : if (!x) return NULL;
1765 2016 : if (zetan)
1766 : {
1767 : GEN *gptr[2];
1768 14 : if (e && absequaliu(p, 2))/*-1 in Q_2*/
1769 : {
1770 7 : tetpil = avma; x = gcopy(x); *zetan = gneg(*zetan);
1771 : }
1772 14 : gptr[0] = &x; gptr[1] = zetan;
1773 14 : gerepilemanysp(av,tetpil,gptr,2);
1774 14 : return x;
1775 : }
1776 2002 : return gerepile(av,tetpil,x);
1777 : }
1778 :
1779 : GEN
1780 23626 : sqrtnint(GEN a, long n)
1781 : {
1782 23626 : pari_sp av = avma;
1783 : GEN x, b, q;
1784 : long s, k, e;
1785 23626 : const ulong nm1 = n - 1;
1786 23626 : if (n == 2) return sqrtint(a);
1787 19643 : if (typ(a) != t_INT)
1788 : {
1789 35 : if (typ(a) == t_REAL)
1790 : {
1791 : long e;
1792 14 : switch(signe(a))
1793 : {
1794 0 : case 0: return gen_0;
1795 7 : case -1: pari_err_DOMAIN("sqrtnint", "argument", "<", gen_0,a);
1796 : }
1797 7 : e = expo(a); if (e < 0) return gen_0;
1798 7 : if (nbits2lg(e+1) > lg(a))
1799 0 : a = floorr(sqrtnr(a,n)); /* try to avoid precision loss in truncation */
1800 : else
1801 7 : a = sqrtnint(truncr(a),n);
1802 : }
1803 : else
1804 : {
1805 21 : GEN b = gfloor(a);
1806 21 : if (typ(b) != t_INT) pari_err_TYPE("sqrtint",a);
1807 14 : if (signe(b) < 0) pari_err_DOMAIN("sqrtnint", "argument", "<", gen_0,b);
1808 7 : a = sqrtnint(b, n);
1809 : }
1810 14 : return gerepileuptoint(av, a);
1811 : }
1812 19608 : if (n <= 0) pari_err_DOMAIN("sqrtnint", "n", "<=", gen_0, stoi(n));
1813 19601 : if (n == 1) return icopy(a);
1814 17459 : s = signe(a);
1815 17459 : if (s < 0) pari_err_DOMAIN("sqrtnint", "x", "<", gen_0, a);
1816 17459 : if (!s) return gen_0;
1817 17382 : if (lgefint(a) == 3) return utoi(usqrtn(itou(a), n));
1818 11606 : e = expi(a); k = e/(2*n);
1819 11606 : if (k == 0)
1820 : {
1821 : long flag;
1822 291 : if (n > e) return gc_const(av, gen_1);
1823 291 : flag = cmpii(a, powuu(3, n)); set_avma(av);
1824 291 : return (flag < 0) ? gen_2: stoi(3);
1825 : }
1826 11315 : if (e < n*BITS_IN_LONG - 1)
1827 : {
1828 : ulong xs, qs;
1829 4181 : b = itor(a, (2*e < n*BITS_IN_LONG)? DEFAULTPREC: MEDDEFAULTPREC);
1830 4181 : x = mpexp(divru(logr_abs(b), n));
1831 4181 : xs = itou(floorr(x)) + 1; /* >= a^(1/n) */
1832 : for(;;) {
1833 8184 : q = divii(a, powuu(xs, nm1));
1834 8184 : if (lgefint(q) > 3) break;
1835 8177 : qs = itou(q); if (qs >= xs) break;
1836 4003 : xs -= (xs - qs + nm1)/n;
1837 : }
1838 4181 : return utoi(xs);
1839 : }
1840 7134 : b = addui(1, shifti(a, -n*k));
1841 7134 : x = shifti(addui(1, sqrtnint(b, n)), k);
1842 7134 : q = divii(a, powiu(x, nm1));
1843 15994 : while (cmpii(q, x) < 0) /* a priori one iteration, no GC necessary */
1844 : {
1845 8860 : x = subii(x, divis(addui(nm1, subii(x, q)), n));
1846 8860 : q = divii(a, powiu(x, nm1));
1847 : }
1848 7134 : return gerepileuptoleaf(av, x);
1849 : }
1850 :
1851 : ulong
1852 7659 : usqrtn(ulong a, ulong n)
1853 : {
1854 : ulong x, s, q;
1855 7659 : const ulong nm1 = n - 1;
1856 7659 : if (!n) pari_err_DOMAIN("sqrtnint", "n", "=", gen_0, utoi(n));
1857 7659 : if (n == 1 || a == 0) return a;
1858 7659 : s = 1 + expu(a)/n; x = 1UL << s;
1859 7659 : q = (nm1*s >= BITS_IN_LONG)? 0: a >> (nm1*s);
1860 19624 : while (q < x) {
1861 : ulong X;
1862 11965 : x -= (x - q + nm1)/n;
1863 11965 : X = upowuu(x, nm1);
1864 11965 : q = X? a/X: 0;
1865 : }
1866 7659 : return x;
1867 : }
1868 :
1869 : static ulong
1870 843606 : cubic_prec_mask(long n)
1871 : {
1872 843606 : long a = n, i;
1873 843606 : ulong mask = 0;
1874 843606 : for(i = 1;; i++, mask *= 3)
1875 4032485 : {
1876 4876091 : long c = a%3;
1877 4876091 : if (c) mask += 3 - c;
1878 4876091 : a = (a+2)/3;
1879 4876091 : if (a==1) return mask + upowuu(3, i);
1880 : }
1881 : }
1882 :
1883 : /* cubic Newton iteration, |a|^(1/n), assuming a != 0 */
1884 : GEN
1885 1094813 : sqrtnr_abs(GEN a, long n)
1886 : {
1887 : pari_sp av;
1888 : GEN x, b;
1889 : long eextra, eold, n1, n2, prec, B, v;
1890 : ulong mask;
1891 :
1892 1094813 : if (n == 1) return mpabs(a);
1893 1094069 : if (n == 2) return sqrtr_abs(a);
1894 :
1895 1009663 : prec = realprec(a); v = expo(a) / n; av = avma;
1896 1009663 : if (v) a = shiftr(a, -n*v);
1897 1009668 : b = rtor(a, DEFAULTPREC);
1898 1009672 : x = mpexp(divru(logr_abs(b), n));
1899 1009677 : if (prec == DEFAULTPREC)
1900 : {
1901 204756 : if (v) shiftr_inplace(x, v);
1902 204755 : return gerepileuptoleaf(av, x);
1903 : }
1904 804921 : n1 = n+1;
1905 804921 : n2 = 2*n;
1906 804921 : B = prec2nbits(prec);
1907 804921 : eextra = expu(n)-1;
1908 804921 : mask = cubic_prec_mask(B + 63);
1909 804921 : eold = 1;
1910 : for(;;)
1911 3197646 : { /* reach 64 */
1912 4002567 : long enew = eold * 3;
1913 4002567 : enew -= mask % 3;
1914 4002567 : if (enew > 64) break; /* back up one step */
1915 3197646 : mask /= 3;
1916 3197646 : eold = enew;
1917 : }
1918 : for(;;)
1919 666627 : {
1920 1471548 : long pr, enew = eold * 3;
1921 : GEN y, z;
1922 1471548 : enew -= mask % 3;
1923 1471548 : mask /= 3;
1924 1471548 : pr = nbits2prec(enew + eextra);
1925 1471548 : b = rtor(a, pr); setsigne(b,1);
1926 1471548 : x = rtor(x, pr);
1927 1471548 : y = subrr(powru(x, n), b);
1928 1471548 : z = divrr(y, addrr(mulur(n1, y), mulur(n2, b)));
1929 1471548 : shiftr_inplace(z,1);
1930 1471548 : x = mulrr(x, subsr(1,z));
1931 1471548 : if (mask == 1)
1932 : {
1933 804921 : if (v) shiftr_inplace(x, v);
1934 804921 : return gerepileuptoleaf(av, gprec_wtrunc(x,prec));
1935 : }
1936 666627 : eold = enew;
1937 : }
1938 : }
1939 :
1940 : static void
1941 55733 : shiftc_inplace(GEN z, long d)
1942 : {
1943 55733 : shiftr_inplace(gel(z,1), d);
1944 55733 : shiftr_inplace(gel(z,2), d);
1945 55733 : }
1946 :
1947 : /* exp(2*Pi*I/n), same iteration as sqrtnr_abs, different initial point */
1948 : static GEN
1949 446256 : sqrtnof1(ulong n, long prec)
1950 : {
1951 : pari_sp av;
1952 : GEN x;
1953 : long eold, n1, n2, B;
1954 : ulong mask;
1955 :
1956 446256 : B = prec2nbits(prec);
1957 446256 : n1 = n+1;
1958 446256 : n2 = 2*n; av = avma;
1959 :
1960 446256 : x = expIr(divru(Pi2n(1, LOWDEFAULTPREC), n));
1961 446256 : if (prec == LOWDEFAULTPREC) return gerepileupto(av, x);
1962 38685 : mask = cubic_prec_mask(B + BITS_IN_LONG-1);
1963 38685 : eold = 1;
1964 : for(;;)
1965 151164 : { /* reach BITS_IN_LONG */
1966 189849 : long enew = eold * 3;
1967 189849 : enew -= mask % 3;
1968 189849 : if (enew > BITS_IN_LONG) break; /* back up one step */
1969 151164 : mask /= 3;
1970 151164 : eold = enew;
1971 : }
1972 : for(;;)
1973 17048 : {
1974 55733 : long pr, enew = eold * 3;
1975 : GEN y, z;
1976 55733 : enew -= mask % 3;
1977 55733 : mask /= 3;
1978 55733 : pr = nbits2prec(enew);
1979 55733 : x = cxtofp(x, pr);
1980 55733 : y = gsub(gpowgs(x, n), gen_1);
1981 55733 : z = gdiv(y, gaddgs(gmulsg(n1, y), n2));
1982 55733 : shiftc_inplace(z,1);
1983 55733 : x = gmul(x, gsubsg(1, z));
1984 55733 : if (mask == 1) return gerepilecopy(av, gprec_w(x,prec));
1985 17048 : eold = enew;
1986 : }
1987 : }
1988 :
1989 : /* exp(2iPi/d) */
1990 : GEN
1991 1177199 : rootsof1u_cx(ulong n, long prec)
1992 : {
1993 1177199 : switch(n)
1994 : {
1995 13006 : case 1: return gen_1;
1996 2954 : case 2: return gen_m1;
1997 257107 : case 4: return gen_I();
1998 9727 : case 3: case 6: case 12:
1999 : {
2000 9727 : pari_sp av = avma;
2001 9727 : GEN a = (n == 3)? mkfrac(gen_m1,gen_2): ghalf;
2002 9727 : GEN sq3 = sqrtr_abs(utor(3, prec));
2003 9727 : shiftr_inplace(sq3, -1);
2004 9727 : a = (n == 12)? mkcomplex(sq3, a): mkcomplex(a, sq3);
2005 9727 : return gerepilecopy(av, a);
2006 : }
2007 448157 : case 8:
2008 : {
2009 448157 : pari_sp av = avma;
2010 448157 : GEN sq2 = sqrtr_abs(utor(2, prec));
2011 448156 : shiftr_inplace(sq2,-1);
2012 448157 : return gerepilecopy(av, mkcomplex(sq2, sq2));
2013 : }
2014 : }
2015 446248 : return sqrtnof1(n, prec);
2016 : }
2017 : /* e(a/b) */
2018 : GEN
2019 14154 : rootsof1q_cx(long a, long b, long prec)
2020 : {
2021 14154 : long g = cgcd(a,b);
2022 : GEN z;
2023 14154 : if (g != 1) { a /= g; b /= g; }
2024 14154 : if (b < 0) { b = -b; a = -a; }
2025 14154 : z = rootsof1u_cx(b, prec);
2026 14154 : if (a < 0) { z = conj_i(z); a = -a; }
2027 14154 : return gpowgs(z, a);
2028 : }
2029 :
2030 : /* initializes powers of e(a/b) */
2031 : GEN
2032 14987 : rootsof1powinit(long a, long b, long prec)
2033 : {
2034 14987 : long g = cgcd(a,b);
2035 14987 : if (g != 1) { a /= g; b /= g; }
2036 14987 : if (b < 0) { b = -b; a = -a; }
2037 14987 : a %= b; if (a < 0) a += b;
2038 14987 : return mkvec2(grootsof1(b,prec), mkvecsmall2(a,b));
2039 : }
2040 : /* T = rootsof1powinit(a,b); return e(a/b)^c */
2041 : GEN
2042 12516441 : rootsof1pow(GEN T, long c)
2043 : {
2044 12516441 : GEN vz = gel(T,1), ab = gel(T,2);
2045 12516441 : long a = ab[1], b = ab[2]; /* a >= 0, b > 0 */
2046 12516441 : c %= b; if (c < 0) c += b;
2047 12516441 : a = Fl_mul(a, c, b);
2048 12516441 : return gel(vz, a + 1);
2049 : }
2050 :
2051 : /* exp(2iPi/d), assume d a t_INT */
2052 : GEN
2053 4788 : rootsof1_cx(GEN d, long prec)
2054 : {
2055 4788 : if (lgefint(d) == 3) return rootsof1u_cx((ulong)d[2], prec);
2056 0 : return expIr(divri(Pi2n(1,prec), d));
2057 : }
2058 :
2059 : GEN
2060 12469 : gsqrtn(GEN x, GEN n, GEN *zetan, long prec)
2061 : {
2062 : long i, lx, tx;
2063 : pari_sp av;
2064 : GEN y, z;
2065 12469 : if (typ(n)!=t_INT) pari_err_TYPE("sqrtn",n);
2066 12469 : if (!signe(n)) pari_err_DOMAIN("sqrtn", "n", "=", gen_0, n);
2067 12469 : if (is_pm1(n))
2068 : {
2069 70 : if (zetan) *zetan = gen_1;
2070 70 : return (signe(n) > 0)? gcopy(x): ginv(x);
2071 : }
2072 12399 : if (zetan) *zetan = gen_0;
2073 12399 : tx = typ(x);
2074 12399 : if (is_matvec_t(tx))
2075 : {
2076 7 : y = cgetg_copy(x, &lx);
2077 21 : for (i=1; i<lx; i++) gel(y,i) = gsqrtn(gel(x,i),n,NULL,prec);
2078 7 : return y;
2079 : }
2080 12392 : av = avma;
2081 12392 : switch(tx)
2082 : {
2083 182 : case t_INTMOD:
2084 : {
2085 182 : GEN p = gel(x,1), s;
2086 182 : z = gen_0;
2087 182 : y = cgetg(3,t_INTMOD); gel(y,1) = icopy(p);
2088 182 : if (zetan) { z = cgetg(3,t_INTMOD); gel(z,1) = gel(y,1); }
2089 182 : s = Fp_sqrtn(gel(x,2),n,p,zetan);
2090 161 : if (!s) {
2091 35 : if (zetan) return gc_const(av,gen_0);
2092 28 : if (!BPSW_psp(p)) pari_err_PRIME("sqrtn [modulus]",p);
2093 14 : pari_err_SQRTN("gsqrtn",x);
2094 : }
2095 126 : gel(y,2) = s;
2096 126 : if (zetan) { gel(z,2) = *zetan; *zetan = z; }
2097 126 : return y;
2098 : }
2099 :
2100 56 : case t_PADIC:
2101 56 : y = Qp_sqrtn(x,n,zetan);
2102 49 : if (!y) {
2103 7 : if (zetan) return gen_0;
2104 7 : pari_err_SQRTN("gsqrtn",x);
2105 : }
2106 42 : return y;
2107 :
2108 196 : case t_FFELT: return FF_sqrtn(x,n,zetan);
2109 :
2110 11545 : case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
2111 11545 : i = precision(x); if (i) prec = i;
2112 11545 : if (isint1(x))
2113 7 : y = real_1(prec);
2114 11538 : else if (gequal0(x))
2115 : {
2116 : long b;
2117 21 : if (signe(n) < 0) pari_err_INV("gsqrtn",x);
2118 21 : if (isinexactreal(x))
2119 14 : b = sdivsi(gexpo(x), n);
2120 : else
2121 7 : b = -prec2nbits(prec);
2122 21 : if (typ(x) == t_COMPLEX)
2123 : {
2124 7 : y = cgetg(3,t_COMPLEX);
2125 7 : gel(y,1) = gel(y,2) = real_0_bit(b);
2126 : }
2127 : else
2128 14 : y = real_0_bit(b);
2129 : }
2130 : else
2131 : {
2132 11517 : long nn = itos_or_0(n);
2133 11517 : if (tx == t_INT) { x = itor(x,prec); tx = t_REAL; }
2134 11517 : if (nn > 0 && tx == t_REAL && signe(x) > 0)
2135 2062 : y = sqrtnr(x, nn);
2136 : else
2137 9455 : y = gexp(gdiv(glog(x,prec), n), prec);
2138 11517 : y = gerepileupto(av, y);
2139 : }
2140 11545 : if (zetan) *zetan = rootsof1_cx(n, prec);
2141 11545 : return y;
2142 :
2143 7 : case t_QUAD:
2144 7 : return gsqrtn(quadtofp(x, prec), n, zetan, prec);
2145 :
2146 406 : default:
2147 406 : av = avma; if (!(y = toser_i(x))) break;
2148 406 : return gerepileupto(av, ser_powfrac(y, ginv(n), prec));
2149 : }
2150 0 : pari_err_TYPE("sqrtn",x);
2151 : return NULL;/* LCOV_EXCL_LINE */
2152 : }
2153 :
2154 : /********************************************************************/
2155 : /** **/
2156 : /** EXP(X) - 1 **/
2157 : /** **/
2158 : /********************************************************************/
2159 : /* exp(|x|) - 1, assume x != 0.
2160 : * For efficiency, x should be reduced mod log(2): if so, we have a < 0 */
2161 : GEN
2162 12970734 : exp1r_abs(GEN x)
2163 : {
2164 12970734 : long l = realprec(x), a = expo(x), b = prec2nbits(l), L, i, n, m, B;
2165 : GEN y, p2, X;
2166 : pari_sp av;
2167 : double d;
2168 :
2169 12970612 : if (b + a <= 0) return mpabs(x);
2170 :
2171 12957468 : y = cgetr(l); av = avma;
2172 12957431 : B = b/3 + BITS_IN_LONG + (BITS_IN_LONG*BITS_IN_LONG)/ b;
2173 12957431 : d = a/2.; m = (long)(d + sqrt(d*d + B)); /* >= 0 */
2174 12957431 : if (m < (-a) * 0.1) m = 0; /* not worth it */
2175 12957431 : L = l + nbits2extraprec(m);
2176 : /* Multiplication is quadratic in this range (l is small, otherwise we
2177 : * use logAGM + Newton). Set Y = 2^(-e-a) x, compute truncated series
2178 : * sum_{k <= n} Y^k/k!: this costs roughly
2179 : * m b^2 + sum_{k <= n} (k e + BITS_IN_LONG)^2
2180 : * bit operations with n ~ b/e, |x| < 2^(1+a), |Y| < 2^(1-e) , m = e+a and
2181 : * b bits of accuracy needed, so
2182 : * B := (b / 3 + BITS_IN_LONG + BITS_IN_LONG^2 / b) ~ m(m-a)
2183 : * we want b ~ 3 m (m-a) or m~b+a hence
2184 : * m = min( a/2 + sqrt(a^2/4 + B), b + a )
2185 : * NB: e ~ (b/3)^(1/2) as b -> oo
2186 : *
2187 : * Truncate the sum at k = n (>= 1), the remainder is
2188 : * sum_{k >= n+1} Y^k / k! < Y^(n+1) / (n+1)! (1-Y) < Y^(n+1) / n!
2189 : * We want Y^(n+1) / n! <= Y 2^-b, hence -n log_2 |Y| + log_2 n! >= b
2190 : * log n! ~ (n + 1/2) log(n+1) - (n+1) + log(2Pi)/2,
2191 : * error bounded by 1/6(n+1) <= 1/12. Finally, we want
2192 : * n (-1/log(2) -log_2 |Y| + log_2(n+1)) >= b */
2193 12957873 : b += m;
2194 12957873 : d = m-dbllog2(x)-1/M_LN2; /* ~ -log_2 Y - 1/log(2) */
2195 12958646 : n = (long)(b / d);
2196 12958646 : if (n > 1)
2197 12217473 : n = (long)(b / (d + log2((double)n+1))); /* log~constant in small ranges */
2198 29262827 : while (n*(d+log2((double)n+1)) < b) n++; /* expect few corrections */
2199 :
2200 12958646 : X = rtor(x,L); shiftr_inplace(X, -m); setsigne(X, 1);
2201 12958858 : if (n == 1) p2 = X;
2202 : else
2203 : {
2204 12958858 : long s = 0, l1 = nbits2prec((long)(d + n + 16));
2205 12958664 : GEN unr = real_1(L);
2206 : pari_sp av2;
2207 :
2208 12958423 : p2 = cgetr(L); av2 = avma;
2209 167262312 : for (i=n; i>=2; i--, set_avma(av2))
2210 : { /* compute X^(n-1)/n! + ... + X/2 + 1 */
2211 : GEN p1, p3;
2212 154412657 : setprec(X,l1); p3 = divru(X,i);
2213 154616378 : l1 += dvmdsBIL(s - expo(p3), &s); if (l1>L) l1=L;
2214 154647123 : setprec(unr,l1); p1 = addrr_sign(unr,1, i == n? p3: mulrr(p3,p2),1);
2215 154164976 : setprec(p2,l1); affrr(p1,p2); /* p2 <- 1 + (X/i)*p2 */
2216 : }
2217 12956247 : setprec(X,L); p2 = mulrr(X,p2);
2218 : }
2219 :
2220 128302900 : for (i=1; i<=m; i++)
2221 : {
2222 115343538 : if (realprec(p2) > L) setprec(p2,L);
2223 115343538 : p2 = mulrr(p2, addsr(2,p2));
2224 : }
2225 12959362 : affrr_fixlg(p2,y); return gc_const(av,y);
2226 : }
2227 :
2228 : GEN
2229 10581 : mpexpm1(GEN x)
2230 : {
2231 10581 : const long s = 6;
2232 10581 : long l, sx = signe(x);
2233 : GEN y, z;
2234 : pari_sp av;
2235 10581 : if (!sx) return real_0_bit(expo(x));
2236 10574 : l = realprec(x);
2237 10574 : if (l > maxss(EXPNEWTON_LIMIT, (1L<<s) + 2))
2238 : {
2239 6 : long e = expo(x);
2240 6 : if (e < 0) x = rtor(x, l + nbits2extraprec(-e));
2241 6 : return subrs(mpexp(x), 1);
2242 : }
2243 10568 : if (sx > 0) return exp1r_abs(x);
2244 : /* compute exp(x) * (1 - exp(-x)) */
2245 4867 : av = avma; y = exp1r_abs(x);
2246 4867 : z = addsr(1, y); setsigne(z, -1);
2247 4867 : return gerepileupto(av, divrr(y, z));
2248 : }
2249 :
2250 : static GEN serexp(GEN x, long prec);
2251 : GEN
2252 12334 : gexpm1(GEN x, long prec)
2253 : {
2254 12334 : switch(typ(x))
2255 : {
2256 4410 : case t_REAL: return mpexpm1(x);
2257 5880 : case t_COMPLEX: return cxexpm1(x,prec);
2258 14 : case t_PADIC: return gsubgs(Qp_exp(x), 1);
2259 2030 : default:
2260 : {
2261 2030 : pari_sp av = avma;
2262 : long ey;
2263 : GEN y;
2264 2030 : if (!(y = toser_i(x))) break;
2265 2009 : ey = valp(y);
2266 2009 : if (ey < 0) pari_err_DOMAIN("expm1","valuation", "<", gen_0, x);
2267 2009 : if (gequal0(y)) return gcopy(y);
2268 2002 : if (ey)
2269 504 : return gerepileupto(av, gsubgs(serexp(y,prec), 1));
2270 : else
2271 : {
2272 1498 : GEN e1 = gexpm1(gel(y,2), prec), e = gaddgs(e1,1);
2273 1498 : y = gmul(e, serexp(serchop0(y),prec));
2274 1498 : gel(y,2) = e1;
2275 1498 : return gerepilecopy(av, y);
2276 : }
2277 : }
2278 : }
2279 21 : return trans_eval("expm1",gexpm1,x,prec);
2280 : }
2281 : /********************************************************************/
2282 : /** **/
2283 : /** EXP(X) **/
2284 : /** **/
2285 : /********************************************************************/
2286 : static GEN
2287 12908786 : mpexp_basecase(GEN x)
2288 : {
2289 12908786 : pari_sp av = avma;
2290 12908786 : long sh, l = realprec(x);
2291 : GEN y, z;
2292 :
2293 12908786 : y = modlog2(x, &sh);
2294 12908085 : if (!y) { set_avma(av); return real2n(sh, l); }
2295 12908085 : z = addsr(1, exp1r_abs(y));
2296 12907002 : if (signe(y) < 0) z = invr(z);
2297 12907415 : if (sh) {
2298 10526127 : shiftr_inplace(z, sh);
2299 10526132 : if (realprec(z) > l) z = rtor(z, l); /* spurious precision increase */
2300 : }
2301 : #ifdef DEBUG
2302 : {
2303 : GEN t = mplog(z), u = divrr(subrr(x, t),x);
2304 : if (signe(u) && expo(u) > 5-prec2nbits(minss(l,realprec(t))))
2305 : pari_err_BUG("exp");
2306 : }
2307 : #endif
2308 12907846 : return gerepileuptoleaf(av, z); /* NOT affrr, precision often increases */
2309 : }
2310 :
2311 : GEN
2312 12977561 : mpexp(GEN x)
2313 : {
2314 12977561 : const long s = 6; /*Initial steps using basecase*/
2315 12977561 : long i, p, l = realprec(x), sh;
2316 : GEN a, t, z;
2317 : ulong mask;
2318 :
2319 12977561 : if (l <= maxss(EXPNEWTON_LIMIT, (1L<<s) + 2))
2320 : {
2321 12977663 : if (!signe(x)) return mpexp0(x);
2322 12908822 : return mpexp_basecase(x);
2323 : }
2324 11 : z = cgetr(l); /* room for result */
2325 13 : x = modlog2(x, &sh);
2326 13 : if (!x) { set_avma((pari_sp)(z+lg(z))); return real2n(sh, l); }
2327 13 : constpi(l); /* precompute for later logr_abs() */
2328 13 : mask = quadratic_prec_mask(prec2nbits(l)+BITS_IN_LONG);
2329 168 : for(i=0, p=1; i<s+TWOPOTBITS_IN_LONG; i++) { p <<= 1; if (mask & 1) p-=1; mask >>= 1; }
2330 13 : a = mpexp_basecase(rtor(x, nbits2prec(p)));
2331 13 : x = addrs(x,1);
2332 13 : if (realprec(x) < l+EXTRAPRECWORD) x = rtor(x, l+EXTRAPRECWORD);
2333 13 : a = rtor(a, l+EXTRAPRECWORD); /*append 0s */
2334 13 : t = NULL;
2335 : for(;;)
2336 : {
2337 14 : p <<= 1; if (mask & 1) p--;
2338 14 : mask >>= 1;
2339 14 : setprec(x, nbits2prec(p));
2340 14 : setprec(a, nbits2prec(p));
2341 14 : t = mulrr(a, subrr(x, logr_abs(a))); /* a (x - log(a)) */
2342 14 : if (mask == 1) break;
2343 1 : affrr(t, a); set_avma((pari_sp)a);
2344 : }
2345 13 : affrr(t,z);
2346 13 : if (sh) shiftr_inplace(z, sh);
2347 13 : return gc_const((pari_sp)z, z);
2348 : }
2349 :
2350 : /* x != 0; k = ceil(tn / (te-1)), t = p-1 */
2351 : long
2352 84 : Qp_exp_prec(GEN x)
2353 : {
2354 84 : long e = valp(x), n = precp(x);
2355 : ulong a, b, q, r, p, t;
2356 :
2357 84 : if (e < 1) return -1;
2358 63 : if (e > n) return 1;
2359 63 : p = itos_or_0(gel(x,2));
2360 63 : if (!p) return n / e + 1;
2361 63 : if (p == 2) return e < 2? -1: ceildivuu(n, e - 1);
2362 : /* n >= e > 0, n = qe + r */
2363 : /* tn = q (te-1) + rt + q = (q+1)(te-1) - t(e-r) + q + 1 */
2364 63 : t = p - 1;
2365 63 : if (e == 1) return n + ceildivuu(n, t - 1);
2366 0 : q = n / e;
2367 0 : r = n % e; /* k = q + 1 if rt + q < te */
2368 0 : a = umuluu_or_0(e - r, t); if (!a || a > q) return q + 1;
2369 0 : b = umuluu_or_0(e, t); if (!b) return q + 2;
2370 0 : return q + 1 + ceildivuu(q + 1 - a, b - 1);
2371 : }
2372 :
2373 : static GEN
2374 108739 : Qp_exp_safe(GEN x)
2375 : {
2376 108739 : pari_sp av = avma;
2377 108739 : GEN p = gel(x,2), a = gel(x,4), z;
2378 108739 : long d = precp(x), v = valp(x), e = d+v;
2379 108739 : if (gequal0(x)) return gaddgs(x,1);
2380 107234 : if (v < (equaliu(p,2)? 2:1)) return NULL;
2381 107227 : z = Zp_exp(mulii(a,powiu(p,v)), p, e);
2382 107229 : return gerepileupto(av, Z_to_padic(z, p, e));
2383 : }
2384 :
2385 : GEN
2386 108277 : Qp_exp(GEN x)
2387 : {
2388 108277 : GEN y = Qp_exp_safe(x);
2389 108279 : if (!y) pari_err_DOMAIN("gexp(t_PADIC)","argument","",gen_0,x);
2390 108272 : return y;
2391 : }
2392 :
2393 : static GEN
2394 49 : cos_p(GEN x)
2395 : {
2396 : long k;
2397 : pari_sp av;
2398 : GEN x2, y;
2399 :
2400 49 : if (gequal0(x)) return gaddgs(x,1);
2401 28 : k = Qp_exp_prec(x);
2402 28 : if (k < 0) return NULL;
2403 21 : av = avma; x2 = gsqr(x);
2404 21 : if (k & 1) k--;
2405 105 : for (y=gen_1; k; k-=2)
2406 : {
2407 84 : GEN t = gdiv(gmul(y,x2), muluu(k, k-1));
2408 84 : y = gsubsg(1, t);
2409 : }
2410 21 : return gerepileupto(av, y);
2411 : }
2412 : static GEN
2413 63 : sin_p(GEN x)
2414 : {
2415 : long k;
2416 : pari_sp av;
2417 : GEN x2, y;
2418 :
2419 63 : if (gequal0(x)) return gcopy(x);
2420 42 : k = Qp_exp_prec(x);
2421 42 : if (k < 0) return NULL;
2422 28 : av = avma; x2 = gsqr(x);
2423 28 : if (k & 1) k--;
2424 133 : for (y=gen_1; k; k-=2)
2425 : {
2426 105 : GEN t = gdiv(gmul(y,x2), muluu(k, k+1));
2427 105 : y = gsubsg(1, t);
2428 : }
2429 28 : return gerepileupto(av, gmul(y, x));
2430 : }
2431 :
2432 : static GEN
2433 2500573 : cxexp(GEN x, long prec)
2434 : {
2435 2500573 : GEN r, p1, p2, y = cgetg(3,t_COMPLEX);
2436 2500657 : pari_sp av = avma, tetpil;
2437 : long l;
2438 2500657 : l = precision(x); if (l > prec) prec = l;
2439 2500752 : if (gequal0(gel(x,1)))
2440 : {
2441 346027 : gsincos(gel(x,2),&gel(y,2),&gel(y,1),prec);
2442 346043 : return y;
2443 : }
2444 2154708 : r = gexp(gel(x,1),prec);
2445 2154916 : gsincos(gel(x,2),&p2,&p1,prec);
2446 2155087 : tetpil = avma;
2447 2155087 : gel(y,1) = gmul(r,p1);
2448 2154761 : gel(y,2) = gmul(r,p2);
2449 2154799 : gerepilecoeffssp(av,tetpil,y+1,2);
2450 2155175 : return y;
2451 : }
2452 :
2453 : /* given a t_SER x^v s(x), with s(0) != 0, return x^v(s - s(0)), shallow */
2454 : GEN
2455 35427 : serchop0(GEN s)
2456 : {
2457 35427 : long i, l = lg(s);
2458 : GEN y;
2459 35427 : if (l == 2) return s;
2460 35427 : if (l == 3 && isexactzero(gel(s,2))) return s;
2461 35427 : y = cgetg(l, t_SER); y[1] = s[1];
2462 158844 : gel(y,2) = gen_0; for (i=3; i <l; i++) gel(y,i) = gel(s,i);
2463 35427 : return normalizeser(y);
2464 : }
2465 :
2466 : GEN
2467 42 : serchop_i(GEN s, long n)
2468 : {
2469 42 : long i, m, l = lg(s);
2470 : GEN y;
2471 42 : if (l == 2 || (l == 3 && isexactzero(gel(s,2))))
2472 : {
2473 14 : if (valp(s) < n) { s = shallowcopy(s); setvalp(s,n); }
2474 14 : return s;
2475 : }
2476 28 : m = n - valp(s); if (m < 0) return s;
2477 21 : if (l-m <= 2) return zeroser(varn(s), n);
2478 14 : y = cgetg(l-m, t_SER); y[1] = s[1]; setvalp(y, valp(y)+m);
2479 42 : for (i=m+2; i < l; i++) gel(y,i-m) = gel(s,i);
2480 14 : return normalizeser(y);
2481 : }
2482 : GEN
2483 42 : serchop(GEN s, long n)
2484 : {
2485 42 : pari_sp av = avma;
2486 42 : if (typ(s) != t_SER) pari_err_TYPE("serchop",s);
2487 42 : return gerepilecopy(av, serchop_i(s,n));
2488 : }
2489 :
2490 : static GEN
2491 75278 : serexp(GEN x, long prec)
2492 : {
2493 75278 : long i, j, lx, ly, mi, e = valp(x);
2494 : GEN y, xd, yd;
2495 : pari_sp av;
2496 :
2497 75278 : if (e < 0) pari_err_DOMAIN("exp","valuation", "<", gen_0, x);
2498 75271 : if (gequal0(x)) return gaddsg(1,x);
2499 65177 : lx = lg(x);
2500 65177 : if (e)
2501 : {
2502 : GEN X;
2503 51485 : ly = lx+e; y = cgetg(ly,t_SER);
2504 556346 : mi = lx-1; while (mi>=3 && isrationalzero(gel(x,mi))) mi--;
2505 51485 : mi += e-2;
2506 51485 : y[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(x));
2507 : /* zd[i] = coefficient of X^i in z */
2508 51485 : xd = x+2-e; yd = y+2; ly -= 2;
2509 51485 : X = gel(xd,e); if (e != 1) X = gmulgu(X, e); /* left on stack */
2510 51485 : X = isint1(X)? NULL: X;
2511 51485 : gel(yd,0) = gen_1;
2512 51856 : for (i = 1; i < e; i++) gel(yd,i) = gen_0;
2513 646912 : for ( ; i < ly; i++)
2514 : {
2515 595427 : GEN t = gel(yd,i-e);
2516 595427 : long J = minss(i, mi);
2517 595427 : av = avma; if (X) t = gmul(t, X);
2518 2556918 : for (j = e + 1; j <= J; j++)
2519 1961491 : t = gadd(t, gmulgu(gmul(gel(xd,j),gel(yd,i-j)), j));
2520 595427 : gel(yd,i) = gerepileupto(av, gdivgu(t, i));
2521 : }
2522 51485 : return y;
2523 : }
2524 13692 : av = avma;
2525 13692 : return gerepileupto(av, gmul(gexp(gel(x,2),prec), serexp(serchop0(x),prec)));
2526 : }
2527 :
2528 : static GEN
2529 1468553 : expQ(GEN x, long prec)
2530 : {
2531 1468553 : GEN p, q, z, z0 = NULL;
2532 : pari_sp av;
2533 1468553 : long n, nmax, s, e, b = prec2nbits(prec);
2534 : double ex;
2535 : struct abpq_res R;
2536 : struct abpq S;
2537 :
2538 1468552 : if (typ(x) == t_INT)
2539 : {
2540 24663 : if (!signe(x)) return real_1(prec);
2541 24606 : p = x; q = gen_1;
2542 24606 : e = expi(p);
2543 24606 : if (e > b) return mpexp(itor(x, prec));
2544 : }
2545 : else
2546 : {
2547 1443889 : long ep, eq, B = usqrt(b) / 2;
2548 1443890 : p = gel(x,1); ep = expi(p);
2549 1443890 : q = gel(x,2); eq = expi(q);
2550 1443890 : if (ep > B || eq > B) return mpexp(fractor(x, prec));
2551 14637 : e = ep - eq;
2552 14637 : if (e < -3) prec += nbits2extraprec(-e); /* see addrr 'extend' rule */
2553 : }
2554 39243 : if (e > 2) { z0 = cgetr(prec); prec += EXTRAPRECWORD; b += BITS_IN_LONG; }
2555 39243 : z = cgetr(prec); av = avma;
2556 39240 : if (e > 0)
2557 : { /* simplify x/2^e = p / (q * 2^e) */
2558 2478 : long v = minss(e, vali(p));
2559 2478 : if (v) p = shifti(p, -v);
2560 2478 : if (e - v) q = shifti(q, e - v);
2561 : }
2562 39240 : s = signe(p);
2563 39240 : if (s < 0) p = negi(p);
2564 39239 : ex = exp2(dbllog2(x) - e) * 2.718281828; /* exp(1) * x / 2^e, x / 2^e < 2 */
2565 39238 : nmax = (long)(1 + exp(dbllambertW0(M_LN2 * b / ex)) * ex);
2566 39241 : abpq_init(&S, nmax);
2567 39266 : S.a[0] = S.b[0] = S.p[0] = S.q[0] = gen_1;
2568 3369920 : for (n = 1; n <= nmax; n++)
2569 : {
2570 3330689 : S.a[n] = gen_1;
2571 3330689 : S.b[n] = gen_1;
2572 3330689 : S.p[n] = p;
2573 3330689 : S.q[n] = muliu(q, n);
2574 : }
2575 39231 : abpq_sum(&R, 0, nmax, &S);
2576 39248 : if (s > 0) rdiviiz(R.T, R.Q, z); else rdiviiz(R.Q, R.T, z);
2577 39251 : if (e > 0)
2578 : {
2579 17136 : q = z; while (e--) q = sqrr(q);
2580 2478 : if (z0) { affrr(q, z0); z = z0; } else affrr(q,z);
2581 : }
2582 39251 : return gc_const(av,z);
2583 : }
2584 :
2585 : GEN
2586 11844416 : gexp(GEN x, long prec)
2587 : {
2588 11844416 : switch(typ(x))
2589 : {
2590 1468554 : case t_INT: case t_FRAC: return expQ(x, prec);
2591 7144183 : case t_REAL: return mpexp(x);
2592 2500581 : case t_COMPLEX: return cxexp(x,prec);
2593 70 : case t_PADIC: return Qp_exp(x);
2594 731028 : default:
2595 : {
2596 731028 : pari_sp av = avma;
2597 : GEN y;
2598 731028 : if (!(y = toser_i(x))) break;
2599 59584 : return gerepileupto(av, serexp(y,prec));
2600 : }
2601 : }
2602 672125 : return trans_eval("exp",gexp,x,prec);
2603 : }
2604 :
2605 : /********************************************************************/
2606 : /** **/
2607 : /** AGM(X, Y) **/
2608 : /** **/
2609 : /********************************************************************/
2610 : static int
2611 20246005 : agmr_gap(GEN a, GEN b, long L)
2612 : {
2613 20246005 : GEN d = subrr(b, a);
2614 20245684 : return (signe(d) && expo(d) - expo(b) >= L);
2615 : }
2616 : /* assume x > 0 */
2617 : static GEN
2618 1390878 : agm1r_abs(GEN x)
2619 : {
2620 1390878 : long l = realprec(x), L = 5-prec2nbits(l);
2621 1390878 : GEN a1, b1, y = cgetr(l);
2622 1390879 : pari_sp av = avma;
2623 :
2624 1390879 : a1 = addrr(real_1(l), x); shiftr_inplace(a1, -1);
2625 1390878 : b1 = sqrtr_abs(x);
2626 20246038 : while (agmr_gap(a1,b1,L))
2627 : {
2628 18854897 : GEN a = a1;
2629 18854897 : a1 = addrr(a,b1); shiftr_inplace(a1, -1);
2630 18854950 : b1 = sqrtr_abs(mulrr(a,b1));
2631 : }
2632 1390789 : affrr_fixlg(a1,y); return gc_const(av,y);
2633 : }
2634 :
2635 : struct agmcx_gap_t { long L, ex, cnt; };
2636 :
2637 : static void
2638 508287 : agmcx_init(GEN x, long *prec, struct agmcx_gap_t *S)
2639 : {
2640 508287 : long l = precision(x);
2641 508287 : if (l) *prec = l;
2642 508287 : S->L = 1-prec2nbits(*prec);
2643 508287 : S->cnt = 0;
2644 508287 : S->ex = LONG_MAX;
2645 508287 : }
2646 :
2647 : static long
2648 508287 : agmcx_a_b(GEN x, GEN *a1, GEN *b1, long prec)
2649 : {
2650 508287 : long rotate = 0;
2651 508287 : if (gsigne(real_i(x))<0)
2652 : { /* Rotate by +/-Pi/2, so that the choice of the principal square
2653 : * root gives the optimal AGM. So a1 = +/-I*a1, b1=sqrt(-x). */
2654 1484 : if (gsigne(imag_i(x))<0) { *a1=mulcxI(*a1); rotate=-1; }
2655 980 : else { *a1=mulcxmI(*a1); rotate=1; }
2656 1484 : x = gneg(x);
2657 : }
2658 508287 : *b1 = gsqrt(x, prec);
2659 508287 : return rotate;
2660 : }
2661 : /* return 0 if we must stop the AGM loop (a=b or a ~ b), 1 otherwise */
2662 : static int
2663 8703619 : agmcx_gap(GEN a, GEN b, struct agmcx_gap_t *S)
2664 : {
2665 8703619 : GEN d = gsub(b, a);
2666 8703619 : long ex = S->ex;
2667 8703619 : S->ex = gexpo(d);
2668 8703619 : if (gequal0(d) || S->ex - gexpo(b) < S->L) return 0;
2669 : /* if (S->ex >= ex) we're no longer making progress; twice in a row */
2670 8348623 : if (S->ex < ex) S->cnt = 0;
2671 : else
2672 306737 : if (S->cnt++) return 0;
2673 8195332 : return 1;
2674 : }
2675 : static GEN
2676 508238 : agm1cx(GEN x, long prec)
2677 : {
2678 : struct agmcx_gap_t S;
2679 : GEN a1, b1;
2680 508238 : pari_sp av = avma;
2681 : long rotate;
2682 508238 : agmcx_init(x, &prec, &S);
2683 508238 : a1 = gtofp(gmul2n(gadd(real_1(prec), x), -1), prec);
2684 508238 : rotate = agmcx_a_b(x, &a1, &b1, prec);
2685 8703297 : while (agmcx_gap(a1,b1,&S))
2686 : {
2687 8195059 : GEN a = a1;
2688 8195059 : a1 = gmul2n(gadd(a,b1),-1);
2689 8195059 : b1 = gsqrt(gmul(a,b1), prec);
2690 : }
2691 508238 : if (rotate) a1 = rotate>0 ? mulcxI(a1):mulcxmI(a1);
2692 508238 : return gerepilecopy(av,a1);
2693 : }
2694 :
2695 : GEN
2696 49 : zellagmcx(GEN a0, GEN b0, GEN r, GEN t, long prec)
2697 : {
2698 : struct agmcx_gap_t S;
2699 49 : pari_sp av = avma;
2700 49 : GEN x = gdiv(a0, b0), a1, b1;
2701 : long rotate;
2702 49 : agmcx_init(x, &prec, &S);
2703 49 : a1 = gtofp(gmul2n(gadd(real_1(prec), x), -1), prec);
2704 49 : r = gsqrt(gdiv(gmul(a1,gaddgs(r, 1)),gadd(r, x)), prec);
2705 49 : t = gmul(r, t);
2706 49 : rotate = agmcx_a_b(x, &a1, &b1, prec);
2707 322 : while (agmcx_gap(a1,b1,&S))
2708 : {
2709 273 : GEN a = a1, b = b1;
2710 273 : a1 = gmul2n(gadd(a,b),-1);
2711 273 : b1 = gsqrt(gmul(a,b), prec);
2712 273 : r = gsqrt(gdiv(gmul(a1,gaddgs(r, 1)),gadd(gmul(b, r), a )), prec);
2713 273 : t = gmul(r, t);
2714 : }
2715 49 : if (rotate) a1 = rotate>0 ? mulcxI(a1):mulcxmI(a1);
2716 49 : a1 = gmul(a1, b0);
2717 49 : t = gatan(gdiv(a1,t), prec);
2718 : /* send t to the fundamental domain if necessary */
2719 49 : if (gsigne(real_i(t))<0) t = gadd(t, mppi(prec));
2720 49 : return gerepileupto(av,gdiv(t,a1));
2721 : }
2722 :
2723 : static long
2724 49 : ser_cmp_expo(GEN A, GEN B)
2725 : {
2726 49 : long e = -(long)HIGHEXPOBIT, d = valp(B) - valp(A);
2727 49 : long i, la = lg(A), v = varn(B);
2728 9849 : for (i = 2; i < la; i++)
2729 : {
2730 9800 : GEN a = gel(A,i), b;
2731 : long ei;
2732 9800 : if (isexactzero(a)) continue;
2733 9800 : b = polcoef_i(B, i-2 + d, v);
2734 9800 : ei = gexpo(a);
2735 9800 : if (!isexactzero(b)) ei -= gexpo(b);
2736 9800 : e = maxss(e, ei);
2737 : }
2738 49 : return e;
2739 : }
2740 :
2741 : static GEN
2742 21 : ser_agm1(GEN y, long prec)
2743 : {
2744 21 : GEN a1 = y, b1 = gen_1;
2745 21 : long l = lg(y)-2, l2 = 6-prec2nbits(prec), eold = LONG_MAX;
2746 : for(;;)
2747 84 : {
2748 105 : GEN a = a1, p1;
2749 105 : a1 = gmul2n(gadd(a,b1),-1);
2750 105 : b1 = gsqrt(gmul(a,b1), prec);
2751 105 : p1 = gsub(b1,a1);
2752 105 : if (isinexactreal(p1))
2753 : {
2754 49 : long e = ser_cmp_expo(p1, b1);
2755 49 : if (e < l2 || e >= eold) break;
2756 42 : eold = e;
2757 : }
2758 56 : else if (valp(p1)-valp(b1) >= l || gequal0(p1)) break;
2759 : }
2760 21 : return a1;
2761 : }
2762 :
2763 : /* agm(1,x) */
2764 : static GEN
2765 17007 : agm1(GEN x, long prec)
2766 : {
2767 : GEN y;
2768 : pari_sp av;
2769 :
2770 17007 : if (gequal0(x)) return gcopy(x);
2771 17007 : switch(typ(x))
2772 : {
2773 28 : case t_INT:
2774 28 : if (!is_pm1(x)) break;
2775 21 : return (signe(x) > 0)? real_1(prec): real_0(prec);
2776 :
2777 11631 : case t_REAL: return signe(x) > 0? agm1r_abs(x): agm1cx(x, prec);
2778 :
2779 5208 : case t_COMPLEX:
2780 5208 : if (gequal0(gel(x,2))) return agm1(gel(x,1), prec);
2781 5201 : return agm1cx(x, prec);
2782 :
2783 14 : case t_PADIC:
2784 : {
2785 14 : GEN a1 = x, b1 = gen_1;
2786 14 : long l = precp(x);
2787 14 : av = avma;
2788 : for(;;)
2789 14 : {
2790 28 : GEN a = a1, p1;
2791 : long ep;
2792 28 : a1 = gmul2n(gadd(a,b1),-1);
2793 28 : a = gmul(a,b1);
2794 28 : b1 = Qp_sqrt(a); if (!b1) pari_err_SQRTN("Qp_sqrt",a);
2795 21 : p1 = gsub(b1,a1); ep = valp(p1)-valp(b1);
2796 21 : if (ep<=0) { b1 = gneg_i(b1); p1 = gsub(b1,a1); ep=valp(p1)-valp(b1); }
2797 21 : if (ep >= l || gequal0(p1)) return gerepilecopy(av,a1);
2798 : }
2799 : }
2800 :
2801 126 : default:
2802 126 : av = avma; if (!(y = toser_i(x))) break;
2803 21 : return gerepilecopy(av, ser_agm1(y, prec));
2804 : }
2805 112 : return trans_eval("agm",agm1,x,prec);
2806 : }
2807 :
2808 : GEN
2809 16853 : agm(GEN x, GEN y, long prec)
2810 : {
2811 : pari_sp av;
2812 16853 : if (is_matvec_t(typ(y)))
2813 : {
2814 14 : if (is_matvec_t(typ(x))) pari_err_TYPE2("agm",x,y);
2815 7 : swap(x, y);
2816 : }
2817 16846 : if (gequal0(y)) return gcopy(y);
2818 16846 : av = avma;
2819 16846 : return gerepileupto(av, gmul(y, agm1(gdiv(x,y), prec)));
2820 : }
2821 :
2822 : /* b2 != 0 */
2823 : static GEN
2824 35 : ellK_i(GEN b2, long prec)
2825 35 : { return gdiv(Pi2n(-1, prec), agm1(gsqrt(b2, prec), prec)); }
2826 : GEN
2827 28 : ellK(GEN k, long prec)
2828 : {
2829 28 : pari_sp av = avma;
2830 28 : GEN k2 = gsqr(k), b2 = gsubsg(1, k2);
2831 28 : if (gequal0(b2)) pari_err_DOMAIN("ellK", "k^2", "=", gen_1, k2);
2832 21 : return gerepileupto(av, ellK_i(b2, prec));
2833 : }
2834 :
2835 : static int
2836 84 : magm_gap(GEN a, GEN b, long L)
2837 : {
2838 84 : GEN d = gsub(b, a);
2839 84 : return !gequal0(d) && gexpo(d) - gexpo(b) >= L;
2840 : }
2841 :
2842 : /* http://www.ams.org/notices/201208/rtx120801094p.pdf
2843 : * An Eloquent Formula for the Perimeter of an Ellipse
2844 : * Semjon Adlaj, Notices of the AMS */
2845 : static GEN
2846 14 : magm(GEN a, GEN b, long prec)
2847 : {
2848 14 : long L = -prec2nbits(prec) + 16;
2849 14 : GEN c = gen_0;
2850 84 : while (magm_gap(a, b, L))
2851 : {
2852 70 : GEN u = gsqrt(gmul(gsub(a, c), gsub(b, c)), prec);
2853 70 : a = gmul2n(gadd(a, b), -1);
2854 70 : b = gadd(c, u); c = gsub(c, u);
2855 : }
2856 14 : return gmul2n(gadd(a, b), -1);
2857 : }
2858 :
2859 : GEN
2860 21 : ellE(GEN k, long prec)
2861 : {
2862 21 : pari_sp av = avma;
2863 21 : GEN b2 = gsubsg(1, gsqr(k));
2864 21 : if (gequal0(b2)) { set_avma(av); return real_1(prec); }
2865 14 : return gerepileupto(av, gmul(ellK_i(b2, prec), magm(gen_1, b2, prec)));
2866 : }
2867 :
2868 : /********************************************************************/
2869 : /** **/
2870 : /** LOG(X) **/
2871 : /** **/
2872 : /********************************************************************/
2873 : /* log(2) = 18*atanh(1/26)-2*atanh(1/4801)+8*atanh(1/8749)
2874 : * faster than 10*atanh(1/17)+4*atanh(13/499) for all precisions,
2875 : * and than Pi/2M(1,4/2^n) ~ n log(2) for bitprec at least up to 10^8 */
2876 : static GEN
2877 38445 : log2_split(long prec)
2878 : {
2879 38445 : GEN u = atanhuu(1, 26, prec);
2880 38444 : GEN v = atanhuu(1, 4801, prec);
2881 38448 : GEN w = atanhuu(1, 8749, prec);
2882 38442 : shiftr_inplace(v, 1); setsigne(v, -1);
2883 38446 : shiftr_inplace(w, 3);
2884 38447 : return addrr(mulur(18, u), addrr(v, w));
2885 : }
2886 : GEN
2887 19811327 : constlog2(long prec)
2888 : {
2889 : pari_sp av;
2890 : GEN tmp;
2891 19811327 : if (glog2 && realprec(glog2) >= prec) return glog2;
2892 :
2893 37964 : tmp = cgetr_block(prec);
2894 38445 : av = avma;
2895 38445 : affrr(log2_split(prec+EXTRAPRECWORD), tmp);
2896 38431 : swap_clone(&glog2,tmp);
2897 38451 : return gc_const(av,glog2);
2898 : }
2899 :
2900 : GEN
2901 19811141 : mplog2(long prec) { return rtor(constlog2(prec), prec); }
2902 :
2903 : /* dont check that q != 2^expo(q), done in logr_abs */
2904 : static GEN
2905 1379290 : logagmr_abs(GEN q)
2906 : {
2907 1379290 : long prec = realprec(q), e = expo(q), lim;
2908 1379290 : GEN z = cgetr(prec), y, Q, _4ovQ;
2909 1379292 : pari_sp av = avma;
2910 :
2911 1379292 : incrprec(prec);
2912 1379292 : lim = prec2nbits(prec) >> 1;
2913 1379292 : Q = rtor(q,prec);
2914 1379293 : shiftr_inplace(Q,lim-e); setsigne(Q,1);
2915 :
2916 1379292 : _4ovQ = invr(Q); shiftr_inplace(_4ovQ, 2); /* 4/Q */
2917 : /* Pi / 2agm(1, 4/Q) ~ log(Q), q = Q * 2^(e-lim) */
2918 1379290 : y = divrr(Pi2n(-1, prec), agm1r_abs(_4ovQ));
2919 1379290 : y = addrr(y, mulsr(e - lim, mplog2(prec)));
2920 1379291 : affrr_fixlg(y, z); return gc_const(av,z);
2921 : }
2922 :
2923 : /* sum_{k >= 0} y^(2k+1) / (2k+1), y close to 0 */
2924 : static GEN
2925 7430071 : logr_aux(GEN y)
2926 : {
2927 7430071 : long k, L = realprec(y); /* should be ~ l+1 - (k-2) */
2928 : /* log(x) = log(1+y) - log(1-y) = 2 sum_{k odd} y^k / k
2929 : * Truncate the sum at k = 2n+1, the remainder is
2930 : * 2 sum_{k >= 2n+3} y^k / k < 2y^(2n+3) / (2n+3)(1-y) < y^(2n+3)
2931 : * We want y^(2n+3) < y 2^(-prec2nbits(L)), hence
2932 : * n+1 > -prec2nbits(L) /-log_2(y^2) */
2933 7430071 : double d = -2*dbllog2r(y); /* ~ -log_2(y^2) */
2934 7430135 : k = (long)(2*(prec2nbits(L) / d));
2935 7430198 : k |= 1;
2936 7430198 : if (k >= 3)
2937 : {
2938 7411590 : GEN T, S = cgetr(L), y2 = sqrr(y), unr = real_1(L);
2939 7411575 : pari_sp av = avma;
2940 7411575 : long s = 0, incs = (long)d, l1 = nbits2prec((long)d);
2941 7411631 : setprec(S, l1);
2942 7411586 : setprec(unr,l1); affrr(divru(unr,k), S);
2943 7411270 : for (k -= 2;; k -= 2) /* k = 2n+1, ..., 1 */
2944 : { /* S = y^(2n+1-k)/(2n+1) + ... + 1 / k */
2945 132922321 : setprec(y2, l1); T = mulrr(S,y2);
2946 132900866 : if (k == 1) break;
2947 :
2948 125489574 : l1 += dvmdsBIL(s + incs, &s); if (l1>L) l1=L;
2949 125446312 : setprec(S, l1);
2950 125508772 : setprec(unr,l1);
2951 125532005 : affrr(addrr(divru(unr, k), T), S); set_avma(av);
2952 : }
2953 : /* k = 1 special-cased for eficiency */
2954 7411292 : y = mulrr(y, addsr(1,T)); /* = log(X)/2 */
2955 : }
2956 7430079 : return y;
2957 : }
2958 : /*return log(|x|), assuming x != 0 */
2959 : GEN
2960 9268728 : logr_abs(GEN X)
2961 : {
2962 9268728 : long EX, L, m, k, a, b, l = realprec(X);
2963 : GEN z, x, y;
2964 : ulong u;
2965 : double d;
2966 :
2967 : /* Assuming 1 < x < 2, we want delta = x-1, 1-x/2, 1-1/x, or 2/x-1 small.
2968 : * We have 2/x-1 > 1-x/2, 1-1/x < x-1. So one should be choosing between
2969 : * 1-1/x and 1-x/2 ( crossover sqrt(2), worse ~ 0.29 ). To avoid an inverse,
2970 : * we choose between x-1 and 1-x/2 ( crossover 4/3, worse ~ 0.33 ) */
2971 9268728 : EX = expo(X);
2972 9268728 : u = uel(X,2);
2973 9268728 : k = 2;
2974 9268728 : if (u > (~0UL / 3) * 2) { /* choose 1-x/2 */
2975 5206192 : EX++; u = ~u;
2976 5265928 : while (!u && ++k < l) { u = uel(X,k); u = ~u; }
2977 : } else { /* choose x - 1 */
2978 4062536 : u &= ~HIGHBIT; /* u - HIGHBIT, assuming HIGHBIT set */
2979 5212015 : while (!u && ++k < l) u = uel(X,k);
2980 : }
2981 9268728 : if (k == l) return EX? mulsr(EX, mplog2(l)): real_0(l);
2982 8809716 : a = prec2nbits(k) + bfffo(u); /* ~ -log2 |1-x| */
2983 8809854 : L = l+1;
2984 8809854 : b = prec2nbits(L - (k-2)); /* take loss of accuracy into account */
2985 15730722 : if (b > 24*a*log2(L) &&
2986 8300009 : prec2nbits(l) > prec2nbits(LOGAGM_LIMIT)) return logagmr_abs(X);
2987 :
2988 7430713 : z = cgetr(EX? l: l - (k-2));
2989 :
2990 : /* Multiplication is quadratic in this range (l is small, otherwise we
2991 : * use AGM). Set Y = x^(1/2^m), y = (Y - 1) / (Y + 1) and compute truncated
2992 : * series sum y^(2k+1)/(2k+1): the costs is less than
2993 : * m b^2 + sum_{k <= n} ((2k+1) e + BITS_IN_LONG)^2
2994 : * bit operations with |x-1| < 2^(1-a), |Y| < 2^(1-e) , m = e-a and b bits of
2995 : * accuracy needed (+ BITS_IN_LONG since bit accuracies increase by
2996 : * increments of BITS_IN_LONG), so
2997 : * 4n^3/3 e^2 + n^2 2e BITS_IN_LONG+ n BITS_IN_LONG ~ m b^2, with n ~ b/2e
2998 : * or b/6e + BITS_IN_LONG/2e + BITS_IN_LONG/2be ~ m
2999 : * B := (b / 6 + BITS_IN_LONG/2 + BITS_IN_LONG^2 / 2b) ~ m(m+a)
3000 : * m = min( -a/2 + sqrt(a^2/4 + B), b - a )
3001 : * NB: e ~ (b/6)^(1/2) as b -> oo
3002 : * Instead of the above pessimistic estimate for the cost of the sum, use
3003 : * optimistic estimate (BITS_IN_LONG -> 0) */
3004 7430602 : d = -a/2.; m = (long)(d + sqrt(d*d + b/6)); /* >= 0 */
3005 :
3006 7430602 : if (m > b-a) m = b-a;
3007 7430602 : if (m < 0.2*a) m = 0; else L += nbits2extraprec(m);
3008 7430616 : x = rtor(X,L);
3009 7430707 : setsigne(x,1); shiftr_inplace(x,-EX);
3010 : /* 2/3 < x < 4/3 */
3011 43480035 : for (k=1; k<=m; k++) x = sqrtr_abs(x);
3012 :
3013 7430539 : y = divrr(subrs(x,1), addrs(x,1)); /* = (x-1) / (x+1), close to 0 */
3014 7430072 : y = logr_aux(y); /* log(1+y) - log(1-y) = log(x) */
3015 7430001 : shiftr_inplace(y, m + 1);
3016 7429857 : if (EX) y = addrr(y, mulsr(EX, mplog2(l+1)));
3017 7429884 : affrr_fixlg(y, z); return gc_const((pari_sp)z, z);
3018 : }
3019 :
3020 : /* assume Im(q) != 0 and precision(q) >= prec. Compute log(q) with accuracy
3021 : * prec [disregard input accuracy] */
3022 : GEN
3023 502995 : logagmcx(GEN q, long prec)
3024 : {
3025 502995 : GEN z = cgetc(prec), y, Q, a, b;
3026 : long lim, e, ea, eb;
3027 502995 : pari_sp av = avma;
3028 502995 : int neg = 0;
3029 :
3030 502995 : incrprec(prec);
3031 502995 : if (gsigne(gel(q,1)) < 0) { q = gneg(q); neg = 1; }
3032 502995 : lim = prec2nbits(prec) >> 1;
3033 502995 : Q = gtofp(q, prec);
3034 502995 : a = gel(Q,1);
3035 502995 : b = gel(Q,2);
3036 502995 : if (gequal0(a)) {
3037 0 : affrr_fixlg(logr_abs(b), gel(z,1));
3038 0 : y = Pi2n(-1, prec);
3039 0 : if (signe(b) < 0) setsigne(y, -1);
3040 0 : affrr_fixlg(y, gel(z,2)); return gc_const(av,z);
3041 : }
3042 502995 : ea = expo(a);
3043 502995 : eb = expo(b);
3044 502995 : e = ea <= eb ? lim - eb : lim - ea;
3045 502995 : shiftr_inplace(a, e);
3046 502995 : shiftr_inplace(b, e);
3047 :
3048 : /* Pi / 2agm(1, 4/Q) ~ log(Q), q = Q * 2^e */
3049 502995 : y = gdiv(Pi2n(-1, prec), agm1cx( gdivsg(4, Q), prec ));
3050 502995 : a = gel(y,1);
3051 502995 : b = gel(y,2);
3052 502995 : a = addrr(a, mulsr(-e, mplog2(prec)));
3053 502995 : if (realprec(a) <= LOWDEFAULTPREC) a = real_0_bit(expo(a));
3054 752194 : if (neg) b = gsigne(b) <= 0? gadd(b, mppi(prec))
3055 249199 : : gsub(b, mppi(prec));
3056 502995 : affrr_fixlg(a, gel(z,1));
3057 502995 : affrr_fixlg(b, gel(z,2)); return gc_const(av,z);
3058 : }
3059 :
3060 : GEN
3061 142959 : mplog(GEN x)
3062 : {
3063 142959 : if (signe(x)<=0) pari_err_DOMAIN("mplog", "argument", "<=", gen_0, x);
3064 142959 : return logr_abs(x);
3065 : }
3066 :
3067 : /* pe = p^e, p prime, 0 < x < pe a t_INT coprime to p. Return the (p-1)-th
3068 : * root of 1 in (Z/pe)^* congruent to x mod p, resp x mod 4 if p = 2.
3069 : * Simplified form of Zp_sqrtnlift: 1/(p-1) is trivial to compute */
3070 : GEN
3071 9932 : Zp_teichmuller(GEN x, GEN p, long e, GEN pe)
3072 : {
3073 : GEN q, z, p1;
3074 : pari_sp av;
3075 : ulong mask;
3076 9932 : if (absequaliu(p,2)) return (mod4(x) & 2)? subiu(pe,1): gen_1;
3077 9393 : if (e == 1) return icopy(x);
3078 9393 : av = avma;
3079 9393 : p1 = subiu(p, 1);
3080 9393 : mask = quadratic_prec_mask(e);
3081 9393 : q = p; z = remii(x, p);
3082 31189 : while (mask > 1)
3083 : { /* Newton iteration solving z^{1 - p} = 1, z = x (mod p) */
3084 21796 : GEN w, t, qold = q;
3085 21796 : if (mask <= 3) /* last iteration */
3086 9393 : q = pe;
3087 : else
3088 : {
3089 12403 : q = sqri(q);
3090 12403 : if (mask & 1) q = diviiexact(q, p);
3091 : }
3092 21796 : mask >>= 1;
3093 : /* q <= qold^2 */
3094 21796 : if (lgefint(q) == 3)
3095 : {
3096 21648 : ulong Z = uel(z,2), Q = uel(q,2), P1 = uel(p1,2);
3097 21648 : ulong W = (Q-1) / P1; /* -1/(p-1) + O(qold) */
3098 21648 : ulong T = Fl_mul(W, Fl_powu(Z,P1,Q) - 1, Q);
3099 21648 : Z = Fl_mul(Z, 1 + T, Q);
3100 21648 : z = utoi(Z);
3101 : }
3102 : else
3103 : {
3104 148 : w = diviiexact(subiu(qold,1),p1); /* -1/(p-1) + O(qold) */
3105 148 : t = Fp_mul(w, subiu(Fp_pow(z,p1,q), 1), q);
3106 148 : z = Fp_mul(z, addui(1,t), q);
3107 : }
3108 : }
3109 9393 : return gerepileuptoint(av, z);
3110 : }
3111 :
3112 : GEN
3113 1225 : teichmullerinit(long p, long n)
3114 : {
3115 : GEN t, pn, g, v;
3116 : ulong gp, tp;
3117 : long a, m;
3118 :
3119 1225 : if (p == 2) return mkvec(gen_1);
3120 1225 : if (!uisprime(p)) pari_err_PRIME("teichmullerinit",utoipos(p));
3121 :
3122 1225 : m = p >> 1; /* (p-1)/2 */
3123 1225 : tp= gp= pgener_Fl(p); /* order (p-1), gp^m = -1 */
3124 1225 : pn = powuu(p, n);
3125 1225 : v = cgetg(p, t_VEC);
3126 1225 : t = g = Zp_teichmuller(utoipos(gp), utoipos(p), n, pn);
3127 1225 : gel(v, 1) = gen_1;
3128 1225 : gel(v, p-1) = subiu(pn,1);
3129 3031 : for (a = 1; a < m; a++)
3130 : {
3131 1806 : gel(v, tp) = t;
3132 1806 : gel(v, p - tp) = Fp_neg(t, pn); /* g^(m+a) = -g^a */
3133 1806 : if (a < m-1)
3134 : {
3135 1029 : t = Fp_mul(t, g, pn); /* g^(a+1) */
3136 1029 : tp = Fl_mul(tp, gp, p); /* t mod p */
3137 : }
3138 : }
3139 1225 : return v;
3140 : }
3141 :
3142 : /* tab from teichmullerinit or NULL */
3143 : GEN
3144 4977 : teichmuller(GEN x, GEN tab)
3145 : {
3146 : GEN p, q, y, z;
3147 4977 : long n, tx = typ(x);
3148 :
3149 4977 : if (!tab)
3150 : {
3151 4865 : if (tx == t_VEC && lg(x) == 3)
3152 : {
3153 7 : p = gel(x,1);
3154 7 : q = gel(x,2);
3155 7 : if (typ(p) == t_INT && typ(q) == t_INT)
3156 7 : return teichmullerinit(itos(p), itos(q));
3157 : }
3158 : }
3159 112 : else if (typ(tab) != t_VEC) pari_err_TYPE("teichmuller",tab);
3160 4970 : if (tx!=t_PADIC) pari_err_TYPE("teichmuller",x);
3161 4970 : z = gel(x,4);
3162 4970 : if (!signe(z)) return gcopy(x);
3163 4970 : p = gel(x,2);
3164 4970 : q = gel(x,3);
3165 4970 : n = precp(x);
3166 4970 : y = cgetg(5,t_PADIC);
3167 4970 : y[1] = evalprecp(n) | _evalvalp(0);
3168 4970 : gel(y,2) = icopy(p);
3169 4970 : gel(y,3) = icopy(q);
3170 4970 : if (tab)
3171 : {
3172 112 : ulong pp = itou_or_0(p);
3173 112 : if (lg(tab) != (long)pp) pari_err_TYPE("teichmuller",tab);
3174 112 : z = gel(tab, umodiu(z, pp));
3175 112 : if (typ(z) != t_INT) pari_err_TYPE("teichmuller",tab);
3176 112 : z = remii(z, q);
3177 : }
3178 : else
3179 4858 : z = Zp_teichmuller(z, p, n, q);
3180 4970 : gel(y,4) = z;
3181 4970 : return y;
3182 : }
3183 : GEN
3184 4739 : teich(GEN x) { return teichmuller(x, NULL); }
3185 :
3186 : GEN
3187 13629399 : glog(GEN x, long prec)
3188 : {
3189 : pari_sp av, tetpil;
3190 : GEN y, p1;
3191 : long l;
3192 :
3193 13629399 : switch(typ(x))
3194 : {
3195 7451383 : case t_REAL:
3196 7451383 : if (signe(x) >= 0)
3197 : {
3198 6428780 : if (!signe(x)) pari_err_DOMAIN("log", "argument", "=", gen_0, x);
3199 6428766 : return logr_abs(x);
3200 : }
3201 1022603 : retmkcomplex(logr_abs(x), mppi(realprec(x)));
3202 :
3203 518918 : case t_FRAC:
3204 : {
3205 : GEN a, b;
3206 : long e1, e2;
3207 518918 : av = avma;
3208 518918 : a = gel(x,1);
3209 518918 : b = gel(x,2);
3210 518918 : e1 = expi(subii(a,b)); e2 = expi(b);
3211 518918 : if (e2 > e1) prec += nbits2nlong(e2 - e1);
3212 518918 : x = fractor(x, prec);
3213 518918 : return gerepileupto(av, glog(x, prec));
3214 : }
3215 4075596 : case t_COMPLEX:
3216 4075596 : if (ismpzero(gel(x,2))) return glog(gel(x,1), prec);
3217 4065371 : l = precision(x); if (l > prec) prec = l;
3218 4065388 : if (ismpzero(gel(x,1)))
3219 : {
3220 71830 : GEN a = gel(x,2), b;
3221 71830 : av = avma; b = Pi2n(-1,prec);
3222 71823 : if (gsigne(a) < 0) { setsigne(b, -1); a = gabs(a,prec); }
3223 71823 : a = isint1(a) ? gen_0: glog(a,prec);
3224 71823 : return gerepilecopy(av, mkcomplex(a, b));
3225 : }
3226 3993561 : if (prec >= LOGAGMCX_LIMIT) return logagmcx(x, prec);
3227 3490762 : y = cgetg(3,t_COMPLEX);
3228 3490755 : gel(y,2) = garg(x,prec);
3229 3490806 : av = avma; p1 = glog(cxnorm(x),prec); tetpil = avma;
3230 3490807 : gel(y,1) = gerepile(av,tetpil,gmul2n(p1,-1)); return y;
3231 :
3232 322 : case t_PADIC: return Qp_log(x);
3233 1583180 : default:
3234 1583180 : av = avma; if (!(y = toser_i(x))) break;
3235 139 : if (!signe(y)) pari_err_DOMAIN("log", "argument", "=", gen_0, x);
3236 139 : if (valp(y)) pari_err_DOMAIN("log", "series valuation", "!=", gen_0, x);
3237 132 : p1 = integser(gdiv(derivser(y), y)); /* log(y)' = y'/y */
3238 133 : if (!gequal1(gel(y,2))) p1 = gadd(p1, glog(gel(y,2),prec));
3239 133 : return gerepileupto(av, p1);
3240 : }
3241 1583433 : return trans_eval("log",glog,x,prec);
3242 : }
3243 :
3244 : static GEN
3245 63 : mplog1p(GEN x)
3246 : {
3247 : long ex, a, b, l, L;
3248 63 : if (!signe(x)) return rcopy(x);
3249 63 : ex = expo(x); if (ex >= -3) return glog(addrs(x,1), 0);
3250 42 : a = -ex;
3251 42 : b = realprec(x); L = b+1;
3252 42 : if (b > a*log2(L) && prec2nbits(b) > prec2nbits(LOGAGM_LIMIT))
3253 : {
3254 0 : x = addrs(x,1); l = b + nbits2extraprec(a);
3255 0 : if (realprec(x) < l) x = rtor(x,l);
3256 0 : return logagmr_abs(x);
3257 : }
3258 42 : x = rtor(x, L);
3259 42 : x = logr_aux(divrr(x, addrs(x,2)));
3260 42 : if (realprec(x) > b) fixlg(x, b);
3261 42 : shiftr_inplace(x,1); return x;
3262 : }
3263 :
3264 : static GEN log1p_i(GEN x, long prec);
3265 : static GEN
3266 14 : cxlog1p(GEN x, long prec)
3267 : {
3268 : pari_sp av;
3269 14 : GEN z, a, b = gel(x,2);
3270 : long l;
3271 14 : if (ismpzero(b)) return log1p_i(gel(x,1), prec);
3272 14 : l = precision(x); if (l > prec) prec = l;
3273 14 : if (prec >= LOGAGMCX_LIMIT) return logagmcx(gaddgs(x,1), prec);
3274 14 : a = gel(x,1);
3275 14 : z = cgetg(3,t_COMPLEX); av = avma;
3276 14 : a = gadd(gadd(gmul2n(a,1), gsqr(a)), gsqr(b));
3277 14 : a = log1p_i(a, prec); shiftr_inplace(a,-1);
3278 14 : gel(z,1) = gerepileupto(av, a);
3279 14 : gel(z,2) = garg(gaddgs(x,1),prec); return z;
3280 : }
3281 : static GEN
3282 133 : log1p_i(GEN x, long prec)
3283 : {
3284 133 : switch(typ(x))
3285 : {
3286 63 : case t_REAL: return mplog1p(x);
3287 14 : case t_COMPLEX: return cxlog1p(x, prec);
3288 7 : case t_PADIC: return Qp_log(gaddgs(x,1));
3289 49 : default:
3290 : {
3291 : long ey;
3292 : GEN y;
3293 49 : if (!(y = toser_i(x))) break;
3294 21 : ey = valp(y);
3295 21 : if (ey < 0) pari_err_DOMAIN("log1p","valuation", "<", gen_0, x);
3296 21 : if (gequal0(y)) return gcopy(y);
3297 14 : if (ey)
3298 7 : return glog(gaddgs(y,1),prec);
3299 : else
3300 : {
3301 7 : GEN a = gel(y,2), a1 = gaddgs(a,1);
3302 7 : y = gdiv(y, a1); gel(y,2) = gen_1;
3303 7 : return gadd(glog1p(a,prec), glog(y, prec));
3304 : }
3305 : }
3306 : }
3307 28 : return trans_eval("log1p",glog1p,x,prec);
3308 : }
3309 : GEN
3310 119 : glog1p(GEN x, long prec)
3311 : {
3312 119 : pari_sp av = avma;
3313 119 : return gerepileupto(av, log1p_i(x, prec));
3314 : }
3315 : /********************************************************************/
3316 : /** **/
3317 : /** SINE, COSINE **/
3318 : /** **/
3319 : /********************************************************************/
3320 :
3321 : /* Reduce x0 mod Pi/2 to x in [-Pi/4, Pi/4]. Return cos(x)-1 */
3322 : static GEN
3323 9494781 : mpcosm1(GEN x, long *ptmod8)
3324 : {
3325 9494781 : long a = expo(x), l = realprec(x), b, L, i, n, m, B;
3326 : GEN y, u, x2;
3327 : double d;
3328 :
3329 9494781 : n = 0;
3330 9494781 : if (a >= 0)
3331 : {
3332 : long p;
3333 : GEN q;
3334 7698319 : if (a > 30)
3335 : {
3336 7 : GEN z, P = Pi2n(-2, nbits2prec(a + 32));
3337 7 : z = addrr(x,P); /* = x + Pi/4 */
3338 7 : if (expo(z) >= bit_prec(z) + 3) pari_err_PREC("mpcosm1");
3339 7 : shiftr_inplace(P, 1);
3340 7 : q = floorr(divrr(z, P)); /* round ( x / (Pi/2) ) */
3341 7 : p = l+EXTRAPRECWORD; x = rtor(x,p);
3342 : } else {
3343 7698312 : q = stoi((long)floor(rtodbl(x) / (M_PI/2) + 0.5));
3344 7698274 : p = l;
3345 : }
3346 7698498 : if (signe(q))
3347 : {
3348 7698285 : GEN y = subrr(x, mulir(q, Pi2n(-1,p))); /* x mod Pi/2 */
3349 7697949 : long b = expo(y);
3350 7697949 : if (a - b < 7) x = y;
3351 : else
3352 : {
3353 3841679 : p += nbits2extraprec(a-b); x = rtor(x, p);
3354 3841701 : x = subrr(x, mulir(q, Pi2n(-1,p)));
3355 : }
3356 7697891 : a = b;
3357 7697891 : if (!signe(x) && a >= 0) pari_err_PREC("mpcosm1");
3358 7697891 : n = Mod4(q);
3359 : }
3360 : }
3361 : /* a < 0 */
3362 9494682 : b = signe(x); *ptmod8 = (b < 0)? 4 + n: n;
3363 9494682 : if (!b) return real_0_bit(expo(x)*2 - 1);
3364 :
3365 9494682 : b = prec2nbits(l);
3366 9494650 : if (b + 2*a <= 0) {
3367 732342 : y = sqrr(x); shiftr_inplace(y, -1); setsigne(y, -1);
3368 732336 : return y;
3369 : }
3370 :
3371 8762308 : y = cgetr(l);
3372 8762285 : B = b/6 + BITS_IN_LONG/2 + (BITS_IN_LONG*BITS_IN_LONG/2)/ b;
3373 8762285 : d = a/2.; m = (long)(d + sqrt(d*d + B)); /* >= 0 ,*/
3374 8762285 : if (m < (-a) * 0.1) m = 0; /* not worth it */
3375 8762285 : L = l + nbits2extraprec(m);
3376 :
3377 8762344 : b += m;
3378 8762344 : d = 2.0 * (m-dbllog2r(x)-1/M_LN2); /* ~ 2( - log_2 Y - 1/log(2) ) */
3379 8762374 : n = (long)(b / d);
3380 8762374 : if (n > 1)
3381 8704560 : n = (long)(b / (d + log2((double)n+1))); /* log~constant in small ranges */
3382 18599099 : while (n*(d+log2((double)n+1)) < b) n++; /* expect few corrections */
3383 :
3384 : /* Multiplication is quadratic in this range (l is small, otherwise we
3385 : * use logAGM + Newton). Set Y = 2^(-e-a) x, compute truncated series
3386 : * sum Y^2k/(2k)!: this costs roughly
3387 : * m b^2 + sum_{k <= n} (2k e + BITS_IN_LONG)^2
3388 : * ~ (b/2e) b^2 / 3 + m b^2
3389 : * bit operations with n ~ b/2e, |x| < 2^(1+a), |Y| < 2^(1-e) , m = e+a and
3390 : * b bits of accuracy needed, so
3391 : * B := (b / 6 + BITS_IN_LONG/2 + BITS_IN_LONG^2 / 2b) ~ m(m-a)
3392 : * we want b ~ 6 m (m-a) or m~b+a hence
3393 : * m = min( a/2 + sqrt(a^2/4 + b/6), b/2 + a )
3394 : * NB: e ~ (b/6)^(1/2) or b/2.
3395 : *
3396 : * Truncate the sum at k = n (>= 1), the remainder is
3397 : * < sum_{k >= n+1} Y^2k / 2k! < Y^(2n+2) / (2n+2)!(1-Y^2) < Y^(2n+2)/(2n+1)!
3398 : * We want ... <= Y^2 2^-b, hence -2n log_2 |Y| + log_2 (2n+1)! >= b
3399 : * log n! ~ (n + 1/2) log(n+1) - (n+1) + log(2Pi)/2,
3400 : * error bounded by 1/6(n+1) <= 1/12. Finally, we want
3401 : * 2n (-1/log(2) - log_2 |Y| + log_2(2n+2)) >= b */
3402 8762374 : x = rtor(x, L); shiftr_inplace(x, -m); setsigne(x, 1);
3403 8762344 : x2 = sqrr(x);
3404 8762036 : if (n == 1) { u = x2; shiftr_inplace(u, -1); setsigne(u, -1); } /*-Y^2/2*/
3405 : else
3406 : {
3407 8762036 : GEN un = real_1(L);
3408 : pari_sp av;
3409 8762151 : long s = 0, l1 = nbits2prec((long)(d + n + 16));
3410 :
3411 8762129 : u = cgetr(L); av = avma;
3412 113909032 : for (i = n; i >= 2; i--)
3413 : {
3414 : GEN t;
3415 105146765 : setprec(x2,l1); t = divrunextu(x2, 2*i-1);
3416 105143576 : l1 += dvmdsBIL(s - expo(t), &s); if (l1 > L) l1 = L;
3417 105144413 : if (i != n) t = mulrr(t,u);
3418 105143014 : setprec(un,l1); t = addrr_sign(un,1, t,-signe(t));
3419 105127834 : setprec(u,l1); affrr(t,u); set_avma(av);
3420 : }
3421 8762267 : shiftr_inplace(u, -1); togglesign(u); /* u := -u/2 */
3422 8762174 : setprec(x2,L); u = mulrr(x2,u);
3423 : }
3424 : /* Now u = sum {1<= i <=n} (-1)^i x^(2i) / (2i)! ~ cos(x) - 1 */
3425 75389066 : for (i = 1; i <= m; i++)
3426 : { /* u = cos(x)-1 <- cos(2x)-1 = 2cos(x)^2 - 2 = 4u + 2u^2*/
3427 66629544 : GEN q = sqrr(u);
3428 66637262 : shiftr_inplace(u, 1); u = addrr(u, q);
3429 66628593 : shiftr_inplace(u, 1);
3430 66626433 : if ((i & 31) == 0) u = gerepileuptoleaf((pari_sp)y, u);
3431 : }
3432 8759522 : affrr_fixlg(u, y); return y;
3433 : }
3434 :
3435 : /* sqrt (|1 - (1+x)^2|) = sqrt(|x*(x+2)|). Sends cos(x)-1 to |sin(x)| */
3436 : static GEN
3437 7788842 : mpaut(GEN x)
3438 : {
3439 7788842 : GEN t = mulrr(x, addsr(2,x)); /* != 0 */
3440 7788858 : if (!signe(t)) return real_0_bit(expo(t) >> 1);
3441 7788858 : return sqrtr_abs(t);
3442 : }
3443 :
3444 : /********************************************************************/
3445 : /** COSINE **/
3446 : /********************************************************************/
3447 :
3448 : GEN
3449 2744997 : mpcos(GEN x)
3450 : {
3451 : long mod8;
3452 : pari_sp av;
3453 : GEN y, z;
3454 :
3455 2744997 : if (!signe(x)) {
3456 75 : long l = nbits2prec(-expo(x));
3457 75 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
3458 75 : return real_1(l);
3459 : }
3460 2744922 : av = avma; z = mpcosm1(x,&mod8);
3461 2744914 : switch(mod8)
3462 : {
3463 759990 : case 0: case 4: y = addsr(1,z); break;
3464 688609 : case 1: case 7: y = mpaut(z); togglesign(y); break;
3465 682287 : case 2: case 6: y = subsr(-1,z); break;
3466 614028 : default: y = mpaut(z); break; /* case 3: case 5: */
3467 : }
3468 2744941 : return gerepileuptoleaf(av, y);
3469 : }
3470 :
3471 : /* convert INT or FRAC to REAL, which is later reduced mod 2Pi : avoid
3472 : * cancellation */
3473 : static GEN
3474 13253 : tofp_safe(GEN x, long prec)
3475 : {
3476 13253 : return (typ(x) == t_INT || gexpo(x) > 0)? gadd(x, real_0(prec))
3477 26501 : : fractor(x, prec);
3478 : }
3479 :
3480 : GEN
3481 154770 : gcos(GEN x, long prec)
3482 : {
3483 : pari_sp av;
3484 : GEN a, b, u, v, y, u1, v1;
3485 : long i;
3486 :
3487 154770 : switch(typ(x))
3488 : {
3489 153502 : case t_REAL: return mpcos(x);
3490 28 : case t_COMPLEX:
3491 28 : a = gel(x,1);
3492 28 : b = gel(x,2);
3493 28 : if (isintzero(a)) return gcosh(b, prec);
3494 14 : i = precision(x); if (i) prec = i;
3495 14 : y = cgetc(prec); av = avma;
3496 14 : if (typ(b) != t_REAL) b = gtofp(b, prec);
3497 14 : mpsinhcosh(b, &u1, &v1); u1 = mpneg(u1);
3498 14 : if (typ(a) != t_REAL) a = gtofp(a, prec);
3499 14 : mpsincos(a, &u, &v);
3500 14 : affrr_fixlg(gmul(v1,v), gel(y,1));
3501 14 : affrr_fixlg(gmul(u1,u), gel(y,2)); return gc_const(av,y);
3502 :
3503 1156 : case t_INT: case t_FRAC:
3504 1156 : y = cgetr(prec); av = avma;
3505 1156 : affrr_fixlg(mpcos(tofp_safe(x,prec)), y); return gc_const(av,y);
3506 :
3507 49 : case t_PADIC: y = cos_p(x);
3508 49 : if (!y) pari_err_DOMAIN("gcos(t_PADIC)","argument","",gen_0,x);
3509 42 : return y;
3510 :
3511 35 : default:
3512 35 : av = avma; if (!(y = toser_i(x))) break;
3513 28 : if (gequal0(y)) return gerepileupto(av, gaddsg(1,y));
3514 28 : if (valp(y) < 0)
3515 7 : pari_err_DOMAIN("cos","valuation", "<", gen_0, x);
3516 21 : gsincos(y,&u,&v,prec);
3517 21 : return gerepilecopy(av,v);
3518 : }
3519 7 : return trans_eval("cos",gcos,x,prec);
3520 : }
3521 : /********************************************************************/
3522 : /** SINE **/
3523 : /********************************************************************/
3524 :
3525 : GEN
3526 849471 : mpsin(GEN x)
3527 : {
3528 : long mod8;
3529 : pari_sp av;
3530 : GEN y, z;
3531 :
3532 849471 : if (!signe(x)) return real_0_bit(expo(x));
3533 849262 : av = avma; z = mpcosm1(x,&mod8);
3534 849239 : switch(mod8)
3535 : {
3536 316662 : case 0: case 6: y = mpaut(z); break;
3537 134023 : case 1: case 5: y = addsr(1,z); break;
3538 269091 : case 2: case 4: y = mpaut(z); togglesign(y); break;
3539 129463 : default: y = subsr(-1,z); break; /* case 3: case 7: */
3540 : }
3541 849276 : return gerepileuptoleaf(av, y);
3542 : }
3543 :
3544 : GEN
3545 882401 : gsin(GEN x, long prec)
3546 : {
3547 : pari_sp av;
3548 : GEN a, b, u, v, y, v1, u1;
3549 : long i;
3550 :
3551 882401 : switch(typ(x))
3552 : {
3553 844294 : case t_REAL: return mpsin(x);
3554 32711 : case t_COMPLEX:
3555 32711 : a = gel(x,1);
3556 32711 : b = gel(x,2);
3557 32711 : if (isintzero(a)) retmkcomplex(gen_0,gsinh(b,prec));
3558 17402 : i = precision(x); if (i) prec = i;
3559 17402 : y = cgetc(prec); av = avma;
3560 17402 : if (typ(b) != t_REAL) b = gtofp(b, prec);
3561 17402 : mpsinhcosh(b, &u1, &v1);
3562 17402 : if (typ(a) != t_REAL) a = gtofp(a, prec);
3563 17402 : mpsincos(a, &u, &v);
3564 17402 : affrr_fixlg(gmul(v1,u), gel(y,1));
3565 17402 : affrr_fixlg(gmul(u1,v), gel(y,2)); return gc_const(av,y);
3566 :
3567 5118 : case t_INT: case t_FRAC:
3568 5118 : y = cgetr(prec); av = avma;
3569 5118 : affrr_fixlg(mpsin(tofp_safe(x,prec)), y); return gc_const(av,y);
3570 :
3571 49 : case t_PADIC: y = sin_p(x);
3572 49 : if (!y) pari_err_DOMAIN("gsin(t_PADIC)","argument","",gen_0,x);
3573 42 : return y;
3574 :
3575 229 : default:
3576 229 : av = avma; if (!(y = toser_i(x))) break;
3577 224 : if (gequal0(y)) return gerepilecopy(av, y);
3578 224 : if (valp(y) < 0)
3579 7 : pari_err_DOMAIN("sin","valuation", "<", gen_0, x);
3580 217 : gsincos(y,&u,&v,prec);
3581 217 : return gerepilecopy(av,u);
3582 : }
3583 7 : return trans_eval("sin",gsin,x,prec);
3584 : }
3585 : /********************************************************************/
3586 : /** SINE, COSINE together **/
3587 : /********************************************************************/
3588 :
3589 : void
3590 5898083 : mpsincos(GEN x, GEN *s, GEN *c)
3591 : {
3592 : long mod8;
3593 : pari_sp av, tetpil;
3594 : GEN z, *gptr[2];
3595 :
3596 5898083 : if (!signe(x))
3597 : {
3598 3276 : long e = expo(x);
3599 3276 : *s = real_0_bit(e);
3600 3276 : *c = e >= 0? real_0_bit(e): real_1_bit(-e);
3601 3276 : return;
3602 : }
3603 :
3604 5894807 : av = avma; z = mpcosm1(x, &mod8); tetpil = avma;
3605 5894800 : switch(mod8)
3606 : {
3607 1494454 : case 0: *c = addsr( 1,z); *s = mpaut(z); break;
3608 433453 : case 1: *s = addsr( 1,z); *c = mpaut(z); togglesign(*c); break;
3609 754183 : case 2: *c = subsr(-1,z); *s = mpaut(z); togglesign(*s); break;
3610 435340 : case 3: *s = subsr(-1,z); *c = mpaut(z); break;
3611 893790 : case 4: *c = addsr( 1,z); *s = mpaut(z); togglesign(*s); break;
3612 413779 : case 5: *s = addsr( 1,z); *c = mpaut(z); break;
3613 1043696 : case 6: *c = subsr(-1,z); *s = mpaut(z); break;
3614 426190 : case 7: *s = subsr(-1,z); *c = mpaut(z); togglesign(*c); break;
3615 : }
3616 5894845 : gptr[0] = s; gptr[1] = c; gerepilemanysp(av,tetpil,gptr,2);
3617 : }
3618 :
3619 : /* SINE and COSINE - 1 */
3620 : void
3621 5838 : mpsincosm1(GEN x, GEN *s, GEN *c)
3622 : {
3623 : long mod8;
3624 : pari_sp av, tetpil;
3625 : GEN z, *gptr[2];
3626 :
3627 5838 : if (!signe(x))
3628 : {
3629 0 : long e = expo(x);
3630 0 : *s = real_0_bit(e);
3631 0 : *c = real_0_bit(2*e-1);
3632 0 : return;
3633 : }
3634 5838 : av = avma; z = mpcosm1(x,&mod8); tetpil = avma;
3635 5838 : switch(mod8)
3636 : {
3637 4872 : case 0: *c = rcopy(z); *s = mpaut(z); break;
3638 42 : case 1: *s = addsr(1,z); *c = addrs(mpaut(z),1); togglesign(*c); break;
3639 0 : case 2: *c = subsr(-2,z); *s = mpaut(z); togglesign(*s); break;
3640 0 : case 3: *s = subsr(-1,z); *c = subrs(mpaut(z),1); break;
3641 819 : case 4: *c = rcopy(z); *s = mpaut(z); togglesign(*s); break;
3642 91 : case 5: *s = addsr( 1,z); *c = subrs(mpaut(z),1); break;
3643 7 : case 6: *c = subsr(-2,z); *s = mpaut(z); break;
3644 7 : case 7: *s = subsr(-1,z); *c = subsr(-1,mpaut(z)); break;
3645 : }
3646 5838 : gptr[0] = s; gptr[1] = c;
3647 5838 : gerepilemanysp(av,tetpil,gptr,2);
3648 : }
3649 :
3650 : /* return exp(ix), x a t_REAL */
3651 : GEN
3652 683518 : expIr(GEN x)
3653 : {
3654 683518 : pari_sp av = avma;
3655 683518 : GEN v = cgetg(3,t_COMPLEX);
3656 683523 : mpsincos(x, (GEN*)(v+2), (GEN*)(v+1));
3657 683528 : if (!signe(gel(v,2))) return gerepilecopy(av, gel(v,1));
3658 681393 : return v;
3659 : }
3660 :
3661 : /* return exp(ix)-1, x a t_REAL */
3662 : static GEN
3663 5838 : expm1_Ir(GEN x)
3664 : {
3665 5838 : pari_sp av = avma;
3666 5838 : GEN v = cgetg(3,t_COMPLEX);
3667 5838 : mpsincosm1(x, (GEN*)(v+2), (GEN*)(v+1));
3668 5838 : if (!signe(gel(v,2))) return gerepilecopy(av, gel(v,1));
3669 5838 : return v;
3670 : }
3671 :
3672 : /* return exp(z)-1, z complex */
3673 : GEN
3674 5894 : cxexpm1(GEN z, long prec)
3675 : {
3676 5894 : pari_sp av = avma;
3677 5894 : GEN X, Y, x = real_i(z), y = imag_i(z);
3678 5894 : long l = precision(z);
3679 5894 : if (l) prec = l;
3680 5894 : if (typ(x) != t_REAL) x = gtofp(x, prec);
3681 5894 : if (typ(y) != t_REAL) y = gtofp(y, prec);
3682 5894 : if (gequal0(y)) return mpexpm1(x);
3683 5838 : if (gequal0(x)) return expm1_Ir(y);
3684 5705 : X = mpexpm1(x); /* t_REAL */
3685 5705 : Y = expm1_Ir(y);
3686 : /* exp(x+iy) - 1 = (exp(x)-1)(exp(iy)-1) + exp(x)-1 + exp(iy)-1 */
3687 5705 : return gerepileupto(av, gadd(gadd(X,Y), gmul(X,Y)));
3688 : }
3689 :
3690 : void
3691 2510158 : gsincos(GEN x, GEN *s, GEN *c, long prec)
3692 : {
3693 : long i, j, ex, ex2, lx, ly, mi;
3694 : pari_sp av, tetpil;
3695 : GEN y, r, u, v, u1, v1, p1, p2, p3, p4, ps, pc;
3696 : GEN *gptr[4];
3697 :
3698 2510158 : switch(typ(x))
3699 : {
3700 6951 : case t_INT: case t_FRAC:
3701 6951 : *s = cgetr(prec);
3702 6950 : *c = cgetr(prec); av = avma;
3703 6945 : mpsincos(tofp_safe(x, prec), &ps, &pc);
3704 6951 : affrr_fixlg(ps,*s);
3705 2510317 : affrr_fixlg(pc,*c); set_avma(av); return;
3706 :
3707 2498612 : case t_REAL:
3708 2498612 : mpsincos(x,s,c); return;
3709 :
3710 4130 : case t_COMPLEX:
3711 4130 : i = precision(x); if (i) prec = i;
3712 4130 : ps = cgetc(prec); *s = ps;
3713 4130 : pc = cgetc(prec); *c = pc; av = avma;
3714 4130 : r = gexp(gel(x,2),prec);
3715 4130 : v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
3716 4130 : u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
3717 4130 : gsincos(gel(x,1), &u,&v, prec);
3718 4130 : affrr_fixlg(mulrr(v1,u), gel(ps,1));
3719 4130 : affrr_fixlg(mulrr(u1,v), gel(ps,2));
3720 4130 : affrr_fixlg(mulrr(v1,v), gel(pc,1));
3721 4130 : affrr_fixlg(mulrr(u1,u), gel(pc,2)); togglesign(gel(pc,2));
3722 4130 : set_avma(av); return;
3723 :
3724 0 : case t_QUAD:
3725 0 : av = avma; gsincos(quadtofp(x, prec), s, c, prec);
3726 0 : gerepileall(av, 2, s, c); return;
3727 :
3728 465 : default:
3729 465 : av = avma; if (!(y = toser_i(x))) break;
3730 504 : if (gequal0(y)) { *s = gerepilecopy(av,y); *c = gaddsg(1,*s); return; }
3731 :
3732 504 : ex = valp(y); lx = lg(y); ex2 = 2*ex+2;
3733 504 : if (ex < 0) pari_err_DOMAIN("gsincos","valuation", "<", gen_0, x);
3734 504 : if (ex2 > lx)
3735 : {
3736 98 : *s = x == y? gcopy(y): gerepilecopy(av, y); av = avma;
3737 98 : *c = gerepileupto(av, gsubsg(1, gdivgu(gsqr(y),2)));
3738 98 : return;
3739 : }
3740 406 : if (!ex)
3741 : {
3742 105 : gsincos(serchop0(y),&u,&v,prec);
3743 105 : gsincos(gel(y,2),&u1,&v1,prec);
3744 105 : p1 = gmul(v1,v);
3745 105 : p2 = gmul(u1,u);
3746 105 : p3 = gmul(v1,u);
3747 105 : p4 = gmul(u1,v); tetpil = avma;
3748 105 : *c = gsub(p1,p2);
3749 105 : *s = gadd(p3,p4);
3750 105 : gptr[0]=s; gptr[1]=c;
3751 105 : gerepilemanysp(av,tetpil,gptr,2);
3752 105 : return;
3753 : }
3754 :
3755 301 : ly = lx+2*ex;
3756 2842 : mi = lx-1; while (mi>=3 && isrationalzero(gel(y,mi))) mi--;
3757 301 : mi += ex-2;
3758 301 : pc = cgetg(ly,t_SER); *c = pc;
3759 301 : ps = cgetg(lx,t_SER); *s = ps;
3760 301 : pc[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(y));
3761 301 : gel(pc,2) = gen_1; ps[1] = y[1];
3762 609 : for (i=2; i<ex+2; i++) gel(ps,i) = gcopy(gel(y,i));
3763 616 : for (i=3; i< ex2; i++) gel(pc,i) = gen_0;
3764 3339 : for (i=ex2; i<ly; i++)
3765 : {
3766 3038 : long ii = i-ex;
3767 3038 : av = avma; p1 = gen_0;
3768 7028 : for (j=ex; j<=minss(ii-2,mi); j++)
3769 3990 : p1 = gadd(p1, gmulgu(gmul(gel(y,j-ex+2),gel(ps,ii-j)),j));
3770 3038 : gel(pc,i) = gerepileupto(av, gdivgs(p1,2-i));
3771 3038 : if (ii < lx)
3772 : {
3773 2730 : av = avma; p1 = gen_0;
3774 5796 : for (j=ex; j<=minss(i-ex2,mi); j++)
3775 3066 : p1 = gadd(p1,gmulgu(gmul(gel(y,j-ex+2),gel(pc,i-j)),j));
3776 2730 : p1 = gdivgu(p1,i-2);
3777 2730 : gel(ps,ii) = gerepileupto(av, gadd(p1,gel(y,ii)));
3778 : }
3779 : }
3780 301 : return;
3781 : }
3782 0 : pari_err_TYPE("gsincos",x);
3783 : }
3784 :
3785 : /********************************************************************/
3786 : /** **/
3787 : /** SINC **/
3788 : /** **/
3789 : /********************************************************************/
3790 : static GEN
3791 2428754 : mpsinc(GEN x)
3792 : {
3793 2428754 : pari_sp av = avma;
3794 : GEN s, c;
3795 :
3796 2428754 : if (!signe(x)) {
3797 0 : long l = nbits2prec(-expo(x));
3798 0 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
3799 0 : return real_1(l);
3800 : }
3801 :
3802 2428754 : mpsincos(x,&s,&c);
3803 2428754 : return gerepileuptoleaf(av, divrr(s,x));
3804 : }
3805 :
3806 : GEN
3807 2428866 : gsinc(GEN x, long prec)
3808 : {
3809 : pari_sp av;
3810 : GEN r, u, v, y, u1, v1;
3811 : long i;
3812 :
3813 2428866 : switch(typ(x))
3814 : {
3815 2428733 : case t_REAL: return mpsinc(x);
3816 49 : case t_COMPLEX:
3817 49 : if (isintzero(gel(x,1)))
3818 : {
3819 28 : av = avma; x = gel(x,2);
3820 28 : if (gequal0(x)) return gcosh(x,prec);
3821 14 : return gerepileuptoleaf(av,gdiv(gsinh(x,prec),x));
3822 : }
3823 21 : i = precision(x); if (i) prec = i;
3824 21 : y = cgetc(prec); av = avma;
3825 21 : r = gexp(gel(x,2),prec);
3826 21 : v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
3827 21 : u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
3828 21 : gsincos(gel(x,1),&u,&v,prec);
3829 21 : affc_fixlg(gdiv(mkcomplex(gmul(v1,u), gmul(u1,v)), x), y);
3830 21 : return gc_const(av,y);
3831 :
3832 14 : case t_INT:
3833 14 : if (!signe(x)) return real_1(prec); /*fall through*/
3834 : case t_FRAC:
3835 21 : y = cgetr(prec); av = avma;
3836 21 : affrr_fixlg(mpsinc(tofp_safe(x,prec)), y); return gc_const(av,y);
3837 :
3838 21 : case t_PADIC:
3839 21 : if (gequal0(x)) return cvtop(gen_1, gel(x,2), valp(x));
3840 14 : av = avma; y = sin_p(x);
3841 14 : if (!y) pari_err_DOMAIN("gsinc(t_PADIC)","argument","",gen_0,x);
3842 7 : return gerepileupto(av,gdiv(y,x));
3843 :
3844 35 : default:
3845 : {
3846 : long ex;
3847 35 : av = avma; if (!(y = toser_i(x))) break;
3848 35 : if (gequal0(y)) return gerepileupto(av, gaddsg(1,y));
3849 35 : ex = valp(y);
3850 35 : if (ex < 0) pari_err_DOMAIN("sinc","valuation", "<", gen_0, x);
3851 28 : if (ex)
3852 : {
3853 28 : gsincos(y,&u,&v,prec);
3854 28 : y = gerepileupto(av, gdiv(u,y));
3855 28 : if (lg(y) > 2) gel(y,2) = gen_1;
3856 28 : return y;
3857 : }
3858 : else
3859 : {
3860 0 : GEN z0, y0 = gel(y,2), y1 = serchop0(y), y10 = y1;
3861 0 : if (!gequal1(y0)) y10 = gdiv(y10, y0);
3862 0 : gsincos(y1,&u,&v,prec);
3863 0 : z0 = gdiv(gcos(y0,prec), y0);
3864 0 : y = gaddsg(1, y10);
3865 0 : u = gadd(gmul(gsinc(y0, prec),v), gmul(z0, u));
3866 0 : return gerepileupto(av,gdiv(u,y));
3867 : }
3868 : }
3869 : }
3870 0 : return trans_eval("sinc",gsinc,x,prec);
3871 : }
3872 :
3873 : /********************************************************************/
3874 : /** **/
3875 : /** TANGENT and COTANGENT **/
3876 : /** **/
3877 : /********************************************************************/
3878 : static GEN
3879 133 : mptan(GEN x)
3880 : {
3881 133 : pari_sp av = avma;
3882 : GEN s, c;
3883 :
3884 133 : mpsincos(x,&s,&c);
3885 133 : if (!signe(c))
3886 0 : pari_err_DOMAIN("tan", "argument", "=", strtoGENstr("Pi/2 + kPi"),x);
3887 133 : return gerepileuptoleaf(av, divrr(s,c));
3888 : }
3889 :
3890 : /* If exp(-|im(x)|) << 1, avoid overflow in sincos(x) */
3891 : static int
3892 4018 : tan_huge_im(GEN ix, long prec)
3893 : {
3894 4018 : long b, p = precision(ix);
3895 4018 : if (!p) p = prec;
3896 4018 : b = bit_accuracy(p);
3897 4018 : return (gexpo(ix) > b || fabs(gtodouble(ix)) > (M_LN2 / 2) * b);
3898 : }
3899 : /* \pm I */
3900 : static GEN
3901 35 : real_I(long s, long prec)
3902 : {
3903 35 : GEN z = cgetg(3, t_COMPLEX);
3904 35 : gel(z,1) = real_0(prec);
3905 35 : gel(z,2) = s > 0? real_1(prec): real_m1(prec); return z;
3906 : }
3907 :
3908 : GEN
3909 217 : gtan(GEN x, long prec)
3910 : {
3911 : pari_sp av;
3912 : GEN y, s, c;
3913 :
3914 217 : switch(typ(x))
3915 : {
3916 126 : case t_REAL: return mptan(x);
3917 :
3918 42 : case t_COMPLEX: {
3919 42 : if (isintzero(gel(x,1))) retmkcomplex(gen_0,gtanh(gel(x,2),prec));
3920 28 : if (tan_huge_im(gel(x,2), prec)) return real_I(gsigne(gel(x,2)), prec);
3921 14 : av = avma; y = mulcxmI(gtanh(mulcxI(x), prec)); /* tan x = -I th(I x) */
3922 14 : gel(y,1) = gcopy(gel(y,1)); return gerepileupto(av, y);
3923 : }
3924 7 : case t_INT: case t_FRAC:
3925 7 : y = cgetr(prec); av = avma;
3926 7 : affrr_fixlg(mptan(tofp_safe(x,prec)), y); return gc_const(av,y);
3927 :
3928 14 : case t_PADIC:
3929 14 : av = avma;
3930 14 : return gerepileupto(av, gdiv(gsin(x,prec), gcos(x,prec)));
3931 :
3932 28 : default:
3933 28 : av = avma; if (!(y = toser_i(x))) break;
3934 21 : if (gequal0(y)) return gerepilecopy(av, y);
3935 21 : if (valp(y) < 0)
3936 7 : pari_err_DOMAIN("tan","valuation", "<", gen_0, x);
3937 14 : gsincos(y,&s,&c,prec);
3938 14 : return gerepileupto(av, gdiv(s,c));
3939 : }
3940 7 : return trans_eval("tan",gtan,x,prec);
3941 : }
3942 :
3943 : static GEN
3944 63 : mpcotan(GEN x)
3945 : {
3946 63 : pari_sp av=avma, tetpil;
3947 : GEN s,c;
3948 :
3949 63 : mpsincos(x,&s,&c); tetpil=avma;
3950 63 : return gerepile(av,tetpil,divrr(c,s));
3951 : }
3952 :
3953 : GEN
3954 4200 : gcotan(GEN x, long prec)
3955 : {
3956 : pari_sp av;
3957 : GEN y, s, c;
3958 :
3959 4200 : switch(typ(x))
3960 : {
3961 56 : case t_REAL:
3962 56 : return mpcotan(x);
3963 :
3964 4011 : case t_COMPLEX:
3965 4011 : if (isintzero(gel(x,1))) {
3966 21 : GEN z = cgetg(3, t_COMPLEX);
3967 21 : gel(z,1) = gen_0; av = avma;
3968 21 : gel(z,2) = gerepileupto(av, gneg(ginv(gtanh(gel(x,2),prec))));
3969 21 : return z;
3970 : }
3971 3990 : if (tan_huge_im(gel(x,2), prec)) return real_I(-gsigne(gel(x,2)), prec);
3972 3969 : av = avma; gsincos(x,&s,&c,prec);
3973 3969 : return gerepileupto(av, gdiv(c,s));
3974 :
3975 7 : case t_INT: case t_FRAC:
3976 7 : y = cgetr(prec); av = avma;
3977 7 : affrr_fixlg(mpcotan(tofp_safe(x,prec)), y); return gc_const(av,y);
3978 :
3979 14 : case t_PADIC:
3980 14 : av = avma;
3981 14 : return gerepileupto(av, gdiv(gcos(x,prec), gsin(x,prec)));
3982 :
3983 112 : default:
3984 112 : av = avma; if (!(y = toser_i(x))) break;
3985 105 : if (gequal0(y)) pari_err_DOMAIN("cotan", "argument", "=", gen_0, y);
3986 105 : if (valp(y) < 0) pari_err_DOMAIN("cotan","valuation", "<", gen_0, x);
3987 98 : gsincos(y,&s,&c,prec);
3988 98 : return gerepileupto(av, gdiv(c,s));
3989 : }
3990 7 : return trans_eval("cotan",gcotan,x,prec);
3991 : }
|
