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