1 /* $Id: trans1.c 12114 2010-02-03 22:59:09Z bill $
2
3 Copyright (C) 2000 The PARI group.
4
5 This file is part of the PARI/GP package.
6
7 PARI/GP is free software; you can redistribute it and/or modify it under the
8 terms of the GNU General Public License as published by the Free Software
9 Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10 ANY WARRANTY WHATSOEVER.
11
12 Check the License for details. You should have received a copy of it, along
13 with the package; see the file 'COPYING'. If not, write to the Free Software
14 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15
16 /********************************************************************/
17 /** **/
18 /** TRANSCENDENTAL FUNCTIONS **/
19 /** **/
20 /********************************************************************/
21 #include "pari.h"
22 #include "paripriv.h"
23
24 #ifdef LONG_IS_64BIT
25 # define SQRTVERYBIGINT 3037000500 /* ceil(sqrt(VERYBIGINT)) */
26 # define CBRTVERYBIGINT 2097152 /* ceil(cbrt(VERYBIGINT)) */
27 #else
28 # define SQRTVERYBIGINT 46341
29 # define CBRTVERYBIGINT 1291
30 #endif
31
32 static GEN glog2;
33 void
pari_init_floats(void)34 pari_init_floats(void)
35 {
36 geuler = gpi = bernzone = glog2 = NULL;
37 }
38
39 /********************************************************************/
40 /** **/
41 /** PI **/
42 /** **/
43 /********************************************************************/
44 #if 0
45 /* Ramanujan's formula:
46 * ----
47 * 53360 (640320)^(1/2) \ (6n)! (545140134 n + 13591409)
48 * -------------------- = / ------------------------------
49 * Pi ---- (n!)^3 (3n)! (-640320)^(3n)
50 * n>=0
51 */
52 GEN
53 piold(long prec)
54 {
55 const long k1 = 545140134, k2 = 13591409, k3 = 640320;
56 const double alpha2 = 47.11041314/BITS_IN_LONG; /* 3log(k3/12) / log(2^BIL) */
57 GEN p1,p2,p3,tmppi;
58 long l, n, n1;
59 pari_sp av = avma, av2;
60 double alpha;
61
62 tmppi = cgetr(prec);
63 prec++;
64 n = (long)(1 + (prec-2)/alpha2);
65 n1 = 6*n - 1;
66 p2 = addsi(k2, mulss(n,k1));
67 p1 = itor(p2, prec);
68
69 /* initialize mantissa length */
70 if (prec>=4) l=4; else l=prec;
71 setlg(p1,l); alpha = (double)l;
72
73 av2 = avma;
74 while (n)
75 {
76 if (n < CBRTVERYBIGINT) /* p3 = n1(n1-2)(n1-4) / n^3 */
77 p3 = divrs(mulsr(n1-4,mulsr(n1*(n1-2),p1)),n*n*n);
78 else
79 {
80 if (n1 < SQRTVERYBIGINT)
81 p3 = divrs(divrs(mulsr(n1-4,mulsr(n1*(n1-2),p1)),n*n),n);
82 else
83 p3 = divrs(divrs(divrs(mulsr(n1-4,mulsr(n1,mulsr(n1-2,p1))),n),n),n);
84 }
85 p3 = divrs(divrs(p3,100100025), 327843840);
86 subisz(p2,k1,p2);
87 subirz(p2,p3,p1);
88 alpha += alpha2; l = (long)(1+alpha); /* new mantissa length */
89 if (l > prec) l = prec;
90 setlg(p1,l); avma = av2;
91 n--; n1-=6;
92 }
93 p1 = divsr(53360,p1);
94 return gerepileuptoleaf(av, mulrr(p1,sqrtr_abs(stor(k3,prec))));
95 }
96 #endif
97 /* Gauss - Brent-Salamin AGM iteration */
98 void
constpi(long prec)99 constpi(long prec)
100 {
101 GEN A, B, C, tmppi;
102 long i, G;
103 pari_sp av, av2;
104
105 if (gpi && lg(gpi) >= prec) return;
106
107 av = avma; tmppi = newbloc(prec);
108 *tmppi = evaltyp(t_REAL) | evallg(prec);
109 G = - bit_accuracy(prec);
110 prec++;
111
112 A = real_1(prec);
113 B = sqrtr_abs(real2n(1,prec)); setexpo(B, -1); /* = 1/sqrt(2) */
114 C = real2n(-2, prec); av2 = avma;
115 for (i = 0;; i++)
116 {
117 GEN y, a, b, B_A = subrr(B, A);
118 if (expo(B_A) < G) break;
119 a = addrr(A,B); setexpo(a, expo(a)-1);
120 b = sqrtr_abs( mulrr(A, B) );
121 y = gsqr(B_A); setexpo(y, expo(y) + i - 2);
122 affrr(subrr(C, y), C);
123 affrr(a, A);
124 affrr(b, B); avma = av2;
125 }
126 setexpo(C, expo(C)+2);
127 affrr(divrr(gsqr(addrr(A,B)), C), tmppi);
128 if (gpi) gunclone(gpi);
129 avma = av; gpi = tmppi;
130 }
131
132 GEN
mppi(long prec)133 mppi(long prec)
134 {
135 GEN x = cgetr(prec);
136 constpi(prec); affrr(gpi,x); return x;
137 }
138
139 /* Pi * 2^n */
140 GEN
Pi2n(long n,long prec)141 Pi2n(long n, long prec)
142 {
143 GEN x = mppi(prec); setexpo(x, 1+n);
144 return x;
145 }
146
147 /* I * Pi * 2^n */
148 GEN
PiI2n(long n,long prec)149 PiI2n(long n, long prec)
150 {
151 GEN z = cgetg(3,t_COMPLEX);
152 gel(z,1) = gen_0;
153 gel(z,2) = Pi2n(n, prec); return z;
154 }
155
156 /* 2I * Pi */
157 GEN
PiI2(long prec)158 PiI2(long prec) { return PiI2n(1, prec); }
159
160 /********************************************************************/
161 /** **/
162 /** EULER CONSTANT **/
163 /** **/
164 /********************************************************************/
165
166 void
consteuler(long prec)167 consteuler(long prec)
168 {
169 GEN u,v,a,b,tmpeuler;
170 long l, n1, n, k, x;
171 pari_sp av1, av2;
172
173 if (geuler && lg(geuler) >= prec) return;
174
175 av1 = avma; tmpeuler = newbloc(prec);
176 *tmpeuler = evaltyp(t_REAL) | evallg(prec);
177
178 prec++;
179
180 l = prec+1; x = (long) (1 + bit_accuracy_mul(l, LOG2/4));
181 a = stor(x,l); u=logr_abs(a); setsigne(u,-1); affrr(u,a);
182 b = real_1(l);
183 v = real_1(l);
184 n = (long)(1+3.591*x); /* z=3.591: z*[ ln(z)-1 ]=1 */
185 n1 = min(n, SQRTVERYBIGINT);
186 if (x < SQRTVERYBIGINT)
187 {
188 long xx = x*x;
189 av2 = avma;
190 for (k=1; k<n1; k++)
191 {
192 divrsz(mulsr(xx,b),k*k, b);
193 divrsz(addrr(divrs(mulsr(xx,a),k),b),k, a);
194 addrrz(u,a,u);
195 addrrz(v,b,v); avma = av2;
196 }
197 for ( ; k<=n; k++)
198 {
199 divrsz(divrs(mulsr(xx,b),k), k, b);
200 divrsz(addrr(divrs(mulsr(xx,a),k),b),k, a);
201 addrrz(u,a,u);
202 addrrz(v,b,v); avma = av2;
203 }
204 }
205 else
206 {
207 GEN xx = mulss(x,x);
208 av2 = avma;
209 for (k=1; k<n1; k++)
210 {
211 divrsz(mulir(xx,b),k*k, b);
212 divrsz(addrr(divrs(mulir(xx,a),k),b),k, a);
213 addrrz(u,a,u);
214 addrrz(v,b,v); avma = av2;
215 }
216 for ( ; k<=n; k++)
217 {
218 divrsz(divrs(mulir(xx,b),k), k, b);
219 divrsz(addrr(divrs(mulir(xx,a),k),b),k, a);
220 addrrz(u,a,u);
221 addrrz(v,b,v); avma = av2;
222 }
223 }
224 divrrz(u,v,tmpeuler);
225 if (geuler) gunclone(geuler);
226 avma = av1; geuler = tmpeuler;
227 }
228
229 GEN
mpeuler(long prec)230 mpeuler(long prec)
231 {
232 GEN x = cgetr(prec);
233 consteuler(prec); affrr(geuler,x); return x;
234 }
235
236 /********************************************************************/
237 /** **/
238 /** TYPE CONVERSION FOR TRANSCENDENTAL FUNCTIONS **/
239 /** **/
240 /********************************************************************/
241
242 GEN
transc(GEN (* f)(GEN,long),GEN x,long prec)243 transc(GEN (*f)(GEN,long), GEN x, long prec)
244 {
245 pari_sp tetpil, av = avma;
246 GEN p1, y;
247 long lx, i;
248
249 if (prec < 2) pari_err(talker, "incorrect precision in transc");
250 switch(typ(x))
251 {
252 case t_INT:
253 p1 = itor(x, prec); tetpil=avma;
254 return gerepile(av,tetpil,f(p1,prec));
255
256 case t_FRAC:
257 p1 = rdivii(gel(x,1), gel(x,2), prec); tetpil=avma;
258 return gerepile(av,tetpil,f(p1,prec));
259
260 case t_QUAD:
261 p1 = quadtoc(x, prec); tetpil = avma;
262 return gerepile(av,tetpil,f(p1,prec));
263
264 case t_POL: case t_RFRAC:
265 return gerepileupto(av, f(toser_i(x), prec));
266
267 case t_VEC: case t_COL: case t_MAT:
268 lx = lg(x); y = cgetg(lx,typ(x));
269 for (i=1; i<lx; i++) gel(y,i) = f(gel(x,i),prec);
270 return y;
271
272 case t_POLMOD:
273 p1 = cleanroots(gel(x,1),prec); lx = lg(p1);
274 for (i=1; i<lx; i++) gel(p1,i) = poleval(gel(x,2),gel(p1,i));
275 tetpil = avma; y = cgetg(lx,t_COL);
276 for (i=1; i<lx; i++) gel(y,i) = f(gel(p1,i),prec);
277 return gerepile(av,tetpil,y);
278
279 default: pari_err(typeer,"a transcendental function");
280 }
281 return f(x,prec);
282 }
283
284 /*******************************************************************/
285 /* */
286 /* POWERING */
287 /* */
288 /*******************************************************************/
289 static GEN
puiss0(GEN x)290 puiss0(GEN x)
291 {
292 long lx, i;
293 GEN y;
294
295 switch(typ(x))
296 {
297 case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX:
298 case t_PADIC: case t_QUAD:
299 return gen_1;
300
301 case t_INTMOD:
302 y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
303 gel(y,2) = gen_1; return y;
304
305 case t_POLMOD:
306 y = cgetg(3,t_POLMOD); gel(y,1) = gcopy(gel(x,1));
307 gel(y,2) = pol_1[varn(x[1])]; return y;
308
309 case t_POL: case t_SER: case t_RFRAC:
310 return pol_1[gvar(x)];
311
312 case t_MAT:
313 lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
314 if (lx != lg(x[1])) pari_err(mattype1,"gpow");
315 y = matid(lx-1);
316 for (i=1; i<lx; i++) gcoeff(y,i,i) = puiss0(gcoeff(x,i,i));
317 return y;
318 case t_QFR: return qfr_unit(x);
319 case t_QFI: return qfi_unit(x);
320 case t_VECSMALL: return perm_identity(lg(x) - 1);
321 }
322 pari_err(typeer,"gpow");
323 return NULL; /* not reached */
324 }
325
326 static GEN
_sqr(void * data,GEN x)327 _sqr(void *data /* ignored */, GEN x) { (void)data; return gsqr(x); }
328 static GEN
_mul(void * data,GEN x,GEN y)329 _mul(void *data /* ignored */, GEN x, GEN y) { (void)data; return gmul(x,y); }
330 static GEN
_sqri(void * data,GEN x)331 _sqri(void *data /* ignored */, GEN x) { (void)data; return sqri(x); }
332 static GEN
_muli(void * data,GEN x,GEN y)333 _muli(void *data /* ignored */, GEN x, GEN y) { (void)data; return mulii(x,y); }
334
335 /* INTEGER POWERING (a^n for integer a != 0 and integer n > 0)
336 *
337 * Use left shift binary algorithm (RS is wasteful: multiplies big numbers,
338 * with LS one of them is the base, hence small). Sign of result is set
339 * to s (= 1,-1). Makes life easier for caller, which otherwise might do a
340 * setsigne(gen_1 / gen_m1) */
341 static GEN
powiu_sign(GEN a,ulong N,long s)342 powiu_sign(GEN a, ulong N, long s)
343 {
344 pari_sp av;
345 GEN y;
346
347 if (lgefint(a) == 3)
348 { /* easy if |a| < 3 */
349 if (a[2] == 1) return (s>0)? gen_1: gen_m1;
350 if (a[2] == 2) { a = int2u(N); setsigne(a,s); return a; }
351 }
352 if (N == 1) { a = icopy(a); setsigne(a,s); return a; }
353 if (N == 2) return sqri(a);
354 av = avma;
355 y = leftright_pow_u(a, N, NULL, &_sqri, &_muli);
356 setsigne(y,s); return gerepileuptoint(av, y);
357 }
358 /* a^N */
359 GEN
powiu(GEN a,ulong N)360 powiu(GEN a, ulong N)
361 {
362 long s;
363 if (!N) return gen_1;
364 s = signe(a);
365 if (!s) return gen_0;
366 return powiu_sign(a, N, (s < 0 && odd(N))? -1: 1);
367 }
368 GEN
powuu(ulong p,ulong N)369 powuu(ulong p, ulong N)
370 {
371 long P[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
372 if (!N) return gen_1;
373 if (!p) return gen_0;
374 P[2] = p;
375 return powiu_sign(P, N, 1);
376 }
377
378 /* assume p^k is SMALL */
379 ulong
upowuu(ulong p,ulong k)380 upowuu(ulong p, ulong k)
381 {
382 ulong i, pk;
383
384 if (!k) return 1;
385 if (p == 2) return 1UL<<k;
386 pk = p; for (i=2; i<=k; i++) pk *= p;
387 return pk;
388 }
389
390 typedef struct {
391 long prec, a;
392 GEN (*sqr)(GEN);
393 GEN (*mulug)(ulong,GEN);
394 } sr_muldata;
395
396 static GEN
_rpowuu_mul(void * data,GEN x,GEN y)397 _rpowuu_mul(void *data, GEN x, GEN y/*unused*/)
398 {
399 sr_muldata *D = (sr_muldata *)data;
400 (void)y; return D->mulug(D->a, x);
401 }
402
403 static GEN
_rpowuu_sqr(void * data,GEN x)404 _rpowuu_sqr(void *data, GEN x)
405 {
406 sr_muldata *D = (sr_muldata *)data;
407 if (typ(x) == t_INT && lgefint(x) >= D->prec)
408 { /* switch to t_REAL */
409 D->sqr = &gsqr;
410 D->mulug = &mulur; x = itor(x, D->prec);
411 }
412 return D->sqr(x);
413 }
414
415 /* return a^n as a t_REAL of precision prec. Assume a > 0, n > 0 */
416 GEN
rpowuu(ulong a,ulong n,long prec)417 rpowuu(ulong a, ulong n, long prec)
418 {
419 pari_sp av;
420 GEN y;
421 sr_muldata D;
422
423 if (a == 1) return real_1(prec);
424 if (a == 2) return real2n(n, prec);
425 if (n == 1) return stor(a, prec);
426 av = avma;
427 D.sqr = &sqri;
428 D.mulug = &mului;
429 D.prec = prec;
430 D.a = (long)a;
431 y = leftright_pow_u(utoipos(a), n, (void*)&D, &_rpowuu_sqr, &_rpowuu_mul);
432 if (typ(y) == t_INT) y = itor(y, prec);
433 return gerepileuptoleaf(av, y);
434 }
435
436 /* x^(s/2), assume x t_REAL */
437 GEN
powrshalf(GEN x,long s)438 powrshalf(GEN x, long s)
439 {
440 if (s & 1) return sqrtr(gpowgs(x, s));
441 return gpowgs(x, s>>1);
442 }
443 /* x^(n/d), assume x t_REAL, return t_REAL */
444 GEN
powrfrac(GEN x,long n,long d)445 powrfrac(GEN x, long n, long d)
446 {
447 long z;
448 if (!n) return real_1(lg(x));
449 z = cgcd(n, d); if (z > 1) { n /= z; d /= z; }
450 if (d == 1) return gpowgs(x, n);
451 x = gpowgs(x, n);
452 if (d == 2) return sqrtr(x);
453 return sqrtnr(x, d);
454 }
455
456 /* assume x != 0 */
457 static GEN
pow_monome(GEN x,long n)458 pow_monome(GEN x, long n)
459 {
460 long i, d, dx = degpol(x);
461 GEN A, b, y;
462
463 if (n < 0) { n = -n; y = cgetg(3, t_RFRAC); } else y = NULL;
464
465 if (HIGHWORD(dx) || HIGHWORD(n))
466 {
467 LOCAL_HIREMAINDER;
468 d = (long)mulll((ulong)dx, (ulong)n);
469 if (hiremainder || (d &~ LGBITS)) d = LGBITS; /* overflow */
470 d += 2;
471 }
472 else
473 d = dx*n + 2;
474 if ((d + 1) & ~LGBITS) pari_err(talker,"degree overflow in pow_monome");
475 A = cgetg(d+1, t_POL); A[1] = x[1];
476 for (i=2; i < d; i++) gel(A,i) = gen_0;
477 b = gpowgs(gel(x,dx+2), n); /* not memory clean if (n < 0) */
478 if (!y) y = A;
479 else {
480 GEN c = denom(b);
481 gel(y,1) = c; if (c != gen_1) b = gmul(b,c);
482 gel(y,2) = A;
483 }
484 gel(A,d) = b; return y;
485 }
486
487 /* x t_PADIC */
488 static GEN
powps(GEN x,long n)489 powps(GEN x, long n)
490 {
491 long e = n*valp(x), v;
492 GEN t, y, mod, p = gel(x,2);
493 pari_sp av;
494
495 if (!signe(x[4])) {
496 if (n < 0) pari_err(gdiver);
497 return zeropadic(p, e);
498 }
499 v = z_pval(n, p);
500
501 y = cgetg(5,t_PADIC);
502 mod = gel(x,3);
503 if (v == 0) mod = icopy(mod);
504 else
505 {
506 if (precp(x) == 1 && equaliu(p, 2)) v++;
507 mod = mulii(mod, powiu(p,v));
508 mod = gerepileuptoint((pari_sp)y, mod);
509 }
510 y[1] = evalprecp(precp(x) + v) | evalvalp(e);
511 gel(y,2) = icopy(p);
512 gel(y,3) = mod;
513
514 av = avma; t = gel(x,4);
515 if (n < 0) { t = Fp_inv(t, mod); n = -n; }
516 t = Fp_powu(t, n, mod);
517 gel(y,4) = gerepileuptoint(av, t);
518 return y;
519 }
520 /* x t_PADIC */
521 static GEN
powp(GEN x,GEN n)522 powp(GEN x, GEN n)
523 {
524 long v;
525 GEN y, mod, p = gel(x,2);
526
527 if (valp(x)) pari_err(errlg);
528
529 if (!signe(x[4])) {
530 if (signe(n) < 0) pari_err(gdiver);
531 return zeropadic(p, 0);
532 }
533 v = Z_pval(n, p);
534
535 y = cgetg(5,t_PADIC);
536 mod = gel(x,3);
537 if (v == 0) mod = icopy(mod);
538 else
539 {
540 mod = mulii(mod, powiu(p,v));
541 mod = gerepileuptoint((pari_sp)y, mod);
542 }
543 y[1] = evalprecp(precp(x) + v) | evalvalp(0);
544 gel(y,2) = icopy(p);
545 gel(y,3) = mod;
546 gel(y,4) = Fp_pow(gel(x,4), n, mod);
547 return y;
548 }
549
550 GEN
gpowgs(GEN x,long n)551 gpowgs(GEN x, long n)
552 {
553 long m;
554 pari_sp av;
555 GEN y;
556
557 if (n == 0) return puiss0(x);
558 if (n == 1) return gcopy(x);
559 if (n ==-1) return ginv(x);
560 switch(typ(x))
561 {
562 case t_INT:
563 {
564 long sx = signe(x), s;
565 GEN t;
566 if (!sx) {
567 if (n < 0) pari_err(gdiver);
568 return gen_0;
569 }
570 s = (sx < 0 && odd(n))? -1: 1;
571 if (n > 0) return powiu_sign(x, n, s);
572 t = (s > 0)? gen_1: gen_m1;
573 if (is_pm1(x)) return t;
574 /* n < 0, |x| > 1 */
575 y = cgetg(3,t_FRAC);
576 gel(y,1) = t;
577 gel(y,2) = powiu_sign(x, -n, 1); /* force denominator > 0 */
578 return y;
579 }
580 case t_INTMOD:
581 y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
582 gel(y,2) = Fp_pows(gel(x,2), n, gel(x,1));
583 return y;
584 case t_FRAC:
585 {
586 GEN a = gel(x,1), b = gel(x,2);
587 long sx = signe(a), s;
588 if (!sx) {
589 if (n < 0) pari_err(gdiver);
590 return gen_0;
591 }
592 s = (sx < 0 && odd(n))? -1: 1;
593 if (n < 0) {
594 n = -n;
595 if (is_pm1(a)) return powiu_sign(b, n, s); /* +-1/x[2] inverts to t_INT */
596 swap(a, b);
597 }
598 y = cgetg(3, t_FRAC);
599 gel(y,1) = powiu_sign(a, n, s);
600 gel(y,2) = powiu_sign(b, n, 1);
601 return y;
602 }
603 case t_PADIC: return powps(x, n);
604 case t_RFRAC:
605 {
606 av = avma; y = cgetg(3, t_RFRAC); m = labs(n);
607 gel(y,1) = gpowgs(gel(x,1),m);
608 gel(y,2) = gpowgs(gel(x,2),m);
609 if (n < 0) y = ginv(y);
610 return gerepileupto(av,y);
611 }
612 case t_POL:
613 if (ismonome(x)) return pow_monome(x, n);
614 default: {
615 pari_sp av = avma;
616 y = leftright_pow_u(x, (ulong)labs(n), NULL, &_sqr, &_mul);
617 if (n < 0) y = ginv(y);
618 return gerepileupto(av,y);
619 }
620 }
621 }
622
623 /* n a t_INT */
624 GEN
powgi(GEN x,GEN n)625 powgi(GEN x, GEN n)
626 {
627 GEN y;
628
629 if (!is_bigint(n)) return gpowgs(x, itos(n));
630 /* probable overflow for non-modular types (typical exception: (X^0)^N) */
631 switch(typ(x))
632 {
633 case t_INTMOD:
634 y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
635 gel(y,2) = Fp_pow(gel(x,2), n, gel(x,1));
636 return y;
637 case t_PADIC: return powp(x, n);
638
639 case t_INT:
640 if (is_pm1(x)) return (signe(x) < 0 && mpodd(n))? gen_m1: gen_1;
641 if (signe(x)) pari_err(errlg);
642 if (signe(n) < 0) pari_err(gdiver);
643 return gen_0;
644 case t_FRAC:
645 if (signe(x[1])) pari_err(errlg);
646 if (signe(n) < 0) pari_err(gdiver);
647 return gen_0;
648
649 case t_QFR: return qfr_pow(x,n);
650 default: {
651 pari_sp av = avma;
652 y = leftright_pow(x, n, NULL, &_sqr, &_mul);
653 if (signe(n) < 0) y = ginv(y);
654 return gerepileupto(av,y);
655 }
656 }
657 }
658
659 /* we suppose n != 0, valp(x) = 0 and leading-term(x) != 0. Not stack clean */
660 static GEN
ser_pow(GEN x,GEN n,long prec)661 ser_pow(GEN x, GEN n, long prec)
662 {
663 pari_sp av, tetpil;
664 GEN y, p1, p2, lead;
665
666 if (gvar(n) <= varn(x)) return gexp(gmul(n, glog(x,prec)), prec);
667 lead = gel(x,2);
668 if (gcmp1(lead))
669 {
670 GEN X, Y, alp = gadd(n,gen_1); /* will be left on the stack */
671 long lx, mi, i, j, d;
672
673 lx = lg(x);
674 y = cgetg(lx,t_SER);
675 y[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(x));
676 X = x+2;
677 Y = y+2;
678
679 d = mi = lx-3; while (mi>=1 && isexactzero(gel(X,mi))) mi--;
680 gel(Y,0) = gen_1;
681 for (i=1; i<=d; i++)
682 {
683 av = avma; p1 = gen_0;
684 for (j=1; j<=min(i,mi); j++)
685 {
686 p2 = gsubgs(gmulgs(alp,j),i);
687 p1 = gadd(p1, gmul(gmul(p2,gel(X,j)),gel(Y,i-j)));
688 }
689 tetpil = avma; gel(Y,i) = gerepile(av,tetpil, gdivgs(p1,i));
690 }
691 return y;
692 }
693 p1 = gdiv(x,lead);
694 if (typ(p1) != t_SER) pari_err(typeer, "ser_pow");
695 gel(p1,2) = gen_1; /* in case it's inexact */
696 if (typ(n) == t_FRAC && !isinexact(lead) && ispower(lead, gel(n,2), &p2))
697 p2 = powgi(p2, gel(n,1));
698 else
699 p2 = gpow(lead,n, prec);
700 return gmul(p2, gpow(p1, n, prec));
701 }
702
703 static long
val_from_i(GEN E)704 val_from_i(GEN E)
705 {
706 if (is_bigint(E)) pari_err(talker,"valuation overflow in sqrtn");
707 return itos(E);
708 }
709
710 /* return x^q, assume typ(x) = t_SER, typ(q) = t_FRAC and q != 0 */
711 static GEN
ser_powfrac(GEN x,GEN q,long prec)712 ser_powfrac(GEN x, GEN q, long prec)
713 {
714 long e = valp(x);
715 GEN y, E = gmulsg(e, q);
716
717 if (gcmp0(x)) return zeroser(varn(x), val_from_i(gfloor(E)));
718 if (typ(E) != t_INT)
719 pari_err(talker,"%Z should divide valuation (= %ld) in sqrtn",q[2], e);
720 e = val_from_i(E);
721 y = shallowcopy(x); setvalp(y, 0);
722 y = ser_pow(y, q, prec);
723 if (typ(y) == t_SER) /* generic case */
724 y[1] = evalsigne(1) | evalvalp(e) | evalvarn(varn(x));
725 else /* e.g coeffs are POLMODs */
726 y = gmul(y, monomial(gen_1, e, varn(x)));
727 return y;
728 }
729
730 GEN
gpow(GEN x,GEN n,long prec)731 gpow(GEN x, GEN n, long prec)
732 {
733 long i, lx, tx, tn = typ(n);
734 pari_sp av;
735 GEN y;
736
737 if (tn == t_INT) return powgi(x,n);
738 tx = typ(x);
739 if (is_matvec_t(tx))
740 {
741 lx = lg(x); y = cgetg(lx,tx);
742 for (i=1; i<lx; i++) gel(y,i) = gpow(gel(x,i),n,prec);
743 return y;
744 }
745 av = avma;
746 if (tx == t_POL || tx == t_RFRAC) { x = toser_i(x); tx = t_SER; }
747 if (tx == t_SER)
748 {
749 if (tn == t_FRAC) return gerepileupto(av, ser_powfrac(x, n, prec));
750 if (valp(x))
751 pari_err(talker,"gpow: need integer exponent if series valuation != 0");
752 if (lg(x) == 2) return gcopy(x); /* O(1) */
753 return gerepileupto(av, ser_pow(x, n, prec));
754 }
755 if (gcmp0(x))
756 {
757 if (!is_scalar_t(tn) || tn == t_PADIC || tn == t_INTMOD)
758 pari_err(talker,"gpow: 0 to a forbidden power");
759 n = real_i(n);
760 if (gsigne(n) <= 0)
761 pari_err(talker,"gpow: 0 to a non positive exponent");
762 if (!precision(x)) return gcopy(x);
763
764 x = ground(gmulsg(gexpo(x),n));
765 if (is_bigint(x) || (ulong)x[2] >= HIGHEXPOBIT)
766 pari_err(talker,"gpow: underflow or overflow");
767 avma = av; return real_0_bit(itos(x));
768 }
769 if (tn == t_FRAC)
770 {
771 GEN z, d = gel(n,2), a = gel(n,1);
772 if (tx == t_INTMOD)
773 {
774 if (!BSW_psp(gel(x,1))) pari_err(talker,"gpow: modulus %Z is not prime",x[1]);
775 y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
776 av = avma;
777 z = Fp_sqrtn(gel(x,2), d, gel(x,1), NULL);
778 if (!z) pari_err(talker,"gpow: nth-root does not exist");
779 gel(y,2) = gerepileuptoint(av, Fp_pow(z, a, gel(x,1)));
780 return y;
781 }
782 else if (tx == t_PADIC)
783 {
784 z = equaliu(d, 2)? padic_sqrt(x): padic_sqrtn(x, d, NULL);
785 if (!z) pari_err(talker, "gpow: nth-root does not exist");
786 return gerepileupto(av, powgi(z, a));
787 }
788 }
789 i = (long) precision(n); if (i) prec=i;
790 y = gmul(n, glog(x,prec));
791 return gerepileupto(av, gexp(y,prec));
792 }
793
794 /********************************************************************/
795 /** **/
796 /** RACINE CARREE **/
797 /** **/
798 /********************************************************************/
799
800 GEN
sqrtr(GEN x)801 sqrtr(GEN x) {
802 long s = signe(x);
803 GEN y;
804 if (typ(x) != t_REAL) pari_err(typeer,"sqrtr");
805 if (s == 0) return real_0_bit(expo(x) >> 1);
806 if (s >= 0) return sqrtr_abs(x);
807 y = cgetg(3,t_COMPLEX);
808 gel(y,2) = sqrtr_abs(x);
809 gel(y,1) = gen_0; return y;
810 }
811
812 /* assume x unit, precp(x) = pp > 3 */
813 static GEN
sqrt_2adic(GEN x,long pp)814 sqrt_2adic(GEN x, long pp)
815 {
816 GEN z = mod16(x)==1? gen_1: utoipos(3);
817 long zp;
818 pari_sp av, lim;
819
820 if (pp == 4) return z;
821 zp = 3; /* number of correct bits in z (compared to sqrt(x)) */
822
823 av = avma; lim = stack_lim(av,2);
824 for(;;)
825 {
826 GEN mod;
827 zp = (zp<<1) - 1;
828 if (zp > pp) zp = pp;
829 mod = int2n(zp);
830 z = addii(z, resmod2n(mulii(x, Fp_inv(z,mod)), zp));
831 z = shifti(z, -1); /* (z + x/z) / 2 */
832 if (pp == zp) return z;
833
834 if (zp < pp) zp--;
835
836 if (low_stack(lim,stack_lim(av,2)))
837 {
838 if (DEBUGMEM > 1) pari_warn(warnmem,"padic_sqrt");
839 z = gerepileuptoint(av,z);
840 }
841 }
842 }
843
844 /* x unit defined modulo modx = p^pp, p != 2, pp > 0 */
845 static GEN
sqrt_padic(GEN x,GEN modx,long pp,GEN p)846 sqrt_padic(GEN x, GEN modx, long pp, GEN p)
847 {
848 GEN mod, z = Fp_sqrt(x, p);
849 long zp = 1;
850 pari_sp av, lim;
851
852 if (!z) pari_err(sqrter5);
853 if (pp <= zp) return z;
854
855 av = avma; lim = stack_lim(av,2);
856 mod = p;
857 for(;;)
858 { /* could use the hensel_lift_accel idea. Not really worth it */
859 GEN inv2;
860 zp <<= 1;
861 if (zp < pp) mod = sqri(mod); else { zp = pp; mod = modx; }
862 inv2 = shifti(addis(mod,1), -1); /* = (mod + 1)/2 = 1/2 */
863 z = addii(z, remii(mulii(x, Fp_inv(z,mod)), mod));
864 z = mulii(z, inv2);
865 z = modii(z, mod); /* (z + x/z) / 2 */
866 if (pp <= zp) return z;
867
868 if (low_stack(lim,stack_lim(av,2)))
869 {
870 GEN *gptr[2]; gptr[0]=&z; gptr[1]=&mod;
871 if (DEBUGMEM>1) pari_warn(warnmem,"padic_sqrt");
872 gerepilemany(av,gptr,2);
873 }
874 }
875 }
876
877 GEN
padic_sqrt(GEN x)878 padic_sqrt(GEN x)
879 {
880 long pp, e = valp(x);
881 GEN z,y,mod, p = gel(x,2);
882 pari_sp av;
883
884 if (gcmp0(x)) return zeropadic(p, (e+1) >> 1);
885 if (e & 1) pari_err(talker,"odd exponent in p-adic sqrt");
886
887 y = cgetg(5,t_PADIC);
888 pp = precp(x);
889 mod = gel(x,3);
890 x = gel(x,4); /* lift to t_INT */
891 e >>= 1; av = avma;
892 if (equaliu(p,2))
893 {
894 long r = mod8(x);
895 if (pp <= 3)
896 {
897 switch(pp) {
898 case 1: break;
899 case 2: if ((r&3) == 1) break;
900 case 3: if (r == 1) break;
901 default: pari_err(sqrter5);
902 }
903 z = gen_1;
904 pp = 1;
905 }
906 else
907 {
908 if (r != 1) pari_err(sqrter5);
909 z = sqrt_2adic(x, pp);
910 z = gerepileuptoint(av, z);
911 pp--;
912 }
913 mod = int2n(pp);
914 }
915 else /* p != 2 */
916 {
917 z = sqrt_padic(x, mod, pp, p);
918 z = gerepileuptoint(av, z);
919 mod = icopy(mod);
920 }
921 y[1] = evalprecp(pp) | evalvalp(e);
922 gel(y,2) = icopy(p);
923 gel(y,3) = mod;
924 gel(y,4) = z; return y;
925 }
926
927 GEN
gsqrt(GEN x,long prec)928 gsqrt(GEN x, long prec)
929 {
930 pari_sp av;
931 GEN y, p1, p2;
932
933 switch(typ(x))
934 {
935 case t_REAL: return sqrtr(x);
936
937 case t_INTMOD:
938 y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
939 p1 = Fp_sqrt(gel(x,2),gel(y,1));
940 if (!p1) pari_err(sqrter5);
941 gel(y,2) = p1; return y;
942
943 case t_COMPLEX:
944 if (isexactzero(gel(x,2))) return gsqrt(gel(x,1), prec);
945 y = cgetg(3,t_COMPLEX); av = avma;
946
947 p1 = gsqr(gel(x,1));
948 p2 = gsqr(gel(x,2)); p1 = gsqrt(gadd(p1,p2), prec);
949 if (gcmp0(p1)) { gel(y,1) = gel(y,2) = sqrtr(p1); return y; }
950 if (gsigne(gel(x,1)) < 0)
951 {
952 p1 = sqrtr( gmul2n(gsub(p1,gel(x,1)), -1) );
953 if (gsigne(gel(x,2)) < 0) setsigne(p1, -signe(p1));
954 gel(y,2) = gerepileuptoleaf(av, p1); av = avma;
955 gel(y,1) = gerepileuptoleaf(av, gdiv(gel(x,2), gmul2n(p1,1)));
956 } else {
957 p1 = sqrtr( gmul2n(gadd(p1,gel(x,1)), -1) );
958 gel(y,1) = gerepileuptoleaf(av, p1); av = avma;
959 gel(y,2) = gerepileuptoleaf(av, gdiv(gel(x,2), gmul2n(p1,1)));
960 }
961 return y;
962
963 case t_PADIC:
964 return padic_sqrt(x);
965
966 default:
967 av = avma; if (!(y = toser_i(x))) break;
968 return gerepileupto(av, ser_powfrac(y, ghalf, prec));
969 }
970 return transc(gsqrt,x,prec);
971 }
972 /********************************************************************/
973 /** **/
974 /** N-th ROOT **/
975 /** **/
976 /********************************************************************/
977 /* exp(2Ipi/n), assume n positive t_INT */
978 GEN
rootsof1complex(GEN n,long prec)979 rootsof1complex(GEN n, long prec)
980 {
981 pari_sp av = avma;
982 if (is_pm1(n)) return real_1(prec);
983 if (lgefint(n)==3 && n[2]==2) return stor(-1, prec);
984 return gerepileupto(av, exp_Ir( divri(Pi2n(1, prec), n) ));
985 }
986
987 /*Only the O() of y is used*/
988 GEN
rootsof1padic(GEN n,GEN y)989 rootsof1padic(GEN n, GEN y)
990 {
991 pari_sp av0 = avma, av;
992 GEN z, r = cgetp(y);
993
994 av = avma; (void)Fp_sqrtn(gen_1,n,gel(y,2),&z);
995 if (z==gen_0) { avma = av0; return gen_0; }/*should not happen*/
996 z = padicsqrtnlift(gen_1, n, z, gel(y,2), precp(y));
997 affii(z, gel(r,4)); avma = av; return r;
998 }
999
1000 static GEN exp_p(GEN x);
1001 /*compute the p^e th root of x p-adic, assume x != 0 */
1002 GEN
padic_sqrtn_ram(GEN x,long e)1003 padic_sqrtn_ram(GEN x, long e)
1004 {
1005 pari_sp ltop=avma;
1006 GEN a, p = gel(x,2), n = powiu(p,e);
1007 long v = valp(x);
1008 if (v)
1009 {
1010 long z;
1011 v = sdivsi_rem(v, n, &z);
1012 if (z) return NULL;
1013 x = gcopy(x); setvalp(x,0);
1014 }
1015 /*If p=2 -1 is an root of unity in U1,we need an extra check*/
1016 if (lgefint(p)==3 && p[2]==2 && mod8(gel(x,4))!=signe(gel(x,4)))
1017 return NULL;
1018 a = exp_p(gdiv(palog(x), n));
1019 if (!a) return NULL;
1020 /*Here n=p^e and a^n=z*x where z is a (p-1)th-root of unity. Note that
1021 z^p=z; hence for b = a/z, then b^n=x. We say b=x/a^(n-1)*/
1022 a = gdiv(x, powgi(a,addis(n,-1))); if (v) setvalp(a,v);
1023 return gerepileupto(ltop,a);
1024 }
1025
1026 /*compute the nth root of x p-adic p prime with n*/
1027 GEN
padic_sqrtn_unram(GEN x,GEN n,GEN * zetan)1028 padic_sqrtn_unram(GEN x, GEN n, GEN *zetan)
1029 {
1030 pari_sp av;
1031 GEN Z, a, r, p = gel(x,2);
1032 long v = valp(x);
1033 if (v)
1034 {
1035 long z;
1036 v = sdivsi_rem(v,n,&z);
1037 if (z) return NULL;
1038 }
1039 r = cgetp(x); setvalp(r,v);
1040 Z = NULL; /* -Wall */
1041 if (zetan) Z = cgetp(x);
1042 av = avma; a = Fp_sqrtn(gel(x,4), n, p, zetan);
1043 if (!a) return NULL;
1044 affii(padicsqrtnlift(gel(x,4), n, a, p, precp(x)), gel(r,4));
1045 if (zetan)
1046 {
1047 affii(padicsqrtnlift(gen_1, n, *zetan, p, precp(x)), gel(Z,4));
1048 *zetan = Z;
1049 }
1050 avma = av; return r;
1051 }
1052
1053 GEN
padic_sqrtn(GEN x,GEN n,GEN * zetan)1054 padic_sqrtn(GEN x, GEN n, GEN *zetan)
1055 {
1056 pari_sp av = avma, tetpil;
1057 GEN q, p = gel(x,2);
1058 long e;
1059 if (!signe(x[4]))
1060 {
1061 long m = itos(n);
1062 if (zetan) *zetan = gen_1;
1063 return zeropadic(p, (valp(x)+m-1)/m);
1064 }
1065 /* treat the ramified part using logarithms */
1066 e = Z_pvalrem(n, p, &q);
1067 if (e) { x = padic_sqrtn_ram(x,e); if (!x) return NULL; }
1068 if (is_pm1(q))
1069 { /* finished */
1070 if (signe(q) < 0) x = ginv(x);
1071 x = gerepileupto(av, x);
1072 if (zetan)
1073 *zetan = (e && lgefint(p)==3 && p[2]==2)? gen_m1 /*-1 in Q_2*/
1074 : gen_1;
1075 return x;
1076 }
1077 tetpil = avma;
1078 /* use hensel lift for unramified case */
1079 x = padic_sqrtn_unram(x, q, zetan);
1080 if (!x) return NULL;
1081 if (zetan)
1082 {
1083 GEN *gptr[2];
1084 if (e && lgefint(p)==3 && p[2]==2)/*-1 in Q_2*/
1085 {
1086 tetpil = avma; x = gcopy(x); *zetan = gneg(*zetan);
1087 }
1088 gptr[0] = &x; gptr[1] = zetan;
1089 gerepilemanysp(av,tetpil,gptr,2);
1090 return x;
1091 }
1092 return gerepile(av,tetpil,x);
1093 }
1094
1095 /* x^(1/n) */
1096 GEN
sqrtnr(GEN x,long n)1097 sqrtnr(GEN x, long n)
1098 {
1099 if (typ(x) != t_REAL) pari_err(typeer,"sqrtnr");
1100 return mpexp(divrs(mplog(x), n));
1101 }
1102
1103 GEN
gsqrtn(GEN x,GEN n,GEN * zetan,long prec)1104 gsqrtn(GEN x, GEN n, GEN *zetan, long prec)
1105 {
1106 long i, lx, tx;
1107 pari_sp av;
1108 GEN y, z;
1109 if (typ(n)!=t_INT) pari_err(talker,"second arg must be integer in gsqrtn");
1110 if (!signe(n)) pari_err(talker,"1/0 exponent in gsqrtn");
1111 if (is_pm1(n))
1112 {
1113 if (zetan) *zetan = gen_1;
1114 return (signe(n) > 0)? gcopy(x): ginv(x);
1115 }
1116 if (zetan) *zetan = gen_0;
1117 tx = typ(x);
1118 if (is_matvec_t(tx))
1119 {
1120 lx = lg(x); y = cgetg(lx,tx);
1121 for (i=1; i<lx; i++) gel(y,i) = gsqrtn(gel(x,i),n,NULL,prec);
1122 return y;
1123 }
1124 av = avma;
1125 switch(tx)
1126 {
1127 case t_INTMOD:
1128 z = gen_0;
1129 y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
1130 if (zetan) { z = cgetg(3,t_INTMOD); gel(z,1) = gel(y,1); }
1131 gel(y,2) = Fp_sqrtn(gel(x,2),n,gel(x,1),zetan);
1132 if (!y[2]) {
1133 if (zetan) return gen_0;
1134 pari_err(talker,"nth-root does not exist in gsqrtn");
1135 }
1136 if (zetan) { gel(z,2) = *zetan; *zetan = z; }
1137 return y;
1138
1139 case t_PADIC:
1140 y = padic_sqrtn(x,n,zetan);
1141 if (!y) {
1142 if (zetan) return gen_0;
1143 pari_err(talker,"nth-root does not exist in gsqrtn");
1144 }
1145 return y;
1146
1147 case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
1148 i = precision(x); if (i) prec = i;
1149 if (tx==t_INT && is_pm1(x) && signe(x) > 0)
1150 /*speed-up since there is no way to call rootsof1complex from gp*/
1151 y = real_1(prec);
1152 else if (gcmp0(x))
1153 {
1154 if (signe(n) < 0) pari_err(gdiver);
1155 if (isinexactreal(x))
1156 {
1157 long e = gexpo(x), junk;
1158 y = real_0_bit(e < 2? 0: sdivsi_rem(e, n, &junk));
1159 }
1160 else
1161 y = real_0(prec);
1162 }
1163 else
1164 y = gerepileupto(av, gexp(gdiv(glog(x,prec), n), prec));
1165 if (zetan) *zetan = rootsof1complex(n,prec);
1166 return y;
1167
1168 case t_QUAD:
1169 return gsqrtn(quadtoc(x, prec), n, zetan, prec);
1170
1171 default:
1172 av = avma; if (!(y = toser_i(x))) break;
1173 return gerepileupto(av, ser_powfrac(y, ginv(n), prec));
1174 }
1175 pari_err(typeer,"gsqrtn");
1176 return NULL;/* not reached */
1177 }
1178
1179 /********************************************************************/
1180 /** **/
1181 /** EXP(X) - 1 **/
1182 /** **/
1183 /********************************************************************/
1184 /* exp(|x|) - 1, assume x != 0 */
1185 GEN
exp1r_abs(GEN x)1186 exp1r_abs(GEN x)
1187 {
1188 long l = lg(x), l2 = l+1, ex = expo(x), l1, i, n, m, s;
1189 GEN y = cgetr(l), p1, p2, p3, X, unr;
1190 pari_sp av2, av = avma;
1191 double a, b, beta, gama = 2.0 /* optimized for SUN3 */;
1192 /* KB: 3.0 is better for UltraSparc */
1193 beta = 5. + bit_accuracy_mul(l, LOG2);
1194 a = sqrt(beta/(gama*LOG2));
1195 b = (BITS_IN_LONG-1-ex)
1196 + log2(a * (gama/2.718281828459045) / (double)(ulong)x[2]);
1197 if (a >= b)
1198 {
1199 n = (long)(1+a*gama);
1200 m = (long)(1+a-b);
1201 l2 += m>>TWOPOTBITS_IN_LONG;
1202 } else { /* rare ! */
1203 b = -1 - log((double)(ulong)x[2]) + (BITS_IN_LONG-1-ex)*LOG2; /*-1-log(x)*/
1204 n = (long)(1.1 + beta/b);
1205 m = 0;
1206 }
1207 unr=real_1(l2);
1208 p2 =real_1(l2); setlg(p2,3);
1209 X = cgetr(l2); affrr(x, X); setsigne(X, 1);
1210 if (m) setexpo(X, ex-m);
1211
1212 s = 0; l1 = 3; av2 = avma;
1213 for (i=n; i>=2; i--)
1214 { /* compute X^(n-1)/n! + ... + X/2 + 1 */
1215 setlg(X,l1); p3 = divrs(X,i);
1216 s -= expo(p3); p1 = mulrr(p3,p2); setlg(p1,l1);
1217 l1 += s>>TWOPOTBITS_IN_LONG; if (l1>l2) l1=l2;
1218 s &= (BITS_IN_LONG-1);
1219 setlg(unr,l1); p1 = addrr_sign(unr,1, p1,1);
1220 setlg(p2,l1); affrr(p1,p2); avma = av2; /* p2 <- 1 + (X/i)*p2 */
1221 }
1222 setlg(X,l2); p2 = mulrr(X,p2);
1223
1224 for (i=1; i<=m; i++)
1225 {
1226 setlg(p2,l2);
1227 p2 = mulrr(p2, addsr(2,p2));
1228 }
1229 affr_fixlg(p2,y); avma = av; return y;
1230 }
1231
1232 GEN
mpexp1(GEN x)1233 mpexp1(GEN x)
1234 {
1235 long sx = signe(x);
1236 GEN y, z;
1237 pari_sp av;
1238 if (!sx) return real_0_bit(expo(x));
1239 if (sx > 0) return exp1r_abs(x);
1240 /* compute exp(x) * (1 - exp(-x)) */
1241 av = avma; y = exp1r_abs(x);
1242 z = addsr(1, y); setsigne(z, -1);
1243 return gerepileupto(av, divrr(y, z));
1244 }
1245
1246 /********************************************************************/
1247 /** **/
1248 /** EXP(X) **/
1249 /** **/
1250 /********************************************************************/
1251
1252 static GEN
mpexp_basecase(GEN x)1253 mpexp_basecase(GEN x)
1254 {
1255 pari_sp av = avma;
1256 GEN y = addsr(1, exp1r_abs(x));
1257 if (signe(x) < 0) y = ginv(y);
1258 return gerepileupto(av,y);
1259 }
1260
1261 GEN
mpexp(GEN x)1262 mpexp(GEN x)
1263 {
1264 const long s = 6; /*Initial steps using basecase*/
1265 long i, n, mask, p, l, sx = signe(x), sh=0;
1266 GEN a, z;
1267
1268 if (!sx) {
1269 long e = expo(x);
1270 return e >= 0? real_0_bit(e): real_1(nbits2prec(-e));
1271 }
1272
1273 l = lg(x);
1274 if (l <= max(EXPNEWTON_LIMIT, 1<<s)) return mpexp_basecase(x);
1275 z = cgetr(l); /* room for result */
1276 if (expo(x) >= 0)
1277 { /* x>=1 : we do x %= log(2) to keep x small*/
1278 sh = (long) (rtodbl(x)/LOG2);
1279 x = subrr(rtor(x,l+1), mulsr(sh, mplog2(l+1)));
1280 if (!signe(x)) { avma = (pari_sp)(z+l); return real2n(sh, l); }
1281 }
1282 n = hensel_lift_accel(l-1,&mask);
1283 for(i=0, p=1; i<s; i++) { p <<= 1; if (mask&(1<<i)) p--; }
1284 a = mpexp_basecase(rtor(x, p+2));
1285 x = addrs(x,1);
1286 if (lg(x) < l+1) x = rtor(x, l+1);
1287 for(i=s; i<n; i++)
1288 {
1289 p <<= 1; if (mask&(1<<i)) p--;
1290 setlg(x, p+2); a = rtor(a, p+2);
1291 a = mulrr(a, subrr(x, logr_abs(a))); /* a := a (x - log(a)) */
1292 }
1293 affrr(a,z);
1294 if (sh) setexpo(z, expo(z) + sh);
1295 avma = (pari_sp)z; return z;
1296 }
1297
1298 static long
exp_p_prec(GEN x)1299 exp_p_prec(GEN x)
1300 {
1301 long k, e = valp(x), n = e + precp(x);
1302 GEN p = gel(x,2);
1303 int is2 = equaliu(p,2);
1304 if (e < 1 || (e == 1 && is2)) return -1;
1305 if (is2)
1306 {
1307 n--; e--; k = n/e;
1308 if (n%e == 0) k--;
1309 }
1310 else
1311 {
1312 GEN r, t = subis(p, 1);
1313 k = itos(dvmdii(subis(mulis(t,n), 1), subis(mulis(t,e), 1), &r));
1314 if (!signe(r)) k--;
1315 }
1316 return k;
1317 }
1318
1319 static GEN
exp_p(GEN x)1320 exp_p(GEN x)
1321 {
1322 long k;
1323 pari_sp av;
1324 GEN y;
1325
1326 if (gcmp0(x)) return gaddgs(x,1);
1327 k = exp_p_prec(x);
1328 if (k < 0) return NULL;
1329 av = avma;
1330 for (y=gen_1; k; k--) y = gaddsg(1, gdivgs(gmul(y,x), k));
1331 return gerepileupto(av, y);
1332 }
1333 static GEN
cos_p(GEN x)1334 cos_p(GEN x)
1335 {
1336 long k;
1337 pari_sp av;
1338 GEN x2, y;
1339
1340 if (gcmp0(x)) return gaddgs(x,1);
1341 k = exp_p_prec(x);
1342 if (k < 0) return NULL;
1343 av = avma; x2 = gsqr(x);
1344 if (k & 1) k--;
1345 for (y=gen_1; k; k-=2)
1346 {
1347 GEN t = gdiv(gmul(y,x2), mulss(k, k-1));
1348 y = gsubsg(1, t);
1349 }
1350 return gerepileupto(av, y);
1351 }
1352 static GEN
sin_p(GEN x)1353 sin_p(GEN x)
1354 {
1355 long k;
1356 pari_sp av;
1357 GEN x2, y;
1358
1359 if (gcmp0(x)) return gaddgs(x,1);
1360 k = exp_p_prec(x);
1361 if (k < 0) return NULL;
1362 av = avma; x2 = gsqr(x);
1363 if (k & 1) k--;
1364 for (y=gen_1; k; k-=2)
1365 {
1366 GEN t = gdiv(gmul(y,x2), mulss(k, k+1));
1367 y = gsubsg(1, t);
1368 }
1369 return gerepileupto(av, gmul(y, x));
1370 }
1371
1372 static GEN
cxexp(GEN x,long prec)1373 cxexp(GEN x, long prec)
1374 {
1375 GEN r,p1,p2, y = cgetg(3,t_COMPLEX);
1376 pari_sp av = avma, tetpil;
1377 r = gexp(gel(x,1),prec);
1378 if (gcmp0(r)) { gel(y,1) = r; gel(y,2) = r; return y; }
1379 gsincos(gel(x,2),&p2,&p1,prec);
1380 tetpil = avma;
1381 gel(y,1) = gmul(r,p1);
1382 gel(y,2) = gmul(r,p2);
1383 gerepilecoeffssp(av,tetpil,y+1,2);
1384 return y;
1385 }
1386
1387 static GEN
serexp(GEN x,long prec)1388 serexp(GEN x, long prec)
1389 {
1390 pari_sp av;
1391 long i,j,lx,ly,ex,mi;
1392 GEN p1,y,xd,yd;
1393
1394 ex = valp(x);
1395 if (ex < 0) pari_err(negexper,"gexp");
1396 if (gcmp0(x)) return gaddsg(1,x);
1397 lx = lg(x);
1398 if (ex)
1399 {
1400 ly = lx+ex; y = cgetg(ly,t_SER);
1401 mi = lx-1; while (mi>=3 && isexactzero(gel(x,mi))) mi--;
1402 mi += ex-2;
1403 y[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(x));
1404 /* zd[i] = coefficient of X^i in z */
1405 xd = x+2-ex; yd = y+2; ly -= 2;
1406 gel(yd,0) = gen_1;
1407 for (i=1; i<ex; i++) gel(yd,i) = gen_0;
1408 for ( ; i<ly; i++)
1409 {
1410 av = avma; p1 = gen_0;
1411 for (j=ex; j<=min(i,mi); j++)
1412 p1 = gadd(p1, gmulgs(gmul(gel(xd,j),gel(yd,i-j)),j));
1413 gel(yd,i) = gerepileupto(av, gdivgs(p1,i));
1414 }
1415 return y;
1416 }
1417 av = avma; y = cgetg(lx, t_SER);
1418 y[1] = x[1]; gel(y,2) = gen_0;
1419 for (i=3; i <lx; i++) y[i] = x[i];
1420 p1 = gexp(gel(x,2),prec);
1421 y = gmul(p1, serexp(normalize(y),prec));
1422 return gerepileupto(av, y);
1423 }
1424
1425 GEN
gexp(GEN x,long prec)1426 gexp(GEN x, long prec)
1427 {
1428 switch(typ(x))
1429 {
1430 case t_REAL: return mpexp(x);
1431 case t_COMPLEX: return cxexp(x,prec);
1432 case t_PADIC: x = exp_p(x);
1433 if (!x) pari_err(talker,"p-adic argument out of range in gexp");
1434 return x;
1435 case t_INTMOD: pari_err(typeer,"gexp");
1436 default:
1437 {
1438 pari_sp av = avma;
1439 GEN y;
1440 if (!(y = toser_i(x))) break;
1441 return gerepileupto(av, serexp(y,prec));
1442 }
1443 }
1444 return transc(gexp,x,prec);
1445 }
1446
1447 /********************************************************************/
1448 /** **/
1449 /** AGM(X, Y) **/
1450 /** **/
1451 /********************************************************************/
1452 static int
agmr_gap(GEN a,GEN b,long L)1453 agmr_gap(GEN a, GEN b, long L)
1454 {
1455 GEN d = subrr(b, a);
1456 return (signe(d) && expo(d) - expo(b) >= L);
1457 }
1458 /* assume x > 0 */
1459 static GEN
agm1r_abs(GEN x)1460 agm1r_abs(GEN x)
1461 {
1462 long l = lg(x), L = 5-bit_accuracy(l);
1463 GEN a1, b1, y = cgetr(l);
1464 pari_sp av = avma;
1465
1466 a1 = addrr(real_1(l), x); setexpo(a1, expo(a1)-1);
1467 b1 = sqrtr_abs(x);
1468 while (agmr_gap(a1,b1,L))
1469 {
1470 GEN a = a1;
1471 a1 = addrr(a,b1); setexpo(a1, expo(a1)-1);
1472 b1 = sqrtr_abs(mulrr(a,b1));
1473 }
1474 affr_fixlg(a1,y); avma = av; return y;
1475 }
1476
1477 static int
agmcx_gap(GEN a,GEN b,long L)1478 agmcx_gap(GEN a, GEN b, long L)
1479 {
1480 GEN d = gsub(b, a);
1481 return (!gcmp0(d) && gexpo(d) - gexpo(b) >= L);
1482 }
1483 static GEN
agm1cx(GEN x,long prec)1484 agm1cx(GEN x, long prec)
1485 {
1486 GEN a1, b1;
1487 pari_sp av = avma, av2;
1488 long L, l = precision(x); if (!l) l = prec;
1489
1490 L = 5-bit_accuracy(l);
1491 a1 = gtofp(gmul2n(gadd(real_1(l), x), -1), l); /* avoid loss of accuracy */
1492 av2 = avma;
1493 b1 = gsqrt(x, prec);
1494 while (agmcx_gap(a1,b1,L))
1495 {
1496 GEN a = a1;
1497 a1 = gmul2n(gadd(a,b1),-1);
1498 av2 = avma;
1499 b1 = gsqrt(gmul(a,b1), prec);
1500 }
1501 avma = av2; return gerepileupto(av,a1);
1502 }
1503
1504 /* agm(1,x) */
1505 static GEN
agm1(GEN x,long prec)1506 agm1(GEN x, long prec)
1507 {
1508 GEN p1, a, a1, b1, y;
1509 long l, l2, ep;
1510 pari_sp av;
1511
1512 if (gcmp0(x)) return gcopy(x);
1513 switch(typ(x))
1514 {
1515 case t_INT:
1516 if (!is_pm1(x)) break;
1517 return (signe(x) > 0)? real_1(prec): real_0(prec);
1518
1519 case t_REAL: return signe(x) > 0? agm1r_abs(x): agm1cx(x, prec);
1520
1521 case t_COMPLEX:
1522 if (gcmp0(gel(x,2)) && gsigne(gel(x,1)) > 0)
1523 return agm1(gel(x,1), prec);
1524 return agm1cx(x, prec);
1525
1526 case t_PADIC:
1527 av = avma;
1528 a1 = x; b1 = gen_1; l = precp(x);
1529 do
1530 {
1531 a = a1;
1532 a1 = gmul2n(gadd(a,b1),-1);
1533 b1 = padic_sqrt(gmul(a,b1));
1534 p1 = gsub(b1,a1); ep = valp(p1)-valp(b1);
1535 if (ep<=0) { b1 = gneg_i(b1); p1 = gsub(b1,a1); ep=valp(p1)-valp(b1); }
1536 }
1537 while (ep<l && !gcmp0(p1));
1538 return gerepilecopy(av,a1);
1539
1540 default:
1541 av = avma; if (!(y = toser_i(x))) break;
1542 a1 = y; b1 = gen_1; l = lg(y)-2;
1543 l2 = 5-bit_accuracy(prec);
1544 do
1545 {
1546 a = a1;
1547 a1 = gmul2n(gadd(a,b1),-1);
1548 b1 = ser_powfrac(gmul(a,b1), ghalf, prec);
1549 p1 = gsub(b1,a1); ep = valp(p1)-valp(b1);
1550 }
1551 while (ep<l && !gcmp0(p1)
1552 && (!isinexactreal(p1) || gexpo(p1) - gexpo(b1) >= l2));
1553 return gerepilecopy(av,a1);
1554 }
1555 return transc(agm1,x,prec);
1556 }
1557
1558 GEN
agm(GEN x,GEN y,long prec)1559 agm(GEN x, GEN y, long prec)
1560 {
1561 long ty = typ(y);
1562 pari_sp av;
1563
1564 if (is_matvec_t(ty))
1565 {
1566 ty = typ(x);
1567 if (is_matvec_t(ty)) pari_err(talker,"agm of two vector/matrices");
1568 swap(x, y);
1569 }
1570 if (gcmp0(y)) return gcopy(y);
1571 av = avma;
1572 return gerepileupto(av, gmul(y, agm1(gdiv(x,y), prec)));
1573 }
1574
1575 /********************************************************************/
1576 /** **/
1577 /** LOG(X) **/
1578 /** **/
1579 /********************************************************************/
1580 /* cf logagmr_abs(). Compute Pi/2agm(1, 4/2^n) ~ log(2^n) = n log(2) */
1581 GEN
constlog2(long prec)1582 constlog2(long prec)
1583 {
1584 pari_sp av;
1585 long l, n;
1586 GEN y, tmplog2;
1587
1588 if (glog2 && lg(glog2) >= prec) return glog2;
1589
1590 tmplog2 = newbloc(prec);
1591 *tmplog2 = evaltyp(t_REAL) | evallg(prec);
1592 av = avma;
1593 l = prec+1;
1594 n = bit_accuracy(l) >> 1;
1595 y = divrr(Pi2n(-1, l), agm1r_abs( real2n(2 - n, l) ));
1596 affrr(divrs(y,n), tmplog2);
1597 if (glog2) gunclone(glog2);
1598 glog2 = tmplog2; avma = av; return glog2;
1599 }
1600
1601 GEN
mplog2(long prec)1602 mplog2(long prec)
1603 {
1604 GEN x = cgetr(prec);
1605 affrr(constlog2(prec), x); return x;
1606 }
1607
1608 /*return log(|x|), assuming x != 0 */
1609 GEN
logr_abs(GEN X)1610 logr_abs(GEN X)
1611 {
1612 pari_sp ltop, av;
1613 long EX, l1, l2, m, n, k, e, s, l = lg(X);
1614 double a, b;
1615 GEN z, x, y, y2, S, unr;
1616 ulong u, v;
1617
1618 if (l > LOGAGM_LIMIT) return logagmr_abs(X);
1619 EX = expo(X);
1620 if (absrnz_egal2n(X)) return EX? mulsr(EX, mplog2(l)): real_0(l);
1621
1622 av = avma; z = cgetr(l); ltop = avma;
1623 l2 = l+1; x = cgetr(l2); affrr(X,x);
1624 x[1] = evalsigne(1) | evalexpo(0);
1625 /* X = x 2^EX, 1 < x < 2 */
1626 av = avma; l -= 2;
1627 k = 2;
1628 u = ((ulong)x[k]) & (~HIGHBIT); /* x[2] - HIGHBIT, assuming HIGHBIT set */
1629 v = BITS_IN_LONG-1;
1630 while (!u) { v += BITS_IN_LONG; u = (ulong)x[++k]; } /* terminates: x>1 */
1631 a = (double)v - log2((double)u); /* ~ -log2(x - 1) */
1632 b = sqrt((BITS_IN_HALFULONG/3.0) * l);
1633 if (a <= b)
1634 {
1635 n = 1 + (long)(3*b);
1636 m = 1 + (long)(b-a);
1637 if ((ulong)m >= BITS_IN_LONG) { GEN t;
1638 l2 += m>>TWOPOTBITS_IN_LONG;
1639 t = cgetr(l2); affrr(x,t); x = t;
1640 }
1641 for (k=1; k<=m; k++) x = sqrtr_abs(x);
1642 }
1643 else
1644 {
1645 n = 1 + (long)(BITS_IN_HALFULONG*l / a);
1646 m = 0;
1647 }
1648 y = divrr(subrex01(x), addrex01(x)); /* = (x-1) / (x+1) ~ 0 */
1649 y2 = gsqr(y);
1650 /* log(x) = log(1+y) - log(1-y) = 2 \sum_{k odd} y^k / k */
1651 k = 2*n + 1;
1652 unr = real_1(l2); S = x; av = avma;
1653 setlg(S, 3);
1654 setlg(unr,3); affrr(divrs(unr,k), S); /* destroy x, not needed anymore */
1655
1656 s = 0; e = expo(y2); l1 = 3;
1657 for (k -= 2; k > 0; k -= 2) /* k = 2n+1, ..., 1 */
1658 {
1659 GEN T; /* S = y^(2n+1-k)/(2n+1) + ... + 1 / k */
1660 setlg(y2, l1); T = mulrr(S,y2);
1661 setlg(unr,l1);
1662 s -= e; /* >= 0 */
1663 l1 += s>>TWOPOTBITS_IN_LONG; if (l1>l2) l1=l2;
1664 s &= (BITS_IN_LONG-1);
1665 setlg(S, l1);
1666 affrr(addrr(divrs(unr, k), T), S); avma = av;
1667 }
1668 setlg(S, l2); y = mulrr(y,S); /* = log(X)/2 */
1669 setexpo(y, expo(y)+m+1);
1670 if (EX) y = addrr(y, mulsr(EX, mplog2(l2)));
1671 affr_fixlg(y, z); avma = ltop; return z;
1672 }
1673
1674 GEN
logagmr_abs(GEN q)1675 logagmr_abs(GEN q)
1676 {
1677 long prec = lg(q), lim, e = expo(q);
1678 GEN z, y, Q;
1679 pari_sp av;
1680
1681 if (absrnz_egal2n(q)) return e? mulsr(e, mplog2(prec)): real_0(prec);
1682 z = cgetr(prec); av = avma; prec++;
1683 lim = bit_accuracy(prec) >> 1;
1684 Q = cgetr(prec); affrr(q, Q);
1685 Q[1] = evalsigne(1) | evalexpo(lim);
1686
1687 /* Pi / 2agm(1, 4/Q) ~ log(Q), q = Q * 2^(e-lim) */
1688 y = divrr(Pi2n(-1, prec), agm1r_abs( divsr(4, Q) ));
1689 y = addrr(y, mulsr(e - lim, mplog2(prec)));
1690 affr_fixlg(y, z); avma = av; return z;
1691 }
1692
1693 /* assume Im(q) != 0 */
1694 GEN
logagmcx(GEN q,long prec)1695 logagmcx(GEN q, long prec)
1696 {
1697 GEN z, y, Q, a, b;
1698 long lim, e, ea, eb;
1699 pari_sp av;
1700 int neg = 0;
1701
1702 z = cgetc(prec); av = avma; prec++;
1703 if (gsigne(gel(q,1)) < 0) { q = gneg(q); neg = 1; }
1704 lim = bit_accuracy(prec) >> 1;
1705 Q = gtofp(q, prec);
1706 a = gel(Q,1);
1707 b = gel(Q,2);
1708 if (gcmp0(a)) {
1709 affr_fixlg(logr_abs(b), gel(z,1));
1710 y = Pi2n(-1, prec);
1711 if (signe(b) < 0) setsigne(y, -1);
1712 affr_fixlg(y, gel(z,2)); avma = av; return z;
1713 }
1714 ea = expo(a);
1715 eb = expo(b);
1716 if (ea <= eb)
1717 {
1718 setexpo(Q[1], lim - eb + ea);
1719 setexpo(Q[2], lim);
1720 e = lim - eb;
1721 } else {
1722 setexpo(Q[1], lim);
1723 setexpo(Q[2], lim - ea + eb);
1724 e = lim - ea;
1725 }
1726
1727 /* Pi / 2agm(1, 4/Q) ~ log(Q), q = Q * 2^e */
1728 y = gdiv(Pi2n(-1, prec), agm1cx( gdivsg(4, Q), prec ));
1729 a = gel(y,1);
1730 b = gel(y,2);
1731 a = addrr(a, mulsr(-e, mplog2(prec)));
1732 if (neg) b = gsigne(b) <= 0? gadd(b, mppi(prec))
1733 : gsub(b, mppi(prec));
1734 affr_fixlg(a, gel(z,1));
1735 affr_fixlg(b, gel(z,2)); avma = av; return z;
1736 }
1737
1738 GEN
mplog(GEN x)1739 mplog(GEN x)
1740 {
1741 if (signe(x)<=0) pari_err(talker,"non positive argument in mplog");
1742 return logr_abs(x);
1743 }
1744
1745 GEN
teich(GEN x)1746 teich(GEN x)
1747 {
1748 GEN p,q,y,z,aux,p1;
1749 long n, k;
1750 pari_sp av;
1751
1752 if (typ(x)!=t_PADIC) pari_err(talker,"not a p-adic argument in teichmuller");
1753 if (!signe(x[4])) return gcopy(x);
1754 p = gel(x,2);
1755 q = gel(x,3);
1756 z = gel(x,4); y = cgetp(x); av = avma;
1757 if (equaliu(p,2))
1758 z = (mod4(z) & 2)? addsi(-1,q): gen_1;
1759 else
1760 {
1761 p1 = addsi(-1, p);
1762 z = remii(z, p);
1763 aux = diviiexact(addsi(-1,q),p1); n = precp(x);
1764 for (k=1; k<n; k<<=1)
1765 z = modii(mulii(z,addsi(1,mulii(aux,addsi(-1,Fp_pow(z,p1,q))))), q);
1766 }
1767 affii(z, gel(y,4)); avma = av; return y;
1768 }
1769
1770 /* Let x = 1 mod p and y := (x-1)/(x+1) = 0 (p). Then
1771 * log(x) = log(1+y) - log(1-y) = 2 \sum_{k odd} y^k / k.
1772 * palogaux(x) returns the last sum (not multiplied by 2) */
1773 static GEN
palogaux(GEN x)1774 palogaux(GEN x)
1775 {
1776 long k,e,pp;
1777 GEN y,s,y2, p = gel(x,2);
1778
1779 if (equalii(gen_1, gel(x,4)))
1780 {
1781 long v = valp(x)+precp(x);
1782 if (equaliu(p,2)) v--;
1783 return zeropadic(p, v);
1784 }
1785 y = gdiv(gaddgs(x,-1), gaddgs(x,1));
1786 e = valp(y); /* > 0 */
1787 pp = e+precp(y);
1788 if (equaliu(p,2)) pp--;
1789 else
1790 {
1791 GEN p1;
1792 for (p1=utoipos(e); cmpui(pp,p1) > 0; pp++) p1 = mulii(p1, p);
1793 pp -= 2;
1794 }
1795 k = pp/e; if (!odd(k)) k--;
1796 y2 = gsqr(y); s = gdivgs(gen_1,k);
1797 while (k > 2)
1798 {
1799 k -= 2; s = gadd(gmul(y2,s), gdivgs(gen_1,k));
1800 }
1801 return gmul(s,y);
1802 }
1803
1804 GEN
palog(GEN x)1805 palog(GEN x)
1806 {
1807 pari_sp av = avma;
1808 GEN y, p = gel(x,2);
1809
1810 if (!signe(x[4])) pari_err(talker,"zero argument in palog");
1811 if (equaliu(p,2))
1812 {
1813 y = gsqr(x); setvalp(y,0);
1814 y = palogaux(y);
1815 }
1816 else
1817 { /* compute log(x^(p-1)) / (p-1) */
1818 GEN mod = gel(x,3), p1 = subis(p,1);
1819 y = cgetp(x);
1820 gel(y,4) = Fp_pow(gel(x,4), p1, mod);
1821 p1 = diviiexact(subis(mod,1), p1); /* 1/(1-p) */
1822 y = gmul(palogaux(y), mulis(p1, -2));
1823 }
1824 return gerepileupto(av,y);
1825 }
1826
1827 GEN
log0(GEN x,long flag,long prec)1828 log0(GEN x, long flag,long prec)
1829 {
1830 switch(flag)
1831 {
1832 case 0:
1833 case 1: return glog(x,prec);
1834 default: pari_err(flagerr,"log");
1835 }
1836 return NULL; /* not reached */
1837 }
1838
1839 GEN
glog(GEN x,long prec)1840 glog(GEN x, long prec)
1841 {
1842 pari_sp av, tetpil;
1843 GEN y, p1;
1844
1845 switch(typ(x))
1846 {
1847 case t_REAL:
1848 if (signe(x) >= 0)
1849 {
1850 if (!signe(x)) pari_err(talker,"zero argument in mplog");
1851 return logr_abs(x);
1852 }
1853 y = cgetg(3,t_COMPLEX);
1854 gel(y,1) = logr_abs(x);
1855 gel(y,2) = mppi(lg(x)); return y;
1856
1857 case t_COMPLEX:
1858 if (gcmp0(gel(x,2))) return glog(gel(x,1), prec);
1859 if (prec > LOGAGMCX_LIMIT) return logagmcx(x, prec);
1860 y = cgetg(3,t_COMPLEX);
1861 gel(y,2) = garg(x,prec);
1862 av = avma; p1 = glog(cxnorm(x),prec); tetpil = avma;
1863 gel(y,1) = gerepile(av,tetpil,gmul2n(p1,-1)); return y;
1864
1865 case t_PADIC:
1866 return palog(x);
1867
1868 case t_INTMOD: pari_err(typeer,"glog");
1869 default:
1870 av = avma; if (!(y = toser_i(x))) break;
1871 if (valp(y) || gcmp0(y)) pari_err(talker,"log is not meromorphic at 0");
1872 p1 = integ(gdiv(derivser(y), y), varn(y)); /* log(y)' = y'/y */
1873 if (!gcmp1(gel(y,2))) p1 = gadd(p1, glog(gel(y,2),prec));
1874 return gerepileupto(av, p1);
1875 }
1876 return transc(glog,x,prec);
1877 }
1878 /********************************************************************/
1879 /** **/
1880 /** SINE, COSINE **/
1881 /** **/
1882 /********************************************************************/
1883
1884 /* Reduce x0 mod Pi/2 to x in [-Pi/4, Pi/4]. Return cos(x)-1 */
1885 static GEN
mpsc1(GEN x,long * ptmod8)1886 mpsc1(GEN x, long *ptmod8)
1887 {
1888 long e = expo(x), l = lg(x), l1, l2, i, n, m, s;
1889 pari_sp av;
1890 double beta, a, b, d;
1891 GEN y, unr, p2, p1, x2;
1892
1893 n = 0;
1894 if (e >= 0)
1895 {
1896 GEN q, z, pitemp = mppi(DEFAULTPREC + (e >> TWOPOTBITS_IN_LONG));
1897 setexpo(pitemp,-1);
1898 z = addrr(x,pitemp); /* = x + Pi/4 */
1899 if (expo(z) >= bit_accuracy(min(l, lg(z))) + 3) pari_err(precer,"mpsc1");
1900 setexpo(pitemp, 0);
1901 q = floorr( divrr(z,pitemp) ); /* round ( x / (Pi/2) ) */
1902 if (signe(q))
1903 {
1904 x = subrr(x, mulir(q, Pi2n(-1, l+1))); /* x mod Pi/2 */
1905 e = expo(x);
1906 n = mod4(q); if (n && signe(q) < 0) n = 4 - n;
1907 }
1908 }
1909 s = signe(x); *ptmod8 = (s < 0)? 4 + n: n;
1910 if (!s) return real_0_bit(-bit_accuracy(l)<<1);
1911
1912 l = lg(x); l2 = l+1; y = cgetr(l);
1913 beta = 5. + bit_accuracy_mul(l2,LOG2);
1914 a = sqrt(beta / LOG2);
1915 d = a + 1/LOG2 - log2(a / (double)(ulong)x[2]) - (BITS_IN_LONG-1-e);
1916 if (d >= 0)
1917 {
1918 n = (long)((1+a) / 2.0);
1919 m = (long)(1+d);
1920 l2 += m>>TWOPOTBITS_IN_LONG;
1921 } else { /* rare ! */
1922 b = -1 - log((double)(ulong)x[2]) + (BITS_IN_LONG-1-e)*LOG2; /*-1-log(x)*/
1923 n = (long)(1 + beta/(2.0*b));
1924 m = 0;
1925 }
1926 unr= real_1(l2);
1927 p2 = real_1(l2);
1928 x2 = cgetr(l2); av = avma;
1929 affrr(gsqr(x), x2);
1930 if (m) setexpo(x2, expo(x2) - (m<<1));
1931
1932 setlg(x2, 3); p1 = divrs(x2, 2*n+1);
1933 s = -expo(p1);
1934 l1 = 3 + (s>>TWOPOTBITS_IN_LONG);
1935 setlg(p2,l1);
1936 s = 0;
1937 for (i=n; i>=2; i--)
1938 {
1939 setlg(x2,l1); p1 = divrsns(x2, 2*i-1);
1940 s -= expo(p1); p1 = mulrr(p1,p2);
1941 l1 += s>>TWOPOTBITS_IN_LONG; if (l1>l2) l1=l2;
1942 s &= (BITS_IN_LONG-1);
1943 setlg(unr,l1); p1 = addrr_sign(unr,1, p1,-signe(p1));
1944 setlg(p2,l1); affrr(p1,p2); avma = av;
1945 }
1946 p2[1] = evalsigne(-signe(p2)) | evalexpo(expo(p2)-1); /* p2 := -p2/2 */
1947 setlg(p2,l2);
1948 setlg(x2,l2); p2 = mulrr(x2,p2);
1949 /* Now p2 = sum {1<= i <=n} (-1)^i x^(2i) / (2i)! ~ cos(x) - 1 */
1950 for (i=1; i<=m; i++)
1951 { /* p2 = cos(x)-1 --> cos(2x)-1 */
1952 setlg(p2,l2);
1953 p2 = mulrr(p2, addsr(2,p2));
1954 setexpo(p2, expo(p2)+1);
1955 }
1956 affr_fixlg(p2,y); return y;
1957 }
1958
1959 /* sqrt (|1 - (1+x)^2|) = sqrt(|x*(x+2)|). Sends cos(x)-1 to |sin(x)| */
1960 static GEN
mpaut(GEN x)1961 mpaut(GEN x)
1962 {
1963 pari_sp av = avma;
1964 GEN t = mulrr(x, addsr(2,x)); /* != 0 */
1965 if (!signe(t)) return real_0_bit(expo(t) >> 1);
1966 return gerepileuptoleaf(av, sqrtr_abs(t));
1967 }
1968
1969 /********************************************************************/
1970 /** COSINE **/
1971 /********************************************************************/
1972
1973 GEN
mpcos(GEN x)1974 mpcos(GEN x)
1975 {
1976 long mod8;
1977 pari_sp av;
1978 GEN y,p1;
1979
1980 if (!signe(x)) return real_1(3 + ((-expo(x)) >> TWOPOTBITS_IN_LONG));
1981
1982 av = avma; p1 = mpsc1(x,&mod8);
1983 switch(mod8)
1984 {
1985 case 0: case 4: y = addsr(1,p1); break;
1986 case 1: case 7: y = mpaut(p1); setsigne(y,-signe(y)); break;
1987 case 2: case 6: y = subsr(-1,p1); break;
1988 default: y = mpaut(p1); break; /* case 3: case 5: */
1989 }
1990 return gerepileuptoleaf(av, y);
1991 }
1992
1993 /* convert INT or FRAC to REAL, which is later reduced mod 2Pi : avoid
1994 * cancellation */
1995 static GEN
tofp_safe(GEN x,long prec)1996 tofp_safe(GEN x, long prec)
1997 {
1998 return (typ(x) == t_INT || gexpo(x) > 0)? gadd(x, real_0(prec))
1999 : fractor(x, prec);
2000 }
2001
2002 GEN
gcos(GEN x,long prec)2003 gcos(GEN x, long prec)
2004 {
2005 pari_sp av;
2006 GEN r, u, v, y, u1, v1;
2007 long i;
2008
2009 switch(typ(x))
2010 {
2011 case t_REAL:
2012 return mpcos(x);
2013
2014 case t_COMPLEX:
2015 i = precision(x); if (!i) i = prec;
2016 y = cgetc(i); av = avma;
2017 r = gexp(gel(x,2),prec);
2018 v1 = gmul2n(addrr(ginv(r),r), -1); /* = cos(I*Im(x)) */
2019 u1 = subrr(v1, r); /* = - I*sin(I*Im(x)) */
2020 gsincos(gel(x,1),&u,&v,prec);
2021 affr_fixlg(gmul(v1,v), gel(y,1));
2022 affr_fixlg(gmul(u1,u), gel(y,2)); return y;
2023
2024 case t_INT: case t_FRAC:
2025 y = cgetr(prec); av = avma;
2026 affr_fixlg(mpcos(tofp_safe(x,prec)), y); avma = av; return y;
2027
2028 case t_INTMOD: pari_err(typeer,"gcos");
2029
2030 case t_PADIC: x = cos_p(x);
2031 if (!x) pari_err(talker,"p-adic argument out of range in gcos");
2032 return x;
2033
2034 default:
2035 av = avma; if (!(y = toser_i(x))) break;
2036 if (gcmp0(y)) return gaddsg(1,y);
2037 if (valp(y) < 0) pari_err(negexper,"gcos");
2038 gsincos(y,&u,&v,prec);
2039 return gerepilecopy(av,v);
2040 }
2041 return transc(gcos,x,prec);
2042 }
2043 /********************************************************************/
2044 /** SINE **/
2045 /********************************************************************/
2046
2047 GEN
mpsin(GEN x)2048 mpsin(GEN x)
2049 {
2050 long mod8;
2051 pari_sp av;
2052 GEN y,p1;
2053
2054 if (!signe(x)) return real_0_bit(expo(x));
2055
2056 av = avma; p1 = mpsc1(x,&mod8);
2057 switch(mod8)
2058 {
2059 case 0: case 6: y=mpaut(p1); break;
2060 case 1: case 5: y=addsr(1,p1); break;
2061 case 2: case 4: y=mpaut(p1); setsigne(y,-signe(y)); break;
2062 default: y=subsr(-1,p1); break; /* case 3: case 7: */
2063 }
2064 return gerepileuptoleaf(av, y);
2065 }
2066
2067 GEN
gsin(GEN x,long prec)2068 gsin(GEN x, long prec)
2069 {
2070 pari_sp av;
2071 GEN r, u, v, y, v1, u1;
2072 long i;
2073
2074 switch(typ(x))
2075 {
2076 case t_REAL:
2077 return mpsin(x);
2078
2079 case t_COMPLEX:
2080 i = precision(x); if (!i) i = prec;
2081 y = cgetc(i); av = avma;
2082 r = gexp(gel(x,2),prec);
2083 v1 = gmul2n(addrr(ginv(r),r), -1); /* = cos(I*Im(x)) */
2084 u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
2085 gsincos(gel(x,1),&u,&v,prec);
2086 affr_fixlg(gmul(v1,u), gel(y,1));
2087 affr_fixlg(gmul(u1,v), gel(y,2)); return y;
2088
2089 case t_INT: case t_FRAC:
2090 y = cgetr(prec); av = avma;
2091 affr_fixlg(mpsin(tofp_safe(x,prec)), y); avma = av; return y;
2092
2093 case t_INTMOD: pari_err(typeer,"gsin");
2094
2095 case t_PADIC: x = sin_p(x);
2096 if (!x) pari_err(talker,"p-adic argument out of range in gsin");
2097 return x;
2098
2099 default:
2100 av = avma; if (!(y = toser_i(x))) break;
2101 if (gcmp0(y)) return gcopy(y);
2102 if (valp(y) < 0) pari_err(negexper,"gsin");
2103 gsincos(y,&u,&v,prec);
2104 return gerepilecopy(av,u);
2105 }
2106 return transc(gsin,x,prec);
2107 }
2108 /********************************************************************/
2109 /** SINE, COSINE together **/
2110 /********************************************************************/
2111
2112 void
mpsincos(GEN x,GEN * s,GEN * c)2113 mpsincos(GEN x, GEN *s, GEN *c)
2114 {
2115 long mod8;
2116 pari_sp av, tetpil;
2117 GEN p1, *gptr[2];
2118
2119 if (!signe(x))
2120 {
2121 long e = expo(x);
2122 *s = real_0_bit(e);
2123 *c = e >= 0? real_0_bit(e): real_1(nbits2prec(-e));
2124 return;
2125 }
2126
2127 av=avma; p1=mpsc1(x,&mod8); tetpil=avma;
2128 switch(mod8)
2129 {
2130 case 0: *c=addsr( 1,p1); *s=mpaut(p1); break;
2131 case 1: *s=addsr( 1,p1); *c=mpaut(p1); setsigne(*c,-signe(*c)); break;
2132 case 2: *c=subsr(-1,p1); *s=mpaut(p1); setsigne(*s,-signe(*s)); break;
2133 case 3: *s=subsr(-1,p1); *c=mpaut(p1); break;
2134 case 4: *c=addsr( 1,p1); *s=mpaut(p1); setsigne(*s,-signe(*s)); break;
2135 case 5: *s=addsr( 1,p1); *c=mpaut(p1); break;
2136 case 6: *c=subsr(-1,p1); *s=mpaut(p1); break;
2137 case 7: *s=subsr(-1,p1); *c=mpaut(p1); setsigne(*c,-signe(*c)); break;
2138 }
2139 gptr[0]=s; gptr[1]=c;
2140 gerepilemanysp(av,tetpil,gptr,2);
2141 }
2142
2143 /* return exp(ix), x a t_REAL */
2144 GEN
exp_Ir(GEN x)2145 exp_Ir(GEN x)
2146 {
2147 pari_sp av = avma;
2148 GEN v = cgetg(3,t_COMPLEX);
2149 mpsincos(x, (GEN*)(v+2), (GEN*)(v+1));
2150 if (!signe(x)) return gerepilecopy(av, gel(v,1));
2151 return v;
2152 }
2153
2154 void
gsincos(GEN x,GEN * s,GEN * c,long prec)2155 gsincos(GEN x, GEN *s, GEN *c, long prec)
2156 {
2157 long ii, i, j, ex, ex2, lx, ly, mi;
2158 pari_sp av, tetpil;
2159 GEN y, r, u, v, u1, v1, p1, p2, p3, p4, ps, pc;
2160 GEN *gptr[4];
2161
2162 switch(typ(x))
2163 {
2164 case t_INT: case t_FRAC:
2165 *s = cgetr(prec);
2166 *c = cgetr(prec); av = avma;
2167 mpsincos(tofp_safe(x, prec), &ps, &pc);
2168 affr_fixlg(ps,*s);
2169 affr_fixlg(pc,*c); avma = av; return;
2170
2171 case t_REAL:
2172 mpsincos(x,s,c); return;
2173
2174 case t_COMPLEX:
2175 i = precision(x); if (!i) i = prec;
2176 ps = cgetc(i); *s = ps;
2177 pc = cgetc(i); *c = pc; av = avma;
2178 r = gexp(gel(x,2),prec);
2179 v1 = gmul2n(addrr(ginv(r),r), -1); /* = cos(I*Im(x)) */
2180 u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
2181 gsincos(gel(x,1), &u,&v, prec);
2182 affr_fixlg(mulrr(v1,u), gel(ps,1));
2183 affr_fixlg(mulrr(u1,v), gel(ps,2));
2184 affr_fixlg(mulrr(v1,v), gel(pc,1));
2185 affr_fixlg(mulrr(mpneg(u1),u),gel(pc,2)); return;
2186
2187 case t_QUAD:
2188 av = avma; gsincos(quadtoc(x, prec), s, c, prec);
2189 gerepileall(av, 2, s, c); return;
2190
2191 default:
2192 av = avma; if (!(y = toser_i(x))) break;
2193 if (gcmp0(y)) { *c = gaddsg(1,y); *s = gcopy(y); return; }
2194
2195 ex = valp(y); lx = lg(y); ex2 = 2*ex+2;
2196 if (ex < 0) pari_err(talker,"non zero exponent in gsincos");
2197 if (ex2 > lx)
2198 {
2199 *s = x == y? gcopy(y): gerepilecopy(av, y); av = avma;
2200 *c = gerepileupto(av, gsubsg(1, gdivgs(gsqr(y),2)));
2201 return;
2202 }
2203 if (!ex)
2204 {
2205 p1 = shallowcopy(y); gel(p1,2) = gen_0;
2206 gsincos(normalize(p1),&u,&v,prec);
2207 gsincos(gel(y,2),&u1,&v1,prec);
2208 p1 = gmul(v1,v);
2209 p2 = gmul(u1,u);
2210 p3 = gmul(v1,u);
2211 p4 = gmul(u1,v); tetpil = avma;
2212 *c = gsub(p1,p2);
2213 *s = gadd(p3,p4);
2214 gptr[0]=s; gptr[1]=c;
2215 gerepilemanysp(av,tetpil,gptr,2);
2216 return;
2217 }
2218
2219 ly = lx+2*ex;
2220 mi = lx-1; while (mi>=3 && isexactzero(gel(y,mi))) mi--;
2221 mi += ex-2;
2222 pc = cgetg(ly,t_SER); *c = pc;
2223 ps = cgetg(lx,t_SER); *s = ps;
2224 pc[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(y));
2225 gel(pc,2) = gen_1; ps[1] = y[1];
2226 for (i=2; i<ex+2; i++) gel(ps,i) = gcopy(gel(y,i));
2227 for (i=3; i< ex2; i++) gel(pc,i) = gen_0;
2228 for (i=ex2; i<ly; i++)
2229 {
2230 ii = i-ex; av = avma; p1 = gen_0;
2231 for (j=ex; j<=min(ii-2,mi); j++)
2232 p1 = gadd(p1, gmulgs(gmul(gel(y,j-ex+2),gel(ps,ii-j)),j));
2233 gel(pc,i) = gerepileupto(av, gdivgs(p1,2-i));
2234 if (ii < lx)
2235 {
2236 av = avma; p1 = gen_0;
2237 for (j=ex; j<=min(i-ex2,mi); j++)
2238 p1 = gadd(p1,gmulgs(gmul(gel(y,j-ex+2),gel(pc,i-j)),j));
2239 p1 = gdivgs(p1,i-2);
2240 gel(ps,i-ex) = gerepileupto(av, gadd(p1,gel(y,i-ex)));
2241 }
2242 }
2243 return;
2244 }
2245 pari_err(typeer,"gsincos");
2246 }
2247
2248 /********************************************************************/
2249 /** **/
2250 /** TANGENT and COTANGENT **/
2251 /** **/
2252 /********************************************************************/
2253 static GEN
mptan(GEN x)2254 mptan(GEN x)
2255 {
2256 pari_sp av = avma;
2257 GEN s, c;
2258
2259 mpsincos(x,&s,&c);
2260 if (!signe(c)) pari_err(talker, "can't compute tan(Pi/2 + kPi)");
2261 return gerepileuptoleaf(av, divrr(s,c));
2262 }
2263
2264 GEN
gtan(GEN x,long prec)2265 gtan(GEN x, long prec)
2266 {
2267 pari_sp av;
2268 GEN y, s, c;
2269
2270 switch(typ(x))
2271 {
2272 case t_REAL: return mptan(x);
2273
2274 case t_COMPLEX:
2275 av = avma; gsincos(x,&s,&c,prec);
2276 return gerepileupto(av, gdiv(s,c));
2277
2278 case t_INT: case t_FRAC:
2279 y = cgetr(prec); av = avma;
2280 affr_fixlg(mptan(tofp_safe(x,prec)), y); avma = av; return y;
2281
2282 case t_PADIC:
2283 av = avma;
2284 return gerepileupto(av, gdiv(gsin(x,prec), gcos(x,prec)));
2285
2286 case t_INTMOD: pari_err(typeer,"gtan");
2287
2288 default:
2289 av = avma; if (!(y = toser_i(x))) break;
2290 if (gcmp0(y)) return gcopy(y);
2291 if (valp(y) < 0) pari_err(negexper,"gtan");
2292 gsincos(y,&s,&c,prec);
2293 return gerepileupto(av, gdiv(s,c));
2294 }
2295 return transc(gtan,x,prec);
2296 }
2297
2298 static GEN
mpcotan(GEN x)2299 mpcotan(GEN x)
2300 {
2301 pari_sp av=avma, tetpil;
2302 GEN s,c;
2303
2304 mpsincos(x,&s,&c); tetpil=avma;
2305 return gerepile(av,tetpil,divrr(c,s));
2306 }
2307
2308 GEN
gcotan(GEN x,long prec)2309 gcotan(GEN x, long prec)
2310 {
2311 pari_sp av;
2312 GEN y, s, c;
2313
2314 switch(typ(x))
2315 {
2316 case t_REAL:
2317 return mpcotan(x);
2318
2319 case t_COMPLEX:
2320 av = avma; gsincos(x,&s,&c,prec);
2321 return gerepileupto(av, gdiv(c,s));
2322
2323 case t_INT: case t_FRAC:
2324 y = cgetr(prec); av = avma;
2325 affr_fixlg(mpcotan(tofp_safe(x,prec)), y); avma = av; return y;
2326
2327 case t_PADIC:
2328 av = avma;
2329 return gerepileupto(av, gdiv(gcos(x,prec), gsin(x,prec)));
2330
2331 case t_INTMOD: pari_err(typeer,"gcotan");
2332
2333 default:
2334 av = avma; if (!(y = toser_i(x))) break;
2335 if (gcmp0(y)) pari_err(talker,"0 argument in cotan");
2336 if (valp(y) < 0) pari_err(negexper,"cotan"); /* fall through */
2337 gsincos(y,&s,&c,prec);
2338 return gerepileupto(av, gdiv(c,s));
2339 }
2340 return transc(gcotan,x,prec);
2341 }
2342