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 318195 : pari_init_floats(void)
34 : {
35 318195 : gcatalan = geuler = gpi = zetazone = bernzone = glog2 = eulerzone = NULL;
36 318195 : }
37 :
38 : void
39 316309 : pari_close_floats(void)
40 : {
41 316309 : guncloneNULL(gcatalan);
42 315243 : guncloneNULL(geuler);
43 313842 : guncloneNULL(gpi);
44 313562 : guncloneNULL(glog2);
45 313288 : guncloneNULL(zetazone);
46 313093 : guncloneNULL_deep(bernzone);
47 313009 : guncloneNULL_deep(eulerzone);
48 312883 : }
49 :
50 : /********************************************************************/
51 : /** GENERIC BINARY SPLITTING **/
52 : /** (Haible, Papanikolaou) **/
53 : /********************************************************************/
54 : void
55 4388379 : abpq_init(struct abpq *A, long n)
56 : {
57 4388379 : A->a = (GEN*)new_chunk(n+1);
58 4399135 : A->b = (GEN*)new_chunk(n+1);
59 4405022 : A->p = (GEN*)new_chunk(n+1);
60 4405173 : A->q = (GEN*)new_chunk(n+1);
61 4406081 : }
62 : static GEN
63 121901576 : mulii3(GEN a, GEN b, GEN c) { return mulii(mulii(a,b),c); }
64 :
65 : /* T_{n1,n1+1} */
66 : static GEN
67 28503327 : T2(struct abpq *A, long n1)
68 : {
69 28503327 : GEN u = mulii3(A->a[n1], A->b[n1+1], A->q[n1+1]);
70 28552371 : GEN v = mulii3(A->b[n1], A->a[n1+1], A->p[n1+1]);
71 28559298 : 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 53484171 : 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 53484171 : switch(n2 - n1)
83 : {
84 : GEN b, q;
85 60 : case 1:
86 60 : r->P = A->p[n1];
87 60 : r->Q = A->q[n1];
88 60 : r->B = A->b[n1];
89 60 : r->T = mulii(A->a[n1], A->p[n1]);
90 28693552 : return;
91 17718113 : case 2:
92 17718113 : r->P = mulii(A->p[n1], A->p[n1+1]);
93 17561898 : r->Q = mulii(A->q[n1], A->q[n1+1]);
94 17537806 : r->B = mulii(A->b[n1], A->b[n1+1]);
95 17545221 : av = avma;
96 17545221 : r->T = gerepileuptoint(av, T2(A, n1));
97 17585192 : return;
98 :
99 11319597 : case 3:
100 11319597 : q = mulii(A->q[n1+1], A->q[n1+2]);
101 11258251 : b = mulii(A->b[n1+1], A->b[n1+2]);
102 11249274 : r->P = mulii3(A->p[n1], A->p[n1+1], A->p[n1+2]);
103 11249525 : r->Q = mulii(A->q[n1], q);
104 11243104 : r->B = mulii(A->b[n1], b);
105 11252102 : av = avma;
106 11252102 : u1 = mulii3(b, q, A->a[n1]);
107 11231566 : u2 = mulii(A->b[n1], T2(A, n1+1));
108 11249043 : r->T = gerepileuptoint(av, mulii(A->p[n1], addii(u1, u2)));
109 11108300 : return;
110 : }
111 :
112 24446401 : av = avma;
113 24446401 : n = (n1 + n2) >> 1;
114 24446401 : abpq_sum(&L, n1, n, A);
115 24498086 : abpq_sum(&R, n, n2, A);
116 :
117 24532696 : r->P = mulii(L.P, R.P);
118 24316624 : r->Q = mulii(L.Q, R.Q);
119 24273767 : r->B = mulii(L.B, R.B);
120 24308931 : u1 = mulii3(R.B, R.Q, L.T);
121 24232100 : u2 = mulii3(L.B, L.P, R.T);
122 24253265 : r->T = addii(u1,u2);
123 24275668 : set_avma(av);
124 24258817 : r->P = icopy(r->P);
125 24558992 : r->Q = icopy(r->Q);
126 24597295 : r->B = icopy(r->B);
127 24597196 : 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 63833 : swap_clone(GEN *old, GEN c)
138 63833 : { 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 31395 : pi_ramanujan(long prec)
149 : {
150 31395 : const ulong B = 545140134, A = 13591409, C = 640320;
151 31395 : 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 31395 : nmax = (long)(1 + prec2nbits(prec)/alpha2);
158 : #ifdef LONG_IS_64BIT
159 30870 : D = utoipos(10939058860032000UL); /* C^3/24 */
160 : #else
161 526 : D = uutoi(2546948UL,495419392UL);
162 : #endif
163 31390 : abpq_init(&S, nmax);
164 31392 : S.a[0] = utoipos(A);
165 31389 : S.b[0] = S.p[0] = S.q[0] = gen_1;
166 282586 : for (n = 1; n <= nmax; n++)
167 : {
168 251228 : S.a[n] = addiu(muluu(B, n), A);
169 251145 : S.b[n] = gen_1;
170 251145 : S.p[n] = mulis(muluu(6*n-5, 2*n-1), 1-6*n);
171 251253 : S.q[n] = mulii(sqru(n), muliu(D,n));
172 : }
173 31358 : abpq_sum(&R, 0, nmax, &S); prec2 = prec+EXTRAPREC64;
174 31396 : u = itor(muliu(R.Q,C/12), prec2);
175 31393 : 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 44940643 : constpi(long prec)
213 : {
214 : pari_sp av;
215 : GEN tmp;
216 44940643 : if (gpi && realprec(gpi) >= prec) return gpi;
217 :
218 31210 : av = avma;
219 31210 : tmp = gclone(pi_ramanujan(prec));
220 31406 : swap_clone(&gpi,tmp);
221 31405 : return gc_const(av, gpi);
222 : }
223 :
224 : GEN
225 44940364 : mppi(long prec) { return rtor(constpi(prec), prec); }
226 :
227 : /* Pi * 2^n */
228 : GEN
229 30337532 : Pi2n(long n, long prec)
230 : {
231 30337532 : GEN x = mppi(prec); shiftr_inplace(x, n);
232 30337007 : return x;
233 : }
234 :
235 : /* I * Pi * 2^n */
236 : GEN
237 268094 : PiI2n(long n, long prec) { retmkcomplex(gen_0, Pi2n(n, prec)); }
238 :
239 : /* 2I * Pi */
240 : GEN
241 265497 : PiI2(long prec) { return PiI2n(1, prec); }
242 :
243 : /********************************************************************/
244 : /** **/
245 : /** EULER CONSTANT **/
246 : /** **/
247 : /********************************************************************/
248 :
249 : GEN
250 57495 : 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 57495 : if (geuler && realprec(geuler) >= prec) return geuler;
257 :
258 504 : av1 = avma; tmpeuler = cgetr_block(prec);
259 :
260 504 : incrprec(prec);
261 :
262 504 : l = prec+EXTRAPREC64; x = (long) (1 + prec2nbits_mul(l, M_LN2/4));
263 504 : a = utor(x,l); u=logr_abs(a); setsigne(u,-1); affrr(u,a);
264 504 : b = real_1(l);
265 504 : v = real_1(l);
266 504 : n = (long)(1+3.591*x); /* z=3.591: z*[ ln(z)-1 ]=1 */
267 504 : n1 = minss(n, SQRTVERYBIGINT);
268 504 : if (x < SQRTVERYBIGINT)
269 : {
270 504 : ulong xx = x*x;
271 504 : av2 = avma;
272 166955 : for (k=1; k<n1; k++)
273 : {
274 166451 : affrr(divru(mulur(xx,b),k*k), b);
275 166470 : affrr(divru(addrr(divru(mulur(xx,a),k),b),k), a);
276 166461 : affrr(addrr(u,a), u);
277 166439 : affrr(addrr(v,b), v); set_avma(av2);
278 : }
279 1008 : for ( ; k<=n; k++)
280 : {
281 504 : affrr(divru(divru(mulur(xx,b),k),k), b);
282 504 : affrr(divru(addrr(divru(mulur(xx,a),k),b),k), a);
283 504 : affrr(addrr(u,a), u);
284 504 : 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 504 : divrrz(u,v,tmpeuler);
307 504 : swap_clone(&geuler,tmpeuler);
308 504 : return gc_const(av1, geuler);
309 : }
310 :
311 : GEN
312 57495 : 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 2013657 : transvec(GEN (*f)(GEN,long), GEN x, long prec)
364 6551710 : { 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 3867303 : trans_eval(const char *fun, GEN (*f)(GEN,long), GEN x, long prec)
371 : {
372 3867303 : pari_sp av = avma;
373 3867303 : if (prec < LOWDEFAULTPREC) pari_err_BUG("trans_eval [prec < 3]");
374 3867318 : switch(typ(x))
375 : {
376 1605606 : case t_INT: x = f(itor(x,prec),prec); break;
377 247998 : 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 2013644 : case t_VEC:
381 : case t_COL:
382 2013644 : case t_MAT: return transvec(f, x, prec);
383 49 : default: pari_err_TYPE(fun,x);
384 : return NULL;/*LCOV_EXCL_LINE*/
385 : }
386 1853609 : return gerepileupto(av, x);
387 : }
388 :
389 : GEN
390 1883 : trans_evalgen(const char *fun, void *E, GEN (*f)(void*,GEN,long),
391 : GEN x, long prec)
392 : {
393 1883 : pari_sp av = avma;
394 1883 : if (prec < LOWDEFAULTPREC) pari_err_BUG("trans_eval [prec < 3]");
395 1883 : switch(typ(x))
396 : {
397 273 : 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 1589 : 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 140090 : mpexp0(GEN x)
418 : {
419 140090 : long e = expo(x);
420 140090 : return e >= 0? real_0_bit(e): real_1_bit(-e);
421 : }
422 : static GEN
423 21055 : powr0(GEN x)
424 21055 : { return signe(x)? real_1(realprec(x)): mpexp0(x); }
425 :
426 : /* assume typ(x) = t_VEC */
427 : static int
428 49 : is_ext_qfr(GEN x)
429 35 : { return lg(x) == 3 && typ(gel(x,1)) == t_QFB && !qfb_is_qfi(gel(x,1))
430 84 : && typ(gel(x,2)) == t_REAL; }
431 :
432 : /* x t_POL or t_SER, return scalarpol(Rg_get_1(x)) */
433 : static GEN
434 370355 : scalarpol_get_1(GEN x)
435 : {
436 370355 : GEN y = cgetg(3,t_POL);
437 370355 : y[1] = evalvarn(varn(x)) | evalsigne(1);
438 370355 : gel(y,2) = Rg_get_1(x); return y;
439 : }
440 : /* to be called by the generic function gpowgs(x,s) when s = 0 */
441 : static GEN
442 2277859 : gpowg0(GEN x)
443 : {
444 : long lx, i;
445 : GEN y;
446 :
447 2277859 : switch(typ(x))
448 : {
449 1862195 : case t_INT: case t_REAL: case t_FRAC: case t_PADIC:
450 1862195 : return gen_1;
451 :
452 7 : case t_QUAD: x++; /*fall through*/
453 39247 : case t_COMPLEX: {
454 39247 : pari_sp av = avma;
455 39247 : GEN a = gpowg0(gel(x,1));
456 39247 : GEN b = gpowg0(gel(x,2));
457 39247 : if (a == gen_1) return b;
458 14 : if (b == gen_1) return a;
459 7 : return gerepileupto(av, gmul(a,b));
460 : }
461 133 : case t_INTMOD:
462 133 : y = cgetg(3,t_INTMOD);
463 133 : gel(y,1) = icopy(gel(x,1));
464 133 : gel(y,2) = is_pm1(gel(x,1))? gen_0: gen_1;
465 133 : return y;
466 :
467 5754 : case t_FFELT: return FF_1(x);
468 :
469 973 : case t_POLMOD:
470 973 : retmkpolmod(scalarpol_get_1(gel(x,1)), gcopy(gel(x,1)));
471 :
472 7 : case t_RFRAC:
473 7 : return scalarpol_get_1(gel(x,2));
474 369375 : case t_POL: case t_SER:
475 369375 : return scalarpol_get_1(x);
476 :
477 84 : case t_MAT:
478 84 : lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
479 77 : if (lx != lgcols(x)) pari_err_DIM("gpow");
480 77 : y = matid(lx-1);
481 252 : for (i=1; i<lx; i++) gcoeff(y,i,i) = gpowg0(gcoeff(x,i,i));
482 77 : return y;
483 21 : case t_VEC: if (!is_ext_qfr(x)) break;
484 : /* fall through handle extended t_QFB */
485 28 : case t_QFB: return qfbpow(x, gen_0);
486 49 : case t_VECSMALL: return identity_perm(lg(x) - 1);
487 : }
488 14 : pari_err_TYPE("gpow",x);
489 : return NULL; /* LCOV_EXCL_LINE */
490 : }
491 :
492 : static GEN
493 5711668 : _sqr(void *data /* ignored */, GEN x) { (void)data; return gsqr(x); }
494 : static GEN
495 3789296 : _mul(void *data /* ignored */, GEN x, GEN y) { (void)data; return gmul(x,y); }
496 : static GEN
497 334594 : _one(void *x) { return gpowg0((GEN) x); }
498 : static GEN
499 81908612 : _sqri(void *data /* ignored */, GEN x) { (void)data; return sqri(x); }
500 : static GEN
501 30001981 : _muli(void *data /* ignored */, GEN x, GEN y) { (void)data; return mulii(x,y); }
502 : static GEN
503 16130440 : _sqrr(void *data /* ignored */, GEN x) { (void)data; return sqrr(x); }
504 : static GEN
505 7053524 : _mulr(void *data /* ignored */, GEN x, GEN y) { (void)data; return mulrr(x,y); }
506 : static GEN
507 13755 : _oner(void *data /* prec */) { return real_1( *(long*) data); }
508 :
509 : /* INTEGER POWERING (a^n for integer a != 0 and integer n > 0)
510 : *
511 : * Use left shift binary algorithm (RS is wasteful: multiplies big numbers,
512 : * with LS one of them is the base, hence small). Sign of result is set
513 : * to s (= 1,-1). Makes life easier for caller, which otherwise might do a
514 : * setsigne(gen_1 / gen_m1) */
515 : static GEN
516 96924313 : powiu_sign(GEN a, ulong N, long s)
517 : {
518 : pari_sp av;
519 : GEN y;
520 :
521 96924313 : if (lgefint(a) == 3)
522 : { /* easy if |a| < 3 */
523 95401345 : ulong q = a[2];
524 95401345 : if (q == 1) return (s>0)? gen_1: gen_m1;
525 86228027 : if (q == 2) { a = int2u(N); setsigne(a,s); return a; }
526 63622204 : q = upowuu(q, N);
527 63625143 : if (q) return s>0? utoipos(q): utoineg(q);
528 : }
529 32647957 : if (N <= 2) {
530 1809826 : if (N == 2) return sqri(a);
531 20517 : a = icopy(a); setsigne(a,s); return a;
532 : }
533 30838131 : av = avma;
534 30838131 : y = gen_powu_i(a, N, NULL, &_sqri, &_muli);
535 30838443 : setsigne(y,s); return gerepileuptoint(av, y);
536 : }
537 : /* a^n */
538 : GEN
539 97083526 : powiu(GEN a, ulong n)
540 : {
541 : long s;
542 97083526 : if (!n) return gen_1;
543 95824108 : s = signe(a);
544 95824108 : if (!s) return gen_0;
545 95749204 : return powiu_sign(a, n, (s < 0 && odd(n))? -1: 1);
546 : }
547 : GEN
548 21105922 : powis(GEN a, long n)
549 : {
550 : long s;
551 : GEN t, y;
552 21105922 : if (n >= 0) return powiu(a, n);
553 619173 : s = signe(a);
554 619173 : if (!s) pari_err_INV("powis",gen_0);
555 619173 : t = (s < 0 && odd(n))? gen_m1: gen_1;
556 619173 : if (is_pm1(a)) return t;
557 : /* n < 0, |a| > 1 */
558 616772 : y = cgetg(3,t_FRAC);
559 616771 : gel(y,1) = t;
560 616771 : gel(y,2) = powiu_sign(a, -n, 1); /* force denominator > 0 */
561 616767 : return y;
562 : }
563 : GEN
564 46240404 : powuu(ulong p, ulong N)
565 : {
566 : pari_sp av;
567 : ulong pN;
568 : GEN y;
569 46240404 : if (!p) return gen_0;
570 46240327 : if (N <= 2)
571 : {
572 40264392 : if (N == 2) return sqru(p);
573 37981148 : if (N == 1) return utoipos(p);
574 5087359 : return gen_1;
575 : }
576 5975935 : pN = upowuu(p, N);
577 5975953 : if (pN) return utoipos(pN);
578 993345 : if (p == 2) return int2u(N);
579 979928 : av = avma;
580 979928 : y = gen_powu_i(utoipos(p), N, NULL, &_sqri, &_muli);
581 979926 : return gerepileuptoint(av, y);
582 : }
583 :
584 : /* return 0 if overflow */
585 : static ulong
586 19147322 : usqru(ulong p) { return p & HIGHMASK? 0: p*p; }
587 : ulong
588 101184906 : upowuu(ulong p, ulong k)
589 : {
590 : #ifdef LONG_IS_64BIT
591 86862779 : const ulong CUTOFF3 = 2642245;
592 86862779 : const ulong CUTOFF4 = 65535;
593 86862779 : const ulong CUTOFF5 = 7131;
594 86862779 : const ulong CUTOFF6 = 1625;
595 86862779 : const ulong CUTOFF7 = 565;
596 86862779 : const ulong CUTOFF8 = 255;
597 86862779 : const ulong CUTOFF9 = 138;
598 86862779 : const ulong CUTOFF10 = 84;
599 86862779 : const ulong CUTOFF11 = 56;
600 86862779 : const ulong CUTOFF12 = 40;
601 86862779 : const ulong CUTOFF13 = 30;
602 86862779 : const ulong CUTOFF14 = 23;
603 86862779 : const ulong CUTOFF15 = 19;
604 86862779 : const ulong CUTOFF16 = 15;
605 86862779 : const ulong CUTOFF17 = 13;
606 86862779 : const ulong CUTOFF18 = 11;
607 86862779 : const ulong CUTOFF19 = 10;
608 86862779 : const ulong CUTOFF20 = 9;
609 : #else
610 14322127 : const ulong CUTOFF3 = 1625;
611 14322127 : const ulong CUTOFF4 = 255;
612 14322127 : const ulong CUTOFF5 = 84;
613 14322127 : const ulong CUTOFF6 = 40;
614 14322127 : const ulong CUTOFF7 = 23;
615 14322127 : const ulong CUTOFF8 = 15;
616 14322127 : const ulong CUTOFF9 = 11;
617 14322127 : const ulong CUTOFF10 = 9;
618 14322127 : const ulong CUTOFF11 = 7;
619 14322127 : const ulong CUTOFF12 = 6;
620 14322127 : const ulong CUTOFF13 = 5;
621 14322127 : const ulong CUTOFF14 = 4;
622 14322127 : const ulong CUTOFF15 = 4;
623 14322127 : const ulong CUTOFF16 = 3;
624 14322127 : const ulong CUTOFF17 = 3;
625 14322127 : const ulong CUTOFF18 = 3;
626 14322127 : const ulong CUTOFF19 = 3;
627 14322127 : const ulong CUTOFF20 = 3;
628 : #endif
629 :
630 101184906 : if (p <= 2)
631 : {
632 9631824 : if (p < 2) return p;
633 9085074 : return k < BITS_IN_LONG? 1UL<<k: 0;
634 : }
635 91553082 : switch(k)
636 : {
637 : ulong p2, p3, p4, p5, p8;
638 8692307 : case 0: return 1;
639 20414805 : case 1: return p;
640 19147313 : case 2: return usqru(p);
641 3790662 : case 3: if (p > CUTOFF3) return 0; return p*p*p;
642 11367997 : case 4: if (p > CUTOFF4) return 0; p2=p*p; return p2*p2;
643 2209768 : case 5: if (p > CUTOFF5) return 0; p2=p*p; return p2*p2*p;
644 6979168 : case 6: if (p > CUTOFF6) return 0; p2=p*p; return p2*p2*p2;
645 376301 : case 7: if (p > CUTOFF7) return 0; p2=p*p; return p2*p2*p2*p;
646 456124 : case 8: if (p > CUTOFF8) return 0; p2=p*p; p4=p2*p2; return p4*p4;
647 396101 : case 9: if (p > CUTOFF9) return 0; p2=p*p; p4=p2*p2; return p4*p4*p;
648 4914760 : case 10: if (p > CUTOFF10)return 0; p2=p*p; p4=p2*p2; return p4*p4*p2;
649 170327 : case 11: if (p > CUTOFF11)return 0; p2=p*p; p4=p2*p2; return p4*p4*p2*p;
650 4792761 : case 12: if (p > CUTOFF12)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4;
651 96977 : case 13: if (p > CUTOFF13)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4*p;
652 4746601 : case 14: if (p > CUTOFF14)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4*p2;
653 159805 : case 15: if (p > CUTOFF15)return 0;
654 100065 : p2=p*p; p3=p2*p; p5=p3*p2; return p5*p5*p5;
655 104239 : case 16: if (p > CUTOFF16)return 0;
656 52261 : p2=p*p; p4=p2*p2; p8=p4*p4; return p8*p8;
657 79533 : case 17: if (p > CUTOFF17)return 0;
658 41875 : p2=p*p; p4=p2*p2; p8=p4*p4; return p*p8*p8;
659 69151 : case 18: if (p > CUTOFF18)return 0;
660 38791 : p2=p*p; p4=p2*p2; p8=p4*p4; return p2*p8*p8;
661 818921 : case 19: if (p > CUTOFF19)return 0;
662 765067 : p2=p*p; p4=p2*p2; p8=p4*p4; return p*p2*p8*p8;
663 63733 : case 20: if (p > CUTOFF20)return 0;
664 21661 : p2=p*p; p4=p2*p2; p8=p4*p4; return p4*p8*p8;
665 : }
666 : #ifdef LONG_IS_64BIT
667 1485785 : switch(p)
668 : {
669 221633 : case 3: if (k > 40) return 0;
670 135638 : break;
671 17028 : case 4: if (k > 31) return 0;
672 774 : return 1UL<<(2*k);
673 637067 : case 5: if (k > 27) return 0;
674 20044 : break;
675 49650 : case 6: if (k > 24) return 0;
676 9180 : break;
677 55577 : case 7: if (k > 22) return 0;
678 2732 : break;
679 504830 : default: return 0;
680 : }
681 : /* no overflow */
682 : {
683 167594 : ulong q = upowuu(p, k >> 1);
684 167594 : q *= q ;
685 167594 : return odd(k)? q*p: q;
686 : }
687 : #else
688 219943 : return 0;
689 : #endif
690 : }
691 :
692 : GEN
693 12017 : upowers(ulong x, long n)
694 : {
695 : long i;
696 12017 : GEN p = cgetg(n + 2, t_VECSMALL);
697 12017 : uel(p,1) = 1; if (n == 0) return p;
698 12017 : uel(p,2) = x;
699 91465 : for (i = 3; i <= n; i++)
700 79448 : uel(p,i) = uel(p,i-1)*x;
701 12017 : return p;
702 : }
703 :
704 : typedef struct {
705 : long prec, a;
706 : GEN (*sqr)(GEN);
707 : GEN (*mulug)(ulong,GEN);
708 : } sr_muldata;
709 :
710 : static GEN
711 1597247 : _rpowuu_sqr(void *data, GEN x)
712 : {
713 1597247 : sr_muldata *D = (sr_muldata *)data;
714 1597247 : if (typ(x) == t_INT && lg2prec(lgefint(x)) >= D->prec)
715 : { /* switch to t_REAL */
716 157825 : D->sqr = &sqrr;
717 157825 : D->mulug = &mulur; x = itor(x, D->prec);
718 : }
719 1597247 : return D->sqr(x);
720 : }
721 :
722 : static GEN
723 624920 : _rpowuu_msqr(void *data, GEN x)
724 : {
725 624920 : GEN x2 = _rpowuu_sqr(data, x);
726 624920 : sr_muldata *D = (sr_muldata *)data;
727 624920 : return D->mulug(D->a, x2);
728 : }
729 :
730 : /* return a^n as a t_REAL of precision prec. Assume a > 0, n > 0 */
731 : GEN
732 428872 : rpowuu(ulong a, ulong n, long prec)
733 : {
734 : pari_sp av;
735 : GEN y, z;
736 : sr_muldata D;
737 :
738 428872 : if (a == 1) return real_1(prec);
739 428872 : if (a == 2) return real2n(n, prec);
740 428872 : if (n == 1) return utor(a, prec);
741 423934 : z = cgetr(prec);
742 423934 : av = avma;
743 423934 : D.sqr = &sqri;
744 423934 : D.mulug = &mului;
745 423934 : D.prec = prec;
746 423934 : D.a = (long)a;
747 423934 : y = gen_powu_fold_i(utoipos(a), n, (void*)&D, &_rpowuu_sqr, &_rpowuu_msqr);
748 423934 : mpaff(y, z); return gc_const(av,z);
749 : }
750 :
751 : GEN
752 5095438 : powrs(GEN x, long n)
753 : {
754 5095438 : pari_sp av = avma;
755 : GEN y;
756 5095438 : if (!n) return powr0(x);
757 5095438 : y = gen_powu_i(x, (ulong)labs(n), NULL, &_sqrr, &_mulr);
758 5096038 : if (n < 0) y = invr(y);
759 5095764 : return gerepileuptoleaf(av,y);
760 : }
761 : GEN
762 6008876 : powru(GEN x, ulong n)
763 : {
764 6008876 : pari_sp av = avma;
765 : GEN y;
766 6008876 : if (!n) return powr0(x);
767 5988332 : y = gen_powu_i(x, n, NULL, &_sqrr, &_mulr);
768 5988211 : return gerepileuptoleaf(av,y);
769 : }
770 :
771 : GEN
772 13755 : powersr(GEN x, long n)
773 : {
774 13755 : long prec = realprec(x);
775 13755 : return gen_powers(x, n, 1, &prec, &_sqrr, &_mulr, &_oner);
776 : }
777 :
778 : /* x^(s/2), assume x t_REAL */
779 : GEN
780 0 : powrshalf(GEN x, long s)
781 : {
782 0 : if (s & 1) return sqrtr(powrs(x, s));
783 0 : return powrs(x, s>>1);
784 : }
785 : /* x^(s/2), assume x t_REAL */
786 : GEN
787 117813 : powruhalf(GEN x, ulong s)
788 : {
789 117813 : if (s & 1) return sqrtr(powru(x, s));
790 7085 : return powru(x, s>>1);
791 : }
792 : /* x^(n/d), assume x t_REAL, return t_REAL */
793 : GEN
794 511 : powrfrac(GEN x, long n, long d)
795 : {
796 : long z;
797 511 : if (!n) return powr0(x);
798 0 : z = cgcd(n, d); if (z > 1) { n /= z; d /= z; }
799 0 : if (d == 1) return powrs(x, n);
800 0 : x = powrs(x, n);
801 0 : if (d == 2) return sqrtr(x);
802 0 : return sqrtnr(x, d);
803 : }
804 :
805 : /* assume x != 0 */
806 : static GEN
807 624699 : pow_monome(GEN x, long n)
808 : {
809 624699 : long i, d, dx = degpol(x);
810 : GEN A, b, y;
811 :
812 624699 : if (n < 0) { n = -n; y = cgetg(3, t_RFRAC); } else y = NULL;
813 :
814 624699 : if (HIGHWORD(dx) || HIGHWORD(n))
815 8 : {
816 : LOCAL_HIREMAINDER;
817 9 : d = (long)mulll((ulong)dx, (ulong)n);
818 9 : if (hiremainder || (d &~ LGBITS)) d = LGBITS; /* overflow */
819 9 : d += 2;
820 : }
821 : else
822 624690 : d = dx*n + 2;
823 624699 : if ((d + 1) & ~LGBITS) pari_err(e_OVERFLOW,"pow_monome [degree]");
824 624692 : A = cgetg(d+1, t_POL); A[1] = x[1];
825 6051402 : for (i=2; i < d; i++) gel(A,i) = gen_0;
826 624691 : b = gpowgs(gel(x,dx+2), n); /* not memory clean if (n < 0) */
827 624692 : if (!y) y = A;
828 : else {
829 20482 : GEN c = denom_i(b);
830 20482 : gel(y,1) = c; if (c != gen_1) b = gmul(b,c);
831 20482 : gel(y,2) = A;
832 : }
833 624692 : gel(A,d) = b; return y;
834 : }
835 :
836 : /* x t_PADIC */
837 : static GEN
838 1305486 : powps(GEN x, long n)
839 : {
840 1305486 : long e = n*valp(x), v;
841 1305486 : GEN t, y, mod, p = gel(x,2);
842 : pari_sp av;
843 :
844 1305486 : if (!signe(gel(x,4))) {
845 84 : if (n < 0) pari_err_INV("powps",x);
846 77 : return zeropadic(p, e);
847 : }
848 1305402 : v = z_pval(n, p);
849 :
850 1305403 : y = cgetg(5,t_PADIC);
851 1305402 : mod = gel(x,3);
852 1305402 : if (v == 0) mod = icopy(mod);
853 : else
854 : {
855 86632 : if (precp(x) == 1 && absequaliu(p, 2)) v++;
856 86632 : mod = mulii(mod, powiu(p,v));
857 86632 : mod = gerepileuptoint((pari_sp)y, mod);
858 : }
859 1305399 : y[1] = evalprecp(precp(x) + v) | evalvalp(e);
860 1305398 : gel(y,2) = icopy(p);
861 1305398 : gel(y,3) = mod;
862 :
863 1305398 : av = avma; t = gel(x,4);
864 1305398 : if (n < 0) { t = Fp_inv(t, mod); n = -n; }
865 1305398 : t = Fp_powu(t, n, mod);
866 1305405 : gel(y,4) = gerepileuptoint(av, t);
867 1305402 : return y;
868 : }
869 : /* x t_PADIC */
870 : static GEN
871 161 : powp(GEN x, GEN n)
872 : {
873 : long v;
874 161 : GEN y, mod, p = gel(x,2);
875 :
876 161 : if (valp(x)) pari_err_OVERFLOW("valp()");
877 :
878 161 : if (!signe(gel(x,4))) {
879 14 : if (signe(n) < 0) pari_err_INV("powp",x);
880 7 : return zeropadic(p, 0);
881 : }
882 147 : v = Z_pval(n, p);
883 :
884 147 : y = cgetg(5,t_PADIC);
885 147 : mod = gel(x,3);
886 147 : if (v == 0) mod = icopy(mod);
887 : else
888 : {
889 70 : mod = mulii(mod, powiu(p,v));
890 70 : mod = gerepileuptoint((pari_sp)y, mod);
891 : }
892 147 : y[1] = evalprecp(precp(x) + v) | _evalvalp(0);
893 147 : gel(y,2) = icopy(p);
894 147 : gel(y,3) = mod;
895 147 : gel(y,4) = Fp_pow(gel(x,4), n, mod);
896 147 : return y;
897 : }
898 : static GEN
899 23753 : pow_polmod(GEN x, GEN n)
900 : {
901 23753 : GEN z = cgetg(3, t_POLMOD), a = gel(x,2), T = gel(x,1);
902 23753 : gel(z,1) = gcopy(T);
903 23753 : if (typ(a) != t_POL || varn(a) != varn(T) || lg(a) <= 3)
904 1269 : a = powgi(a, n);
905 : else {
906 22484 : pari_sp av = avma;
907 22484 : GEN p = NULL;
908 22484 : if (RgX_is_FpX(T, &p) && RgX_is_FpX(a, &p) && p)
909 : {
910 8771 : T = RgX_to_FpX(T, p); a = RgX_to_FpX(a, p);
911 8771 : if (lgefint(p) == 3)
912 : {
913 8764 : ulong pp = p[2];
914 8764 : a = Flxq_pow(ZX_to_Flx(a, pp), n, ZX_to_Flx(T, pp), pp);
915 8764 : a = Flx_to_ZX(a);
916 : }
917 : else
918 7 : a = FpXQ_pow(a, n, T, p);
919 8771 : a = FpX_to_mod(a, p);
920 8771 : a = gerepileupto(av, a);
921 : }
922 : else
923 : {
924 13713 : set_avma(av);
925 13713 : a = RgXQ_pow(a, n, gel(z,1));
926 : }
927 : }
928 23753 : gel(z,2) = a; return z;
929 : }
930 :
931 : GEN
932 114842406 : gpowgs(GEN x, long n)
933 : {
934 : long m;
935 : pari_sp av;
936 : GEN y;
937 :
938 114842406 : if (n == 0) return gpowg0(x);
939 112977820 : if (n == 1)
940 : {
941 73370493 : long t = typ(x);
942 73370493 : if (is_scalar_t(t)) return gcopy(x);
943 719568 : switch(t)
944 : {
945 666517 : case t_POL: case t_SER: case t_RFRAC: case t_MAT: case t_VECSMALL:
946 666517 : return gcopy(x);
947 21 : case t_VEC: if (!is_ext_qfr(x)) break;
948 : /* fall through handle extended t_QFB */
949 53037 : case t_QFB: return qfbred(x);
950 : }
951 14 : pari_err_TYPE("gpow", x);
952 : }
953 39607445 : if (n ==-1) return ginv(x);
954 32265534 : switch(typ(x))
955 : {
956 20947687 : case t_INT: return powis(x,n);
957 5086578 : case t_REAL: return powrs(x,n);
958 29264 : case t_INTMOD:
959 29264 : y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
960 29264 : gel(y,2) = Fp_pows(gel(x,2), n, gel(x,1));
961 29264 : return y;
962 278956 : case t_FRAC:
963 : {
964 278956 : GEN a = gel(x,1), b = gel(x,2);
965 278956 : long s = (signe(a) < 0 && odd(n))? -1: 1;
966 278956 : if (n < 0) {
967 3038 : n = -n;
968 3038 : if (is_pm1(a)) return powiu_sign(b, n, s); /* +-1/x[2] inverts to t_INT */
969 2821 : swap(a, b);
970 : }
971 278739 : y = cgetg(3, t_FRAC);
972 278739 : gel(y,1) = powiu_sign(a, n, s);
973 278739 : gel(y,2) = powiu_sign(b, n, 1);
974 278739 : return y;
975 : }
976 1305486 : case t_PADIC: return powps(x, n);
977 249144 : case t_RFRAC:
978 : {
979 249144 : av = avma; y = cgetg(3, t_RFRAC); m = labs(n);
980 249144 : gel(y,1) = gpowgs(gel(x,1),m);
981 249144 : gel(y,2) = gpowgs(gel(x,2),m);
982 249144 : if (n < 0) y = ginv(y);
983 249144 : return gerepileupto(av,y);
984 : }
985 23746 : case t_POLMOD: {
986 23746 : long N[] = {evaltyp(t_INT) | _evallg(3),0,0};
987 23746 : affsi(n,N); return pow_polmod(x, N);
988 : }
989 7 : case t_VEC: if (!is_ext_qfr(x)) pari_err_TYPE("gpow", x);
990 : /* fall through handle extended t_QFB */
991 1293693 : case t_QFB: return qfbpows(x, n);
992 1183415 : case t_POL:
993 1183415 : if (RgX_is_monomial(x)) return pow_monome(x, n);
994 : default: {
995 2426281 : pari_sp av = avma;
996 2426281 : y = gen_powu_i(x, (ulong)labs(n), NULL, &_sqr, &_mul);
997 2426287 : if (n < 0) y = ginv(y);
998 2426295 : return gerepileupto(av,y);
999 : }
1000 : }
1001 : }
1002 :
1003 : /* n a t_INT */
1004 : GEN
1005 103366176 : powgi(GEN x, GEN n)
1006 : {
1007 : GEN y;
1008 :
1009 103366176 : if (!is_bigint(n)) return gpowgs(x, itos(n));
1010 : /* probable overflow for nonmodular types (typical exception: (X^0)^N) */
1011 25579 : switch(typ(x))
1012 : {
1013 25261 : case t_INTMOD:
1014 25261 : y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
1015 25261 : gel(y,2) = Fp_pow(gel(x,2), n, gel(x,1));
1016 25262 : return y;
1017 101 : case t_FFELT: return FF_pow(x,n);
1018 161 : case t_PADIC: return powp(x, n);
1019 :
1020 35 : case t_INT:
1021 35 : if (is_pm1(x)) return (signe(x) < 0 && mpodd(n))? gen_m1: gen_1;
1022 14 : if (signe(x)) pari_err_OVERFLOW("lg()");
1023 7 : if (signe(n) < 0) pari_err_INV("powgi",gen_0);
1024 7 : return gen_0;
1025 7 : case t_FRAC:
1026 7 : pari_err_OVERFLOW("lg()");
1027 :
1028 0 : case t_VEC: if (!is_ext_qfr(x)) pari_err_TYPE("gpow",x);
1029 : /* fall through handle extended t_QFB */
1030 12 : case t_QFB: return qfbpow(x, n);
1031 7 : case t_POLMOD: return pow_polmod(x, n);
1032 7 : default: {
1033 7 : pari_sp av = avma;
1034 7 : y = gen_pow_i(x, n, NULL, &_sqr, &_mul);
1035 7 : if (signe(n) < 0) return gerepileupto(av, ginv(y));
1036 7 : return gerepilecopy(av,y);
1037 : }
1038 : }
1039 : }
1040 :
1041 : /* Assume x = 1 + O(t), n a scalar. Return x^n */
1042 : static GEN
1043 7854 : ser_pow_1(GEN x, GEN n)
1044 : {
1045 : long lx, mi, i, j, d;
1046 7854 : GEN y = cgetg_copy(x, &lx), X = x+2, Y = y + 2;
1047 7854 : y[1] = evalsigne(1) | _evalvalser(0) | evalvarn(varn(x));
1048 74179 : d = mi = lx-3; while (mi>=1 && isrationalzero(gel(X,mi))) mi--;
1049 7854 : gel(Y,0) = gen_1;
1050 110383 : for (i=1; i<=d; i++)
1051 : {
1052 102529 : pari_sp av = avma;
1053 102529 : GEN s = gen_0;
1054 487704 : for (j=1; j<=minss(i,mi); j++)
1055 : {
1056 385175 : GEN t = gsubgs(gmulgu(n,j),i-j);
1057 385175 : s = gadd(s, gmul(gmul(t, gel(X,j)), gel(Y,i-j)));
1058 : }
1059 102529 : gel(Y,i) = gerepileupto(av, gdivgu(s,i));
1060 : }
1061 7854 : return y;
1062 : }
1063 :
1064 : /* we suppose n != 0, valser(x) = 0 and leading-term(x) != 0. Not stack clean */
1065 : static GEN
1066 7959 : ser_pow(GEN x, GEN n, long prec)
1067 : {
1068 : GEN y, c, lead;
1069 7959 : if (varncmp(gvar(n), varn(x)) <= 0) return gexp(gmul(n, glog(x,prec)), prec);
1070 7854 : lead = gel(x,2);
1071 7854 : if (gequal1(lead)) return ser_pow_1(x, n);
1072 7469 : x = ser_normalize(x);
1073 7469 : if (typ(n) == t_FRAC && !isinexact(lead) && ispower(lead, gel(n,2), &c))
1074 112 : c = powgi(c, gel(n,1));
1075 : else
1076 7357 : c = gpow(lead,n, prec);
1077 7469 : y = gmul(c, ser_pow_1(x, n));
1078 : /* gpow(t_POLMOD,n) can be a t_COL [conjvec] */
1079 7469 : if (typ(y) != t_SER) pari_err_TYPE("gpow", y);
1080 7469 : return y;
1081 : }
1082 :
1083 : static long
1084 7868 : val_from_i(GEN E)
1085 : {
1086 7868 : if (is_bigint(E)) pari_err_OVERFLOW("sqrtn [valuation]");
1087 7861 : return itos(E);
1088 : }
1089 :
1090 : /* return x^q, assume typ(x) = t_SER, typ(q) = t_INT/t_FRAC and q != 0 */
1091 : static GEN
1092 7875 : ser_powfrac(GEN x, GEN q, long prec)
1093 : {
1094 7875 : GEN y, E = gmulsg(valser(x), q);
1095 : long e;
1096 :
1097 7875 : if (!signe(x))
1098 : {
1099 21 : if (gsigne(q) < 0) pari_err_INV("gpow", x);
1100 21 : return zeroser(varn(x), val_from_i(gfloor(E)));
1101 : }
1102 7854 : if (typ(E) != t_INT)
1103 7 : pari_err_DOMAIN("sqrtn", "valuation", "!=", mkintmod(gen_0, gel(q,2)), x);
1104 7847 : e = val_from_i(E);
1105 7847 : y = leafcopy(x); setvalser(y, 0);
1106 7847 : y = ser_pow(y, q, prec);
1107 7847 : setvalser(y, e); return y;
1108 : }
1109 :
1110 : static GEN
1111 126 : gpow0(GEN x, GEN n, long prec)
1112 : {
1113 126 : pari_sp av = avma;
1114 : long i, lx;
1115 : GEN y;
1116 126 : switch(typ(n))
1117 : {
1118 84 : case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD:
1119 84 : break;
1120 35 : case t_VEC: case t_COL: case t_MAT:
1121 35 : y = cgetg_copy(n, &lx);
1122 105 : for (i=1; i<lx; i++) gel(y,i) = gpow0(x,gel(n,i),prec);
1123 35 : return y;
1124 7 : default: pari_err_TYPE("gpow(0,n)", n);
1125 : }
1126 84 : n = real_i(n);
1127 84 : if (gsigne(n) <= 0) pari_err_DOMAIN("gpow(0,n)", "n", "<=", gen_0, n);
1128 77 : if (!precision(x)) return gcopy(x);
1129 :
1130 14 : x = ground(gmulsg(gexpo(x),n));
1131 14 : if (is_bigint(x) || uel(x,2) >= HIGHEXPOBIT)
1132 7 : pari_err_OVERFLOW("gpow");
1133 7 : set_avma(av); return real_0_bit(itos(x));
1134 : }
1135 :
1136 : /* centermod(x, log(2)), set *sh to the quotient */
1137 : static GEN
1138 19238982 : modlog2(GEN x, long *sh)
1139 : {
1140 19238982 : double d = rtodbl(x), qd = (fabs(d) + M_LN2/2)/M_LN2;
1141 : long q;
1142 19238943 : if (dblexpo(qd) >= BITS_IN_LONG-1) pari_err_OVERFLOW("expo()");
1143 19238963 : q = d < 0 ? - (long) qd: (long) qd;
1144 19238963 : *sh = q;
1145 19238963 : if (q) {
1146 15579443 : long l = realprec(x) + EXTRAPRECWORD;
1147 15579443 : x = subrr(rtor(x,l), mulsr(q, mplog2(l)));
1148 15578913 : if (!signe(x)) return NULL;
1149 : }
1150 19238433 : return x;
1151 : }
1152 :
1153 : /* x^n, n a t_FRAC */
1154 : static GEN
1155 9404288 : powfrac(GEN x, GEN n, long prec)
1156 : {
1157 9404288 : GEN a = gel(n,1), d = gel(n,2);
1158 9404288 : long D = itos_or_0(d);
1159 9403778 : if (D == 2)
1160 : {
1161 7808192 : GEN y = gsqrt(x,prec);
1162 7808165 : if (!equali1(a)) y = gmul(y, powgi(x, shifti(subiu(a,1), -1)));
1163 7808985 : return y;
1164 : }
1165 1595586 : if (D && is_real_t(typ(x)) && gsigne(x) > 0)
1166 : { /* x^n = x^q * x^(r/D) */
1167 1591469 : GEN z, r, q = truedvmdis(a, D, &r);
1168 1591473 : if (typ(x) == t_REAL)
1169 : {
1170 171745 : z = sqrtnr(x, D);
1171 171745 : if (!equali1(r)) z = powgi(z, r);
1172 171745 : if (signe(q)) z = gmul(z, powgi(x, q));
1173 : }
1174 : else
1175 : {
1176 1419728 : GEN X = x;
1177 1419728 : x = gtofp(x, prec + nbits2extraprec(expi(r)));
1178 1419728 : z = sqrtnr(x, D);
1179 1419728 : if (!equali1(r)) z = powgi(z, r);
1180 1419728 : if (signe(q))
1181 : {
1182 16816 : long e = typ(X)==t_INT? expi(X): maxuu(expi(gel(X,1)), expi(gel(X,2)));
1183 16816 : z = gmul(z, powgi(cmpiu(muliu(q,e), realprec(x)) > 0? x: X, q));
1184 : }
1185 : }
1186 1591473 : return z;
1187 : }
1188 4117 : return NULL;
1189 : }
1190 :
1191 : /* n = a+ib, x > 0 real, ex ~ |log2(x)|; return precision at which
1192 : * log(x) must be computed to evaluate x^n */
1193 : long
1194 187542 : powcx_prec(long ex, GEN n, long prec)
1195 : {
1196 187542 : GEN a = gel(n,1), b = gel(n,2);
1197 187542 : long e = (ex < 2)? 0: expu(ex);
1198 187542 : e += gexpo_safe(is_rational_t(typ(a))? b: n);
1199 187542 : return e > 2? prec + nbits2extraprec(e): prec;
1200 : }
1201 : GEN
1202 4460934 : powcx(GEN x, GEN logx, GEN n, long prec)
1203 : {
1204 4460934 : GEN sxb, cxb, xa, a = gel(n,1), xb = gmul(gel(n,2), logx);
1205 4461832 : long sh, p = realprec(logx);
1206 4461832 : switch(typ(a))
1207 : {
1208 49969 : case t_INT: xa = powgi(x, a); break;
1209 4318372 : case t_FRAC: xa = powfrac(x, a, prec);
1210 4318252 : if (xa) break;
1211 : default:
1212 93497 : xa = modlog2(gmul(gel(n,1), logx), &sh);
1213 93508 : if (!xa) xa = real2n(sh, prec);
1214 : else
1215 : {
1216 93508 : if (signe(xa) && realprec(xa) > prec) setprec(xa, prec);
1217 93508 : xa = mpexp(xa); shiftr_inplace(xa, sh);
1218 : }
1219 : }
1220 4461695 : if (typ(xb) != t_REAL) return xa;
1221 4461695 : if (gexpo(xb) > 30)
1222 : {
1223 4128728 : GEN q, P = Pi2n(-2, p), z = addrr(xb,P); /* = x + Pi/4 */
1224 4127588 : shiftr_inplace(P, 1);
1225 4125927 : q = floorr(divrr(z, P)); /* round ( x / (Pi/2) ) */
1226 4102512 : xb = subrr(xb, mulir(q, P)); /* x mod Pi/2 */
1227 4126662 : sh = Mod4(q);
1228 : }
1229 : else
1230 : {
1231 332760 : long q = floor(rtodbl(xb) / (M_PI/2) + 0.5);
1232 332766 : if (q) xb = subrr(xb, mulsr(q, Pi2n(-1,p))); /* x mod Pi/2 */
1233 332766 : sh = q & 3;
1234 : }
1235 4458743 : if (signe(xb) && realprec(xb) > prec) setprec(xb, prec);
1236 4458743 : mpsincos(xb, &sxb, &cxb);
1237 4464235 : return gmul(xa, mulcxpowIs(mkcomplex(cxb, sxb), sh));
1238 : }
1239 :
1240 : GEN
1241 21380092 : gpow(GEN x, GEN n, long prec)
1242 : {
1243 21380092 : long prec0, i, lx, tx, tn = typ(n);
1244 : pari_sp av;
1245 : GEN y;
1246 :
1247 21380092 : if (tn == t_INT) return powgi(x,n);
1248 5874666 : tx = typ(x);
1249 5874666 : if (is_matvec_t(tx))
1250 : {
1251 49 : y = cgetg_copy(x, &lx);
1252 133 : for (i=1; i<lx; i++) gel(y,i) = gpow(gel(x,i),n,prec);
1253 49 : return y;
1254 : }
1255 5874645 : av = avma;
1256 5874645 : switch (tx)
1257 : {
1258 28 : case t_POL: case t_RFRAC: x = toser_i(x); /* fall through */
1259 7560 : case t_SER:
1260 7560 : if (tn == t_FRAC) return gerepileupto(av, ser_powfrac(x, n, prec));
1261 140 : if (valser(x))
1262 21 : pari_err_DOMAIN("gpow [irrational exponent]",
1263 : "valuation", "!=", gen_0, x);
1264 119 : if (lg(x) == 2) return gerepilecopy(av, x); /* O(1) */
1265 112 : return gerepileupto(av, ser_pow(x, n, prec));
1266 : }
1267 5867086 : if (gequal0(x)) return gpow0(x, n, prec);
1268 5867029 : if (tn == t_FRAC)
1269 : {
1270 5089440 : GEN p, z, a = gel(n,1), d = gel(n,2);
1271 5089440 : switch (tx)
1272 : {
1273 1474789 : case t_INT:
1274 1474789 : if (signe(x) < 0)
1275 : {
1276 42 : if (equaliu(d, 2) && Z_issquareall(negi(x), &z))
1277 : {
1278 21 : z = powgi(z, a);
1279 21 : if (Mod4(a) == 3) z = gneg(z);
1280 5085065 : return gerepilecopy(av, mkcomplex(gen_0, z));
1281 : }
1282 21 : break;
1283 : }
1284 1474747 : if (ispower(x, d, &z)) return powgi(z, a);
1285 1473011 : break;
1286 69827 : case t_FRAC:
1287 69827 : if (signe(gel(x,1)) < 0)
1288 : {
1289 28 : if (equaliu(d, 2) && ispower(absfrac(x), d, &z))
1290 7 : return gerepilecopy(av, mkcomplex(gen_0, powgi(z, a)));
1291 21 : break;
1292 : }
1293 69799 : if (ispower(x, d, &z)) return powgi(z, a);
1294 68427 : break;
1295 :
1296 21 : case t_INTMOD:
1297 21 : p = gel(x,1);
1298 21 : if (!BPSW_psp(p)) pari_err_PRIME("gpow",p);
1299 14 : y = cgetg(3,t_INTMOD); gel(y,1) = icopy(p);
1300 14 : av = avma;
1301 14 : z = Fp_sqrtn(gel(x,2), d, p, NULL);
1302 14 : if (!z) pari_err_SQRTN("gpow",x);
1303 7 : gel(y,2) = gerepileuptoint(av, Fp_pow(z, a, p));
1304 7 : return y;
1305 :
1306 14 : case t_PADIC:
1307 14 : z = Qp_sqrtn(x, d, NULL); if (!z) pari_err_SQRTN("gpow",x);
1308 7 : return gerepileupto(av, powgi(z, a));
1309 :
1310 21 : case t_FFELT:
1311 21 : return gerepileupto(av,FF_pow(FF_sqrtn(x,d,NULL),a));
1312 : }
1313 5086248 : z = powfrac(x, n, prec);
1314 5086137 : if (z) return gerepileupto(av, z);
1315 : }
1316 781699 : if (tn == t_COMPLEX && is_real_t(typ(x)) && gsigne(x) > 0)
1317 : {
1318 175943 : long p = powcx_prec(fabs(dbllog2(x)), n, prec);
1319 175943 : return gerepileupto(av, powcx(x, glog(x, p), n, prec));
1320 : }
1321 605756 : if (tn == t_PADIC) x = gcvtop(x, gel(n,2), precp(n));
1322 605756 : i = precision(n);
1323 605756 : if (i) prec = i;
1324 605756 : prec0 = prec;
1325 605756 : if (!gprecision(x))
1326 : {
1327 39401 : long e = gexpo_safe(n); /* avoided if n = 0 or gexpo not defined */
1328 39401 : if (e > 2) prec += nbits2extraprec(e);
1329 : }
1330 605756 : y = gmul(n, glog(x,prec));
1331 605728 : y = gexp(y,prec);
1332 605728 : if (prec0 == prec) return gerepileupto(av, y);
1333 29246 : return gerepilecopy(av, gprec_wtrunc(y,prec0));
1334 : }
1335 : GEN
1336 11591 : powPis(GEN s, long prec)
1337 : {
1338 11591 : pari_sp av = avma;
1339 : GEN x;
1340 11591 : if (typ(s) != t_COMPLEX) return gpow(mppi(prec), s, prec);
1341 490 : x = mppi(powcx_prec(1, s, prec));
1342 490 : return gerepileupto(av, powcx(x, logr_abs(x), s, prec));
1343 : }
1344 : GEN
1345 7903 : pow2Pis(GEN s, long prec)
1346 : {
1347 7903 : pari_sp av = avma;
1348 : GEN x;
1349 7903 : if (typ(s) != t_COMPLEX) return gpow(Pi2n(1,prec), s, prec);
1350 1876 : x = Pi2n(1, powcx_prec(2, s, prec));
1351 1876 : return gerepileupto(av, powcx(x, logr_abs(x), s, prec));
1352 : }
1353 :
1354 : GEN
1355 200632 : gpowers0(GEN x, long n, GEN x0)
1356 : {
1357 : long i, l;
1358 : GEN V;
1359 200632 : if (!x0) return gpowers(x,n);
1360 186142 : if (n < 0) return cgetg(1,t_VEC);
1361 186142 : l = n+2; V = cgetg(l, t_VEC); gel(V,1) = gcopy(x0);
1362 7453114 : for (i = 2; i < l; i++) gel(V,i) = gmul(gel(V,i-1),x);
1363 186183 : return V;
1364 : }
1365 :
1366 : GEN
1367 334601 : gpowers(GEN x, long n)
1368 : {
1369 334601 : if (n < 0) return cgetg(1,t_VEC);
1370 334594 : return gen_powers(x, n, 0, (void*)x, &_sqr, &_mul, &_one);
1371 : }
1372 :
1373 : /* return [q^1,q^4,...,q^{n^2}] */
1374 : GEN
1375 37394 : gsqrpowers(GEN q, long n)
1376 : {
1377 37394 : pari_sp av = avma;
1378 37394 : GEN L = gpowers0(gsqr(q), n, q); /* L[i] = q^(2i - 1), i <= n+1 */
1379 37394 : GEN v = cgetg(n+1, t_VEC);
1380 : long i;
1381 37394 : gel(v, 1) = gcopy(q);
1382 6717337 : for (i = 2; i <= n ; ++i) gel(v, i) = q = gmul(q, gel(L,i)); /* q^(i^2) */
1383 37394 : return gerepileupto(av, v);
1384 : }
1385 :
1386 : /* 4 | N. returns a vector RU which contains exp(2*i*k*Pi/N), k=0..N-1 */
1387 : static GEN
1388 958008 : grootsof1_4(long N, long prec)
1389 : {
1390 958008 : GEN z, RU = cgetg(N+1,t_COL), *v = ((GEN*)RU) + 1;
1391 958007 : long i, N2 = (N>>1), N4 = (N>>2), N8 = (N>>3);
1392 : /* z^N2 = -1, z^N4 = I; if z^k = a+I*b, then z^(N4-k) = I*conj(z) = b+a*I */
1393 :
1394 958007 : v[0] = gen_1; v[1] = z = rootsof1u_cx(N, prec);
1395 958009 : if (odd(N4)) N8++;
1396 1064822 : for (i=1; i<N8; i++)
1397 : {
1398 106813 : GEN t = v[i];
1399 106813 : v[i+1] = gmul(z, t);
1400 106813 : v[N4-i] = mkcomplex(gel(t,2), gel(t,1));
1401 : }
1402 2435843 : for (i=0; i<N4; i++) v[i+N4] = mulcxI(v[i]);
1403 3913671 : for (i=0; i<N2; i++) v[i+N2] = gneg(v[i]);
1404 958006 : return RU;
1405 : }
1406 :
1407 : /* as above, N arbitrary */
1408 : GEN
1409 1119673 : grootsof1(long N, long prec)
1410 : {
1411 : GEN z, RU, *v;
1412 : long i, k;
1413 :
1414 1119673 : if (N <= 0) pari_err_DOMAIN("rootsof1", "N", "<=", gen_0, stoi(N));
1415 1119659 : if ((N & 3) == 0) return grootsof1_4(N, prec);
1416 161650 : if (N <= 2) return N == 1? mkcol(gen_1): mkcol2(gen_1, gen_m1);
1417 43746 : k = (N+1)>>1;
1418 43746 : RU = cgetg(N+1,t_COL);
1419 43746 : v = ((GEN*)RU) + 1;
1420 43746 : v[0] = gen_1; v[1] = z = rootsof1u_cx(N, prec);
1421 103748 : for (i=2; i<k; i++) v[i] = gmul(z, v[i-1]);
1422 43746 : if (!odd(N)) v[i++] = gen_m1; /*avoid loss of accuracy*/
1423 147494 : for ( ; i<N; i++) v[i] = gconj(v[N-i]);
1424 43746 : return RU;
1425 : }
1426 :
1427 : /********************************************************************/
1428 : /** **/
1429 : /** RACINE CARREE **/
1430 : /** **/
1431 : /********************************************************************/
1432 : /* assume x unit, e = precp(x) */
1433 : GEN
1434 144690 : Z2_sqrt(GEN x, long e)
1435 : {
1436 144690 : ulong r = signe(x)>=0?mod16(x):16-mod16(x);
1437 : GEN z;
1438 : long ez;
1439 : pari_sp av;
1440 :
1441 144690 : switch(e)
1442 : {
1443 7 : case 1: return gen_1;
1444 161 : case 2: return (r & 3UL) == 1? gen_1: NULL;
1445 28 : case 3: return (r & 7UL) == 1? gen_1: NULL;
1446 71064 : case 4: if (r == 1) return gen_1;
1447 35133 : else return (r == 9)? utoipos(3): NULL;
1448 73430 : default: if ((r&7UL) != 1) return NULL;
1449 : }
1450 73430 : av = avma; z = (r==1)? gen_1: utoipos(3);
1451 73430 : ez = 3; /* number of correct bits in z (compared to sqrt(x)) */
1452 : for(;;)
1453 47978 : {
1454 : GEN mod;
1455 121408 : ez = (ez<<1) - 1;
1456 121408 : if (ez > e) ez = e;
1457 121408 : mod = int2n(ez);
1458 121408 : z = addii(z, remi2n(mulii(x, Fp_inv(z,mod)), ez));
1459 121408 : z = shifti(z, -1); /* (z + x/z) / 2 */
1460 121408 : if (e == ez) return gerepileuptoint(av, z);
1461 47978 : if (ez < e) ez--;
1462 47978 : if (gc_needed(av,2))
1463 : {
1464 0 : if (DEBUGMEM > 1) pari_warn(warnmem,"Qp_sqrt");
1465 0 : z = gerepileuptoint(av,z);
1466 : }
1467 : }
1468 : }
1469 :
1470 : /* x unit defined modulo p^e, e > 0 */
1471 : GEN
1472 1897 : Qp_sqrt(GEN x)
1473 : {
1474 1897 : long pp, e = valp(x);
1475 1897 : GEN z,y,mod, p = gel(x,2);
1476 :
1477 1897 : if (gequal0(x)) return zeropadic(p, (e+1) >> 1);
1478 1883 : if (e & 1) return NULL;
1479 :
1480 1869 : y = cgetg(5,t_PADIC);
1481 1869 : pp = precp(x);
1482 1869 : mod = gel(x,3);
1483 1869 : z = gel(x,4); /* lift to t_INT */
1484 1869 : e >>= 1;
1485 1869 : z = Zp_sqrt(z, p, pp);
1486 1869 : if (!z) return NULL;
1487 1806 : if (absequaliu(p,2))
1488 : {
1489 805 : pp = (pp <= 3) ? 1 : pp-1;
1490 805 : mod = int2n(pp);
1491 : }
1492 1001 : else mod = icopy(mod);
1493 1806 : y[1] = evalprecp(pp) | evalvalp(e);
1494 1806 : gel(y,2) = icopy(p);
1495 1806 : gel(y,3) = mod;
1496 1806 : gel(y,4) = z; return y;
1497 : }
1498 :
1499 : GEN
1500 420 : Zn_sqrt(GEN d, GEN fn)
1501 : {
1502 420 : pari_sp ltop = avma, btop;
1503 420 : GEN b = gen_0, m = gen_1;
1504 : long j, np;
1505 420 : if (typ(d) != t_INT) pari_err_TYPE("Zn_sqrt",d);
1506 420 : if (typ(fn) == t_INT)
1507 0 : fn = absZ_factor(fn);
1508 420 : else if (!is_Z_factorpos(fn))
1509 0 : pari_err_TYPE("Zn_sqrt",fn);
1510 420 : np = nbrows(fn);
1511 420 : btop = avma;
1512 1680 : for (j = 1; j <= np; ++j)
1513 : {
1514 : GEN bp, mp, pr, r;
1515 1260 : GEN p = gcoeff(fn, j, 1);
1516 1260 : long e = itos(gcoeff(fn, j, 2));
1517 1260 : long v = Z_pvalrem(d,p,&r);
1518 1260 : if (v >= e) bp =gen_0;
1519 : else
1520 : {
1521 1134 : if (odd(v)) return NULL;
1522 1134 : bp = Zp_sqrt(r, p, e-v);
1523 1134 : if (!bp) return NULL;
1524 1134 : if (v) bp = mulii(bp, powiu(p, v>>1L));
1525 : }
1526 1260 : mp = powiu(p, e);
1527 1260 : pr = mulii(m, mp);
1528 1260 : b = Z_chinese_coprime(b, bp, m, mp, pr);
1529 1260 : m = pr;
1530 1260 : if (gc_needed(btop, 1))
1531 0 : gerepileall(btop, 2, &b, &m);
1532 : }
1533 420 : return gerepileupto(ltop, b);
1534 : }
1535 :
1536 : static GEN
1537 18669 : sqrt_ser(GEN b, long prec)
1538 : {
1539 18669 : long e = valser(b), vx = varn(b), lx, lold, j;
1540 : ulong mask;
1541 : GEN a, x, lta, ltx;
1542 :
1543 18669 : if (!signe(b)) return zeroser(vx, e>>1);
1544 18669 : a = leafcopy(b);
1545 18669 : x = cgetg_copy(b, &lx);
1546 18669 : if (e & 1)
1547 14 : pari_err_DOMAIN("sqrtn", "valuation", "!=", mkintmod(gen_0, gen_2), b);
1548 18655 : a[1] = x[1] = evalsigne(1) | evalvarn(0) | _evalvalser(0);
1549 18655 : lta = gel(a,2);
1550 18655 : if (gequal1(lta)) ltx = lta;
1551 14833 : else if (!issquareall(lta,<x)) ltx = gsqrt(lta,prec);
1552 18648 : gel(x,2) = ltx;
1553 315399 : for (j = 3; j < lx; j++) gel(x,j) = gen_0;
1554 18648 : setlg(x,3);
1555 18648 : mask = quadratic_prec_mask(lx - 2);
1556 18648 : lold = 1;
1557 96295 : while (mask > 1)
1558 : {
1559 77647 : GEN y, x2 = gmul2n(x,1);
1560 77647 : long l = lold << 1, lx;
1561 :
1562 77647 : if (mask & 1) l--;
1563 77647 : mask >>= 1;
1564 77647 : setlg(a, l + 2);
1565 77647 : setlg(x, l + 2);
1566 77647 : y = sqr_ser_part(x, lold, l-1) - lold;
1567 374398 : for (j = lold+2; j < l+2; j++) gel(y,j) = gsub(gel(y,j), gel(a,j));
1568 77647 : y += lold; setvalser(y, lold);
1569 77647 : y = normalizeser(y);
1570 77647 : y = gsub(x, gdiv(y, x2)); /* = gmul2n(gsub(x, gdiv(a,x)), -1); */
1571 77647 : lx = minss(l+2, lg(y));
1572 374391 : for (j = lold+2; j < lx; j++) gel(x,j) = gel(y,j);
1573 77647 : lold = l;
1574 : }
1575 18648 : x[1] = evalsigne(1) | evalvarn(vx) | _evalvalser(e >> 1);
1576 18648 : return x;
1577 : }
1578 :
1579 : GEN
1580 57421438 : gsqrt(GEN x, long prec)
1581 : {
1582 : pari_sp av;
1583 : GEN y;
1584 :
1585 57421438 : switch(typ(x))
1586 : {
1587 4488781 : case t_INT:
1588 4488781 : if (!signe(x)) return real_0(prec); /* no loss of accuracy */
1589 4488711 : x = itor(x,prec); /* fall through */
1590 51392083 : case t_REAL: return sqrtr(x);
1591 :
1592 35 : case t_INTMOD:
1593 : {
1594 35 : GEN p = gel(x,1), a;
1595 35 : y = cgetg(3,t_INTMOD); gel(y,1) = icopy(p);
1596 35 : a = Fp_sqrt(gel(x,2),p);
1597 21 : if (!a)
1598 : {
1599 7 : if (!BPSW_psp(p)) pari_err_PRIME("sqrt [modulus]",p);
1600 7 : pari_err_SQRTN("gsqrt",x);
1601 : }
1602 14 : gel(y,2) = a; return y;
1603 : }
1604 :
1605 5761868 : case t_COMPLEX:
1606 : { /* (u+iv)^2 = a+ib <=> u^2+v^2 = sqrt(a^2+b^2), u^2-v^2=a, 2uv=b */
1607 5761868 : GEN a = gel(x,1), b = gel(x,2), r, u, v;
1608 5761868 : if (isrationalzero(b)) return gsqrt(a, prec);
1609 5761868 : y = cgetg(3,t_COMPLEX); av = avma;
1610 :
1611 5761867 : r = cxnorm(x);
1612 5761856 : if (typ(r) == t_INTMOD || typ(r) == t_PADIC)
1613 0 : pari_err_IMPL("sqrt(complex of t_INTMODs)");
1614 5761856 : r = gsqrt(r, prec); /* t_REAL, |a+Ib| */
1615 5761863 : if (!signe(r))
1616 67 : u = v = gerepileuptoleaf(av, sqrtr(r));
1617 5761796 : else if (gsigne(a) < 0)
1618 : {
1619 : /* v > 0 since r > 0, a < 0, rounding errors can't make the sum of two
1620 : * positive numbers = 0 */
1621 183140 : v = sqrtr( gmul2n(gsub(r,a), -1) );
1622 183139 : if (gsigne(b) < 0) togglesign(v);
1623 183139 : v = gerepileuptoleaf(av, v); av = avma;
1624 : /* v = 0 is impossible */
1625 183139 : u = gerepileuptoleaf(av, gdiv(b, shiftr(v,1)));
1626 : } else {
1627 5578656 : u = sqrtr( gmul2n(gadd(r,a), -1) );
1628 5578659 : u = gerepileuptoleaf(av, u); av = avma;
1629 5578660 : if (!signe(u)) /* possible if a = 0.0, e.g. sqrt(0.e-10+1e-10*I) */
1630 7 : v = u;
1631 : else
1632 5578653 : v = gerepileuptoleaf(av, gdiv(b, shiftr(u,1)));
1633 : }
1634 5761865 : gel(y,1) = u;
1635 5761865 : gel(y,2) = v; return y;
1636 : }
1637 :
1638 63 : case t_PADIC:
1639 63 : y = Qp_sqrt(x);
1640 63 : if (!y) pari_err_SQRTN("Qp_sqrt",x);
1641 42 : return y;
1642 :
1643 161 : case t_FFELT: return FF_sqrt(x);
1644 :
1645 264938 : default:
1646 264938 : av = avma; if (!(y = toser_i(x))) break;
1647 18669 : return gerepilecopy(av, sqrt_ser(y, prec));
1648 : }
1649 246269 : return trans_eval("sqrt",gsqrt,x,prec);
1650 : }
1651 : /********************************************************************/
1652 : /** **/
1653 : /** N-th ROOT **/
1654 : /** **/
1655 : /********************************************************************/
1656 :
1657 : static GEN
1658 303679 : Z_to_padic(GEN a, GEN p, long e)
1659 : {
1660 303679 : if (signe(a)==0)
1661 1036 : return zeropadic(p, e);
1662 : else
1663 : {
1664 302643 : GEN z = cgetg(5, t_PADIC);
1665 302643 : long v = Z_pvalrem(a, p, &a), d = e - v;
1666 302643 : z[1] = evalprecp(d) | evalvalp(v);
1667 302643 : gel(z,2) = icopy(p);
1668 302643 : gel(z,3) = powiu(p, d);
1669 302642 : gel(z,4) = a;
1670 302642 : return z;
1671 : }
1672 : }
1673 :
1674 : GEN
1675 196485 : Qp_log(GEN x)
1676 : {
1677 196485 : pari_sp av = avma;
1678 196485 : GEN y, p = gel(x,2), a = gel(x,4);
1679 196485 : long e = precp(x);
1680 :
1681 196485 : if (!signe(a)) pari_err_DOMAIN("Qp_log", "argument", "=", gen_0, x);
1682 196464 : if (absequaliu(p,2) || equali1(modii(a, p)))
1683 76986 : y = Zp_log(a, p, e);
1684 : else
1685 : { /* compute log(x^(p-1)) / (p-1) */
1686 119478 : GEN q = gel(x,3), t = subiu(p, 1);
1687 119478 : a = Fp_pow(a, t, q);
1688 119478 : y = Fp_mul(Zp_log(a, p, e), diviiexact(subsi(1, q), t), q);
1689 : }
1690 196464 : return gerepileupto(av, Z_to_padic(y, p, e));
1691 : }
1692 :
1693 : static GEN Qp_exp_safe(GEN x);
1694 :
1695 : /*compute the p^e th root of x p-adic, assume x != 0 */
1696 : static GEN
1697 854 : Qp_sqrtn_ram(GEN x, long e)
1698 : {
1699 854 : pari_sp ltop=avma;
1700 854 : GEN a, p = gel(x,2), n = powiu(p,e);
1701 854 : long v = valp(x), va;
1702 854 : if (v)
1703 : {
1704 : long z;
1705 161 : v = sdivsi_rem(v, n, &z);
1706 161 : if (z) return NULL;
1707 91 : x = leafcopy(x);
1708 91 : setvalp(x,0);
1709 : }
1710 : /*If p = 2, -1 is a root of 1 in U1: need extra check*/
1711 784 : if (absequaliu(p, 2) && mod8(gel(x,4)) != 1) return NULL;
1712 749 : a = Qp_log(x);
1713 749 : va = valp(a) - e;
1714 749 : if (va <= 0)
1715 : {
1716 287 : if (signe(gel(a,4))) return NULL;
1717 : /* all accuracy lost */
1718 119 : a = cvtop(remii(gel(x,4),p), p, 1);
1719 : }
1720 : else
1721 : {
1722 462 : setvalp(a, va); /* divide by p^e */
1723 462 : a = Qp_exp_safe(a);
1724 462 : if (!a) return NULL;
1725 : /* n=p^e and a^n=z*x where z is a (p-1)th-root of 1.
1726 : * Since z^n=z, we have (a/z)^n = x. */
1727 462 : a = gdiv(x, powgi(a,subiu(n,1))); /* = a/z = x/a^(n-1)*/
1728 462 : if (v) setvalp(a,v);
1729 : }
1730 581 : return gerepileupto(ltop,a);
1731 : }
1732 :
1733 : /*compute the nth root of x p-adic p prime with n*/
1734 : static GEN
1735 2037 : Qp_sqrtn_unram(GEN x, GEN n, GEN *zetan)
1736 : {
1737 : pari_sp av;
1738 2037 : GEN Z, a, r, p = gel(x,2);
1739 2037 : long v = valp(x);
1740 2037 : if (v)
1741 : {
1742 : long z;
1743 84 : v = sdivsi_rem(v,n,&z);
1744 84 : if (z) return NULL;
1745 : }
1746 2030 : r = cgetp(x); setvalp(r,v);
1747 2030 : Z = NULL; /* -Wall */
1748 2030 : if (zetan) Z = cgetp(x);
1749 2030 : av = avma; a = Fp_sqrtn(gel(x,4), n, p, zetan);
1750 2030 : if (!a) return NULL;
1751 2016 : affii(Zp_sqrtnlift(gel(x,4), n, a, p, precp(x)), gel(r,4));
1752 2016 : if (zetan)
1753 : {
1754 14 : affii(Zp_sqrtnlift(gen_1, n, *zetan, p, precp(x)), gel(Z,4));
1755 14 : *zetan = Z;
1756 : }
1757 2016 : return gc_const(av,r);
1758 : }
1759 :
1760 : GEN
1761 2604 : Qp_sqrtn(GEN x, GEN n, GEN *zetan)
1762 : {
1763 : pari_sp av, tetpil;
1764 : GEN q, p;
1765 : long e;
1766 2604 : if (absequaliu(n, 2))
1767 : {
1768 70 : if (zetan) *zetan = gen_m1;
1769 70 : if (signe(n) < 0) x = ginv(x);
1770 63 : return Qp_sqrt(x);
1771 : }
1772 2534 : av = avma; p = gel(x,2);
1773 2534 : if (!signe(gel(x,4)))
1774 : {
1775 203 : if (signe(n) < 0) pari_err_INV("Qp_sqrtn", x);
1776 203 : q = divii(addis(n, valp(x)-1), n);
1777 203 : if (zetan) *zetan = gen_1;
1778 203 : set_avma(av); return zeropadic(p, itos(q));
1779 : }
1780 : /* treat the ramified part using logarithms */
1781 2331 : e = Z_pvalrem(n, p, &q);
1782 2331 : if (e) { x = Qp_sqrtn_ram(x,e); if (!x) return NULL; }
1783 2058 : if (is_pm1(q))
1784 : { /* finished */
1785 21 : if (signe(q) < 0) x = ginv(x);
1786 21 : x = gerepileupto(av, x);
1787 21 : if (zetan)
1788 28 : *zetan = (e && absequaliu(p, 2))? gen_m1 /*-1 in Q_2*/
1789 28 : : gen_1;
1790 21 : return x;
1791 : }
1792 2037 : tetpil = avma;
1793 : /* use hensel lift for unramified case */
1794 2037 : x = Qp_sqrtn_unram(x, q, zetan);
1795 2037 : if (!x) return NULL;
1796 2016 : if (zetan)
1797 : {
1798 : GEN *gptr[2];
1799 14 : if (e && absequaliu(p, 2))/*-1 in Q_2*/
1800 : {
1801 7 : tetpil = avma; x = gcopy(x); *zetan = gneg(*zetan);
1802 : }
1803 14 : gptr[0] = &x; gptr[1] = zetan;
1804 14 : gerepilemanysp(av,tetpil,gptr,2);
1805 14 : return x;
1806 : }
1807 2002 : return gerepile(av,tetpil,x);
1808 : }
1809 :
1810 : GEN
1811 24347 : sqrtnint(GEN a, long n)
1812 : {
1813 24347 : pari_sp av = avma;
1814 : GEN x, b, q;
1815 : long s, k, e;
1816 24347 : const ulong nm1 = n - 1;
1817 24347 : if (n == 2) return sqrtint(a);
1818 20105 : if (typ(a) != t_INT)
1819 : {
1820 35 : if (typ(a) == t_REAL)
1821 : {
1822 : long e;
1823 14 : switch(signe(a))
1824 : {
1825 0 : case 0: return gen_0;
1826 7 : case -1: pari_err_DOMAIN("sqrtnint", "argument", "<", gen_0,a);
1827 : }
1828 7 : e = expo(a); if (e < 0) return gen_0;
1829 7 : if (nbits2lg(e+1) > lg(a))
1830 0 : a = floorr(sqrtnr(a,n)); /* try to avoid precision loss in truncation */
1831 : else
1832 7 : a = sqrtnint(truncr(a),n);
1833 : }
1834 : else
1835 : {
1836 21 : GEN b = gfloor(a);
1837 21 : if (typ(b) != t_INT) pari_err_TYPE("sqrtint",a);
1838 14 : if (signe(b) < 0) pari_err_DOMAIN("sqrtnint", "argument", "<", gen_0,b);
1839 7 : a = sqrtnint(b, n);
1840 : }
1841 14 : return gerepileuptoint(av, a);
1842 : }
1843 20070 : if (n <= 0) pari_err_DOMAIN("sqrtnint", "n", "<=", gen_0, stoi(n));
1844 20063 : if (n == 1) return icopy(a);
1845 17907 : s = signe(a);
1846 17907 : if (s < 0) pari_err_DOMAIN("sqrtnint", "x", "<", gen_0, a);
1847 17907 : if (!s) return gen_0;
1848 17830 : if (lgefint(a) == 3) return utoi(usqrtn(itou(a), n));
1849 11606 : e = expi(a); k = e/(2*n);
1850 11606 : if (k == 0)
1851 : {
1852 : long flag;
1853 291 : if (n > e) return gc_const(av, gen_1);
1854 291 : flag = cmpii(a, powuu(3, n)); set_avma(av);
1855 291 : return (flag < 0) ? gen_2: stoi(3);
1856 : }
1857 11315 : if (e < n*BITS_IN_LONG - 1)
1858 : {
1859 : ulong xs, qs;
1860 4181 : b = itor(a, (2*e < n*BITS_IN_LONG)? DEFAULTPREC: MEDDEFAULTPREC);
1861 4181 : x = mpexp(divru(logr_abs(b), n));
1862 4181 : xs = itou(floorr(x)) + 1; /* >= a^(1/n) */
1863 : for(;;) {
1864 8184 : q = divii(a, powuu(xs, nm1));
1865 8184 : if (lgefint(q) > 3) break;
1866 8177 : qs = itou(q); if (qs >= xs) break;
1867 4003 : xs -= (xs - qs + nm1)/n;
1868 : }
1869 4181 : return utoi(xs);
1870 : }
1871 7134 : b = addui(1, shifti(a, -n*k));
1872 7134 : x = shifti(addui(1, sqrtnint(b, n)), k);
1873 7134 : q = divii(a, powiu(x, nm1));
1874 15994 : while (cmpii(q, x) < 0) /* a priori one iteration, no GC necessary */
1875 : {
1876 8860 : x = subii(x, divis(addui(nm1, subii(x, q)), n));
1877 8860 : q = divii(a, powiu(x, nm1));
1878 : }
1879 7134 : return gerepileuptoleaf(av, x);
1880 : }
1881 :
1882 : ulong
1883 8107 : usqrtn(ulong a, ulong n)
1884 : {
1885 : ulong x, s, q;
1886 8107 : const ulong nm1 = n - 1;
1887 8107 : if (!n) pari_err_DOMAIN("sqrtnint", "n", "=", gen_0, utoi(n));
1888 8107 : if (n == 1 || a == 0) return a;
1889 8107 : s = 1 + expu(a)/n; x = 1UL << s;
1890 8107 : q = (nm1*s >= BITS_IN_LONG)? 0: a >> (nm1*s);
1891 21045 : while (q < x) {
1892 : ulong X;
1893 12938 : x -= (x - q + nm1)/n;
1894 12938 : X = upowuu(x, nm1);
1895 12938 : q = X? a/X: 0;
1896 : }
1897 8107 : return x;
1898 : }
1899 :
1900 : static ulong
1901 1725058 : cubic_prec_mask(long n)
1902 : {
1903 1725058 : long a = n, i;
1904 1725058 : ulong mask = 0;
1905 1725058 : for(i = 1;; i++, mask *= 3)
1906 8206721 : {
1907 9931779 : long c = a%3;
1908 9931779 : if (c) mask += 3 - c;
1909 9931779 : a = (a+2)/3;
1910 9931779 : if (a==1) return mask + upowuu(3, i);
1911 : }
1912 : }
1913 :
1914 : /* cubic Newton iteration, |a|^(1/n), assuming a != 0 */
1915 : GEN
1916 2743657 : sqrtnr_abs(GEN a, long n)
1917 : {
1918 : pari_sp av;
1919 : GEN x, b;
1920 : long eextra, eold, n1, n2, prec, B, v;
1921 : ulong mask;
1922 :
1923 2743657 : if (n == 1) return mpabs(a);
1924 2743013 : if (n == 2) return sqrtr_abs(a);
1925 :
1926 2446260 : prec = realprec(a); v = expo(a) / n; av = avma;
1927 2446260 : if (v) a = shiftr(a, -n*v);
1928 2446275 : b = rtor(a, DEFAULTPREC);
1929 2446296 : x = mpexp(divru(logr_abs(b), n));
1930 2446299 : if (prec == DEFAULTPREC)
1931 : {
1932 758651 : if (v) shiftr_inplace(x, v);
1933 758651 : return gerepileuptoleaf(av, x);
1934 : }
1935 1687648 : n1 = n+1;
1936 1687648 : n2 = 2*n;
1937 1687648 : B = prec2nbits(prec);
1938 1687648 : eextra = expu(n)-1;
1939 1687648 : mask = cubic_prec_mask(B + 63);
1940 1687648 : eold = 1;
1941 : for(;;)
1942 6728678 : { /* reach 64 */
1943 8416326 : long enew = eold * 3;
1944 8416326 : enew -= mask % 3;
1945 8416326 : if (enew > 64) break; /* back up one step */
1946 6728678 : mask /= 3;
1947 6728678 : eold = enew;
1948 : }
1949 : for(;;)
1950 1315442 : {
1951 3003090 : long pr, enew = eold * 3;
1952 : GEN y, z;
1953 3003090 : enew -= mask % 3;
1954 3003090 : mask /= 3;
1955 3003090 : pr = nbits2prec(enew + eextra);
1956 3003090 : b = rtor(a, pr); setsigne(b,1);
1957 3003090 : x = rtor(x, pr);
1958 3003090 : y = subrr(powru(x, n), b);
1959 3003090 : z = divrr(y, addrr(mulur(n1, y), mulur(n2, b)));
1960 3003090 : shiftr_inplace(z,1);
1961 3003090 : x = mulrr(x, subsr(1,z));
1962 3003090 : if (mask == 1)
1963 : {
1964 1687648 : if (v) shiftr_inplace(x, v);
1965 1687648 : return gerepileuptoleaf(av, gprec_wtrunc(x,prec));
1966 : }
1967 1315442 : eold = enew;
1968 : }
1969 : }
1970 :
1971 : static void
1972 53896 : shiftc_inplace(GEN z, long d)
1973 : {
1974 53896 : shiftr_inplace(gel(z,1), d);
1975 53896 : shiftr_inplace(gel(z,2), d);
1976 53896 : }
1977 :
1978 : /* exp(2*Pi*I/n), same iteration as sqrtnr_abs, different initial point */
1979 : static GEN
1980 543311 : sqrtnof1(ulong n, long prec)
1981 : {
1982 : pari_sp av;
1983 : GEN x;
1984 : long eold, n1, n2, B;
1985 : ulong mask;
1986 :
1987 543311 : B = prec2nbits(prec);
1988 543311 : n1 = n+1;
1989 543311 : n2 = 2*n; av = avma;
1990 :
1991 543311 : x = expIr(divru(Pi2n(1, LOWDEFAULTPREC), n));
1992 543311 : if (prec == LOWDEFAULTPREC) return gerepileupto(av, x);
1993 37410 : mask = cubic_prec_mask(B + BITS_IN_LONG-1);
1994 37410 : eold = 1;
1995 : for(;;)
1996 146115 : { /* reach BITS_IN_LONG */
1997 183525 : long enew = eold * 3;
1998 183525 : enew -= mask % 3;
1999 183525 : if (enew > BITS_IN_LONG) break; /* back up one step */
2000 146115 : mask /= 3;
2001 146115 : eold = enew;
2002 : }
2003 : for(;;)
2004 16486 : {
2005 53896 : long pr, enew = eold * 3;
2006 : GEN y, z;
2007 53896 : enew -= mask % 3;
2008 53896 : mask /= 3;
2009 53896 : pr = nbits2prec(enew);
2010 53896 : x = cxtofp(x, pr);
2011 53896 : y = gsub(gpowgs(x, n), gen_1);
2012 53896 : z = gdiv(y, gaddgs(gmulsg(n1, y), n2));
2013 53896 : shiftc_inplace(z,1);
2014 53896 : x = gmul(x, gsubsg(1, z));
2015 53896 : if (mask == 1) return gerepilecopy(av, gprec_w(x,prec));
2016 16486 : eold = enew;
2017 : }
2018 : }
2019 :
2020 : /* exp(2iPi/d) */
2021 : GEN
2022 2060133 : rootsof1u_cx(ulong n, long prec)
2023 : {
2024 2060133 : switch(n)
2025 : {
2026 13006 : case 1: return gen_1;
2027 2884 : case 2: return gen_m1;
2028 652924 : case 4: return gen_I();
2029 40092 : case 3: case 6: case 12:
2030 : {
2031 40092 : pari_sp av = avma;
2032 40092 : GEN a = (n == 3)? mkfrac(gen_m1,gen_2): ghalf;
2033 40092 : GEN sq3 = sqrtr_abs(utor(3, prec));
2034 40092 : shiftr_inplace(sq3, -1);
2035 40092 : a = (n == 12)? mkcomplex(sq3, a): mkcomplex(a, sq3);
2036 40092 : return gerepilecopy(av, a);
2037 : }
2038 807923 : case 8:
2039 : {
2040 807923 : pari_sp av = avma;
2041 807923 : GEN sq2 = sqrtr_abs(utor(2, prec));
2042 807898 : shiftr_inplace(sq2,-1);
2043 807906 : return gerepilecopy(av, mkcomplex(sq2, sq2));
2044 : }
2045 : }
2046 543304 : return sqrtnof1(n, prec);
2047 : }
2048 : /* e(a/b) */
2049 : GEN
2050 14154 : rootsof1q_cx(long a, long b, long prec)
2051 : {
2052 14154 : long g = cgcd(a,b);
2053 : GEN z;
2054 14154 : if (g != 1) { a /= g; b /= g; }
2055 14154 : if (b < 0) { b = -b; a = -a; }
2056 14154 : z = rootsof1u_cx(b, prec);
2057 14154 : if (a < 0) { z = conj_i(z); a = -a; }
2058 14154 : return gpowgs(z, a);
2059 : }
2060 :
2061 : /* initializes powers of e(a/b) */
2062 : GEN
2063 14987 : rootsof1powinit(long a, long b, long prec)
2064 : {
2065 14987 : long g = cgcd(a,b);
2066 14987 : if (g != 1) { a /= g; b /= g; }
2067 14987 : if (b < 0) { b = -b; a = -a; }
2068 14987 : a %= b; if (a < 0) a += b;
2069 14987 : return mkvec2(grootsof1(b,prec), mkvecsmall2(a,b));
2070 : }
2071 : /* T = rootsof1powinit(a,b); return e(a/b)^c */
2072 : GEN
2073 12516441 : rootsof1pow(GEN T, long c)
2074 : {
2075 12516441 : GEN vz = gel(T,1), ab = gel(T,2);
2076 12516441 : long a = ab[1], b = ab[2]; /* a >= 0, b > 0 */
2077 12516441 : c %= b; if (c < 0) c += b;
2078 12516441 : a = Fl_mul(a, c, b);
2079 12516441 : return gel(vz, a + 1);
2080 : }
2081 :
2082 : /* exp(2iPi/d), assume d a t_INT */
2083 : GEN
2084 4480 : rootsof1_cx(GEN d, long prec)
2085 : {
2086 4480 : if (lgefint(d) == 3) return rootsof1u_cx((ulong)d[2], prec);
2087 0 : return expIr(divri(Pi2n(1,prec), d));
2088 : }
2089 :
2090 : GEN
2091 41974 : gsqrtn(GEN x, GEN n, GEN *zetan, long prec)
2092 : {
2093 : long i, lx, tx;
2094 : pari_sp av;
2095 : GEN y, z;
2096 41974 : if (typ(n)!=t_INT) pari_err_TYPE("sqrtn",n);
2097 41974 : if (!signe(n)) pari_err_DOMAIN("sqrtn", "n", "=", gen_0, n);
2098 41974 : if (is_pm1(n))
2099 : {
2100 70 : if (zetan) *zetan = gen_1;
2101 70 : return (signe(n) > 0)? gcopy(x): ginv(x);
2102 : }
2103 41904 : if (zetan) *zetan = gen_0;
2104 41904 : tx = typ(x);
2105 41904 : if (is_matvec_t(tx))
2106 : {
2107 7 : y = cgetg_copy(x, &lx);
2108 21 : for (i=1; i<lx; i++) gel(y,i) = gsqrtn(gel(x,i),n,NULL,prec);
2109 7 : return y;
2110 : }
2111 41897 : av = avma;
2112 41897 : switch(tx)
2113 : {
2114 182 : case t_INTMOD:
2115 : {
2116 182 : GEN p = gel(x,1), s;
2117 182 : z = gen_0;
2118 182 : y = cgetg(3,t_INTMOD); gel(y,1) = icopy(p);
2119 182 : if (zetan) { z = cgetg(3,t_INTMOD); gel(z,1) = gel(y,1); }
2120 182 : s = Fp_sqrtn(gel(x,2),n,p,zetan);
2121 161 : if (!s) {
2122 35 : if (zetan) return gc_const(av,gen_0);
2123 28 : if (!BPSW_psp(p)) pari_err_PRIME("sqrtn [modulus]",p);
2124 14 : pari_err_SQRTN("gsqrtn",x);
2125 : }
2126 126 : gel(y,2) = s;
2127 126 : if (zetan) { gel(z,2) = *zetan; *zetan = z; }
2128 126 : return y;
2129 : }
2130 :
2131 56 : case t_PADIC:
2132 56 : y = Qp_sqrtn(x,n,zetan);
2133 49 : if (!y) {
2134 7 : if (zetan) return gen_0;
2135 7 : pari_err_SQRTN("gsqrtn",x);
2136 : }
2137 42 : return y;
2138 :
2139 196 : case t_FFELT: return FF_sqrtn(x,n,zetan);
2140 :
2141 41001 : case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
2142 41001 : i = precision(x); if (i) prec = i;
2143 41001 : if (isint1(x))
2144 7 : y = real_1(prec);
2145 40994 : else if (gequal0(x))
2146 : {
2147 : long b;
2148 21 : if (signe(n) < 0) pari_err_INV("gsqrtn",x);
2149 21 : if (isinexactreal(x))
2150 14 : b = sdivsi(gexpo(x), n);
2151 : else
2152 7 : b = -prec2nbits(prec);
2153 21 : if (typ(x) == t_COMPLEX)
2154 : {
2155 7 : y = cgetg(3,t_COMPLEX);
2156 7 : gel(y,1) = gel(y,2) = real_0_bit(b);
2157 : }
2158 : else
2159 14 : y = real_0_bit(b);
2160 : }
2161 : else
2162 : {
2163 40973 : long nn = itos_or_0(n);
2164 40973 : if (tx == t_INT) { x = itor(x,prec); tx = t_REAL; }
2165 40973 : if (nn > 0 && tx == t_REAL && signe(x) > 0)
2166 31023 : y = sqrtnr(x, nn);
2167 : else
2168 9950 : y = gexp(gdiv(glog(x,prec), n), prec);
2169 40973 : y = gerepileupto(av, y);
2170 : }
2171 41001 : if (zetan) *zetan = rootsof1_cx(n, prec);
2172 41001 : return y;
2173 :
2174 7 : case t_QUAD:
2175 7 : return gsqrtn(quadtofp(x, prec), n, zetan, prec);
2176 :
2177 455 : default:
2178 455 : av = avma; if (!(y = toser_i(x))) break;
2179 455 : return gerepileupto(av, ser_powfrac(y, ginv(n), prec));
2180 : }
2181 0 : pari_err_TYPE("sqrtn",x);
2182 : return NULL;/* LCOV_EXCL_LINE */
2183 : }
2184 :
2185 : /********************************************************************/
2186 : /** **/
2187 : /** EXP(X) - 1 **/
2188 : /** **/
2189 : /********************************************************************/
2190 : /* exp(|x|) - 1, assume x != 0.
2191 : * For efficiency, x should be reduced mod log(2): if so, we have a < 0 */
2192 : GEN
2193 19207884 : exp1r_abs(GEN x)
2194 : {
2195 19207884 : long l = realprec(x), a = expo(x), b = prec2nbits(l), L, i, n, m, B;
2196 : GEN y, p2, X;
2197 : pari_sp av;
2198 : double d;
2199 :
2200 19207643 : if (b + a <= 0) return mpabs(x);
2201 :
2202 19192258 : y = cgetr(l); av = avma;
2203 19191395 : B = b/3 + BITS_IN_LONG + (BITS_IN_LONG*BITS_IN_LONG)/ b;
2204 19191395 : d = a/2.; m = (long)(d + sqrt(d*d + B)); /* >= 0 */
2205 19191395 : if (m < (-a) * 0.1) m = 0; /* not worth it */
2206 : /* Multiplication is quadratic in this range (l is small, otherwise we
2207 : * use logAGM + Newton). Set Y = 2^(-e-a) x, compute truncated series
2208 : * sum_{k <= n} Y^k/k!: this costs roughly
2209 : * m b^2 + sum_{k <= n} (k e + BITS_IN_LONG)^2
2210 : * bit operations with n ~ b/e, |x| < 2^(1+a), |Y| < 2^(1-e) , m = e+a and
2211 : * b bits of accuracy needed, so
2212 : * B := (b / 3 + BITS_IN_LONG + BITS_IN_LONG^2 / b) ~ m(m-a)
2213 : * we want b ~ 3 m (m-a) or m~b+a hence
2214 : * m = min( a/2 + sqrt(a^2/4 + B), b + a )
2215 : * NB: e ~ (b/3)^(1/2) as b -> oo
2216 : *
2217 : * Truncate the sum at k = n (>= 1), the remainder is
2218 : * sum_{k >= n+1} Y^k / k! < Y^(n+1) / (n+1)! (1-Y) < Y^(n+1) / n!
2219 : * We want Y^(n+1) / n! <= Y 2^-b, hence -n log_2 |Y| + log_2 n! >= b
2220 : * log n! ~ (n + 1/2) log(n+1) - (n+1) + log(2Pi)/2,
2221 : * error bounded by 1/6(n+1) <= 1/12. Finally, we want
2222 : * n (-1/log(2) -log_2 |Y| + log_2(n+1)) >= b */
2223 19191395 : d = m-dbllog2(x)-1/M_LN2; /* ~ -log_2 Y - 1/log(2) */
2224 19192504 : while (d <= 0) { d++; m++; } /* d < 0 can occur from expm1 */
2225 19192498 : L = l + nbits2extraprec(m);
2226 19192352 : b += m;
2227 19192352 : n = (long)(b / d); /* > 0 */
2228 19192352 : if (n == 1)
2229 743920 : n = (long)(b / (d + log2((double)n+1))); /* log ~ const in small ranges */
2230 20552005 : while (n*(d+log2((double)n+1)) < b) n++; /* expect few corrections */
2231 :
2232 19192352 : X = rtor(x,L); shiftr_inplace(X, -m); setsigne(X, 1);
2233 19193402 : if (n == 1) p2 = X;
2234 : else
2235 : {
2236 19193402 : long s = 0, l1 = nbits2prec((long)(d + n + 16));
2237 19193366 : GEN unr = real_1(L);
2238 : pari_sp av2;
2239 :
2240 19192935 : p2 = cgetr(L); av2 = avma;
2241 352881528 : for (i=n; i>=2; i--, set_avma(av2))
2242 : { /* compute X^(n-1)/n! + ... + X/2 + 1 */
2243 : GEN p1, p3;
2244 333784561 : setprec(X,l1); p3 = divru(X,i);
2245 334258301 : l1 += nbits2extraprec(dvmdsBIL(s - expo(p3), &s)<<TWOPOTBITS_IN_LONG);
2246 334164701 : if (l1>L) l1=L;
2247 334164701 : setprec(unr,l1); p1 = addrr_sign(unr,1, i == n? p3: mulrr(p3,p2),1);
2248 333526998 : setprec(p2,l1); affrr(p1,p2); /* p2 <- 1 + (X/i)*p2 */
2249 : }
2250 19191616 : setprec(X,L); p2 = mulrr(X,p2);
2251 : }
2252 :
2253 19193686 : B = prec2nbits(L);
2254 205301736 : for (i = 1; i <= m; i++)
2255 : {
2256 186108570 : if (realprec(p2) > L) setprec(p2,L);
2257 186108570 : if (expo(p2) < -B)
2258 0 : shiftr_inplace(p2, 1); /* 2 + p2 ~ 2 and may blow up accuracy */
2259 : else
2260 186108570 : p2 = mulrr(p2, addsr(2,p2));
2261 : }
2262 19193166 : affrr_fixlg(p2,y); return gc_const(av,y);
2263 : }
2264 :
2265 : GEN
2266 10948 : mpexpm1(GEN x)
2267 : {
2268 10948 : const long s = 6;
2269 10948 : long B, l, sx = signe(x);
2270 : GEN y, z;
2271 : pari_sp av;
2272 10948 : if (!sx) return real_0_bit(expo(x));
2273 10941 : l = realprec(x);
2274 10941 : if (l > maxss(EXPNEWTON_LIMIT, BITS_IN_LONG<<s))
2275 : {
2276 6 : long e = expo(x);
2277 6 : if (e < 0) x = rtor(x, l + nbits2extraprec(-e));
2278 6 : return subrs(mpexp(x), 1);
2279 : }
2280 10935 : if (sx > 0) return exp1r_abs(x);
2281 5142 : B = prec2nbits(l);
2282 5142 : if (cmpsr(-B, x) > 0) return real_m1(l);
2283 : /* compute exp(x) * (1 - exp(-x)) */
2284 5135 : av = avma; y = exp1r_abs(x); /* > 0 */
2285 5135 : if (expo(y) >= -B) { z = addsr(1, y); y = divrr(y, z); }
2286 5135 : setsigne(y, -1);
2287 5135 : return gerepileuptoleaf(av, y);
2288 : }
2289 :
2290 : static GEN serexp(GEN x, long prec);
2291 : GEN
2292 12763 : gexpm1(GEN x, long prec)
2293 : {
2294 12763 : switch(typ(x))
2295 : {
2296 4219 : case t_REAL: return mpexpm1(x);
2297 6438 : case t_COMPLEX: return cxexpm1(x,prec);
2298 14 : case t_PADIC: return gsubgs(Qp_exp(x), 1);
2299 2092 : default:
2300 : {
2301 2092 : pari_sp av = avma;
2302 : long ey;
2303 : GEN y;
2304 2092 : if (!(y = toser_i(x))) break;
2305 2071 : ey = valser(y);
2306 2071 : if (ey < 0) pari_err_DOMAIN("expm1","valuation", "<", gen_0, x);
2307 2071 : if (gequal0(y)) return gcopy(y);
2308 2064 : if (ey)
2309 504 : return gerepileupto(av, gsubgs(serexp(y,prec), 1));
2310 : else
2311 : {
2312 1560 : GEN e1 = gexpm1(gel(y,2), prec), e = gaddgs(e1,1);
2313 1560 : y = gmul(e, serexp(serchop0(y),prec));
2314 1560 : gel(y,2) = e1;
2315 1560 : return gerepilecopy(av, y);
2316 : }
2317 : }
2318 : }
2319 21 : return trans_eval("expm1",gexpm1,x,prec);
2320 : }
2321 : /********************************************************************/
2322 : /** **/
2323 : /** EXP(X) **/
2324 : /** **/
2325 : /********************************************************************/
2326 : static GEN
2327 19145205 : mpexp_basecase(GEN x)
2328 : {
2329 19145205 : pari_sp av = avma;
2330 19145205 : long sh, l = realprec(x);
2331 : GEN y, z;
2332 :
2333 19145205 : y = modlog2(x, &sh);
2334 19144822 : if (!y) { set_avma(av); return real2n(sh, l); }
2335 19144822 : z = addsr(1, exp1r_abs(y));
2336 19144383 : if (signe(y) < 0) z = invr(z);
2337 19144480 : if (sh) {
2338 15486185 : shiftr_inplace(z, sh);
2339 15486095 : if (realprec(z) > l) z = rtor(z, l); /* spurious precision increase */
2340 : }
2341 : #ifdef DEBUG
2342 : {
2343 : GEN t = mplog(z), u = divrr(subrr(x, t),x);
2344 : if (signe(u) && expo(u) > 5-prec2nbits(minss(l,realprec(t))))
2345 : pari_err_BUG("exp");
2346 : }
2347 : #endif
2348 19144701 : return gerepileuptoleaf(av, z); /* NOT affrr, precision often increases */
2349 : }
2350 :
2351 : GEN
2352 19285075 : mpexp(GEN x)
2353 : {
2354 19285075 : const long s = 6; /*Initial steps using basecase*/
2355 19285075 : long i, p, l = realprec(x), sh;
2356 : GEN a, t, z;
2357 : ulong mask;
2358 :
2359 19285075 : if (l <= maxss(EXPNEWTON_LIMIT, (BITS_IN_LONG<<s) + 2))
2360 : {
2361 19285205 : if (!signe(x)) return mpexp0(x);
2362 19145115 : return mpexp_basecase(x);
2363 : }
2364 11 : z = cgetr(l); /* room for result */
2365 13 : x = modlog2(x, &sh);
2366 13 : if (!x) { set_avma((pari_sp)(z+lg(z))); return real2n(sh, l); }
2367 13 : constpi(l); /* precompute for later logr_abs() */
2368 13 : mask = quadratic_prec_mask(prec2nbits(l)+BITS_IN_LONG);
2369 168 : for(i=0, p=1; i<s+TWOPOTBITS_IN_LONG; i++) { p <<= 1; if (mask & 1) p-=1; mask >>= 1; }
2370 13 : a = mpexp_basecase(rtor(x, nbits2prec(p)));
2371 13 : x = addrs(x,1);
2372 13 : if (realprec(x) < l+EXTRAPREC64) x = rtor(x, l+EXTRAPREC64);
2373 13 : a = rtor(a, l+EXTRAPREC64); /*append 0s */
2374 13 : t = NULL;
2375 : for(;;)
2376 : {
2377 14 : p <<= 1; if (mask & 1) p--;
2378 14 : mask >>= 1;
2379 14 : setprec(x, nbits2prec(p));
2380 14 : setprec(a, nbits2prec(p));
2381 14 : t = mulrr(a, subrr(x, logr_abs(a))); /* a (x - log(a)) */
2382 14 : if (mask == 1) break;
2383 1 : affrr(t, a); set_avma((pari_sp)a);
2384 : }
2385 13 : affrr(t,z);
2386 13 : if (sh) shiftr_inplace(z, sh);
2387 13 : return gc_const((pari_sp)z, z);
2388 : }
2389 :
2390 : /* x != 0; k = ceil(tn / (te-1)), t = p-1 */
2391 : long
2392 98 : Qp_exp_prec(GEN x)
2393 : {
2394 98 : long e = valp(x), n = precp(x);
2395 : ulong a, b, q, r, p, t;
2396 :
2397 98 : if (e < 1) return -1;
2398 77 : if (e > n) return 1;
2399 77 : p = itos_or_0(gel(x,2));
2400 77 : if (!p) return n / e + 1;
2401 77 : if (p == 2) return e < 2? -1: ceildivuu(n, e - 1);
2402 : /* n >= e > 0, n = qe + r */
2403 : /* tn = q (te-1) + rt + q = (q+1)(te-1) - t(e-r) + q + 1 */
2404 63 : t = p - 1;
2405 63 : if (e == 1) return n + ceildivuu(n, t - 1);
2406 0 : q = n / e;
2407 0 : r = n % e; /* k = q + 1 if rt + q < te */
2408 0 : a = umuluu_or_0(e - r, t); if (!a || a > q) return q + 1;
2409 0 : b = umuluu_or_0(e, t); if (!b) return q + 2;
2410 0 : return q + 1 + ceildivuu(q + 1 - a, b - 1);
2411 : }
2412 :
2413 : static GEN
2414 108727 : Qp_exp_safe(GEN x)
2415 : {
2416 108727 : pari_sp av = avma;
2417 108727 : GEN p = gel(x,2), a = gel(x,4), z;
2418 108727 : long d = precp(x), v = valp(x), e = d+v;
2419 108727 : if (gequal0(x)) return gaddgs(x,1);
2420 107222 : if (v < (equaliu(p,2)? 2:1)) return NULL;
2421 107215 : z = Zp_exp(mulii(a,powiu(p,v)), p, e);
2422 107215 : return gerepileupto(av, Z_to_padic(z, p, e));
2423 : }
2424 :
2425 : GEN
2426 108265 : Qp_exp(GEN x)
2427 : {
2428 108265 : GEN y = Qp_exp_safe(x);
2429 108263 : if (!y) pari_err_DOMAIN("gexp(t_PADIC)","argument","",gen_0,x);
2430 108256 : return y;
2431 : }
2432 :
2433 : static GEN
2434 49 : cos_p(GEN x)
2435 : {
2436 : long k;
2437 : pari_sp av;
2438 : GEN x2, y;
2439 :
2440 49 : if (gequal0(x)) return gaddgs(x,1);
2441 28 : k = Qp_exp_prec(x);
2442 28 : if (k < 0) return NULL;
2443 21 : av = avma; x2 = gsqr(x);
2444 21 : if (k & 1) k--;
2445 105 : for (y=gen_1; k; k-=2)
2446 : {
2447 84 : GEN t = gdiv(gmul(y,x2), muluu(k, k-1));
2448 84 : y = gsubsg(1, t);
2449 : }
2450 21 : return gerepileupto(av, y);
2451 : }
2452 : static GEN
2453 63 : sin_p(GEN x)
2454 : {
2455 : long k;
2456 : pari_sp av;
2457 : GEN x2, y;
2458 :
2459 63 : if (gequal0(x)) return gcopy(x);
2460 42 : k = Qp_exp_prec(x);
2461 42 : if (k < 0) return NULL;
2462 28 : av = avma; x2 = gsqr(x);
2463 28 : if (k & 1) k--;
2464 133 : for (y=gen_1; k; k-=2)
2465 : {
2466 105 : GEN t = gdiv(gmul(y,x2), muluu(k, k+1));
2467 105 : y = gsubsg(1, t);
2468 : }
2469 28 : return gerepileupto(av, gmul(y, x));
2470 : }
2471 :
2472 : static GEN
2473 4166463 : cxexp(GEN x, long prec)
2474 : {
2475 4166463 : GEN r, p1, p2, y = cgetg(3,t_COMPLEX);
2476 4166425 : pari_sp av = avma, tetpil;
2477 : long l;
2478 4166425 : l = precision(x); if (l > prec) prec = l;
2479 4166519 : if (gequal0(gel(x,1)))
2480 : {
2481 347938 : gsincos(gel(x,2),&gel(y,2),&gel(y,1),prec);
2482 347955 : return y;
2483 : }
2484 3818565 : r = gexp(gel(x,1),prec);
2485 3818718 : gsincos(gel(x,2),&p2,&p1,prec);
2486 3818908 : tetpil = avma;
2487 3818908 : gel(y,1) = gmul(r,p1);
2488 3818787 : gel(y,2) = gmul(r,p2);
2489 3818798 : gerepilecoeffssp(av,tetpil,y+1,2);
2490 3818868 : return y;
2491 : }
2492 :
2493 : /* given a t_SER x^v s(x), with s(0) != 0, return x^v(s - s(0)), shallow */
2494 : GEN
2495 37377 : serchop0(GEN s)
2496 : {
2497 37377 : long i, l = lg(s);
2498 : GEN y;
2499 37377 : if (l == 2) return s;
2500 37377 : if (l == 3 && isexactzero(gel(s,2))) return s;
2501 37377 : y = cgetg(l, t_SER); y[1] = s[1];
2502 164953 : gel(y,2) = gen_0; for (i=3; i <l; i++) gel(y,i) = gel(s,i);
2503 37377 : return normalizeser(y);
2504 : }
2505 :
2506 : GEN
2507 42 : serchop_i(GEN s, long n)
2508 : {
2509 42 : long i, m, l = lg(s);
2510 : GEN y;
2511 42 : if (l == 2 || (l == 3 && isexactzero(gel(s,2))))
2512 : {
2513 14 : if (valser(s) < n) { s = shallowcopy(s); setvalser(s,n); }
2514 14 : return s;
2515 : }
2516 28 : m = n - valser(s); if (m < 0) return s;
2517 21 : if (l-m <= 2) return zeroser(varn(s), n);
2518 14 : y = cgetg(l-m, t_SER); y[1] = s[1]; setvalser(y, valser(y)+m);
2519 42 : for (i=m+2; i < l; i++) gel(y,i-m) = gel(s,i);
2520 14 : return normalizeser(y);
2521 : }
2522 : GEN
2523 42 : serchop(GEN s, long n)
2524 : {
2525 42 : pari_sp av = avma;
2526 42 : if (typ(s) != t_SER) pari_err_TYPE("serchop",s);
2527 42 : return gerepilecopy(av, serchop_i(s,n));
2528 : }
2529 :
2530 : static GEN
2531 82962 : serexp(GEN x, long prec)
2532 : {
2533 82962 : long i, j, lx, ly, mi, e = valser(x);
2534 : GEN y, xd, yd;
2535 : pari_sp av;
2536 :
2537 82962 : if (e < 0) pari_err_DOMAIN("exp","valuation", "<", gen_0, x);
2538 82955 : if (gequal0(x)) return gaddsg(1,x);
2539 70224 : lx = lg(x);
2540 70224 : if (e)
2541 : {
2542 : GEN X;
2543 55518 : ly = lx+e; y = cgetg(ly,t_SER);
2544 567057 : mi = lx-1; while (mi>=3 && isrationalzero(gel(x,mi))) mi--;
2545 55518 : mi += e-2;
2546 55518 : y[1] = evalsigne(1) | _evalvalser(0) | evalvarn(varn(x));
2547 : /* zd[i] = coefficient of X^i in z */
2548 55518 : xd = x+2-e; yd = y+2; ly -= 2;
2549 55518 : X = gel(xd,e); if (e != 1) X = gmulgu(X, e); /* left on stack */
2550 55518 : X = isint1(X)? NULL: X;
2551 55518 : gel(yd,0) = gen_1;
2552 55889 : for (i = 1; i < e; i++) gel(yd,i) = gen_0;
2553 664652 : for ( ; i < ly; i++)
2554 : {
2555 609134 : GEN t = gel(yd,i-e);
2556 609134 : long J = minss(i, mi);
2557 609134 : av = avma; if (X) t = gmul(t, X);
2558 2579410 : for (j = e + 1; j <= J; j++)
2559 1970276 : t = gadd(t, gmulgu(gmul(gel(xd,j),gel(yd,i-j)), j));
2560 609134 : gel(yd,i) = gerepileupto(av, gdivgu(t, i));
2561 : }
2562 55518 : return y;
2563 : }
2564 14706 : av = avma;
2565 14706 : return gerepileupto(av, gmul(gexp(gel(x,2),prec), serexp(serchop0(x),prec)));
2566 : }
2567 :
2568 : static GEN
2569 1468531 : expQ(GEN x, long prec)
2570 : {
2571 1468531 : GEN p, q, z, z0 = NULL;
2572 : pari_sp av;
2573 1468531 : long n, nmax, s, e, b = prec2nbits(prec);
2574 : double ex;
2575 : struct abpq_res R;
2576 : struct abpq S;
2577 :
2578 1468530 : if (typ(x) == t_INT)
2579 : {
2580 24640 : if (!signe(x)) return real_1(prec);
2581 24618 : p = x; q = gen_1;
2582 24618 : e = expi(p);
2583 24619 : if (e > b) return mpexp(itor(x, prec));
2584 : }
2585 : else
2586 : {
2587 1443890 : long ep, eq, B = usqrt(b) / 2;
2588 1443890 : p = gel(x,1); ep = expi(p);
2589 1443890 : q = gel(x,2); eq = expi(q);
2590 1443890 : if (ep > B || eq > B) return mpexp(fractor(x, prec));
2591 14637 : e = ep - eq;
2592 14637 : if (e < -3) prec += nbits2extraprec(-e); /* see addrr 'extend' rule */
2593 : }
2594 39256 : if (e > 2) { z0 = cgetr(prec); prec += EXTRAPREC64; b += BITS_IN_LONG; }
2595 39256 : z = cgetr(prec); av = avma;
2596 39252 : if (e > 0)
2597 : { /* simplify x/2^e = p / (q * 2^e) */
2598 2478 : long v = minss(e, vali(p));
2599 2478 : if (v) p = shifti(p, -v);
2600 2478 : if (e - v) q = shifti(q, e - v);
2601 : }
2602 39252 : s = signe(p);
2603 39252 : if (s < 0) p = negi(p);
2604 39255 : ex = exp2(dbllog2(x) - e) * 2.718281828; /* exp(1) * x / 2^e, x / 2^e < 2 */
2605 39261 : nmax = (long)(1 + exp(dbllambertW0(M_LN2 * b / ex)) * ex);
2606 39262 : abpq_init(&S, nmax);
2607 39283 : S.a[0] = S.b[0] = S.p[0] = S.q[0] = gen_1;
2608 3369690 : for (n = 1; n <= nmax; n++)
2609 : {
2610 3330442 : S.a[n] = gen_1;
2611 3330442 : S.b[n] = gen_1;
2612 3330442 : S.p[n] = p;
2613 3330442 : S.q[n] = muliu(q, n);
2614 : }
2615 39248 : abpq_sum(&R, 0, nmax, &S);
2616 39258 : if (s > 0) rdiviiz(R.T, R.Q, z); else rdiviiz(R.Q, R.T, z);
2617 39260 : if (e > 0)
2618 : {
2619 17136 : q = z; while (e--) q = sqrr(q);
2620 2478 : if (z0) { affrr(q, z0); z = z0; } else affrr(q,z);
2621 : }
2622 39260 : return gc_const(av,z);
2623 : }
2624 :
2625 : GEN
2626 17409353 : gexp(GEN x, long prec)
2627 : {
2628 17409353 : switch(typ(x))
2629 : {
2630 1468531 : case t_INT: case t_FRAC: return expQ(x, prec);
2631 10425919 : case t_REAL: return mpexp(x);
2632 4166418 : case t_COMPLEX: return cxexp(x,prec);
2633 70 : case t_PADIC: return Qp_exp(x);
2634 1348415 : default:
2635 : {
2636 1348415 : pari_sp av = avma;
2637 : GEN y;
2638 1348415 : if (!(y = toser_i(x))) break;
2639 66192 : return gerepileupto(av, serexp(y,prec));
2640 : }
2641 : }
2642 1282935 : return trans_eval("exp",gexp,x,prec);
2643 : }
2644 :
2645 : /********************************************************************/
2646 : /** **/
2647 : /** AGM(X, Y) **/
2648 : /** **/
2649 : /********************************************************************/
2650 : static int
2651 13787254 : agmr_gap(GEN a, GEN b, long L)
2652 : {
2653 13787254 : GEN d = subrr(b, a);
2654 13787228 : return (signe(d) && expo(d) - expo(b) >= L);
2655 : }
2656 : /* assume x > 0 */
2657 : static GEN
2658 940495 : agm1r_abs(GEN x)
2659 : {
2660 940495 : long l = realprec(x), L = 5-prec2nbits(l);
2661 940495 : GEN a1, b1, y = cgetr(l);
2662 940493 : pari_sp av = avma;
2663 :
2664 940493 : a1 = addrr(real_1(l), x); shiftr_inplace(a1, -1);
2665 940494 : b1 = sqrtr_abs(x);
2666 13787276 : while (agmr_gap(a1,b1,L))
2667 : {
2668 12846770 : GEN a = a1;
2669 12846770 : a1 = addrr(a,b1); shiftr_inplace(a1, -1);
2670 12846838 : b1 = sqrtr_abs(mulrr(a,b1));
2671 : }
2672 940457 : affrr_fixlg(a1,y); return gc_const(av,y);
2673 : }
2674 :
2675 : struct agmcx_gap_t { long L, ex, cnt; };
2676 :
2677 : static void
2678 324421 : agmcx_init(GEN x, long *prec, struct agmcx_gap_t *S)
2679 : {
2680 324421 : long l = precision(x);
2681 324421 : if (l) *prec = l;
2682 324421 : S->L = 1-prec2nbits(*prec);
2683 324421 : S->cnt = 0;
2684 324421 : S->ex = LONG_MAX;
2685 324421 : }
2686 :
2687 : static long
2688 324421 : agmcx_a_b(GEN x, GEN *a1, GEN *b1, long prec)
2689 : {
2690 324421 : long rotate = 0;
2691 324421 : if (gsigne(real_i(x))<0)
2692 : { /* Rotate by +/-Pi/2, so that the choice of the principal square
2693 : * root gives the optimal AGM. So a1 = +/-I*a1, b1=sqrt(-x). */
2694 11655 : if (gsigne(imag_i(x))<0) { *a1=mulcxI(*a1); rotate=-1; }
2695 11137 : else { *a1=mulcxmI(*a1); rotate=1; }
2696 11655 : x = gneg(x);
2697 : }
2698 324421 : *b1 = gsqrt(x, prec);
2699 324421 : return rotate;
2700 : }
2701 : /* return 0 if we must stop the AGM loop (a=b or a ~ b), 1 otherwise */
2702 : static int
2703 4854196 : agmcx_gap(GEN a, GEN b, struct agmcx_gap_t *S)
2704 : {
2705 4854196 : GEN d = gsub(b, a);
2706 4854196 : long ex = S->ex;
2707 4854196 : S->ex = gexpo(d);
2708 4854196 : if (gequal0(d) || S->ex - gexpo(b) < S->L) return 0;
2709 : /* if (S->ex >= ex) we're no longer making progress; twice in a row */
2710 4620980 : if (S->ex < ex) S->cnt = 0;
2711 : else
2712 182939 : if (S->cnt++) return 0;
2713 4529775 : return 1;
2714 : }
2715 : static GEN
2716 295672 : agm1cx(GEN x, long prec)
2717 : {
2718 : struct agmcx_gap_t S;
2719 : GEN a1, b1;
2720 295672 : pari_sp av = avma;
2721 : long rotate;
2722 295672 : agmcx_init(x, &prec, &S);
2723 295672 : a1 = gtofp(gmul2n(gadd(real_1(prec), x), -1), prec);
2724 295672 : rotate = agmcx_a_b(x, &a1, &b1, prec);
2725 4672181 : while (agmcx_gap(a1,b1,&S))
2726 : {
2727 4376509 : GEN a = a1;
2728 4376509 : a1 = gmul2n(gadd(a,b1),-1);
2729 4376509 : b1 = gsqrt(gmul(a,b1), prec);
2730 : }
2731 295672 : if (rotate) a1 = rotate>0 ? mulcxI(a1):mulcxmI(a1);
2732 295672 : return gerepilecopy(av,a1);
2733 : }
2734 :
2735 : GEN
2736 28749 : zellagmcx(GEN a0, GEN b0, GEN r, GEN t, long prec)
2737 : {
2738 : struct agmcx_gap_t S;
2739 28749 : pari_sp av = avma;
2740 28749 : GEN x = gdiv(a0, b0), a1, b1;
2741 : long rotate;
2742 28749 : agmcx_init(x, &prec, &S);
2743 28749 : a1 = gtofp(gmul2n(gadd(real_1(prec), x), -1), prec);
2744 28749 : r = gsqrt(gdiv(gmul(a1,gaddgs(r, 1)),gadd(r, x)), prec);
2745 28749 : t = gmul(r, t);
2746 28749 : rotate = agmcx_a_b(x, &a1, &b1, prec);
2747 182015 : while (agmcx_gap(a1,b1,&S))
2748 : {
2749 153266 : GEN a = a1, b = b1;
2750 153266 : a1 = gmul2n(gadd(a,b),-1);
2751 153266 : b1 = gsqrt(gmul(a,b), prec);
2752 153266 : r = gsqrt(gdiv(gmul(a1,gaddgs(r, 1)),gadd(gmul(b, r), a )), prec);
2753 153266 : t = gmul(r, t);
2754 : }
2755 28749 : if (rotate) a1 = rotate>0 ? mulcxI(a1):mulcxmI(a1);
2756 28749 : a1 = gmul(a1, b0);
2757 28749 : t = gatan(gdiv(a1,t), prec);
2758 : /* send t to the fundamental domain if necessary */
2759 28749 : if (gsigne(real_i(t))<0) t = gadd(t, mppi(prec));
2760 28749 : return gerepileupto(av,gdiv(t,a1));
2761 : }
2762 :
2763 : static long
2764 49 : ser_cmp_expo(GEN A, GEN B)
2765 : {
2766 49 : long e = -(long)HIGHEXPOBIT, d = valser(B) - valser(A);
2767 49 : long i, la = lg(A), v = varn(B);
2768 9849 : for (i = 2; i < la; i++)
2769 : {
2770 9800 : GEN a = gel(A,i), b;
2771 : long ei;
2772 9800 : if (isexactzero(a)) continue;
2773 9800 : b = polcoef_i(B, i-2 + d, v);
2774 9800 : ei = gexpo(a);
2775 9800 : if (!isexactzero(b)) ei -= gexpo(b);
2776 9800 : e = maxss(e, ei);
2777 : }
2778 49 : return e;
2779 : }
2780 :
2781 : static GEN
2782 21 : ser_agm1(GEN y, long prec)
2783 : {
2784 21 : GEN a1 = y, b1 = gen_1;
2785 21 : long l = lg(y)-2, l2 = 6-prec2nbits(prec), eold = LONG_MAX;
2786 : for(;;)
2787 84 : {
2788 105 : GEN a = a1, p1;
2789 105 : a1 = gmul2n(gadd(a,b1),-1);
2790 105 : b1 = gsqrt(gmul(a,b1), prec);
2791 105 : p1 = gsub(b1,a1);
2792 105 : if (isinexactreal(p1))
2793 : {
2794 49 : long e = ser_cmp_expo(p1, b1);
2795 49 : if (e < l2 || e >= eold) break;
2796 42 : eold = e;
2797 : }
2798 56 : else if (valser(p1)-valser(b1) >= l || gequal0(p1)) break;
2799 : }
2800 21 : return a1;
2801 : }
2802 :
2803 : /* agm(1,x) */
2804 : static GEN
2805 75159 : agm1(GEN x, long prec)
2806 : {
2807 : GEN y;
2808 : pari_sp av;
2809 :
2810 75159 : if (gequal0(x)) return gcopy(x);
2811 75159 : switch(typ(x))
2812 : {
2813 28 : case t_INT:
2814 28 : if (!is_pm1(x)) break;
2815 21 : return (signe(x) > 0)? real_1(prec): real_0(prec);
2816 :
2817 37835 : case t_REAL: return signe(x) > 0? agm1r_abs(x): agm1cx(x, prec);
2818 :
2819 37156 : case t_COMPLEX:
2820 37156 : if (gequal0(gel(x,2))) return agm1(gel(x,1), prec);
2821 37121 : return agm1cx(x, prec);
2822 :
2823 14 : case t_PADIC:
2824 : {
2825 14 : GEN a1 = x, b1 = gen_1;
2826 14 : long l = precp(x);
2827 14 : av = avma;
2828 : for(;;)
2829 14 : {
2830 28 : GEN a = a1, p1;
2831 : long ep;
2832 28 : a1 = gmul2n(gadd(a,b1),-1);
2833 28 : a = gmul(a,b1);
2834 28 : b1 = Qp_sqrt(a); if (!b1) pari_err_SQRTN("Qp_sqrt",a);
2835 21 : p1 = gsub(b1,a1); ep = valp(p1)-valp(b1);
2836 21 : if (ep<=0) { b1 = gneg_i(b1); p1 = gsub(b1,a1); ep=valp(p1)-valp(b1); }
2837 21 : if (ep >= l || gequal0(p1)) return gerepilecopy(av,a1);
2838 : }
2839 : }
2840 :
2841 126 : default:
2842 126 : av = avma; if (!(y = toser_i(x))) break;
2843 21 : return gerepilecopy(av, ser_agm1(y, prec));
2844 : }
2845 112 : return trans_eval("agm",agm1,x,prec);
2846 : }
2847 :
2848 : GEN
2849 74977 : agm(GEN x, GEN y, long prec)
2850 : {
2851 : pari_sp av;
2852 74977 : if (is_matvec_t(typ(y)))
2853 : {
2854 14 : if (is_matvec_t(typ(x))) pari_err_TYPE2("agm",x,y);
2855 7 : swap(x, y);
2856 : }
2857 74970 : if (gequal0(y)) return gcopy(y);
2858 74970 : av = avma;
2859 74970 : return gerepileupto(av, gmul(y, agm1(gdiv(x,y), prec)));
2860 : }
2861 :
2862 : /* b2 != 0 */
2863 : static GEN
2864 35 : ellK_i(GEN b2, long prec)
2865 35 : { return gdiv(Pi2n(-1, prec), agm1(gsqrt(b2, prec), prec)); }
2866 : GEN
2867 28 : ellK(GEN k, long prec)
2868 : {
2869 28 : pari_sp av = avma;
2870 28 : GEN k2 = gsqr(k), b2 = gsubsg(1, k2);
2871 28 : if (gequal0(b2)) pari_err_DOMAIN("ellK", "k^2", "=", gen_1, k2);
2872 21 : return gerepileupto(av, ellK_i(b2, prec));
2873 : }
2874 :
2875 : static int
2876 84 : magm_gap(GEN a, GEN b, long L)
2877 : {
2878 84 : GEN d = gsub(b, a);
2879 84 : return !gequal0(d) && gexpo(d) - gexpo(b) >= L;
2880 : }
2881 :
2882 : /* http://www.ams.org/notices/201208/rtx120801094p.pdf
2883 : * An Eloquent Formula for the Perimeter of an Ellipse
2884 : * Semjon Adlaj, Notices of the AMS */
2885 : static GEN
2886 14 : magm(GEN a, GEN b, long prec)
2887 : {
2888 14 : long L = -prec2nbits(prec) + 16;
2889 14 : GEN c = gen_0;
2890 84 : while (magm_gap(a, b, L))
2891 : {
2892 70 : GEN u = gsqrt(gmul(gsub(a, c), gsub(b, c)), prec);
2893 70 : a = gmul2n(gadd(a, b), -1);
2894 70 : b = gadd(c, u); c = gsub(c, u);
2895 : }
2896 14 : return gmul2n(gadd(a, b), -1);
2897 : }
2898 :
2899 : GEN
2900 21 : ellE(GEN k, long prec)
2901 : {
2902 21 : pari_sp av = avma;
2903 21 : GEN b2 = gsubsg(1, gsqr(k));
2904 21 : if (gequal0(b2)) { set_avma(av); return real_1(prec); }
2905 14 : return gerepileupto(av, gmul(ellK_i(b2, prec), magm(gen_1, b2, prec)));
2906 : }
2907 :
2908 : /********************************************************************/
2909 : /** **/
2910 : /** LOG(X) **/
2911 : /** **/
2912 : /********************************************************************/
2913 : /* log(2) = 18*atanh(1/26)-2*atanh(1/4801)+8*atanh(1/8749)
2914 : * faster than 10*atanh(1/17)+4*atanh(13/499) for all precisions,
2915 : * and than Pi/2M(1,4/2^n) ~ n log(2) for bitprec at least up to 10^8 */
2916 : static GEN
2917 31911 : log2_split(long prec)
2918 : {
2919 31911 : GEN u = atanhuu(1, 26, prec);
2920 31915 : GEN v = atanhuu(1, 4801, prec);
2921 31910 : GEN w = atanhuu(1, 8749, prec);
2922 31908 : shiftr_inplace(v, 1); setsigne(v, -1);
2923 31910 : shiftr_inplace(w, 3);
2924 31908 : return addrr(mulur(18, u), addrr(v, w));
2925 : }
2926 : GEN
2927 28476474 : constlog2(long prec)
2928 : {
2929 : pari_sp av;
2930 : GEN tmp;
2931 28476474 : if (glog2 && realprec(glog2) >= prec) return glog2;
2932 :
2933 31814 : tmp = cgetr_block(prec);
2934 31911 : av = avma;
2935 31911 : affrr(log2_split(prec+EXTRAPREC64), tmp);
2936 31917 : swap_clone(&glog2,tmp);
2937 31922 : return gc_const(av,glog2);
2938 : }
2939 :
2940 : GEN
2941 28476463 : mplog2(long prec) { return rtor(constlog2(prec), prec); }
2942 :
2943 : /* dont check that q != 2^expo(q), done in logr_abs */
2944 : static GEN
2945 902700 : logagmr_abs(GEN q)
2946 : {
2947 902700 : long prec = realprec(q), e = expo(q), lim;
2948 902700 : GEN z = cgetr(prec), y, Q, _4ovQ;
2949 902698 : pari_sp av = avma;
2950 :
2951 902698 : incrprec(prec);
2952 902698 : lim = prec2nbits(prec) >> 1;
2953 902698 : Q = rtor(q,prec);
2954 902701 : shiftr_inplace(Q,lim-e); setsigne(Q,1);
2955 :
2956 902701 : _4ovQ = invr(Q); shiftr_inplace(_4ovQ, 2); /* 4/Q */
2957 : /* Pi / 2agm(1, 4/Q) ~ log(Q), q = Q * 2^(e-lim) */
2958 902702 : y = divrr(Pi2n(-1, prec), agm1r_abs(_4ovQ));
2959 902701 : y = addrr(y, mulsr(e - lim, mplog2(prec)));
2960 902702 : affrr_fixlg(y, z); return gc_const(av,z);
2961 : }
2962 :
2963 : /* sum_{k >= 0} y^(2k+1) / (2k+1), y close to 0 */
2964 : static GEN
2965 11774289 : logr_aux(GEN y)
2966 : {
2967 11774289 : long k, L = realprec(y); /* should be ~ l+1 - (k-2) */
2968 : /* log(x) = log(1+y) - log(1-y) = 2 sum_{k odd} y^k / k
2969 : * Truncate the sum at k = 2n+1, the remainder is
2970 : * 2 sum_{k >= 2n+3} y^k / k < 2y^(2n+3) / (2n+3)(1-y) < y^(2n+3)
2971 : * We want y^(2n+3) < y 2^(-prec2nbits(L)), hence
2972 : * n+1 > -prec2nbits(L) /-log_2(y^2) */
2973 11774289 : double d = -2*dbllog2r(y); /* ~ -log_2(y^2) */
2974 11774239 : k = (long)(2*(prec2nbits(L) / d));
2975 11774181 : k |= 1;
2976 11774181 : if (k >= 3)
2977 : {
2978 11751768 : GEN T, S = cgetr(L), y2 = sqrr(y), unr = real_1(L);
2979 11752034 : pari_sp av = avma;
2980 11752034 : long s = 0, incs = (long)d, l1 = nbits2prec((long)d);
2981 11752039 : setprec(S, l1);
2982 11752009 : setprec(unr,l1); affrr(divru(unr,k), S);
2983 208089960 : for (k -= 2;; k -= 2) /* k = 2n+1, ..., 1 */
2984 : { /* S = y^(2n+1-k)/(2n+1) + ... + 1 / k */
2985 208089960 : setprec(y2, l1); T = mulrr(S,y2);
2986 208215168 : if (k == 1) break;
2987 :
2988 196463054 : l1 += nbits2extraprec(dvmdsBIL(s + incs, &s)<<TWOPOTBITS_IN_LONG);
2989 196442027 : if (l1>L) l1=L;
2990 196442027 : setprec(S, l1);
2991 196432512 : setprec(unr,l1);
2992 196411657 : affrr(addrr(divru(unr, k), T), S); set_avma(av);
2993 : }
2994 : /* k = 1 special-cased for eficiency */
2995 11752114 : y = mulrr(y, addsr(1,T)); /* = log(X)/2 */
2996 : }
2997 11774529 : return y;
2998 : }
2999 : /*return log(|x|), assuming x != 0 */
3000 : GEN
3001 13472518 : logr_abs(GEN X)
3002 : {
3003 13472518 : long EX, L, m, k, a, b, l = lg(X), p = realprec(X);
3004 : GEN z, x, y;
3005 : ulong u;
3006 : double d;
3007 :
3008 : /* Assuming 1 < x < 2, we want delta = x-1, 1-x/2, 1-1/x, or 2/x-1 small.
3009 : * We have 2/x-1 > 1-x/2, 1-1/x < x-1. So one should be choosing between
3010 : * 1-1/x and 1-x/2 ( crossover sqrt(2), worse ~ 0.29 ). To avoid an inverse,
3011 : * we choose between x-1 and 1-x/2 ( crossover 4/3, worse ~ 0.33 ) */
3012 13472518 : EX = expo(X);
3013 13472518 : u = uel(X,2);
3014 13472518 : k = 2;
3015 13472518 : if (u > (~0UL / 3) * 2) { /* choose 1-x/2 */
3016 7641038 : EX++; u = ~u;
3017 7711434 : while (!u && ++k < l) { u = uel(X,k); u = ~u; }
3018 : } else { /* choose x - 1 */
3019 5831480 : u &= ~HIGHBIT; /* u - HIGHBIT, assuming HIGHBIT set */
3020 7142840 : while (!u && ++k < l) u = uel(X,k);
3021 : }
3022 13472518 : if (k == l) return EX? mulsr(EX, mplog2(p)): real_0(p);
3023 12677079 : a = bit_accuracy(k) + bfffo(u); /* ~ -log2 |1-x| */
3024 12677139 : L = p+EXTRAPRECWORD;
3025 12677139 : b = prec2nbits(L - (bit_accuracy(k))); /* take loss of accuracy into account */
3026 12677143 : if (b > 24*a*log2(prec2lg(L)) && p > LOGAGM_LIMIT) return logagmr_abs(X);
3027 :
3028 11774509 : z = cgetr(EX? p: p - bit_accuracy(k));
3029 :
3030 : /* Multiplication is quadratic in this range (l is small, otherwise we
3031 : * use AGM). Set Y = x^(1/2^m), y = (Y - 1) / (Y + 1) and compute truncated
3032 : * series sum y^(2k+1)/(2k+1): the costs is less than
3033 : * m b^2 + sum_{k <= n} ((2k+1) e + BITS_IN_LONG)^2
3034 : * bit operations with |x-1| < 2^(1-a), |Y| < 2^(1-e) , m = e-a and b bits of
3035 : * accuracy needed (+ BITS_IN_LONG since bit accuracies increase by
3036 : * increments of BITS_IN_LONG), so
3037 : * 4n^3/3 e^2 + n^2 2e BITS_IN_LONG+ n BITS_IN_LONG ~ m b^2, with n ~ b/2e
3038 : * or b/6e + BITS_IN_LONG/2e + BITS_IN_LONG/2be ~ m
3039 : * B := (b / 6 + BITS_IN_LONG/2 + BITS_IN_LONG^2 / 2b) ~ m(m+a)
3040 : * m = min( -a/2 + sqrt(a^2/4 + B), b - a )
3041 : * NB: e ~ (b/6)^(1/2) as b -> oo
3042 : * Instead of the above pessimistic estimate for the cost of the sum, use
3043 : * optimistic estimate (BITS_IN_LONG -> 0) */
3044 11774495 : d = -a/2.; m = (long)(d + sqrt(d*d + b/6)); /* >= 0 */
3045 :
3046 11774495 : if (m > b-a) m = b-a;
3047 11774495 : if (m < 0.2*a) m = 0; else L += nbits2extraprec(m);
3048 11774481 : x = rtor(X,L);
3049 11774452 : setsigne(x,1); shiftr_inplace(x,-EX);
3050 : /* 2/3 < x < 4/3 */
3051 68136943 : for (k=1; k<=m; k++) x = sqrtr_abs(x);
3052 :
3053 11774373 : y = divrr(subrs(x,1), addrs(x,1)); /* = (x-1) / (x+1), close to 0 */
3054 11774260 : y = logr_aux(y); /* log(1+y) - log(1-y) = log(x) */
3055 11774461 : shiftr_inplace(y, m + 1);
3056 11774358 : if (EX) y = addrr(y, mulsr(EX, mplog2(p+EXTRAPRECWORD)));
3057 11774107 : affrr_fixlg(y, z); return gc_const((pari_sp)z, z);
3058 : }
3059 :
3060 : /* assume Im(q) != 0 and precision(q) >= prec. Compute log(q) with accuracy
3061 : * prec [disregard input accuracy] */
3062 : GEN
3063 258509 : logagmcx(GEN q, long prec)
3064 : {
3065 258509 : GEN z = cgetc(prec), y, Q, a, b;
3066 : long lim, e, ea, eb;
3067 258509 : pari_sp av = avma;
3068 258509 : int neg = 0;
3069 :
3070 258509 : incrprec(prec);
3071 258509 : if (gsigne(gel(q,1)) < 0) { q = gneg(q); neg = 1; }
3072 258509 : lim = prec2nbits(prec) >> 1;
3073 258509 : Q = gtofp(q, prec);
3074 258509 : a = gel(Q,1);
3075 258509 : b = gel(Q,2);
3076 258509 : if (gequal0(a)) {
3077 0 : affrr_fixlg(logr_abs(b), gel(z,1));
3078 0 : y = Pi2n(-1, prec);
3079 0 : if (signe(b) < 0) setsigne(y, -1);
3080 0 : affrr_fixlg(y, gel(z,2)); return gc_const(av,z);
3081 : }
3082 258509 : ea = expo(a);
3083 258509 : eb = expo(b);
3084 258509 : e = ea <= eb ? lim - eb : lim - ea;
3085 258509 : shiftr_inplace(a, e);
3086 258509 : shiftr_inplace(b, e);
3087 :
3088 : /* Pi / 2agm(1, 4/Q) ~ log(Q), q = Q * 2^e */
3089 258509 : y = gdiv(Pi2n(-1, prec), agm1cx( gdivsg(4, Q), prec ));
3090 258509 : a = gel(y,1);
3091 258509 : b = gel(y,2);
3092 258509 : a = addrr(a, mulsr(-e, mplog2(prec)));
3093 258509 : if (realprec(a) <= LOWDEFAULTPREC) a = real_0_bit(expo(a));
3094 389517 : if (neg) b = gsigne(b) <= 0? gadd(b, mppi(prec))
3095 131008 : : gsub(b, mppi(prec));
3096 258509 : affrr_fixlg(a, gel(z,1));
3097 258509 : affrr_fixlg(b, gel(z,2)); return gc_const(av,z);
3098 : }
3099 :
3100 : GEN
3101 203444 : mplog(GEN x)
3102 : {
3103 203444 : if (signe(x)<=0) pari_err_DOMAIN("mplog", "argument", "<=", gen_0, x);
3104 203444 : return logr_abs(x);
3105 : }
3106 :
3107 : /* pe = p^e, p prime, 0 < x < pe a t_INT coprime to p. Return the (p-1)-th
3108 : * root of 1 in (Z/pe)^* congruent to x mod p, resp x mod 4 if p = 2.
3109 : * Simplified form of Zp_sqrtnlift: 1/(p-1) is trivial to compute */
3110 : GEN
3111 9933 : Zp_teichmuller(GEN x, GEN p, long e, GEN pe)
3112 : {
3113 : GEN q, z, p1;
3114 : pari_sp av;
3115 : ulong mask;
3116 9933 : if (absequaliu(p,2)) return (mod4(x) & 2)? subiu(pe,1): gen_1;
3117 9394 : if (e == 1) return icopy(x);
3118 9394 : av = avma;
3119 9394 : p1 = subiu(p, 1);
3120 9394 : mask = quadratic_prec_mask(e);
3121 9394 : q = p; z = remii(x, p);
3122 31192 : while (mask > 1)
3123 : { /* Newton iteration solving z^{1 - p} = 1, z = x (mod p) */
3124 21798 : GEN w, t, qold = q;
3125 21798 : if (mask <= 3) /* last iteration */
3126 9394 : q = pe;
3127 : else
3128 : {
3129 12404 : q = sqri(q);
3130 12404 : if (mask & 1) q = diviiexact(q, p);
3131 : }
3132 21798 : mask >>= 1;
3133 : /* q <= qold^2 */
3134 21798 : if (lgefint(q) == 3)
3135 : {
3136 21650 : ulong Z = uel(z,2), Q = uel(q,2), P1 = uel(p1,2);
3137 21650 : ulong W = (Q-1) / P1; /* -1/(p-1) + O(qold) */
3138 21650 : ulong T = Fl_mul(W, Fl_powu(Z,P1,Q) - 1, Q);
3139 21650 : Z = Fl_mul(Z, 1 + T, Q);
3140 21650 : z = utoi(Z);
3141 : }
3142 : else
3143 : {
3144 148 : w = diviiexact(subiu(qold,1),p1); /* -1/(p-1) + O(qold) */
3145 148 : t = Fp_mul(w, subiu(Fp_pow(z,p1,q), 1), q);
3146 148 : z = Fp_mul(z, addui(1,t), q);
3147 : }
3148 : }
3149 9394 : return gerepileuptoint(av, z);
3150 : }
3151 :
3152 : GEN
3153 1225 : teichmullerinit(long p, long n)
3154 : {
3155 : GEN t, pn, g, v;
3156 : ulong gp, tp;
3157 : long a, m;
3158 :
3159 1225 : if (p == 2) return mkvec(gen_1);
3160 1225 : if (!uisprime(p)) pari_err_PRIME("teichmullerinit",utoipos(p));
3161 :
3162 1225 : m = p >> 1; /* (p-1)/2 */
3163 1225 : tp= gp= pgener_Fl(p); /* order (p-1), gp^m = -1 */
3164 1225 : pn = powuu(p, n);
3165 1225 : v = cgetg(p, t_VEC);
3166 1225 : t = g = Zp_teichmuller(utoipos(gp), utoipos(p), n, pn);
3167 1225 : gel(v, 1) = gen_1;
3168 1225 : gel(v, p-1) = subiu(pn,1);
3169 3031 : for (a = 1; a < m; a++)
3170 : {
3171 1806 : gel(v, tp) = t;
3172 1806 : gel(v, p - tp) = Fp_neg(t, pn); /* g^(m+a) = -g^a */
3173 1806 : if (a < m-1)
3174 : {
3175 1029 : t = Fp_mul(t, g, pn); /* g^(a+1) */
3176 1029 : tp = Fl_mul(tp, gp, p); /* t mod p */
3177 : }
3178 : }
3179 1225 : return v;
3180 : }
3181 :
3182 : /* tab from teichmullerinit or NULL */
3183 : GEN
3184 4977 : teichmuller(GEN x, GEN tab)
3185 : {
3186 : GEN p, q, y, z;
3187 4977 : long n, tx = typ(x);
3188 :
3189 4977 : if (!tab)
3190 : {
3191 4865 : if (tx == t_VEC && lg(x) == 3)
3192 : {
3193 7 : p = gel(x,1);
3194 7 : q = gel(x,2);
3195 7 : if (typ(p) == t_INT && typ(q) == t_INT)
3196 7 : return teichmullerinit(itos(p), itos(q));
3197 : }
3198 : }
3199 112 : else if (typ(tab) != t_VEC) pari_err_TYPE("teichmuller",tab);
3200 4970 : if (tx!=t_PADIC) pari_err_TYPE("teichmuller",x);
3201 4970 : z = gel(x,4);
3202 4970 : if (!signe(z)) return gcopy(x);
3203 4970 : p = gel(x,2);
3204 4970 : q = gel(x,3);
3205 4970 : n = precp(x);
3206 4970 : y = cgetg(5,t_PADIC);
3207 4970 : y[1] = evalprecp(n) | _evalvalp(0);
3208 4970 : gel(y,2) = icopy(p);
3209 4970 : gel(y,3) = icopy(q);
3210 4970 : if (tab)
3211 : {
3212 112 : ulong pp = itou_or_0(p);
3213 112 : if (lg(tab) != (long)pp) pari_err_TYPE("teichmuller",tab);
3214 112 : z = gel(tab, umodiu(z, pp));
3215 112 : if (typ(z) != t_INT) pari_err_TYPE("teichmuller",tab);
3216 112 : z = remii(z, q);
3217 : }
3218 : else
3219 4858 : z = Zp_teichmuller(z, p, n, q);
3220 4970 : gel(y,4) = z;
3221 4970 : return y;
3222 : }
3223 : GEN
3224 4739 : teich(GEN x) { return teichmuller(x, NULL); }
3225 :
3226 : GEN
3227 17321704 : glog(GEN x, long prec)
3228 : {
3229 : pari_sp av, tetpil;
3230 : GEN y, p1;
3231 : long l;
3232 :
3233 17321704 : switch(typ(x))
3234 : {
3235 10069479 : case t_REAL:
3236 10069479 : if (signe(x) >= 0)
3237 : {
3238 8384030 : if (!signe(x)) pari_err_DOMAIN("log", "argument", "=", gen_0, x);
3239 8384016 : return logr_abs(x);
3240 : }
3241 1685449 : retmkcomplex(logr_abs(x), mppi(realprec(x)));
3242 :
3243 514428 : case t_FRAC:
3244 : {
3245 : GEN a, b;
3246 : long e1, e2;
3247 514428 : av = avma;
3248 514428 : a = gel(x,1);
3249 514428 : b = gel(x,2);
3250 514428 : e1 = expi(subii(a,b)); e2 = expi(b);
3251 514428 : if (e2 > e1) prec += nbits2extraprec(e2 - e1);
3252 514428 : x = fractor(x, prec);
3253 514428 : return gerepileupto(av, glog(x, prec));
3254 : }
3255 4402669 : case t_COMPLEX:
3256 4402669 : if (ismpzero(gel(x,2))) return glog(gel(x,1), prec);
3257 4391723 : l = precision(x); if (l > prec) prec = l;
3258 4391740 : if (ismpzero(gel(x,1)))
3259 : {
3260 72875 : GEN a = gel(x,2), b;
3261 72875 : av = avma; b = Pi2n(-1,prec);
3262 72876 : if (gsigne(a) < 0) { setsigne(b, -1); a = gabs(a,prec); }
3263 72876 : a = isint1(a) ? gen_0: glog(a,prec);
3264 72876 : return gerepilecopy(av, mkcomplex(a, b));
3265 : }
3266 4318871 : if (prec >= LOGAGMCX_LIMIT) return logagmcx(x, prec);
3267 4060546 : y = cgetg(3,t_COMPLEX);
3268 4060547 : gel(y,2) = garg(x,prec);
3269 4060519 : av = avma; p1 = glog(cxnorm(x),prec); tetpil = avma;
3270 4060536 : gel(y,1) = gerepile(av,tetpil,gmul2n(p1,-1)); return y;
3271 :
3272 322 : case t_PADIC: return Qp_log(x);
3273 2334806 : default:
3274 2334806 : av = avma; if (!(y = toser_i(x))) break;
3275 140 : if (!signe(y)) pari_err_DOMAIN("log", "argument", "=", gen_0, x);
3276 140 : if (valser(y)) pari_err_DOMAIN("log", "series valuation", "!=", gen_0, x);
3277 133 : p1 = integser(gdiv(derivser(y), y)); /* log(y)' = y'/y */
3278 133 : if (!gequal1(gel(y,2))) p1 = gadd(p1, glog(gel(y,2),prec));
3279 133 : return gerepileupto(av, p1);
3280 : }
3281 2334966 : return trans_eval("log",glog,x,prec);
3282 : }
3283 :
3284 : static GEN
3285 63 : mplog1p(GEN x)
3286 : {
3287 : long ex, a, b, l, L;
3288 63 : if (!signe(x)) return rcopy(x);
3289 63 : ex = expo(x); if (ex >= -3) return glog(addrs(x,1), 0);
3290 42 : a = -ex;
3291 42 : b = realprec(x); L = b+1;
3292 42 : if (b > a*log2(L) && b > LOGAGM_LIMIT)
3293 : {
3294 0 : x = addrs(x,1); l = b + nbits2extraprec(a);
3295 0 : if (realprec(x) < l) x = rtor(x,l);
3296 0 : return logagmr_abs(x);
3297 : }
3298 42 : x = rtor(x, L);
3299 42 : x = logr_aux(divrr(x, addrs(x,2)));
3300 42 : if (realprec(x) > b) fixlg(x, b);
3301 42 : shiftr_inplace(x,1); return x;
3302 : }
3303 :
3304 : static GEN log1p_i(GEN x, long prec);
3305 : static GEN
3306 14 : cxlog1p(GEN x, long prec)
3307 : {
3308 : pari_sp av;
3309 14 : GEN z, a, b = gel(x,2);
3310 : long l;
3311 14 : if (ismpzero(b)) return log1p_i(gel(x,1), prec);
3312 14 : l = precision(x); if (l > prec) prec = l;
3313 14 : if (prec >= LOGAGMCX_LIMIT) return logagmcx(gaddgs(x,1), prec);
3314 14 : a = gel(x,1);
3315 14 : z = cgetg(3,t_COMPLEX); av = avma;
3316 14 : a = gadd(gadd(gmul2n(a,1), gsqr(a)), gsqr(b));
3317 14 : a = log1p_i(a, prec); shiftr_inplace(a,-1);
3318 14 : gel(z,1) = gerepileupto(av, a);
3319 14 : gel(z,2) = garg(gaddgs(x,1),prec); return z;
3320 : }
3321 : static GEN
3322 133 : log1p_i(GEN x, long prec)
3323 : {
3324 133 : switch(typ(x))
3325 : {
3326 63 : case t_REAL: return mplog1p(x);
3327 14 : case t_COMPLEX: return cxlog1p(x, prec);
3328 7 : case t_PADIC: return Qp_log(gaddgs(x,1));
3329 49 : default:
3330 : {
3331 : long ey;
3332 : GEN y;
3333 49 : if (!(y = toser_i(x))) break;
3334 21 : ey = valser(y);
3335 21 : if (ey < 0) pari_err_DOMAIN("log1p","valuation", "<", gen_0, x);
3336 21 : if (gequal0(y)) return gcopy(y);
3337 14 : if (ey)
3338 7 : return glog(gaddgs(y,1),prec);
3339 : else
3340 : {
3341 7 : GEN a = gel(y,2), a1 = gaddgs(a,1);
3342 7 : y = gdiv(y, a1); gel(y,2) = gen_1;
3343 7 : return gadd(glog1p(a,prec), glog(y, prec));
3344 : }
3345 : }
3346 : }
3347 28 : return trans_eval("log1p",glog1p,x,prec);
3348 : }
3349 : GEN
3350 119 : glog1p(GEN x, long prec)
3351 : {
3352 119 : pari_sp av = avma;
3353 119 : return gerepileupto(av, log1p_i(x, prec));
3354 : }
3355 : /********************************************************************/
3356 : /** **/
3357 : /** SINE, COSINE **/
3358 : /** **/
3359 : /********************************************************************/
3360 :
3361 : /* Reduce x0 mod Pi/2 to x in [-Pi/4, Pi/4]. Return cos(x)-1 */
3362 : static GEN
3363 16883473 : mpcosm1(GEN x, long *ptmod8)
3364 : {
3365 16883473 : long a = expo(x), l = realprec(x), b, L, i, n, m, B;
3366 : GEN y, u, x2;
3367 : double d;
3368 :
3369 16883473 : n = 0;
3370 16883473 : if (a >= 0)
3371 : {
3372 : long p;
3373 : GEN q;
3374 10711775 : if (a > 30)
3375 : {
3376 406468 : GEN z, P = Pi2n(-2, nbits2prec(a + 32));
3377 406468 : z = addrr(x,P); /* = x + Pi/4 */
3378 406468 : if (expo(z) >= bit_prec(z) + 3) pari_err_PREC("mpcosm1");
3379 406468 : shiftr_inplace(P, 1);
3380 406468 : q = floorr(divrr(z, P)); /* round ( x / (Pi/2) ) */
3381 406468 : p = l+EXTRAPREC64; x = rtor(x,p);
3382 : } else {
3383 10305307 : q = stoi((long)floor(rtodbl(x) / (M_PI/2) + 0.5));
3384 10305296 : p = l;
3385 : }
3386 10716158 : if (signe(q))
3387 : {
3388 10711764 : GEN y = subrr(x, mulir(q, Pi2n(-1,p))); /* x mod Pi/2 */
3389 10711656 : long b = expo(y);
3390 10711656 : if (a - b < 7) x = y;
3391 : else
3392 : {
3393 6445413 : p += nbits2extraprec(a-b); x = rtor(x, p);
3394 6445423 : x = subrr(x, mulir(q, Pi2n(-1,p)));
3395 : }
3396 10711628 : a = b;
3397 10711628 : if (!signe(x) && a >= 0) pari_err_PREC("mpcosm1");
3398 10711628 : n = Mod4(q);
3399 : }
3400 : }
3401 : /* a < 0 */
3402 16887702 : b = signe(x); *ptmod8 = (b < 0)? 4 + n: n;
3403 16887702 : if (!b) return real_0_bit(expo(x)*2 - 1);
3404 :
3405 16887702 : b = prec2nbits(l);
3406 16882032 : if (b + 2*a <= 0) {
3407 1382463 : y = sqrr(x); shiftr_inplace(y, -1); setsigne(y, -1);
3408 1382465 : return y;
3409 : }
3410 :
3411 15499569 : y = cgetr(l);
3412 15501318 : B = b/6 + BITS_IN_LONG/2 + (BITS_IN_LONG*BITS_IN_LONG/2)/ b;
3413 15501318 : d = a/2.; m = (long)(d + sqrt(d*d + B)); /* >= 0 ,*/
3414 15501318 : if (m < (-a) * 0.1) m = 0; /* not worth it */
3415 15501318 : L = l + nbits2extraprec(m);
3416 :
3417 15500911 : b += m;
3418 15500911 : d = 2.0 * (m-dbllog2r(x)-1/M_LN2); /* ~ 2( - log_2 Y - 1/log(2) ) */
3419 15500925 : n = (long)(b / d);
3420 15500925 : if (n > 1)
3421 15443011 : n = (long)(b / (d + log2((double)n+1))); /* log~constant in small ranges */
3422 33340811 : while (n*(d+log2((double)n+1)) < b) n++; /* expect few corrections */
3423 :
3424 : /* Multiplication is quadratic in this range (l is small, otherwise we
3425 : * use logAGM + Newton). Set Y = 2^(-e-a) x, compute truncated series
3426 : * sum Y^2k/(2k)!: this costs roughly
3427 : * m b^2 + sum_{k <= n} (2k e + BITS_IN_LONG)^2
3428 : * ~ (b/2e) b^2 / 3 + m b^2
3429 : * bit operations with n ~ b/2e, |x| < 2^(1+a), |Y| < 2^(1-e) , m = e+a and
3430 : * b bits of accuracy needed, so
3431 : * B := (b / 6 + BITS_IN_LONG/2 + BITS_IN_LONG^2 / 2b) ~ m(m-a)
3432 : * we want b ~ 6 m (m-a) or m~b+a hence
3433 : * m = min( a/2 + sqrt(a^2/4 + b/6), b/2 + a )
3434 : * NB: e ~ (b/6)^(1/2) or b/2.
3435 : *
3436 : * Truncate the sum at k = n (>= 1), the remainder is
3437 : * < sum_{k >= n+1} Y^2k / 2k! < Y^(2n+2) / (2n+2)!(1-Y^2) < Y^(2n+2)/(2n+1)!
3438 : * We want ... <= Y^2 2^-b, hence -2n log_2 |Y| + log_2 (2n+1)! >= b
3439 : * log n! ~ (n + 1/2) log(n+1) - (n+1) + log(2Pi)/2,
3440 : * error bounded by 1/6(n+1) <= 1/12. Finally, we want
3441 : * 2n (-1/log(2) - log_2 |Y| + log_2(2n+2)) >= b */
3442 15500925 : x = rtor(x, L); shiftr_inplace(x, -m); setsigne(x, 1);
3443 15501962 : x2 = sqrr(x);
3444 15505582 : if (n == 1) { u = x2; shiftr_inplace(u, -1); setsigne(u, -1); } /*-Y^2/2*/
3445 : else
3446 : {
3447 15505582 : GEN un = real_1(L);
3448 : pari_sp av;
3449 15503026 : long s = 0, l1 = nbits2prec((long)(d + n + 16));
3450 :
3451 15502763 : u = cgetr(L); av = avma;
3452 221610416 : for (i = n; i >= 2; i--)
3453 : {
3454 : GEN t;
3455 206112172 : setprec(x2,l1); t = divrunextu(x2, 2*i-1);
3456 206923879 : l1 += nbits2extraprec(dvmdsBIL(s - expo(t), &s)<<TWOPOTBITS_IN_LONG);
3457 206720747 : if (l1 > L) l1 = L;
3458 206720747 : if (i != n) t = mulrr(t,u);
3459 207086752 : setprec(un,l1); t = addrr_sign(un,1, t,-signe(t));
3460 206442150 : setprec(u,l1); affrr(t,u); set_avma(av);
3461 : }
3462 15498244 : shiftr_inplace(u, -1); togglesign(u); /* u := -u/2 */
3463 15500495 : setprec(x2,L); u = mulrr(x2,u);
3464 : }
3465 : /* Now u = sum {1<= i <=n} (-1)^i x^(2i) / (2i)! ~ cos(x) - 1 */
3466 135625648 : for (i = 1; i <= m; i++)
3467 : { /* u = cos(x)-1 <- cos(2x)-1 = 2cos(x)^2 - 2 = 4u + 2u^2*/
3468 120160681 : GEN q = sqrr(u);
3469 120385686 : shiftr_inplace(u, 1); u = addrr(u, q);
3470 120230496 : shiftr_inplace(u, 1);
3471 120126731 : if ((i & 31) == 0) u = gerepileuptoleaf((pari_sp)y, u);
3472 : }
3473 15464967 : affrr_fixlg(u, y); return y;
3474 : }
3475 :
3476 : /* sqrt (|1 - (1+x)^2|) = sqrt(|x*(x+2)|). Sends cos(x)-1 to |sin(x)| */
3477 : static GEN
3478 15186361 : mpaut(GEN x)
3479 : {
3480 15186361 : GEN t = mulrr(x, addsr(2,x)); /* != 0 */
3481 15192967 : if (!signe(t)) return real_0_bit(expo(t) >> 1);
3482 15192967 : return sqrtr_abs(t);
3483 : }
3484 :
3485 : /********************************************************************/
3486 : /** COSINE **/
3487 : /********************************************************************/
3488 :
3489 : GEN
3490 2745063 : mpcos(GEN x)
3491 : {
3492 : long mod8;
3493 : pari_sp av;
3494 : GEN y, z;
3495 :
3496 2745063 : if (!signe(x)) {
3497 75 : long l = nbits2prec(-expo(x));
3498 75 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
3499 75 : return real_1(l);
3500 : }
3501 2744988 : av = avma; z = mpcosm1(x,&mod8);
3502 2744973 : switch(mod8)
3503 : {
3504 760057 : case 0: case 4: y = addsr(1,z); break;
3505 688604 : case 1: case 7: y = mpaut(z); togglesign(y); break;
3506 682286 : case 2: case 6: y = subsr(-1,z); break;
3507 614026 : default: y = mpaut(z); break; /* case 3: case 5: */
3508 : }
3509 2745013 : return gerepileuptoleaf(av, y);
3510 : }
3511 :
3512 : /* convert INT or FRAC to REAL, which is later reduced mod 2Pi : avoid
3513 : * cancellation */
3514 : static GEN
3515 13253 : tofp_safe(GEN x, long prec)
3516 : {
3517 13253 : return (typ(x) == t_INT || gexpo(x) > 0)? gadd(x, real_0(prec))
3518 26493 : : fractor(x, prec);
3519 : }
3520 :
3521 : GEN
3522 154855 : gcos(GEN x, long prec)
3523 : {
3524 : pari_sp av;
3525 : GEN a, b, u, v, y, u1, v1;
3526 : long i;
3527 :
3528 154855 : switch(typ(x))
3529 : {
3530 153573 : case t_REAL: return mpcos(x);
3531 42 : case t_COMPLEX:
3532 42 : a = gel(x,1);
3533 42 : b = gel(x,2);
3534 42 : if (isintzero(a)) return gcosh(b, prec);
3535 28 : i = precision(x); if (i) prec = i;
3536 28 : y = cgetc(prec); av = avma;
3537 28 : if (typ(b) != t_REAL) b = gtofp(b, prec);
3538 28 : mpsinhcosh(b, &u1, &v1); u1 = mpneg(u1);
3539 28 : if (typ(a) != t_REAL) a = gtofp(a, prec);
3540 28 : mpsincos(a, &u, &v);
3541 28 : affrr_fixlg(gmul(v1,v), gel(y,1));
3542 28 : affrr_fixlg(gmul(u1,u), gel(y,2)); return gc_const(av,y);
3543 :
3544 1156 : case t_INT: case t_FRAC:
3545 1156 : y = cgetr(prec); av = avma;
3546 1156 : affrr_fixlg(mpcos(tofp_safe(x,prec)), y); return gc_const(av,y);
3547 :
3548 49 : case t_PADIC: y = cos_p(x);
3549 49 : if (!y) pari_err_DOMAIN("gcos(t_PADIC)","argument","",gen_0,x);
3550 42 : return y;
3551 :
3552 35 : default:
3553 35 : av = avma; if (!(y = toser_i(x))) break;
3554 28 : if (gequal0(y)) return gerepileupto(av, gaddsg(1,y));
3555 28 : if (valser(y) < 0)
3556 7 : pari_err_DOMAIN("cos","valuation", "<", gen_0, x);
3557 21 : gsincos(y,&u,&v,prec);
3558 21 : return gerepilecopy(av,v);
3559 : }
3560 7 : return trans_eval("cos",gcos,x,prec);
3561 : }
3562 : /********************************************************************/
3563 : /** SINE **/
3564 : /********************************************************************/
3565 :
3566 : GEN
3567 826141 : mpsin(GEN x)
3568 : {
3569 : long mod8;
3570 : pari_sp av;
3571 : GEN y, z;
3572 :
3573 826141 : if (!signe(x)) return real_0_bit(expo(x));
3574 825932 : av = avma; z = mpcosm1(x,&mod8);
3575 825923 : switch(mod8)
3576 : {
3577 309976 : case 0: case 6: y = mpaut(z); break;
3578 129039 : case 1: case 5: y = addsr(1,z); break;
3579 262405 : case 2: case 4: y = mpaut(z); togglesign(y); break;
3580 124503 : default: y = subsr(-1,z); break; /* case 3: case 7: */
3581 : }
3582 825951 : return gerepileuptoleaf(av, y);
3583 : }
3584 :
3585 : GEN
3586 2290228 : gsin(GEN x, long prec)
3587 : {
3588 : pari_sp av;
3589 : GEN a, b, u, v, y, v1, u1;
3590 : long i;
3591 :
3592 2290228 : switch(typ(x))
3593 : {
3594 820967 : case t_REAL: return mpsin(x);
3595 1463863 : case t_COMPLEX:
3596 1463863 : a = gel(x,1);
3597 1463863 : b = gel(x,2);
3598 1463863 : if (isintzero(a)) retmkcomplex(gen_0,gsinh(b,prec));
3599 1448554 : i = precision(x); if (i) prec = i;
3600 1448554 : y = cgetc(prec); av = avma;
3601 1448554 : if (typ(b) != t_REAL) b = gtofp(b, prec);
3602 1448554 : mpsinhcosh(b, &u1, &v1);
3603 1448554 : if (typ(a) != t_REAL) a = gtofp(a, prec);
3604 1448554 : mpsincos(a, &u, &v);
3605 1448554 : affrr_fixlg(gmul(v1,u), gel(y,1));
3606 1448554 : affrr_fixlg(gmul(u1,v), gel(y,2)); return gc_const(av,y);
3607 :
3608 5118 : case t_INT: case t_FRAC:
3609 5118 : y = cgetr(prec); av = avma;
3610 5118 : affrr_fixlg(mpsin(tofp_safe(x,prec)), y); return gc_const(av,y);
3611 :
3612 49 : case t_PADIC: y = sin_p(x);
3613 49 : if (!y) pari_err_DOMAIN("gsin(t_PADIC)","argument","",gen_0,x);
3614 42 : return y;
3615 :
3616 231 : default:
3617 231 : av = avma; if (!(y = toser_i(x))) break;
3618 224 : if (gequal0(y)) return gerepilecopy(av, y);
3619 224 : if (valser(y) < 0)
3620 7 : pari_err_DOMAIN("sin","valuation", "<", gen_0, x);
3621 217 : gsincos(y,&u,&v,prec);
3622 217 : return gerepilecopy(av,u);
3623 : }
3624 7 : return trans_eval("sin",gsin,x,prec);
3625 : }
3626 : /********************************************************************/
3627 : /** SINE, COSINE together **/
3628 : /********************************************************************/
3629 :
3630 : void
3631 13310559 : mpsincos(GEN x, GEN *s, GEN *c)
3632 : {
3633 : long mod8;
3634 : pari_sp av, tetpil;
3635 : GEN z, *gptr[2];
3636 :
3637 13310559 : if (!signe(x))
3638 : {
3639 3994 : long e = expo(x);
3640 3994 : *s = real_0_bit(e);
3641 3994 : *c = e >= 0? real_0_bit(e): real_1_bit(-e);
3642 3994 : return;
3643 : }
3644 :
3645 13306565 : av = avma; z = mpcosm1(x, &mod8); tetpil = avma;
3646 13308927 : switch(mod8)
3647 : {
3648 4035246 : case 0: *c = addsr( 1,z); *s = mpaut(z); break;
3649 828951 : case 1: *s = addsr( 1,z); *c = mpaut(z); togglesign(*c); break;
3650 1095834 : case 2: *c = subsr(-1,z); *s = mpaut(z); togglesign(*s); break;
3651 633907 : case 3: *s = subsr(-1,z); *c = mpaut(z); break;
3652 3235742 : case 4: *c = addsr( 1,z); *s = mpaut(z); togglesign(*s); break;
3653 792527 : case 5: *s = addsr( 1,z); *c = mpaut(z); break;
3654 2011923 : case 6: *c = subsr(-1,z); *s = mpaut(z); break;
3655 677958 : case 7: *s = subsr(-1,z); *c = mpaut(z); togglesign(*c); break;
3656 : }
3657 13308882 : gptr[0] = s; gptr[1] = c; gerepilemanysp(av,tetpil,gptr,2);
3658 : }
3659 :
3660 : /* SINE and COSINE - 1 */
3661 : void
3662 6396 : mpsincosm1(GEN x, GEN *s, GEN *c)
3663 : {
3664 : long mod8;
3665 : pari_sp av, tetpil;
3666 : GEN z, *gptr[2];
3667 :
3668 6396 : if (!signe(x))
3669 : {
3670 0 : long e = expo(x);
3671 0 : *s = real_0_bit(e);
3672 0 : *c = real_0_bit(2*e-1);
3673 0 : return;
3674 : }
3675 6396 : av = avma; z = mpcosm1(x,&mod8); tetpil = avma;
3676 6396 : switch(mod8)
3677 : {
3678 5430 : case 0: *c = rcopy(z); *s = mpaut(z); break;
3679 42 : case 1: *s = addsr(1,z); *c = addrs(mpaut(z),1); togglesign(*c); break;
3680 0 : case 2: *c = subsr(-2,z); *s = mpaut(z); togglesign(*s); break;
3681 0 : case 3: *s = subsr(-1,z); *c = subrs(mpaut(z),1); break;
3682 819 : case 4: *c = rcopy(z); *s = mpaut(z); togglesign(*s); break;
3683 91 : case 5: *s = addsr( 1,z); *c = subrs(mpaut(z),1); break;
3684 7 : case 6: *c = subsr(-2,z); *s = mpaut(z); break;
3685 7 : case 7: *s = subsr(-1,z); *c = subsr(-1,mpaut(z)); break;
3686 : }
3687 6396 : gptr[0] = s; gptr[1] = c;
3688 6396 : gerepilemanysp(av,tetpil,gptr,2);
3689 : }
3690 :
3691 : /* return exp(ix), x a t_REAL */
3692 : GEN
3693 872303 : expIr(GEN x)
3694 : {
3695 872303 : pari_sp av = avma;
3696 872303 : GEN v = cgetg(3,t_COMPLEX);
3697 872304 : mpsincos(x, (GEN*)(v+2), (GEN*)(v+1));
3698 872304 : if (!signe(gel(v,2))) return gerepilecopy(av, gel(v,1));
3699 869619 : return v;
3700 : }
3701 :
3702 : /* return exp(ix)-1, x a t_REAL */
3703 : static GEN
3704 6396 : expm1_Ir(GEN x)
3705 : {
3706 6396 : pari_sp av = avma;
3707 6396 : GEN v = cgetg(3,t_COMPLEX);
3708 6396 : mpsincosm1(x, (GEN*)(v+2), (GEN*)(v+1));
3709 6396 : if (!signe(gel(v,2))) return gerepilecopy(av, gel(v,1));
3710 6396 : return v;
3711 : }
3712 :
3713 : /* return exp(z)-1, z complex */
3714 : GEN
3715 6452 : cxexpm1(GEN z, long prec)
3716 : {
3717 6452 : pari_sp av = avma;
3718 6452 : GEN X, Y, x = real_i(z), y = imag_i(z);
3719 6452 : long l = precision(z);
3720 6452 : if (l) prec = l;
3721 6452 : if (typ(x) != t_REAL) x = gtofp(x, prec);
3722 6452 : if (typ(y) != t_REAL) y = gtofp(y, prec);
3723 6452 : if (gequal0(y)) return mpexpm1(x);
3724 6396 : if (gequal0(x)) return expm1_Ir(y);
3725 6263 : X = mpexpm1(x); /* t_REAL */
3726 6263 : Y = expm1_Ir(y);
3727 : /* exp(x+iy) - 1 = (exp(x)-1)(exp(iy)-1) + exp(x)-1 + exp(iy)-1 */
3728 6263 : return gerepileupto(av, gadd(gadd(X,Y), gmul(X,Y)));
3729 : }
3730 :
3731 : void
3732 4246200 : gsincos(GEN x, GEN *s, GEN *c, long prec)
3733 : {
3734 : long i, j, ex, ex2, lx, ly, mi;
3735 : pari_sp av, tetpil;
3736 : GEN y, r, u, v, u1, v1, p1, p2, p3, p4, ps, pc;
3737 : GEN *gptr[4];
3738 :
3739 4246200 : switch(typ(x))
3740 : {
3741 6945 : case t_INT: case t_FRAC:
3742 6945 : *s = cgetr(prec);
3743 6943 : *c = cgetr(prec); av = avma;
3744 6939 : mpsincos(tofp_safe(x, prec), &ps, &pc);
3745 6956 : affrr_fixlg(ps,*s);
3746 4246348 : affrr_fixlg(pc,*c); set_avma(av); return;
3747 :
3748 4199468 : case t_REAL:
3749 4199468 : mpsincos(x,s,c); return;
3750 :
3751 39268 : case t_COMPLEX:
3752 39268 : i = precision(x); if (i) prec = i;
3753 39268 : ps = cgetc(prec); *s = ps;
3754 39268 : pc = cgetc(prec); *c = pc; av = avma;
3755 39268 : r = gexp(gel(x,2),prec);
3756 39268 : v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
3757 39268 : u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
3758 39268 : gsincos(gel(x,1), &u,&v, prec);
3759 39268 : affrr_fixlg(mulrr(v1,u), gel(ps,1));
3760 39268 : affrr_fixlg(mulrr(u1,v), gel(ps,2));
3761 39268 : affrr_fixlg(mulrr(v1,v), gel(pc,1));
3762 39268 : affrr_fixlg(mulrr(u1,u), gel(pc,2)); togglesign(gel(pc,2));
3763 39268 : set_avma(av); return;
3764 :
3765 0 : case t_QUAD:
3766 0 : av = avma; gsincos(quadtofp(x, prec), s, c, prec);
3767 0 : gerepileall(av, 2, s, c); return;
3768 :
3769 519 : default:
3770 519 : av = avma; if (!(y = toser_i(x))) break;
3771 524 : if (gequal0(y)) { *s = gerepilecopy(av,y); *c = gaddsg(1,*s); return; }
3772 :
3773 524 : ex = valser(y); lx = lg(y); ex2 = 2*ex+2;
3774 524 : if (ex < 0) pari_err_DOMAIN("gsincos","valuation", "<", gen_0, x);
3775 524 : if (ex2 > lx)
3776 : {
3777 98 : *s = x == y? gcopy(y): gerepilecopy(av, y); av = avma;
3778 98 : *c = gerepileupto(av, gsubsg(1, gdivgu(gsqr(y),2)));
3779 98 : return;
3780 : }
3781 426 : if (!ex)
3782 : {
3783 108 : gsincos(serchop0(y),&u,&v,prec);
3784 108 : gsincos(gel(y,2),&u1,&v1,prec);
3785 108 : p1 = gmul(v1,v);
3786 108 : p2 = gmul(u1,u);
3787 108 : p3 = gmul(v1,u);
3788 108 : p4 = gmul(u1,v); tetpil = avma;
3789 108 : *c = gsub(p1,p2);
3790 108 : *s = gadd(p3,p4);
3791 108 : gptr[0]=s; gptr[1]=c;
3792 108 : gerepilemanysp(av,tetpil,gptr,2);
3793 108 : return;
3794 : }
3795 :
3796 318 : ly = lx+2*ex;
3797 3111 : mi = lx-1; while (mi>=3 && isrationalzero(gel(y,mi))) mi--;
3798 318 : mi += ex-2;
3799 318 : pc = cgetg(ly,t_SER); *c = pc;
3800 318 : ps = cgetg(lx,t_SER); *s = ps;
3801 318 : pc[1] = evalsigne(1) | _evalvalser(0) | evalvarn(varn(y));
3802 318 : gel(pc,2) = gen_1; ps[1] = y[1];
3803 643 : for (i=2; i<ex+2; i++) gel(ps,i) = gcopy(gel(y,i));
3804 650 : for (i=3; i< ex2; i++) gel(pc,i) = gen_0;
3805 3625 : for (i=ex2; i<ly; i++)
3806 : {
3807 3307 : long ii = i-ex;
3808 3307 : av = avma; p1 = gen_0;
3809 7566 : for (j=ex; j<=minss(ii-2,mi); j++)
3810 4259 : p1 = gadd(p1, gmulgu(gmul(gel(y,j-ex+2),gel(ps,ii-j)),j));
3811 3307 : gel(pc,i) = gerepileupto(av, gdivgs(p1,2-i));
3812 3307 : if (ii < lx)
3813 : {
3814 2982 : av = avma; p1 = gen_0;
3815 6283 : for (j=ex; j<=minss(i-ex2,mi); j++)
3816 3301 : p1 = gadd(p1,gmulgu(gmul(gel(y,j-ex+2),gel(pc,i-j)),j));
3817 2982 : p1 = gdivgu(p1,i-2);
3818 2982 : gel(ps,ii) = gerepileupto(av, gadd(p1,gel(y,ii)));
3819 : }
3820 : }
3821 318 : return;
3822 : }
3823 0 : pari_err_TYPE("gsincos",x);
3824 : }
3825 :
3826 : /********************************************************************/
3827 : /** **/
3828 : /** SINC **/
3829 : /** **/
3830 : /********************************************************************/
3831 : static GEN
3832 2319450 : mpsinc(GEN x)
3833 : {
3834 2319450 : pari_sp av = avma;
3835 : GEN s, c;
3836 :
3837 2319450 : if (!signe(x)) {
3838 0 : long l = nbits2prec(-expo(x));
3839 0 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
3840 0 : return real_1(l);
3841 : }
3842 :
3843 2319450 : mpsincos(x,&s,&c);
3844 2319450 : return gerepileuptoleaf(av, divrr(s,x));
3845 : }
3846 :
3847 : GEN
3848 2319562 : gsinc(GEN x, long prec)
3849 : {
3850 : pari_sp av;
3851 : GEN r, u, v, y, u1, v1;
3852 : long i;
3853 :
3854 2319562 : switch(typ(x))
3855 : {
3856 2319429 : case t_REAL: return mpsinc(x);
3857 49 : case t_COMPLEX:
3858 49 : if (isintzero(gel(x,1)))
3859 : {
3860 28 : av = avma; x = gel(x,2);
3861 28 : if (gequal0(x)) return gcosh(x,prec);
3862 14 : return gerepileuptoleaf(av,gdiv(gsinh(x,prec),x));
3863 : }
3864 21 : i = precision(x); if (i) prec = i;
3865 21 : y = cgetc(prec); av = avma;
3866 21 : r = gexp(gel(x,2),prec);
3867 21 : v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
3868 21 : u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
3869 21 : gsincos(gel(x,1),&u,&v,prec);
3870 21 : affc_fixlg(gdiv(mkcomplex(gmul(v1,u), gmul(u1,v)), x), y);
3871 21 : return gc_const(av,y);
3872 :
3873 14 : case t_INT:
3874 14 : if (!signe(x)) return real_1(prec); /*fall through*/
3875 : case t_FRAC:
3876 21 : y = cgetr(prec); av = avma;
3877 21 : affrr_fixlg(mpsinc(tofp_safe(x,prec)), y); return gc_const(av,y);
3878 :
3879 21 : case t_PADIC:
3880 21 : if (gequal0(x)) return cvtop(gen_1, gel(x,2), valp(x));
3881 14 : av = avma; y = sin_p(x);
3882 14 : if (!y) pari_err_DOMAIN("gsinc(t_PADIC)","argument","",gen_0,x);
3883 7 : return gerepileupto(av,gdiv(y,x));
3884 :
3885 35 : default:
3886 : {
3887 : long ex;
3888 35 : av = avma; if (!(y = toser_i(x))) break;
3889 35 : if (gequal0(y)) return gerepileupto(av, gaddsg(1,y));
3890 35 : ex = valser(y);
3891 35 : if (ex < 0) pari_err_DOMAIN("sinc","valuation", "<", gen_0, x);
3892 28 : if (ex)
3893 : {
3894 28 : gsincos(y,&u,&v,prec);
3895 28 : y = gerepileupto(av, gdiv(u,y));
3896 28 : if (lg(y) > 2) gel(y,2) = gen_1;
3897 28 : return y;
3898 : }
3899 : else
3900 : {
3901 0 : GEN z0, y0 = gel(y,2), y1 = serchop0(y), y10 = y1;
3902 0 : if (!gequal1(y0)) y10 = gdiv(y10, y0);
3903 0 : gsincos(y1,&u,&v,prec);
3904 0 : z0 = gdiv(gcos(y0,prec), y0);
3905 0 : y = gaddsg(1, y10);
3906 0 : u = gadd(gmul(gsinc(y0, prec),v), gmul(z0, u));
3907 0 : return gerepileupto(av,gdiv(u,y));
3908 : }
3909 : }
3910 : }
3911 0 : return trans_eval("sinc",gsinc,x,prec);
3912 : }
3913 :
3914 : /********************************************************************/
3915 : /** **/
3916 : /** TANGENT and COTANGENT **/
3917 : /** **/
3918 : /********************************************************************/
3919 : static GEN
3920 133 : mptan(GEN x)
3921 : {
3922 133 : pari_sp av = avma;
3923 : GEN s, c;
3924 :
3925 133 : mpsincos(x,&s,&c);
3926 133 : if (!signe(c))
3927 0 : pari_err_DOMAIN("tan", "argument", "=", strtoGENstr("Pi/2 + kPi"),x);
3928 133 : return gerepileuptoleaf(av, divrr(s,c));
3929 : }
3930 :
3931 : /* If exp(-|im(x)|) << 1, avoid overflow in sincos(x) */
3932 : static int
3933 39284 : tan_huge_im(GEN ix, long prec)
3934 : {
3935 39284 : long b, p = precision(ix);
3936 39284 : if (!p) p = prec;
3937 39284 : b = prec2nbits(p);
3938 39284 : return (gexpo(ix) > b || fabs(gtodouble(ix)) > (M_LN2 / 2) * b);
3939 : }
3940 : /* \pm I */
3941 : static GEN
3942 163 : real_I(long s, long prec)
3943 : {
3944 163 : GEN z = cgetg(3, t_COMPLEX);
3945 163 : gel(z,1) = real_0(prec);
3946 163 : gel(z,2) = s > 0? real_1(prec): real_m1(prec); return z;
3947 : }
3948 :
3949 : GEN
3950 224 : gtan(GEN x, long prec)
3951 : {
3952 : pari_sp av;
3953 : GEN y, s, c;
3954 :
3955 224 : switch(typ(x))
3956 : {
3957 126 : case t_REAL: return mptan(x);
3958 :
3959 42 : case t_COMPLEX: {
3960 42 : if (isintzero(gel(x,1))) retmkcomplex(gen_0,gtanh(gel(x,2),prec));
3961 28 : if (tan_huge_im(gel(x,2), prec)) return real_I(gsigne(gel(x,2)), prec);
3962 14 : av = avma; y = mulcxmI(gtanh(mulcxI(x), prec)); /* tan x = -I th(I x) */
3963 14 : gel(y,1) = gcopy(gel(y,1)); return gerepileupto(av, y);
3964 : }
3965 7 : case t_INT: case t_FRAC:
3966 7 : y = cgetr(prec); av = avma;
3967 7 : affrr_fixlg(mptan(tofp_safe(x,prec)), y); return gc_const(av,y);
3968 :
3969 14 : case t_PADIC:
3970 14 : av = avma;
3971 14 : return gerepileupto(av, gdiv(gsin(x,prec), gcos(x,prec)));
3972 :
3973 35 : default:
3974 35 : av = avma; if (!(y = toser_i(x))) break;
3975 28 : if (gequal0(y)) return gerepilecopy(av, y);
3976 28 : if (valser(y) < 0)
3977 7 : pari_err_DOMAIN("tan","valuation", "<", gen_0, x);
3978 21 : gsincos(y,&s,&c,prec);
3979 21 : return gerepileupto(av, gdiv(s,c));
3980 : }
3981 7 : return trans_eval("tan",gtan,x,prec);
3982 : }
3983 :
3984 : static GEN
3985 70 : mpcotan(GEN x)
3986 : {
3987 70 : pari_sp av=avma, tetpil;
3988 : GEN s,c;
3989 :
3990 70 : mpsincos(x,&s,&c); tetpil=avma;
3991 70 : return gerepile(av,tetpil,divrr(c,s));
3992 : }
3993 :
3994 : GEN
3995 39483 : gcotan(GEN x, long prec)
3996 : {
3997 : pari_sp av;
3998 : GEN y, s, c;
3999 :
4000 39483 : switch(typ(x))
4001 : {
4002 63 : case t_REAL:
4003 63 : return mpcotan(x);
4004 :
4005 39277 : case t_COMPLEX:
4006 39277 : if (isintzero(gel(x,1))) {
4007 21 : GEN z = cgetg(3, t_COMPLEX);
4008 21 : gel(z,1) = gen_0; av = avma;
4009 21 : gel(z,2) = gerepileupto(av, gneg(ginv(gtanh(gel(x,2),prec))));
4010 21 : return z;
4011 : }
4012 39256 : if (tan_huge_im(gel(x,2), prec)) return real_I(-gsigne(gel(x,2)), prec);
4013 39107 : av = avma; gsincos(x,&s,&c,prec);
4014 39107 : return gerepileupto(av, gdiv(c,s));
4015 :
4016 7 : case t_INT: case t_FRAC:
4017 7 : y = cgetr(prec); av = avma;
4018 7 : affrr_fixlg(mpcotan(tofp_safe(x,prec)), y); return gc_const(av,y);
4019 :
4020 14 : case t_PADIC:
4021 14 : av = avma;
4022 14 : return gerepileupto(av, gdiv(gcos(x,prec), gsin(x,prec)));
4023 :
4024 122 : default:
4025 122 : av = avma; if (!(y = toser_i(x))) break;
4026 115 : if (gequal0(y)) pari_err_DOMAIN("cotan", "argument", "=", gen_0, y);
4027 115 : if (valser(y) < 0) pari_err_DOMAIN("cotan","valuation", "<", gen_0, x);
4028 108 : gsincos(y,&s,&c,prec);
4029 108 : return gerepileupto(av, gdiv(c,s));
4030 : }
4031 7 : return trans_eval("cotan",gcotan,x,prec);
4032 : }
|