Annotation of OpenXM_contrib/pari-2.2/src/basemath/trans1.c, Revision 1.2
1.2 ! noro 1: /* $Id: trans1.c,v 1.69 2002/08/08 10:24:22 karim Exp $
1.1 noro 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 FONCTIONS **/
19: /** **/
20: /********************************************************************/
21: #include "pari.h"
22:
23: #ifdef LONG_IS_64BIT
24: # define C31 (9223372036854775808.0) /* VERYBIGINT * 1.0 */
25: # define SQRTVERYBIGINT 3037000500 /* ceil(sqrt(VERYBIGINT)) */
26: # define CBRTVERYBIGINT 2097152 /* ceil(cbrt(VERYBIGINT)) */
27: #else
28: # define C31 (2147483648.0)
29: # define SQRTVERYBIGINT 46341
30: # define CBRTVERYBIGINT 1291
31: #endif
32:
33:
34: /********************************************************************/
35: /** **/
36: /** FONCTION PI **/
37: /** **/
38: /********************************************************************/
39:
40: void
41: constpi(long prec)
42: {
43: GEN p1,p2,p3,tmppi;
1.2 ! noro 44: long l, n, n1;
! 45: gpmem_t av1, av2;
1.1 noro 46: double alpha;
47:
48: #define k1 545140134
49: #define k2 13591409
50: #define k3 640320
51: #define alpha2 (1.4722004/(BYTES_IN_LONG/4)) /* 3*log(k3/12)/(32*log(2)) */
52:
53: if (gpi && lg(gpi) >= prec) return;
54:
55: av1 = avma; tmppi = newbloc(prec);
56: *tmppi = evaltyp(t_REAL) | evallg(prec);
57:
58: prec++;
59:
1.2 ! noro 60: n = (long)(1 + (prec-2)/alpha2);
! 61: n1 = 6*n - 1;
! 62: p2 = addsi(k2, mulss(n,k1));
! 63: p1 = itor(p2, prec);
1.1 noro 64:
65: /* initialisation longueur mantisse */
66: if (prec>=4) l=4; else l=prec;
67: setlg(p1,l); alpha=l;
68:
69: av2 = avma;
70: while (n)
71: {
72: if (n < CBRTVERYBIGINT)
73: p3 = divrs(mulsr(n1-4,mulsr(n1*(n1-2),p1)),n*n*n);
74: else
75: {
76: if (n1 < SQRTVERYBIGINT)
77: p3 = divrs(divrs(mulsr(n1-4,mulsr(n1*(n1-2),p1)),n*n),n);
78: else
79: p3 = divrs(divrs(divrs(mulsr(n1-4,mulsr(n1,mulsr(n1-2,p1))),n),n),n);
80: }
81: p3 = divrs(divrs(p3,100100025), 327843840);
82: subisz(p2,k1,p2); subirz(p2,p3,p1);
83: alpha += alpha2; l = (long)(1+alpha); /* nouvelle longueur mantisse */
84: if (l>prec) l=prec;
85: setlg(p1,l); avma = av2;
86: n--; n1-=6;
87: }
88: p1 = divsr(53360,p1);
89: mulrrz(p1,gsqrt(stoi(k3),prec), tmppi);
90: gunclone(gpi); avma = av1; gpi = tmppi;
91: }
92:
93: GEN
94: mppi(long prec)
95: {
96: GEN x = cgetr(prec);
97: constpi(prec); affrr(gpi,x); return x;
98: }
99:
1.2 ! noro 100: /* Pi * 2^n */
! 101: GEN
! 102: Pi2n(long n, long prec)
! 103: {
! 104: GEN x = mppi(prec); setexpo(x, 1+n);
! 105: return x;
! 106: }
! 107:
! 108: /* I * Pi * 2^n */
! 109: GEN
! 110: PiI2n(long n, long prec)
! 111: {
! 112: GEN z = cgetg(3,t_COMPLEX);
! 113: z[1] = zero;
! 114: z[2] = (long)Pi2n(n, prec); return z;
! 115: }
! 116:
! 117: /* 2I * Pi */
! 118: GEN
! 119: PiI2(long prec) { return PiI2n(1, prec); }
! 120:
1.1 noro 121: /********************************************************************/
122: /** **/
123: /** FONCTION EULER **/
124: /** **/
125: /********************************************************************/
126:
127: void
128: consteuler(long prec)
129: {
130: GEN u,v,a,b,tmpeuler;
1.2 ! noro 131: long l, n, k, x;
! 132: gpmem_t av1, av2;
1.1 noro 133:
134: if (geuler && lg(geuler) >= prec) return;
135:
136: av1 = avma; tmpeuler = newbloc(prec);
137: *tmpeuler = evaltyp(t_REAL) | evallg(prec);
138:
139: prec++;
140:
141: l = prec+1; x = (long) (1 + (bit_accuracy(l) >> 2) * LOG2);
1.2 ! noro 142: a = stor(x,l); u=mplog(a); setsigne(u,-1); affrr(u,a);
! 143: b = realun(l);
! 144: v = realun(l);
1.1 noro 145: n=(long)(1+3.591*x); /* z=3.591: z*[ ln(z)-1 ]=1 */
146: av2 = avma;
147: if (x < SQRTVERYBIGINT)
148: {
149: long xx = x*x;
150: for (k=1; k<=n; k++)
151: {
152: divrsz(mulsr(xx,b),k*k, b);
153: divrsz(addrr(divrs(mulsr(xx,a),k),b),k, a);
154: addrrz(u,a,u); addrrz(v,b,v); avma = av2;
155: }
156: }
157: else
158: {
159: GEN xx = mulss(x,x);
160: for (k=1; k<=n; k++)
161: {
162: divrsz(mulir(xx,b),k*k, b);
163: divrsz(addrr(divrs(mulir(xx,a),k),b),k, a);
164: addrrz(u,a,u); addrrz(v,b,v); avma = av2;
165: }
166: }
167: divrrz(u,v,tmpeuler);
168: gunclone(geuler); avma = av1; geuler = tmpeuler;
169: }
170:
171: GEN
172: mpeuler(long prec)
173: {
174: GEN x = cgetr(prec);
175: consteuler(prec); affrr(geuler,x); return x;
176: }
177:
178: /********************************************************************/
179: /** **/
180: /** CONVERSION DE TYPES POUR FONCTIONS TRANSCENDANTES **/
181: /** **/
182: /********************************************************************/
183:
184: GEN
185: transc(GEN (*f)(GEN,long), GEN x, long prec)
186: {
187: GEN p1,p2,y;
1.2 ! noro 188: long lx, i;
! 189: gpmem_t av, tetpil;
1.1 noro 190:
191: av=avma;
192: switch(typ(x))
193: {
194: case t_INT: case t_FRAC: case t_FRACN:
195: p1=cgetr(prec); gaffect(x,p1); tetpil=avma;
196: return gerepile(av,tetpil,f(p1,prec));
197:
198: case t_COMPLEX: case t_QUAD:
199: av=avma; p1=gmul(x,realun(prec)); tetpil=avma;
200: return gerepile(av,tetpil,f(p1,prec));
201:
202: case t_POL: case t_RFRAC: case t_RFRACN:
203: p1=tayl(x,gvar(x),precdl); tetpil=avma;
204: return gerepile(av,tetpil,f(p1,prec));
205:
206: case t_VEC: case t_COL: case t_MAT:
207: lx=lg(x); y=cgetg(lx,typ(x));
208: for (i=1; i<lx; i++)
209: y[i] = (long) f((GEN)x[i],prec);
210: return y;
211:
212: case t_POLMOD:
213: av=avma; p1=roots((GEN)x[1],prec); lx=lg(p1); p2=cgetg(lx,t_COL);
214: for (i=1; i<lx; i++) p2[i]=lpoleval((GEN)x[2],(GEN)p1[i]);
215: tetpil=avma; y=cgetg(lx,t_COL);
216: for (i=1; i<lx; i++)
217: y[i] = (long)f((GEN)p2[i],prec);
218: return gerepile(av,tetpil,y);
219:
220: default:
221: err(typeer,"a transcendental function");
222: }
223: return f(x,prec);
224: }
225:
226: /*******************************************************************/
227: /* */
228: /* POWERING */
229: /* */
230: /*******************************************************************/
231: extern GEN real_unit_form(GEN x);
232: extern GEN imag_unit_form(GEN x);
233:
234: static GEN
235: puiss0(GEN x)
236: {
237: long lx,i;
238: GEN y;
239:
240: switch(typ(x))
241: {
242: case t_INT: case t_REAL: case t_FRAC: case t_FRACN: case t_PADIC:
243: return gun;
244:
245: case t_INTMOD:
246: y=cgetg(3,t_INTMOD); copyifstack(x[1],y[1]); y[2]=un; break;
247:
248: case t_COMPLEX:
249: y=cgetg(3,t_COMPLEX); y[1]=un; y[2]=zero; break;
250:
251: case t_QUAD:
252: y=cgetg(4,t_QUAD); copyifstack(x[1],y[1]);
253: y[2]=un; y[3]=zero; break;
254:
255: case t_POLMOD:
256: y=cgetg(3,t_POLMOD); copyifstack(x[1],y[1]);
257: y[2]=lpolun[varn(x[1])]; break;
258:
259: case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:
260: return polun[gvar(x)];
261:
262: case t_MAT:
263: lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
264: if (lx != (lg(x[1]))) err(mattype1,"gpowgs");
265: y = idmat(lx-1);
266: for (i=1; i<lx; i++) coeff(y,i,i) = lpuigs(gcoeff(x,i,i),0);
267: break;
268: case t_QFR: return real_unit_form(x);
269: case t_QFI: return imag_unit_form(x);
270: case t_VECSMALL:
271: lx = lg(x);
272: y = cgetg(lx, t_VECSMALL);
273: for (i=1; i<lx; i++) y[i] = i;
274: return y;
275: default: err(typeer,"gpowgs");
276: return NULL; /* not reached */
277: }
278: return y;
279: }
280:
1.2 ! noro 281: static GEN
! 282: _sqr(void *data /* ignored */, GEN x) {
! 283: (void)data; return gsqr(x);
! 284: }
! 285: static GEN
! 286: _mul(void *data /* ignored */, GEN x, GEN y) {
! 287: (void)data; return gmul(x,y);
! 288: }
! 289: static GEN
! 290: _sqri(void *data /* ignored */, GEN x) {
! 291: (void)data; return sqri(x);
! 292: }
! 293: static GEN
! 294: _muli(void *data /* ignored */, GEN x, GEN y) {
! 295: (void)data; return mulii(x,y);
! 296: }
! 297:
1.1 noro 298: /* INTEGER POWERING (a^|n| for integer a and integer |n| > 1)
299: *
300: * Nonpositive powers and exponent one should be handled by gpow() and
301: * gpowgs() in trans1.c, which at the time of this writing is the only place
302: * where the following is [slated to be] used.
303: *
304: * Use left shift binary algorithm (RS is wasteful: multiplies big numbers,
305: * with LS one of them is the base, hence small). If result is nonzero, its
306: * sign is set to s (= +-1) regardless of what it would have been. This makes
307: * life easier for gpow()/gpowgs(), which otherwise might do a
308: * setsigne(gun,-1)... -- GN1998May04
309: */
310: static GEN
311: puissii(GEN a, GEN n, long s)
312: {
1.2 ! noro 313: gpmem_t av;
1.1 noro 314: GEN y;
315:
316: if (!signe(a)) return gzero; /* a==0 */
317: if (lgefint(a)==3)
318: { /* easy if |a| < 3 */
319: if (a[2] == 1) return (s>0)? gun: negi(gun);
320: if (a[2] == 2) { a = shifti(gun, labs(itos(n))); setsigne(a,s); return a; }
321: }
322: if (lgefint(n)==3)
323: { /* or if |n| < 3 */
324: if (n[2] == 1) { a = icopy(a); setsigne(a,s); return a; }
325: if (n[2] == 2) return sqri(a);
326: }
1.2 ! noro 327: av = avma;
! 328: y = leftright_pow(a, n, NULL, &_sqri, &_muli);
! 329: setsigne(y,s); return gerepileuptoint(av,y);
! 330: }
! 331:
! 332: typedef struct {
! 333: long prec, a;
! 334: GEN (*sqr)(GEN);
! 335: GEN (*mulsg)(long,GEN);
! 336: } sr_muldata;
! 337:
! 338: static GEN
! 339: _rpowsi_mul(void *data, GEN x, GEN y/*unused*/)
! 340: {
! 341: sr_muldata *D = (sr_muldata *)data;
! 342: return D->mulsg(D->a, x);
! 343: }
! 344:
! 345: static GEN
! 346: _rpowsi_sqr(void *data, GEN x)
! 347: {
! 348: sr_muldata *D = (sr_muldata *)data;
! 349: if (typ(x) == t_INT && lgefint(x) >= D->prec)
! 350: { /* switch to t_REAL */
! 351: D->sqr = &gsqr;
! 352: D->mulsg = &mulsr;
! 353: x = itor(x, D->prec);
1.1 noro 354: }
1.2 ! noro 355: return D->sqr(x);
1.1 noro 356: }
357:
358: /* return a^n as a t_REAL of precision prec. Adapted from puissii().
359: * Assume a > 0, n > 0 */
360: GEN
361: rpowsi(ulong a, GEN n, long prec)
362: {
1.2 ! noro 363: gpmem_t av;
! 364: GEN y;
! 365: sr_muldata D;
1.1 noro 366:
1.2 ! noro 367: if (a == 1) return realun(prec);
! 368: if (a == 2) return real2n(itos(n), prec);
! 369: if (is_pm1(n)) return stor(a, prec);
! 370: av = avma;
! 371: D.sqr = &sqri;
! 372: D.mulsg = &mulsi;
! 373: D.prec = prec;
! 374: D.a = a;
! 375: y = leftright_pow(stoi(a), n, (void*)&D, &_rpowsi_sqr, &_rpowsi_mul);
! 376: if (typ(y) == t_INT) y = itor(y, prec);
! 377: return gerepileuptoleaf(av, y);
1.1 noro 378: }
379:
380: GEN
381: gpowgs(GEN x, long n)
382: {
1.2 ! noro 383: long m, tx;
! 384: gpmem_t lim, av;
! 385: static long gn[3] = {evaltyp(t_INT)|_evallg(3), 0, 0};
1.1 noro 386: GEN y;
387:
388: if (n == 0) return puiss0(x);
389: if (n == 1) return gcopy(x);
390: if (n ==-1) return ginv(x);
391: if (n>0) { gn[1] = evalsigne( 1) | evallgefint(3); gn[2]= n; }
392: else { gn[1] = evalsigne(-1) | evallgefint(3); gn[2]=-n; }
393: /* If gpowgs were only ever called from gpow, the switch wouldn't be needed.
394: * In fact, under current word and bit field sizes, an integer power with
395: * multiword exponent will always overflow. But it seems easier to call
396: * puissii|powmodulo() with a mock-up GEN as 2nd argument than to write
397: * separate versions for a long exponent. Note that n = HIGHBIT is an
398: * invalid argument. --GN
399: */
400: switch((tx=typ(x)))
401: {
402: case t_INT:
403: {
404: long sx=signe(x), sr = (sx<0 && (n&1))? -1: 1;
405: if (n>0) return puissii(x,(GEN)gn,sr);
406: if (!sx) err(talker, "division by zero in gpowgs");
407: if (is_pm1(x)) return (sr < 0)? icopy(x): gun;
408: /* n<0, |x|>1 */
409: y=cgetg(3,t_FRAC); setsigne(gn,1);
410: y[1]=(sr>0)? un: lnegi(gun);
411: y[2]=(long)puissii(x,(GEN)gn,1); /* force denominator > 0 */
412: return y;
413: }
414: case t_INTMOD:
415: y=cgetg(3,tx); copyifstack(x[1],y[1]);
416: y[2]=(long)powmodulo((GEN)(x[2]),(GEN)gn,(GEN)(x[1]));
417: return y;
418: case t_FRAC: case t_FRACN:
419: {
420: GEN a = (GEN)x[1], b = (GEN)x[2];
421: long sr = (n&1 && (signe(a)!=signe(b))) ? -1 : 1;
422: if (n > 0) { if (!signe(a)) return gzero; }
423: else
424: { /* n < 0 */
425: if (!signe(a)) err(talker, "division by zero fraction in gpowgs");
426: /* +-1/x[2] inverts to an integer */
427: if (is_pm1(a)) return puissii(b,(GEN)gn,sr);
428: y = b; b = a; a = y;
429: }
430: /* HACK: puissii disregards the sign of gn */
431: y = cgetg(3,tx);
432: y[1] = (long)puissii(a,(GEN)gn,sr);
433: y[2] = (long)puissii(b,(GEN)gn,1);
434: return y;
435: }
436: case t_PADIC: case t_POL: case t_POLMOD:
437: return powgi(x,gn);
438: case t_RFRAC: case t_RFRACN:
439: {
440: av=avma; y=cgetg(3,tx); m = labs(n);
441: y[1]=lpuigs((GEN)x[1],m);
442: y[2]=lpuigs((GEN)x[2],m);
443: if (n<0) y=ginv(y); /* let ginv worry about normalizations */
444: return gerepileupto(av,y);
445: }
446: default:
447: m = labs(n);
448: av=avma; y=NULL; lim=stack_lim(av,1);
449: for (; m>1; m>>=1)
450: {
451: if (m&1) y = y? gmul(y,x): x;
452: x=gsqr(x);
453: if (low_stack(lim, stack_lim(av,1)))
454: {
455: GEN *gptr[2]; gptr[0]=&x; gptr[1]=&y;
456: if(DEBUGMEM>1) err(warnmem,"[3]: gpowgs");
457: gerepilemany(av,gptr,y? 2: 1);
458: }
459: }
460: y = y? gmul(y,x): x;
461: if (n<=0) y=ginv(y);
462: return gerepileupto(av,y);
463: }
464: }
465:
466: GEN
467: pow_monome(GEN x, GEN nn)
468: {
1.2 ! noro 469: long n, m, i, j, dd;
! 470: gpmem_t tetpil, av = avma;
1.1 noro 471: GEN p1,y;
472: if (is_bigint(nn)) err(talker,"power overflow in pow_monome");
473: n = itos(nn); m = labs(n);
474: j=lgef(x); dd=(j-3)*m+3;
475: p1=cgetg(dd,t_POL); m = labs(n);
476: p1[1] = evalsigne(1) | evallgef(dd) | evalvarn(varn(x));
477: for (i=2; i<dd-1; i++) p1[i]=zero;
478: p1[i]=lpuigs((GEN)x[j-1],m);
479: if (n>0) return p1;
480:
481: tetpil=avma; y=cgetg(3,t_RFRAC);
482: y[1]=(long)denom((GEN)p1[i]);
483: y[2]=lmul(p1,(GEN)y[1]); return gerepile(av,tetpil,y);
484: }
485:
486: extern GEN powrealform(GEN x, GEN n);
487:
488: /* n is assumed to be an integer */
489: GEN
490: powgi(GEN x, GEN n)
491: {
1.2 ! noro 492: long tx, sn = signe(n);
! 493: GEN y;
1.1 noro 494:
495: if (typ(n) != t_INT) err(talker,"not an integral exponent in powgi");
496: if (!sn) return puiss0(x);
497:
498: switch(tx=typ(x))
499: {/* catch some types here, instead of trying gpowgs() first, because of
500: * the simpler call interface to puissii() and powmodulo() -- even though
501: * for integer/rational bases other than 0,+-1 and non-wordsized
502: * exponents the result will be certain to overflow. --GN
503: */
504: case t_INT:
505: {
506: long sx=signe(x), sr = (sx<0 && mod2(n))? -1: 1;
507: if (sn>0) return puissii(x,n,sr);
508: if (!sx) err(talker, "division by zero in powgi");
509: if (is_pm1(x)) return (sr < 0)? icopy(x): gun;
510: /* n<0, |x|>1 */
511: y=cgetg(3,t_FRAC); setsigne(n,1); /* temporarily replace n by abs(n) */
512: y[1]=(sr>0)? un: lnegi(gun);
513: y[2]=(long)puissii(x,n,1);
514: setsigne(n,-1); return y;
515: }
516: case t_INTMOD:
517: y=cgetg(3,tx); copyifstack(x[1],y[1]);
518: y[2]=(long)powmodulo((GEN)x[2],n,(GEN)x[1]);
519: return y;
520: case t_FRAC: case t_FRACN:
521: {
522: GEN a = (GEN)x[1], b = (GEN)x[2];
523: long sr = (mod2(n) && (signe(a)!=signe(b))) ? -1 : 1;
524: if (sn > 0) { if (!signe(a)) return gzero; }
525: else
526: { /* n < 0 */
527: if (!signe(a)) err(talker, "division by zero fraction in powgi");
528: /* +-1/b inverts to an integer */
529: if (is_pm1(a)) return puissii(b,n,sr);
530: y = b; b = a; a = y;
531: }
532: /* HACK: puissii disregards the sign of n */
533: y = cgetg(3,tx);
534: y[1] = (long)puissii(a,n,sr);
535: y[2] = (long)puissii(b,n,1);
536: return y;
537: }
538: case t_PADIC:
539: {
1.2 ! noro 540: long e = itos(n)*valp(x), v;
! 541: GEN mod, p = (GEN)x[2];
! 542:
1.1 noro 543: if (!signe(x[4]))
544: {
545: if (sn < 0) err(talker, "division by 0 p-adic in powgi");
1.2 ! noro 546: return padiczero(p, e);
1.1 noro 547: }
548: y = cgetg(5,t_PADIC);
1.2 ! noro 549: mod = (GEN)x[3]; v = ggval(n, p);
! 550: if (v == 0) mod = icopy(mod);
1.1 noro 551: else
552: {
1.2 ! noro 553: mod = mulii(mod, gpowgs(p,v));
! 554: mod = gerepileuptoint((gpmem_t)y, mod);
1.1 noro 555: }
1.2 ! noro 556: y[1] = evalprecp(precp(x)+v) | evalvalp(e);
1.1 noro 557: icopyifstack(p, y[2]);
1.2 ! noro 558: y[3] = (long)mod;
! 559: y[4] = (long)powmodulo((GEN)x[4], n, mod);
1.1 noro 560: return y;
561: }
562: case t_QFR:
563: if (signe(x[4])) return powrealform(x,n);
564: case t_POL:
565: if (ismonome(x)) return pow_monome(x,n);
566: default:
1.2 ! noro 567: {
! 568: gpmem_t av = avma;
! 569: y = leftright_pow(x, n, NULL, &_sqr, &_mul);
1.1 noro 570: if (sn < 0) y = ginv(y);
571: return av==avma? gcopy(y): gerepileupto(av,y);
1.2 ! noro 572: }
1.1 noro 573: }
574: }
575:
576: /* we suppose n != 0, and valp(x) = 0 */
577: static GEN
578: ser_pui(GEN x, GEN n, long prec)
579: {
1.2 ! noro 580: gpmem_t av, tetpil;
! 581: GEN y, p1, p2, lead;
1.1 noro 582:
1.2 ! noro 583: if (gvar9(n) <= varn(x)) return gexp(gmul(n, glog(x,prec)), prec);
! 584: lead = (GEN)x[2];
! 585: if (gcmp1(lead))
! 586: {
! 587: GEN X, Y, alp;
! 588: long lx, mi, i, j, d;
! 589:
! 590: av = avma; alp = gclone(gadd(n,gun)); avma = av;
! 591: lx = lg(x);
! 592: y = cgetg(lx,t_SER);
! 593: y[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(x));
! 594: X = x+2;
! 595: Y = y+2;
! 596:
! 597: d = mi = lx-3; while (mi>=1 && gcmp0((GEN)X[mi])) mi--;
! 598: Y[0] = un;
! 599: for (i=1; i<=d; i++)
1.1 noro 600: {
1.2 ! noro 601: av = avma; p1 = gzero;
! 602: for (j=1; j<=min(i,mi); j++)
1.1 noro 603: {
1.2 ! noro 604: p2 = gsubgs(gmulgs(alp,j),i);
! 605: p1 = gadd(p1, gmul(gmul(p2,(GEN)X[j]),(GEN)Y[i-j]));
1.1 noro 606: }
1.2 ! noro 607: tetpil = avma; Y[i] = lpile(av,tetpil, gdivgs(p1,i));
1.1 noro 608: }
1.2 ! noro 609: gunclone(alp); return y;
1.1 noro 610: }
1.2 ! noro 611: p1 = gdiv(x,lead); p1[2] = un; /* in case it's inexact */
! 612: p1 = gpow(p1, n, prec);
! 613: p2 = gpow(lead,n, prec); return gmul(p1, p2);
1.1 noro 614: }
615:
616: GEN
617: gpow(GEN x, GEN n, long prec)
618: {
1.2 ! noro 619: long i, lx, tx;
! 620: gpmem_t av, tetpil;
1.1 noro 621: GEN y;
622:
623: if (typ(n)==t_INT) return powgi(x,n);
624: tx = typ(x);
625: if (is_matvec_t(tx))
626: {
627: lx=lg(x); y=cgetg(lx,tx);
628: for (i=1; i<lx; i++) y[i]=lpui((GEN)x[i],n,prec);
629: return y;
630: }
1.2 ! noro 631: if (tx==t_POL)
! 632: {
! 633: av=avma; y=tayl(x,gvar(x),precdl); tetpil=avma;
! 634: return gerepile(av,tetpil,gpow(y,n,prec));
! 635: }
1.1 noro 636: if (tx==t_SER)
637: {
638: if (valp(x))
639: err(talker,"not an integer exponent for non invertible series in gpow");
640: if (lg(x) == 2) return gcopy(x); /* O(1) */
1.2 ! noro 641: av = avma;
! 642: return gerepileupto(av, ser_pui(x, n, prec));
1.1 noro 643: }
644: av=avma;
645: if (gcmp0(x))
646: {
647: long tn = typ(n);
648: if (!is_scalar_t(tn) || tn == t_PADIC || tn == t_INTMOD)
649: err(talker,"zero to a forbidden power in gpow");
650: n = greal(n);
651: if (gsigne(n) <= 0)
652: err(talker,"zero to a non positive exponent in gpow");
653: if (!precision(x)) return gcopy(x);
654:
655: x = ground(gmulsg(gexpo(x),n));
656: if (is_bigint(x) || (ulong)x[2] >= (ulong)HIGHEXPOBIT)
657: err(talker,"underflow or overflow in gpow");
1.2 ! noro 658: avma = av; return realzero_bit(itos(x));
1.1 noro 659: }
660: if (tx==t_INTMOD && typ(n)==t_FRAC)
661: {
662: GEN p1;
1.2 ! noro 663: if (!BSW_psp((GEN)x[1])) err(talker,"modulus must be prime in gpow");
1.1 noro 664: y=cgetg(3,tx); copyifstack(x[1],y[1]);
665: av=avma;
666: p1=mpsqrtnmod((GEN)x[2],(GEN)n[2],(GEN)x[1],NULL);
667: if(!p1) err(talker,"n-root does not exists in gpow");
668: p1=powmodulo(p1,(GEN)n[1],(GEN)x[1]);
669: y[2]=lpileupto(av,p1);
670: return y;
671: }
672: i = (long) precision(n); if (i) prec=i;
673: y=gmul(n,glog(x,prec)); tetpil=avma;
674: return gerepile(av,tetpil,gexp(y,prec));
675: }
676:
677: /********************************************************************/
678: /** **/
679: /** RACINE CARREE **/
680: /** **/
681: /********************************************************************/
682:
683: GEN
684: mpsqrt(GEN x)
685: {
1.2 ! noro 686: gpmem_t av, av0;
! 687: long l,l0,l1,l2,s,eps,n,i,ex;
1.1 noro 688: double beta;
689: GEN y,p1,p2,p3;
690:
1.2 ! noro 691: if (typ(x) != t_REAL) err(typeer,"mpsqrt");
! 692: s = signe(x); if (s < 0) err(talker,"negative argument in mpsqrt");
! 693: if (!s) return realzero_bit(expo(x) >> 1);
! 694:
! 695: l = lg(x); y = cgetr(l); av0 = avma;
! 696:
! 697: p1 = cgetr(l+1); affrr(x,p1);
! 698: ex = expo(x); eps = ex & 1; ex >>= 1;
! 699: setexpo(p1,eps); /* x = 2^(2 ex) p1 */
! 700:
! 701: n = (long)(2 + log((double) (l-2))/LOG2);
! 702: beta = sqrt((eps+1) * (2 + p1[2]/C31));
! 703: p2 = cgetr(l+1);
! 704: p2[1] = evalexpo(0) | evalsigne(1);
! 705: p2[2] = (long)((beta-2)*C31);
! 706: if (!p2[2]) { p2[2] = (long)HIGHBIT; setexpo(p2,1); }
! 707: for (i=3; i<=l; i++) p2[i] = 0;
1.1 noro 708:
1.2 ! noro 709: p3 = cgetr(l+1); l -= 2; l1 = 1; l2 = 3;
! 710: av = avma;
1.1 noro 711: for (i=1; i<=n; i++)
1.2 ! noro 712: { /* execute p2 := (p2 + p1/p2)/2 */
1.1 noro 713: l0 = l1<<1;
1.2 ! noro 714: if (l0 <= l)
1.1 noro 715: { l2 += l1; l1=l0; }
716: else
717: { l2 += l+1-l1; l1=l+1; }
1.2 ! noro 718: setlg(p1, l1+2);
! 719: setlg(p3, l1+2);
! 720: setlg(p2, l2);
! 721: affrr(divrr(p1,p2), p3);
! 722: affrr(addrr(p2,p3), p2); setexpo(p2, expo(p2)-1);
! 723: avma = av;
! 724: }
! 725: affrr(p2,y); setexpo(y, expo(y)+ex);
! 726: avma = av0; return y;
! 727: }
! 728:
! 729: /* O(p^e) */
! 730: GEN
! 731: padiczero(GEN p, long e)
! 732: {
! 733: GEN y = cgetg(5,t_PADIC);
! 734: y[4] = zero;
! 735: y[3] = un;
! 736: copyifstack(p,y[2]);
! 737: y[1] = evalvalp(e) | evalprecp(0);
! 738: return y;
! 739: }
! 740:
! 741: /* assume x unit, precp(x) = pp > 3 */
! 742: static GEN
! 743: sqrt_2adic(GEN x, long pp)
! 744: {
! 745: GEN z = mod16(x)==1? gun: stoi(3);
! 746: long zp;
! 747: gpmem_t av, lim;
! 748:
! 749: if (pp == 4) return z;
! 750: zp = 3; /* number of correct bits in z (compared to sqrt(x)) */
! 751:
! 752: av = avma; lim = stack_lim(av,2);
! 753: for(;;)
! 754: {
! 755: GEN mod;
! 756: zp = (zp<<1) - 1;
! 757: if (zp > pp) zp = pp;
! 758: mod = shifti(gun, zp);
! 759: z = addii(z, resmod2n(mulii(x, mpinvmod(z,mod)), zp));
! 760: z = shifti(z, -1); /* (z + x/z) / 2 */
! 761: if (pp == zp) return z;
! 762:
! 763: if (zp < pp) zp--;
! 764:
! 765: if (low_stack(lim,stack_lim(av,2)))
! 766: {
! 767: if (DEBUGMEM > 1) err(warnmem,"padic_sqrt");
! 768: z = gerepileuptoint(av,z);
! 769: }
! 770: }
! 771: }
! 772:
! 773: /* x unit defined modulo modx = p^pp, p != 2, pp > 0 */
! 774: static GEN
! 775: sqrt_padic(GEN x, GEN modx, long pp, GEN p)
! 776: {
! 777: GEN mod, z = mpsqrtmod(x, p);
! 778: long zp = 1;
! 779: gpmem_t av, lim;
! 780:
! 781: if (!z) err(sqrter5);
! 782: if (pp <= zp) return z;
! 783:
! 784: av = avma; lim = stack_lim(av,2);
! 785: mod = p;
! 786: for(;;)
! 787: { /* could use the hensel_lift_accel idea. Not really worth it */
! 788: GEN inv2;
! 789: zp <<= 1;
! 790: if (zp < pp) mod = sqri(mod); else { zp = pp; mod = modx; }
! 791: inv2 = shifti(mod, -1); /* = (mod + 1)/2 = 1/2 */
! 792: z = addii(z, resii(mulii(x, mpinvmod(z,mod)), mod));
! 793: z = mulii(z, inv2);
! 794: z = modii(z, mod); /* (z + x/z) / 2 */
! 795: if (pp <= zp) return z;
! 796:
! 797: if (low_stack(lim,stack_lim(av,2)))
! 798: {
! 799: GEN *gptr[2]; gptr[0]=&z; gptr[1]=&mod;
! 800: if (DEBUGMEM>1) err(warnmem,"padic_sqrt");
! 801: gerepilemany(av,gptr,2);
! 802: }
1.1 noro 803: }
804: }
805:
806: static GEN
807: padic_sqrt(GEN x)
808: {
1.2 ! noro 809: long pp, e = valp(x);
! 810: GEN z,y,mod, p = (GEN)x[2];
! 811: gpmem_t av;
1.1 noro 812:
1.2 ! noro 813: if (gcmp0(x)) return padiczero(p, (e+1) >> 1);
1.1 noro 814: if (e & 1) err(sqrter6);
1.2 ! noro 815:
! 816: y = cgetg(5,t_PADIC);
1.1 noro 817: pp = precp(x);
1.2 ! noro 818: mod = (GEN)x[3];
! 819: x = (GEN)x[4]; /* lift to int */
! 820: e >>= 1; av = avma;
! 821: if (egalii(gdeux, p))
! 822: {
! 823: long r = mod8(x);
! 824: if (pp <= 3)
! 825: {
! 826: switch(pp) {
! 827: case 1: break;
! 828: case 2: if ((r&3) == 1) break;
! 829: case 3: if (r == 1) break;
! 830: default: err(sqrter5);
! 831: }
! 832: z = gun;
! 833: pp = 1;
! 834: }
! 835: else
1.1 noro 836: {
1.2 ! noro 837: if (r != 1) err(sqrter5);
! 838: z = sqrt_2adic(x, pp);
! 839: z = gerepileuptoint(av, z);
! 840: pp--;
1.1 noro 841: }
1.2 ! noro 842: mod = shifti(gun, pp);
1.1 noro 843: }
844: else /* p != 2 */
845: {
1.2 ! noro 846: z = sqrt_padic(x, mod, pp, p);
! 847: z = gerepileuptoint(av, z);
! 848: mod = icopy(mod);
! 849: }
! 850: y[1] = evalprecp(pp) | evalvalp(e);
! 851: copyifstack(p,y[2]);
! 852: y[3] = (long)mod;
! 853: y[4] = (long)z; return y;
1.1 noro 854: }
855:
856: GEN
857: gsqrt(GEN x, long prec)
858: {
1.2 ! noro 859: long e;
! 860: gpmem_t av, tetpil;
1.1 noro 861: GEN y,p1,p2;
862:
863: switch(typ(x))
864: {
865: case t_REAL:
866: if (signe(x)>=0) return mpsqrt(x);
867: y=cgetg(3,t_COMPLEX); y[1]=zero;
868: setsigne(x,1); y[2]=lmpsqrt(x); setsigne(x,-1);
869: return y;
870:
871: case t_INTMOD:
872: y=cgetg(3,t_INTMOD); copyifstack(x[1],y[1]);
873: p1 = mpsqrtmod((GEN)x[2],(GEN)y[1]);
874: if (!p1) err(sqrter5);
875: y[2] = (long)p1; return y;
876:
877: case t_COMPLEX:
878: y=cgetg(3,t_COMPLEX); av=avma;
879: if (gcmp0((GEN)x[2]))
880: {
881: long tx=typ(x[1]);
882:
883: if ((is_intreal_t(tx) || is_frac_t(tx)) && gsigne((GEN)x[1]) < 0)
884: {
885: y[1]=zero; p1=gneg_i((GEN)x[1]); tetpil=avma;
886: y[2]=lpile(av,tetpil,gsqrt(p1,prec));
887: return y;
888: }
889: y[1]=lsqrt((GEN)x[1],prec);
890: y[2]=zero; return y;
891: }
892:
893: p1 = gsqr((GEN)x[1]);
894: p2 = gsqr((GEN)x[2]);
895: p1 = gsqrt(gadd(p1,p2),prec);
896: if (gcmp0(p1))
897: {
898: y[1]=lsqrt(p1,prec);
899: y[2]=lcopy((GEN)y[1]);
900: return y;
901: }
902:
903: if (gsigne((GEN)x[1]) < 0)
904: {
905: p1 = gmul2n(gsub(p1,(GEN)x[1]), -1);
906: y[2] = lsqrt(p1,prec);
907: y[1] = ldiv((GEN)x[2],gmul2n((GEN)y[2],1));
908: tetpil=avma;
909: y = (gsigne((GEN)x[2]) > 0)? gcopy(y): gneg(y);
910: return gerepile(av,tetpil,y);
911: }
912:
913: p1 = gmul2n(gadd(p1,(GEN)x[1]), -1); tetpil=avma;
914: y[1] = lpile(av,tetpil,gsqrt(p1,prec));
915: av=avma; p1=gmul2n((GEN)y[1],1); tetpil=avma;
916: y[2] = lpile(av,tetpil,gdiv((GEN)x[2],p1));
917: return y;
918:
919: case t_PADIC:
920: return padic_sqrt(x);
921:
922: case t_SER:
1.2 ! noro 923: e = valp(x);
1.1 noro 924: if (gcmp0(x)) return zeroser(varn(x), (e+1)>>1);
925: if (e & 1) err(sqrter6);
926: av = avma;
1.2 ! noro 927: y = dummycopy(x); setvalp(y, 0);
! 928: y = ser_pui(y, ghalf, prec);
1.1 noro 929: if (typ(y) == t_SER) /* generic case */
930: y[1] = evalsigne(1) | evalvalp(e>>1) | evalvarn(varn(x));
931: else /* e.g coeffs are POLMODs */
1.2 ! noro 932: y = gmul(y, gpowgs(polx[varn(x)], e>>1));
! 933: return gerepileupto(av, y);
1.1 noro 934: }
935: return transc(gsqrt,x,prec);
936: }
937:
938: void
939: gsqrtz(GEN x, GEN y)
940: {
1.2 ! noro 941: long prec = precision(y);
! 942: gpmem_t av=avma;
1.1 noro 943:
944: if (!prec) err(infprecer,"gsqrtz");
945: gaffect(gsqrt(x,prec),y); avma=av;
946: }
947: /********************************************************************/
948: /** **/
949: /** FONCTION RACINE N-IEME **/
950: /** **/
951: /********************************************************************/
952: GEN
953: rootsof1complex(GEN n, long prec)
954: {
1.2 ! noro 955: gpmem_t ltop = avma;
! 956: GEN z;
1.1 noro 957: if (is_pm1(n)) return realun(prec);
1.2 ! noro 958: if (lgefint(n)==3 && n[2]==2) return stor(-1, prec);
! 959:
! 960: z = exp_Ir( divri(Pi2n(1, prec), n) ); /* exp(2Ipi/n) */
1.1 noro 961: return gerepileupto(ltop,z);
962: }
963: /*Only the O() of y is used*/
964: GEN
965: rootsof1padic(GEN n, GEN y)
966: {
1.2 ! noro 967: gpmem_t ltop=avma;
1.1 noro 968: GEN z,r;
969: (void)mpsqrtnmod(gun,n,(GEN)y[2],&z);
970: if (z==gzero){avma=ltop;return gzero;}/*should not happen*/
971: r=cgetg(5,t_PADIC);
972: r[1]=y[1];setvalp(r,0);/*rootsofunity are unramified*/
973: r[2]=licopy((GEN)y[2]);
974: r[3]=licopy((GEN)y[3]);
975: r[4]=(long)padicsqrtnlift(gun,n,z,(GEN)y[2],precp(y));
976: return gerepileupto(ltop,r);
977: }
978: static GEN paexp(GEN x);
979: /*compute the p^e th root of x p-adic*/
1.2 ! noro 980: GEN
! 981: padic_sqrtn_ram(GEN x, long e)
1.1 noro 982: {
1.2 ! noro 983: gpmem_t ltop=avma;
1.1 noro 984: GEN n,a;
985: GEN p=(GEN)x[2];
986: long v=0;
987: n=gpowgs((GEN) x[2],e);
988: if (valp(x))
989: {
990: GEN p1,z;
991: p1=dvmdsi(valp(x),n,&z);
992: if (signe(z))
993: err(talker,"n-root does not exists in gsqrtn");
994: v=itos(p1);
995: x=gcopy(x);setvalp(x,0);
996: }
997: /*If p=2 -1 is an root of unity in U1,we need an extra check*/
998: if (lgefint(p)==3 && p[2]==2 && mod8((GEN)x[4])!=signe((GEN)x[4]))
999: err(talker,"n-root does not exists in gsqrtn");
1000: /*Other "n-root does not exists in gsqrtn" are caught by paexp...*/
1001: a=paexp(gdiv(palog(x),n));
1002: /*Here n=p^e and a^n=z*x where z is a root of unity. note that
1003: z^p=z, so z^n=z. and if b=a/z then b^n=x. We say b=x/(a^(n-1))*/
1004: a=gdiv(x,powgi(a,addis(n,-1)));
1005: if (v) { a=gcopy(a);setvalp(a,v); }
1006: a=gerepileupto(ltop,a);
1007: return a;
1008: }
1009: /*compute the nth root of x p-adic p prime with n*/
1.2 ! noro 1010: GEN
! 1011: padic_sqrtn_unram(GEN x, GEN n, GEN *zetan)
1.1 noro 1012: {
1.2 ! noro 1013: gpmem_t ltop=avma, tetpil;
1.1 noro 1014: GEN a,r;
1015: GEN p=(GEN)x[2];
1016: long v=0;
1017: /*check valuation*/
1018: if (valp(x))
1019: {
1020: GEN p1,z;
1021: p1=dvmdsi(valp(x),n,&z);
1022: if (signe(z))
1023: err(talker,"n-root does not exists in gsqrtn");
1024: v=itos(p1);
1025: }
1026: a=mpsqrtnmod((GEN)x[4],n,p,zetan);
1027: if (!a)
1028: err(talker,"n-root does not exists in gsqrtn");
1029: tetpil=avma;
1030: r=cgetg(5,t_PADIC);
1031: r[1]=x[1];setvalp(r,v);
1032: r[2]=licopy(p);
1033: r[3]=licopy((GEN)x[3]);
1034: r[4]=(long)padicsqrtnlift((GEN)x[4],n,a,p,precp(x));
1035: if (zetan)
1036: {
1037: GEN z,*gptr[2];
1038: z=cgetg(5,t_PADIC);
1039: z[1]=x[1];setvalp(z,0);
1040: z[2]=licopy(p);
1041: z[3]=licopy((GEN)x[3]);
1042: z[4]=(long)padicsqrtnlift(gun,n,*zetan,p,precp(x));
1043: gptr[0]=&r;gptr[1]=&z;
1044: gerepilemanysp(ltop,tetpil,gptr,2);
1045: *zetan=z;
1046: return r;
1047: }
1048: else
1049: return gerepile(ltop,tetpil,r);
1050: }
1.2 ! noro 1051:
! 1052: GEN
! 1053: padic_sqrtn(GEN x, GEN n, GEN *zetan)
1.1 noro 1054: {
1.2 ! noro 1055: gpmem_t ltop=avma, tetpil;
1.1 noro 1056: GEN p=(GEN)x[2];
1057: GEN q;
1058: long e;
1059: GEN *gptr[2];
1060: if (gcmp0(x))
1061: {
1062: long m = itos(n);
1.2 ! noro 1063: return padiczero(p, (valp(x)+m-1)/m);
1.1 noro 1064: }
1065: /*First treat the ramified part using logarithms*/
1066: e=pvaluation(n,p,&q);
1067: tetpil=avma;
1068: if (e)
1069: {
1070: x=padic_sqrtn_ram(x,e);
1071: }
1072: /*finished ?*/
1073: if (is_pm1(q))
1074: {
1075: if (signe(q)<0)
1076: {
1077: tetpil=avma;
1078: x=ginv(x);
1079: }
1080: if (zetan && e && lgefint(p)==3 && p[2]==2)/*-1 in Q_2*/
1081: {
1082: *zetan=negi(gun);
1083: gptr[0]=&x;gptr[1]=zetan;
1084: gerepilemanysp(ltop,tetpil,gptr,2);
1085: return x;
1086: }
1087: if (zetan) *zetan=gun;
1088: return gerepile(ltop,tetpil,x);
1089: }
1090: /*Now we use hensel lift for unramified case. 4x faster.*/
1091: tetpil=avma;
1092: x=padic_sqrtn_unram(x,q,zetan);
1093: if (zetan)
1094: {
1095: if (e && lgefint(p)==3 && p[2]==2)/*-1 in Q_2*/
1096: {
1097: tetpil=avma;
1098: x=gcopy(x);
1099: *zetan=gneg(*zetan);
1100: }
1101: gptr[0]=&x;gptr[1]=zetan;
1102: gerepilemanysp(ltop,tetpil,gptr,2);
1103: return x;
1104: }
1105: return gerepile(ltop,tetpil,x);
1106: }
1107:
1108: GEN
1109: gsqrtn(GEN x, GEN n, GEN *zetan, long prec)
1110: {
1.2 ! noro 1111: long i, lx, tx;
! 1112: gpmem_t av, tetpil;
1.1 noro 1113: long e,m;
1114: GEN y,z;
1115: if (zetan) *zetan=gzero;
1116: if (typ(n)!=t_INT) err(talker,"second arg must be integer in gsqrtn");
1117: if (!signe(n)) err(talker,"1/0 exponent in gsqrtn");
1118: if (is_pm1(n))
1119: {
1120: if (zetan) *zetan=gun;
1121: if (signe(n)>0)
1122: return gcopy(x);
1123: return ginv(x);
1124: }
1125: tx = typ(x);
1126: if (is_matvec_t(tx))
1127: {
1128: lx=lg(x); y=cgetg(lx,tx);
1129: for (i=1; i<lx; i++) y[i]=(long)gsqrtn((GEN)x[i],n,NULL,prec);
1130: return y;
1131: }
1132: switch(tx)
1133: {
1134: case t_POL: case t_RFRAC: case t_RFRACN:
1135: av = avma;
1136: y=tayl(x,gvar(x),precdl); tetpil=avma;
1137: return gerepile(av,tetpil,gsqrtn(y,n,zetan,prec));
1138: case t_SER:
1139: e=valp(x);m=itos(n);
1140: if (gcmp0(x)) return zeroser(varn(x), (e+m-1)/m);
1.2 ! noro 1141: if (e % m) err(talker,"incorrect valuation in gsqrtn");
1.1 noro 1142: av = avma;
1.2 ! noro 1143: y = dummycopy(x); setvalp(y, 0);
! 1144: y = ser_pui(y, ginv(n), prec);
1.1 noro 1145: if (typ(y) == t_SER) /* generic case */
1146: y[1] = evalsigne(1) | evalvalp(e/m) | evalvarn(varn(x));
1147: else /* e.g coeffs are POLMODs */
1.2 ! noro 1148: y = gmul(y, gpowgs(polx[varn(x)], e/m));
! 1149: return gerepileupto(av, y);
1.1 noro 1150: case t_INTMOD:
1151: z=gzero;
1152: /*This is not great, but else it will generate too much trouble*/
1.2 ! noro 1153: if (!BSW_psp((GEN)x[1])) err(talker,"modulus must be prime in gsqrtn");
1.1 noro 1154: if (zetan)
1155: {
1156: z=cgetg(3,tx); copyifstack(x[1],z[1]);
1157: z[2]=lgetg(lgefint(z[1]),t_INT);
1158: }
1159: y=cgetg(3,tx); copyifstack(x[1],y[1]);
1160: y[2]=(long)mpsqrtnmod((GEN)x[2],n,(GEN)x[1],zetan);
1161: if (zetan)
1162: {
1163: cgiv(*zetan);/*Ole!*/
1164: affii(*zetan,(GEN)z[2]);
1165: *zetan=z;
1166: }
1167: if(!y[2]) err(talker,"n-root does not exists in gsqrtn");
1168: return y;
1169:
1170: case t_PADIC:
1171: return padic_sqrtn(x,n,zetan);
1172: case t_INT: case t_FRAC: case t_FRACN: case t_REAL: case t_COMPLEX:
1173: i = (long) precision(n); if (i) prec=i;
1174: if (tx==t_INT && is_pm1(x) && signe(x)>0)
1175: y=gun; /*speed-up since there is no way to call rootsof1complex
1176: directly from gp*/
1.2 ! noro 1177: else if (gcmp0(x))
! 1178: {
! 1179: if (gsigne(n) < 0) err(gdiver2);
! 1180: if (isinexactreal(x))
! 1181: y = realzero_bit( itos( gfloor(gdivsg(gexpo(x), n)) ) );
! 1182: else
! 1183: y = realzero(prec);
! 1184: }
1.1 noro 1185: else
1186: {
1187: av=avma;
1188: y=gmul(ginv(n),glog(x,prec)); tetpil=avma;
1189: y=gerepile(av,tetpil,gexp(y,prec));
1190: }
1191: if (zetan)
1192: *zetan=rootsof1complex(n,prec);
1193: return y;
1194: default:
1195: err(typeer,"gsqrtn");
1196: }
1197: return NULL;/*keep GCC happy*/
1198: }
1199:
1200: /********************************************************************/
1201: /** **/
1202: /** FONCTION EXPONENTIELLE-1 **/
1203: /** **/
1204: /********************************************************************/
1205: #ifdef LONG_IS_64BIT
1206: # define EXMAX 46
1207: #else
1208: # define EXMAX 22
1209: #endif
1210:
1211: GEN
1212: mpexp1(GEN x)
1213: {
1.2 ! noro 1214: long l, l1, l2, i, n, m, ex, s, sx = signe(x);
! 1215: gpmem_t av, av2;
1.1 noro 1216: double a,b,alpha,beta, gama = 2.0 /* optimized for SUN3 */;
1217: /* KB: 3.0 is better for UltraSparc */
1218: GEN y,p1,p2,p3,p4,unr;
1219:
1220: if (typ(x)!=t_REAL) err(typeer,"mpexp1");
1.2 ! noro 1221: if (!sx) return realzero_bit(expo(x));
1.1 noro 1222: l=lg(x); y=cgetr(l); av=avma; /* room for result */
1223:
1224: l2 = l+1; ex = expo(x);
1225: if (ex > EXMAX) err(talker,"exponent too large in exp");
1226: alpha = -1-log(2+x[2]/C31)-ex*LOG2;
1227: beta = 5 + bit_accuracy(l)*LOG2;
1228: a = sqrt(beta/(gama*LOG2));
1229: b = (alpha + 0.5*log(beta*gama/LOG2))/LOG2;
1230: if (a>=b)
1231: {
1232: n=(long)(1+sqrt(beta*gama/LOG2));
1233: m=(long)(1+a-b);
1234: l2 += m>>TWOPOTBITS_IN_LONG;
1235: }
1236: else
1237: {
1238: n=(long)(1+beta/alpha);
1239: m=0;
1240: }
1241: unr=realun(l2);
1242: p2=rcopy(unr); setlg(p2,4);
1243: p4=cgetr(l2); affrr(x,p4); setsigne(p4,1);
1244: if (m) setexpo(p4,ex-m);
1245:
1246: s=0; l1=4; av2 = avma;
1247: for (i=n; i>=2; i--)
1248: {
1249: setlg(p4,l1); p3 = divrs(p4,i);
1250: s -= expo(p3); p1 = mulrr(p3,p2); setlg(p1,l1);
1251: l1 += s>>TWOPOTBITS_IN_LONG; if (l1>l2) l1=l2;
1252: s %= BITS_IN_LONG;
1253: setlg(unr,l1); p1 = addrr(unr,p1);
1254: setlg(p2,l1); affrr(p1,p2); avma = av2;
1255: }
1256: setlg(p2,l2); setlg(p4,l2); p2 = mulrr(p4,p2);
1257:
1258: for (i=1; i<=m; i++)
1259: {
1260: setlg(p2,l2);
1261: p2 = mulrr(addsr(2,p2),p2);
1262: }
1263: if (sx == -1)
1264: {
1265: setlg(unr,l2); p2 = addrr(unr,p2);
1266: setlg(p2,l2); p2 = ginv(p2);
1267: p2 = subrr(p2,unr);
1268: }
1269: affrr(p2,y); avma = av; return y;
1270: }
1271: #undef EXMAX
1272:
1273: /********************************************************************/
1274: /** **/
1275: /** FONCTION EXPONENTIELLE **/
1276: /** **/
1277: /********************************************************************/
1278:
1279: GEN
1280: mpexp(GEN x)
1281: {
1.2 ! noro 1282: long sx = signe(x);
! 1283: gpmem_t av;
1.1 noro 1284: GEN y;
1285:
1286: if (!sx) return addsr(1,x);
1287: if (sx<0) setsigne(x,1);
1288: av = avma; y = addsr(1, mpexp1(x));
1.2 ! noro 1289: if (sx<0) { y = ginv(y); setsigne(x,-1); }
1.1 noro 1290: return gerepileupto(av,y);
1291: }
1292:
1293: static GEN
1294: paexp(GEN x)
1295: {
1.2 ! noro 1296: long k, e = valp(x), pp = precp(x), n = e + pp;
! 1297: gpmem_t av;
1.1 noro 1298: GEN y,r,p1;
1299:
1300: if (gcmp0(x)) return gaddgs(x,1);
1301: if (e<=0 || (!cmpis((GEN)x[2],2) && e==1))
1302: err(talker,"p-adic argument out of range in gexp");
1303: av = avma;
1304: if (egalii(gdeux,(GEN)x[2]))
1305: {
1306: n--; e--; k = n/e;
1307: if (n%e==0) k--;
1308: }
1309: else
1310: {
1311: p1 = subis((GEN)x[2], 1);
1312: k = itos(dvmdii(subis(mulis(p1,n), 1), subis(mulis(p1,e), 1), &r));
1313: if (!signe(r)) k--;
1314: }
1315: for (y=gun; k; k--) y = gaddsg(1, gdivgs(gmul(y,x), k));
1316: return gerepileupto(av, y);
1317: }
1318:
1.2 ! noro 1319: static GEN
! 1320: cxexp(GEN x, long prec)
! 1321: {
! 1322: GEN r,p1,p2, y = cgetg(3,t_COMPLEX);
! 1323: gpmem_t av = avma, tetpil;
! 1324: r=gexp((GEN)x[1],prec);
! 1325: gsincos((GEN)x[2],&p2,&p1,prec);
! 1326: tetpil = avma;
! 1327: y[1] = lmul(r,p1);
! 1328: y[2] = lmul(r,p2);
! 1329: gerepilemanyvec(av,tetpil,y+1,2);
! 1330: return y;
! 1331: }
! 1332:
! 1333: static GEN
! 1334: serexp(GEN x, long prec)
! 1335: {
! 1336: gpmem_t av;
! 1337: long i,j,lx,ly,ex,mi;
! 1338: GEN p1,y,xd,yd;
! 1339:
! 1340: ex = valp(x);
! 1341: if (ex < 0) err(negexper,"gexp");
! 1342: if (gcmp0(x)) return gaddsg(1,x);
! 1343: lx = lg(x);
! 1344: if (ex)
! 1345: {
! 1346: ly = lx+ex; y = cgetg(ly,t_SER);
! 1347: mi = lx-1; while (mi>=3 && gcmp0((GEN)x[mi])) mi--;
! 1348: mi += ex-2;
! 1349: y[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(x));
! 1350: /* zd[i] = coeff of X^i in z */
! 1351: xd = x+2-ex; yd = y+2; ly -= 2;
! 1352: yd[0] = un;
! 1353: for (i=1; i<ex; i++) yd[i]=zero;
! 1354: for ( ; i<ly; i++)
! 1355: {
! 1356: av = avma; p1 = gzero;
! 1357: for (j=ex; j<=min(i,mi); j++)
! 1358: p1 = gadd(p1, gmulgs(gmul((GEN)xd[j],(GEN)yd[i-j]),j));
! 1359: yd[i] = lpileupto(av, gdivgs(p1,i));
! 1360: }
! 1361: return y;
! 1362: }
! 1363: av = avma; y = cgetg(lx, t_SER);
! 1364: y[1] = x[1]; y[2] = zero;
! 1365: for (i=3; i <lx; i++) y[i] = x[i];
! 1366: p1 = gexp((GEN)x[2],prec);
! 1367: y = gmul(p1, serexp(normalize(y),prec));
! 1368: return gerepileupto(av, y);
! 1369: }
! 1370:
1.1 noro 1371: GEN
1372: gexp(GEN x, long prec)
1373: {
1374: switch(typ(x))
1375: {
1.2 ! noro 1376: case t_REAL: return mpexp(x);
! 1377: case t_COMPLEX: return cxexp(x,prec);
! 1378: case t_PADIC: return paexp(x);
! 1379: case t_SER: return serexp(x,prec);
! 1380: case t_INTMOD: err(typeer,"gexp");
1.1 noro 1381: }
1382: return transc(gexp,x,prec);
1383: }
1384:
1385: void
1386: gexpz(GEN x, GEN y)
1387: {
1.2 ! noro 1388: long prec = precision(y);
! 1389: gpmem_t av=avma;
1.1 noro 1390:
1391: if (!prec) err(infprecer,"gexpz");
1392: gaffect(gexp(x,prec),y); avma=av;
1393: }
1394:
1395: /********************************************************************/
1396: /** **/
1397: /** FONCTION LOGARITHME **/
1398: /** **/
1399: /********************************************************************/
1.2 ! noro 1400: /* 2 * atanh(1/3) */
! 1401: GEN
! 1402: constlog2(long prec)
! 1403: {
! 1404: static GEN glog2 = NULL;
! 1405: const long _3 = 3, _9 = _3*_3;
! 1406: gpmem_t av0, av;
! 1407: long k, l, G;
! 1408: GEN s, u, S, U, tmplog2;
! 1409:
! 1410: if (glog2 && lg(glog2) >= prec) return glog2;
! 1411:
! 1412: tmplog2 = newbloc(prec);
! 1413: *tmplog2 = evaltyp(t_REAL) | evallg(prec);
! 1414: av0 = avma;
! 1415: l = prec+1; G = bit_accuracy(l+1);
! 1416:
! 1417: s = S = divrs(realun(l), _3);
! 1418: u = U = mpcopy(s);
! 1419: av = avma;
! 1420: for (k = 3; ; k += 2)
! 1421: {
! 1422: u = divrs(u, _9);
! 1423: if (bit_accuracy(l) - expo(u) > G) {
! 1424: l--; if (l <= 2) break;
! 1425: setlg(U,l);
! 1426: affrr(s,S); s = S;
! 1427: affrr(u,U); u = U; avma = av;
! 1428: }
! 1429: s = addrr(s, divrs(u,k));
! 1430: }
! 1431: setexpo(s, -1); affrr(s, tmplog2);
! 1432: gunclone(glog2); glog2 = tmplog2;
! 1433: avma = av0; return glog2;
! 1434: }
! 1435:
! 1436: GEN
! 1437: mplog2(long prec)
! 1438: {
! 1439: GEN x = cgetr(prec);
! 1440: affrr(constlog2(prec), x); return x;
! 1441: }
1.1 noro 1442:
1443: GEN
1444: mplog(GEN x)
1445: {
1.2 ! noro 1446: gpmem_t ltop, av;
! 1447: long EX,l,l1,l2,m,n,k,ex,s;
! 1448: double alpha,a,b;
! 1449: GEN z,p1,y,y2,p4,p5,unr;
1.1 noro 1450:
1451: if (typ(x)!=t_REAL) err(typeer,"mplog");
1452: if (signe(x)<=0) err(talker,"non positive argument in mplog");
1453:
1.2 ! noro 1454: av = avma;
! 1455: l = lg(x); EX = expo(x);
! 1456: z = cgetr(l); ltop = avma;
1.1 noro 1457:
1.2 ! noro 1458: l2 = l+1; y=p1=cgetr(l2); affrr(x,p1);
! 1459: setexpo(p1, 0);
! 1460: if (gcmp1(p1)) {
! 1461: if (!EX) { avma = av; return realzero(l); }
! 1462: affrr(mulsr(EX, mplog2(l)), z);
! 1463: avma = ltop; return z;
! 1464: }
! 1465: /* 1 < p1 < 2 */
! 1466: av = avma; l -= 2;
1.1 noro 1467: alpha = 1+p1[2]/C31; if (!alpha) alpha = 0.00000001;
1.2 ! noro 1468: a = -log(alpha)/LOG2;
! 1469: b = sqrt(BITS_IN_HALFULONG*l/3.0);
! 1470: if (a <= b)
1.1 noro 1471: {
1.2 ! noro 1472: n = 1 + (long)sqrt(BITS_IN_HALFULONG*3.0*l);
! 1473: m = 1 + (long)(b-a);
1.1 noro 1474: l2 += m>>TWOPOTBITS_IN_LONG;
1.2 ! noro 1475: p4 = cgetr(l2); affrr(p1,p4);
! 1476: p1 = p4; av = avma;
! 1477: for (k=1; k<=m; k++) p1 = mpsqrt(p1);
! 1478: affrr(p1,p4); avma = av;
1.1 noro 1479: }
1480: else
1481: {
1.2 ! noro 1482: n = 1 + (long)(BITS_IN_HALFULONG*l / a);
! 1483: m = 0; p4 = p1;
1.1 noro 1484: }
1.2 ! noro 1485: unr = realun(l2);
! 1486: y = cgetr(l2);
! 1487: y2 = cgetr(l2); av = avma;
! 1488:
! 1489: /* want to compute log(X), X ~ 1 (X = p4) */
! 1490: /* y = (X-1)/(X+1). log(X) = log(1+y) - log(1-y) = 2 \sum_{k odd} y^k / * k */
1.1 noro 1491:
1492: /* affrr needed here instead of setlg since prec may DECREASE */
1493: p1 = cgetr(l2); affrr(subrr(p4,unr), p1);
1494:
1495: p5 = addrr(p4,unr); setlg(p5,l2);
1.2 ! noro 1496: affrr(divrr(p1,p5), y); /* = (X-1) / (X+1) ~ 0 */
! 1497: affrr(gsqr(y), y2); /* = y^2 */
! 1498: k = 2*n + 1;
! 1499: affrr(divrs(unr,k), p4); setlg(p4,4); avma = av;
! 1500:
! 1501: s = 0; ex = expo(y2); l1 = 4;
! 1502: for (k -= 2; k > 0; k -= 2)
! 1503: { /* compute sum_i=0^n y^2i / (2i + 1), k = 2i+1 */
! 1504: setlg(y2, l1); p5 = mulrr(p4,y2);
! 1505: setlg(unr,l1); p1 = divrs(unr, k);
1.1 noro 1506: s -= ex;
1507: l1 += s>>TWOPOTBITS_IN_LONG; if (l1>l2) l1=l2;
1.2 ! noro 1508: s &= (BITS_IN_LONG-1);
! 1509: setlg(p1, l1);
! 1510: setlg(p4, l1);
! 1511: setlg(p5, l1); affrr(addrr(p1,p5), p4); avma=av;
! 1512: }
! 1513: setlg(p4, l2);
! 1514: y = mulrr(y,p4); /* = log(X)/2 */
! 1515: setexpo(y, expo(y)+m+1);
! 1516: if (EX) y = addrr(y, mulsr(EX, mplog2(l2)));
! 1517: affrr(y, z); avma = ltop; return z;
1.1 noro 1518: }
1519:
1520: GEN
1521: teich(GEN x)
1522: {
1.2 ! noro 1523: GEN p,q,y,z,aux,p1;
! 1524: long n, k;
! 1525: gpmem_t av;
1.1 noro 1526:
1527: if (typ(x)!=t_PADIC) err(talker,"not a p-adic argument in teichmuller");
1528: if (!signe(x[4])) return gcopy(x);
1.2 ! noro 1529: p = (GEN)x[2];
! 1530: q = (GEN)x[3];
! 1531: z = (GEN)x[4];
! 1532: y = cgetp(x); av = avma;
! 1533: if (egalii(p, gdeux))
! 1534: z = (mod4(z) & 2)? addsi(-1,q): gun;
! 1535: else
1.1 noro 1536: {
1.2 ! noro 1537: p1 = addsi(-1, p);
! 1538: z = resii(z, p);
! 1539: aux = divii(addsi(-1,q),p1); n = precp(x);
! 1540: for (k=1; k<n; k<<=1)
! 1541: z = modii(mulii(z,addsi(1,mulii(aux,addsi(-1,powmodulo(z,p1,q))))), q);
! 1542: }
! 1543: affii(z, (GEN)y[4]); avma = av; return y;
1.1 noro 1544: }
1545:
1.2 ! noro 1546: /* Let x = 1 mod p and y := (x-1)/(x+1) = 0 (p). Then
! 1547: * log(x) = log(1+y) - log(1-y) = 2 \sum_{k odd} y^k / k.
! 1548: * palogaux(x) returns the last sum (not multiplied by 2) */
1.1 noro 1549: static GEN
1550: palogaux(GEN x)
1551: {
1.2 ! noro 1552: long k,e,pp;
! 1553: GEN y,s,y2, p = (GEN)x[2];
1.1 noro 1554:
1.2 ! noro 1555: if (egalii(gun, (GEN)x[4]))
1.1 noro 1556: {
1.2 ! noro 1557: long v = valp(x)+precp(x);
! 1558: if (egalii(gdeux,p)) v--;
! 1559: return padiczero(p, v);
! 1560: }
! 1561: y = gdiv(gaddgs(x,-1), gaddgs(x,1));
! 1562: e = valp(y); pp = e+precp(y);
! 1563: if (egalii(gdeux,p)) pp--;
1.1 noro 1564: else
1565: {
1566: GEN p1;
1.2 ! noro 1567: for (p1=stoi(e); cmpsi(pp,p1)>0; pp++) p1 = mulii(p1, p);
! 1568: pp -= 2;
1.1 noro 1569: }
1.2 ! noro 1570: k = pp/e; if (!odd(k)) k--;
! 1571: y2 = gsqr(y); s = gdivgs(gun,k);
! 1572: while (k > 2)
1.1 noro 1573: {
1.2 ! noro 1574: k -= 2; s = gadd(gmul(y2,s), gdivgs(gun,k));
1.1 noro 1575: }
1.2 ! noro 1576: return gmul(s,y);
1.1 noro 1577: }
1578:
1579: GEN
1580: palog(GEN x)
1581: {
1.2 ! noro 1582: gpmem_t av = avma;
! 1583: GEN y, p = (GEN)x[2];
1.1 noro 1584:
1585: if (!signe(x[4])) err(talker,"zero argument in palog");
1.2 ! noro 1586: if (egalii(p, gdeux))
1.1 noro 1587: {
1.2 ! noro 1588: y = gsqr(x); setvalp(y,0);
! 1589: y = palogaux(y);
1.1 noro 1590: }
1.2 ! noro 1591: else
! 1592: { /* compute log(x^(p-1)) / (p-1) */
! 1593: GEN mod = (GEN)x[3], p1 = subis(p,1);
! 1594: y = cgetp(x);
! 1595: y[4] = (long)powmodulo((GEN)x[4], p1, mod);
! 1596: p1 = divii(subis(mod,1), p1); /* 1/(1-p) */
! 1597: y = gmul(palogaux(y), mulis(p1, -2));
! 1598: }
! 1599: return gerepileupto(av,y);
1.1 noro 1600: }
1601:
1602: GEN
1603: log0(GEN x, long flag,long prec)
1604: {
1605: switch(flag)
1606: {
1607: case 0: return glog(x,prec);
1608: case 1: return glogagm(x,prec);
1609: default: err(flagerr,"log");
1610: }
1611: return NULL; /* not reached */
1612: }
1613:
1614: GEN
1615: glog(GEN x, long prec)
1616: {
1.2 ! noro 1617: gpmem_t av, tetpil;
1.1 noro 1618: GEN y,p1,p2;
1619:
1620: switch(typ(x))
1621: {
1622: case t_REAL:
1623: if (signe(x)>=0) return mplog(x);
1624: y=cgetg(3,t_COMPLEX); y[2]=lmppi(lg(x));
1625: setsigne(x,1); y[1]=lmplog(x);
1626: setsigne(x,-1); return y;
1627:
1628: case t_COMPLEX:
1629: y=cgetg(3,t_COMPLEX); y[2]=larg(x,prec);
1630: av=avma; p1=glog(gnorm(x),prec); tetpil=avma;
1631: y[1]=lpile(av,tetpil,gmul2n(p1,-1));
1632: return y;
1633:
1634: case t_PADIC:
1635: return palog(x);
1636:
1637: case t_SER:
1.2 ! noro 1638: av=avma;
! 1639: if (valp(x) || gcmp0(x)) err(talker,"log is not analytic at 0");
1.1 noro 1640: p1=gdiv(derivser(x),x);
1641: tetpil=avma; p1=integ(p1,varn(x));
1642: if (gcmp1((GEN)x[2])) return gerepile(av,tetpil,p1);
1643:
1644: p2=glog((GEN)x[2],prec); tetpil=avma;
1645: return gerepile(av,tetpil,gadd(p1,p2));
1646:
1647: case t_INTMOD:
1648: err(typeer,"glog");
1649: }
1650: return transc(glog,x,prec);
1651: }
1652:
1653: void
1654: glogz(GEN x, GEN y)
1655: {
1.2 ! noro 1656: long prec = precision(y);
! 1657: gpmem_t av=avma;
1.1 noro 1658:
1659: if (!prec) err(infprecer,"glogz");
1660: gaffect(glog(x,prec),y); avma=av;
1661: }
1662:
1663: /********************************************************************/
1664: /** **/
1.2 ! noro 1665: /** SINE, COSINE **/
1.1 noro 1666: /** **/
1667: /********************************************************************/
1668:
1.2 ! noro 1669: /* Reduce x0 mod Pi/2 to x in [-Pi/4, Pi/4]. Return cos(x)-1 */
1.1 noro 1670: static GEN
1.2 ! noro 1671: mpsc1(GEN x0, long *ptmod8)
1.1 noro 1672: {
1.2 ! noro 1673: const long mmax = 23169; /* largest m such that (2m+2)(2m+1) < 2^31 */
1.1 noro 1674: /* on a 64-bit machine with true 128 bit/64 bit division, one could
1.2 ! noro 1675: * take mmax=1518500248; on the alpha it does not seem worthwhile */
! 1676: long l, l0, l1, l2, l4, i, n, m, s, t;
! 1677: gpmem_t av;
1.1 noro 1678: double alpha,beta,a,b,c,d;
1.2 ! noro 1679: GEN y,p1,p2,p3,x2,pitemp, x = x0;
! 1680:
! 1681: if (typ(x) != t_REAL) err(typeer,"mpsc1");
! 1682: l = lg(x); y = cgetr(l); av = avma;
! 1683:
! 1684: l++; pitemp = mppi(l);
! 1685: setexpo(pitemp,-1);
! 1686: p1 = addrr(x,pitemp); /* = x + Pi/4 */
! 1687: l0 = min(l, lg(p1));
! 1688: if (expo(p1) >= bit_accuracy(l0) + 3) err(precer,"mpsc1");
1.1 noro 1689:
1.2 ! noro 1690: setexpo(pitemp, 0);
! 1691: p1 = mpent( divrr(p1,pitemp) ); /* round ( x / (Pi/2) ) */
! 1692: if (signe(p1)) x = subrr(x, mulir(p1,pitemp)); /* x mod Pi/2 */
! 1693: p2 = cgetr(l); affrr(x, p2); x = p2;
! 1694:
! 1695: *ptmod8 = (signe(x) < 0)? 4: 0;
! 1696: if (signe(p1))
1.1 noro 1697: {
1.2 ! noro 1698: long k = mod4(p1);
! 1699: if (signe(p1) < 0 && k) k = 4-k;
! 1700: /* x = x0 - k*Pi/2 (mod 2Pi) */
! 1701: *ptmod8 += k;
1.1 noro 1702: }
1.2 ! noro 1703:
! 1704: if (gcmp0(x)) alpha = 1000000.0;
1.1 noro 1705: else
1706: {
1.2 ! noro 1707: long e = expo(x);
! 1708: alpha = -1 - ((e<-1022)? e*LOG2: log(fabs(rtodbl(x))));
1.1 noro 1709: }
1710: beta = 5 + bit_accuracy(l)*LOG2;
1.2 ! noro 1711: a = 0.5 / LOG2;
! 1712: b = 0.5 * a;
! 1713: c = a + sqrt((beta+b) / LOG2);
! 1714: d = ((beta/c) - alpha - log(c)) / LOG2;
! 1715: if (d >= 0)
! 1716: {
! 1717: m = (long)(1+d);
! 1718: n = (long)((1+c) / 2.0);
! 1719: setexpo(x, expo(x)-m);
1.1 noro 1720: }
1.2 ! noro 1721: else { m = 0; n = (long)((2+beta/alpha) / 2.0); }
1.1 noro 1722: l2=l+1+(m>>TWOPOTBITS_IN_LONG);
1.2 ! noro 1723: p2 = realun(l2);
! 1724: x2 = cgetr(l2); av = avma;
! 1725: affrr(gsqr(x), x2);
! 1726:
! 1727: setlg(x2, 3);
! 1728: if (n > mmax)
! 1729: p3 = divrs(divrs(x2, 2*n+2), 2*n+1);
1.1 noro 1730: else
1.2 ! noro 1731: p3 = divrs(x2, (2*n+2)*(2*n+1));
! 1732: s = -expo(p3);
! 1733: l4 = l1 = 3 + (s>>TWOPOTBITS_IN_LONG);
! 1734: if (l4<=l2) { setlg(p2,l4); setlg(x2,l4); }
! 1735: s = 0;
! 1736: for (i=n; i>=2; i--)
1.1 noro 1737: {
1.2 ! noro 1738: if (i > mmax)
! 1739: p3 = divrs(divrs(x2, 2*i), 2*i-1);
! 1740: else
! 1741: p3 = divrs(x2, 2*i*(2*i-1));
! 1742: s -= expo(p3);
! 1743: t = s & (BITS_IN_LONG-1); l0 = (s>>TWOPOTBITS_IN_LONG);
1.1 noro 1744: if (t) l0++;
1745: l1 += l0; if (l1>l2) { l0 += l2-l1; l1=l2; }
1746: l4 += l0;
1747: p3 = mulrr(p3,p2);
1.2 ! noro 1748: if (l4<=l2) { setlg(p2,l4); setlg(x2,l4); }
! 1749: subsrz(1,p3, p2); avma = av;
1.1 noro 1750: }
1.2 ! noro 1751: setlg(p2,l2); setlg(x2,l2);
! 1752: setexpo(p2, expo(p2)-1); setsigne(p2, -signe(p2));
! 1753: p2 = mulrr(x2,p2);
! 1754: /* Now p2 = sum {1<= i <=n} (-1)^i x^(2i) / (2i)! ~ cos(x) - 1 */
1.1 noro 1755: for (i=1; i<=m; i++)
1.2 ! noro 1756: { /* p2 = cos(x)-1 --> cos(2x)-1 */
! 1757: p2 = mulrr(p2, addsr(2,p2)); setexpo(p2, expo(p2)+1);
1.1 noro 1758: }
1759: affrr(p2,y); avma=av; return y;
1760: }
1761:
1.2 ! noro 1762: /* sqrt (1 - (1+x)^2) = sqrt(-x*(x+2)). Sends cos(x)-1 to |sin(x)| */
1.1 noro 1763: static GEN
1764: mpaut(GEN x)
1765: {
1.2 ! noro 1766: gpmem_t av = avma;
! 1767: GEN p1 = mulrr(x, addsr(2,x));
! 1768: setsigne(p1, -signe(p1));
1.1 noro 1769: return gerepileuptoleaf(av, mpsqrt(p1));
1770: }
1771:
1772: /********************************************************************/
1.2 ! noro 1773: /** COSINE **/
1.1 noro 1774: /********************************************************************/
1775:
1776: static GEN
1777: mpcos(GEN x)
1778: {
1.2 ! noro 1779: long mod8;
! 1780: gpmem_t av;
1.1 noro 1781: GEN y,p1;
1782:
1783: if (typ(x)!=t_REAL) err(typeer,"mpcos");
1784: if (!signe(x)) return addsr(1,x);
1785:
1.2 ! noro 1786: av = avma; p1 = mpsc1(x,&mod8);
1.1 noro 1787: switch(mod8)
1788: {
1789: case 0: case 4:
1790: y=addsr(1,p1); break;
1791: case 1: case 7:
1792: y=mpaut(p1); setsigne(y,-signe(y)); break;
1793: case 2: case 6:
1794: y=subsr(-1,p1); break;
1795: default: /* case 3: case 5: */
1796: y=mpaut(p1); break;
1797: }
1.2 ! noro 1798: return gerepileuptoleaf(av, y);
1.1 noro 1799: }
1800:
1801: GEN
1802: gcos(GEN x, long prec)
1803: {
1.2 ! noro 1804: gpmem_t av, tetpil;
1.1 noro 1805: GEN r,u,v,y,p1,p2;
1806:
1807: switch(typ(x))
1808: {
1809: case t_REAL:
1810: return mpcos(x);
1811:
1812: case t_COMPLEX:
1813: y=cgetg(3,t_COMPLEX); av=avma;
1814: r=gexp((GEN)x[2],prec); p1=ginv(r);
1815: p2=gmul2n(gadd(p1,r),-1);
1816: p1=gsub(p2,r);
1817: gsincos((GEN)x[1],&u,&v,prec);
1818: tetpil=avma;
1819: y[1]=lmul(p2,v); y[2]=lmul(p1,u);
1820: gerepilemanyvec(av,tetpil,y+1,2);
1821: return y;
1822:
1823: case t_SER:
1824: if (gcmp0(x)) return gaddsg(1,x);
1825: if (valp(x)<0) err(negexper,"gcos");
1826:
1827: av=avma; gsincos(x,&u,&v,prec);
1828: return gerepilecopy(av,v);
1829:
1830: case t_INTMOD: case t_PADIC:
1831: err(typeer,"gcos");
1832: }
1833: return transc(gcos,x,prec);
1834: }
1835:
1836: void
1837: gcosz(GEN x, GEN y)
1838: {
1.2 ! noro 1839: long prec = precision(y);
! 1840: gpmem_t av = avma;
1.1 noro 1841:
1842: if (!prec) err(infprecer,"gcosz");
1843: gaffect(gcos(x,prec),y); avma=av;
1844: }
1845:
1846: /********************************************************************/
1.2 ! noro 1847: /** SINE **/
1.1 noro 1848: /********************************************************************/
1849:
1850: GEN
1851: mpsin(GEN x)
1852: {
1.2 ! noro 1853: long mod8;
! 1854: gpmem_t av;
1.1 noro 1855: GEN y,p1;
1856:
1857: if (typ(x)!=t_REAL) err(typeer,"mpsin");
1.2 ! noro 1858: if (!signe(x)) return realzero_bit(expo(x));
1.1 noro 1859:
1.2 ! noro 1860: av = avma; p1 = mpsc1(x,&mod8);
1.1 noro 1861: switch(mod8)
1862: {
1863: case 0: case 6:
1864: y=mpaut(p1); break;
1865: case 1: case 5:
1866: y=addsr(1,p1); break;
1867: case 2: case 4:
1868: y=mpaut(p1); setsigne(y,-signe(y)); break;
1869: default: /* case 3: case 7: */
1870: y=subsr(-1,p1); break;
1871: }
1.2 ! noro 1872: return gerepileuptoleaf(av, y);
1.1 noro 1873: }
1874:
1875: GEN
1876: gsin(GEN x, long prec)
1877: {
1.2 ! noro 1878: gpmem_t av, tetpil;
1.1 noro 1879: GEN r,u,v,y,p1,p2;
1880:
1881: switch(typ(x))
1882: {
1883: case t_REAL:
1884: return mpsin(x);
1885:
1886: case t_COMPLEX:
1887: y=cgetg(3,t_COMPLEX); av=avma;
1888: r=gexp((GEN)x[2],prec); p1=ginv(r);
1889: p2=gmul2n(gadd(p1,r),-1);
1890: p1=gsub(p2,p1);
1891: gsincos((GEN)x[1],&u,&v,prec);
1892: tetpil=avma;
1.2 ! noro 1893: y[1]=lmul(p2,u);
! 1894: y[2]=lmul(p1,v);
1.1 noro 1895: gerepilemanyvec(av,tetpil,y+1,2);
1896: return y;
1897:
1898: case t_SER:
1899: if (gcmp0(x)) return gcopy(x);
1900: if (valp(x)<0) err(negexper,"gsin");
1901:
1902: av=avma; gsincos(x,&u,&v,prec);
1903: return gerepilecopy(av,u);
1904:
1905: case t_INTMOD: case t_PADIC:
1906: err(typeer,"gsin");
1907: }
1908: return transc(gsin,x,prec);
1909: }
1910:
1911: void
1912: gsinz(GEN x, GEN y)
1913: {
1.2 ! noro 1914: long prec = precision(y);
! 1915: gpmem_t av=avma;
1.1 noro 1916:
1917: if (!prec) err(infprecer,"gsinz");
1918: gaffect(gsin(x,prec),y); avma=av;
1919: }
1920:
1921: /********************************************************************/
1.2 ! noro 1922: /** SINE, COSINE together **/
1.1 noro 1923: /********************************************************************/
1924:
1925: void
1926: mpsincos(GEN x, GEN *s, GEN *c)
1927: {
1.2 ! noro 1928: long mod8;
! 1929: gpmem_t av, tetpil;
1.1 noro 1930: GEN p1, *gptr[2];
1931:
1932: if (typ(x)!=t_REAL) err(typeer,"mpsincos");
1933: if (!signe(x))
1934: {
1.2 ! noro 1935: *s = realzero_bit(expo(x));
! 1936: *c = addsr(1,x); return;
1.1 noro 1937: }
1938:
1939: av=avma; p1=mpsc1(x,&mod8); tetpil=avma;
1940: switch(mod8)
1941: {
1942: case 0: *c=addsr( 1,p1); *s=mpaut(p1); break;
1943: case 1: *s=addsr( 1,p1); *c=mpaut(p1); setsigne(*c,-signe(*c)); break;
1944: case 2: *c=subsr(-1,p1); *s=mpaut(p1); setsigne(*s,-signe(*s)); break;
1945: case 3: *s=subsr(-1,p1); *c=mpaut(p1); break;
1946: case 4: *c=addsr( 1,p1); *s=mpaut(p1); setsigne(*s,-signe(*s)); break;
1947: case 5: *s=addsr( 1,p1); *c=mpaut(p1); break;
1948: case 6: *c=subsr(-1,p1); *s=mpaut(p1); break;
1949: case 7: *s=subsr(-1,p1); *c=mpaut(p1); setsigne(*c,-signe(*c)); break;
1950: }
1951: gptr[0]=s; gptr[1]=c;
1952: gerepilemanysp(av,tetpil,gptr,2);
1953: }
1954:
1.2 ! noro 1955: /* return exp(ix), x a t_REAL */
! 1956: GEN
! 1957: exp_Ir(GEN x)
! 1958: {
! 1959: gpmem_t av = avma;
! 1960: GEN v = cgetg(3,t_COMPLEX);
! 1961: mpsincos(x, (GEN*)(v+2), (GEN*)(v+1));
! 1962: if (!signe(x)) return gerepilecopy(av, (GEN)v[1]);
! 1963: return v;
! 1964: }
! 1965:
1.1 noro 1966: void
1967: gsincos(GEN x, GEN *s, GEN *c, long prec)
1968: {
1.2 ! noro 1969: long ii, i, j, ex, ex2, lx, ly, mi;
! 1970: gpmem_t av, tetpil;
1.1 noro 1971: GEN r,u,v,u1,v1,p1,p2,p3,p4,ps,pc;
1972: GEN *gptr[4];
1973:
1974: switch(typ(x))
1975: {
1976: case t_INT: case t_FRAC: case t_FRACN:
1977: av=avma; p1=cgetr(prec); gaffect(x,p1); tetpil=avma;
1978: mpsincos(p1,s,c); gptr[0]=s; gptr[1]=c;
1979: gerepilemanysp(av,tetpil,gptr,2);
1980: return;
1981:
1982: case t_REAL:
1983: mpsincos(x,s,c); return;
1984:
1985: case t_COMPLEX:
1986: ps=cgetg(3,t_COMPLEX); pc=cgetg(3,t_COMPLEX);
1987: *s=ps; *c=pc; av=avma;
1988: r=gexp((GEN)x[2],prec); p1=ginv(r);
1989: v1=gmul2n(gadd(p1,r),-1);
1990: u1=gsub(v1,p1); r=gsub(v1,r);/*u1=I*sin(I*Im(x));v1=cos(I*Im(x));r=-u1*/
1991: gsincos((GEN)x[1],&u,&v,prec);
1992: tetpil=avma;
1993: p1=gmul(v1,u); p2=gmul(u1,v);
1994: p3=gmul(v1,v); p4=gmul(r,u);
1995: gptr[0]=&p1; gptr[1]=&p2; gptr[2]=&p3; gptr[3]=&p4;
1996: gerepilemanysp(av,tetpil,gptr,4);
1.2 ! noro 1997: ps[1]=(long)p1; pc[1]=(long)p3;
! 1998: ps[2]=(long)p2; pc[2]=(long)p4;
1.1 noro 1999: return;
2000:
2001: case t_SER:
2002: if (gcmp0(x)) { *c=gaddsg(1,x); *s=gcopy(x); return; }
2003:
2004: ex=valp(x); lx=lg(x); ex2=2*ex+2;
2005: if (ex < 0) err(talker,"non zero exponent in gsincos");
2006: if (ex2 > lx)
2007: {
2008: *s=gcopy(x); av=avma; p1=gdivgs(gsqr(x),2);
2009: tetpil=avma; *c=gerepile(av,tetpil,gsubsg(1,p1));
2010: return;
2011: }
2012: if (!ex)
2013: {
2014: av=avma; p1=gcopy(x); p1[2]=zero;
2015: gsincos(normalize(p1),&u,&v,prec);
2016: gsincos((GEN)x[2],&u1,&v1,prec);
2017: p1=gmul(v1,v); p2=gmul(u1,u); p3=gmul(v1,u);
2018: p4=gmul(u1,v); tetpil=avma;
2019: *c=gsub(p1,p2); *s=gadd(p3,p4);
2020: gptr[0]=s; gptr[1]=c;
2021: gerepilemanysp(av,tetpil,gptr,2);
2022: return;
2023: }
2024:
2025: ly=lx+2*ex;
1.2 ! noro 2026: mi = lx-1; while (mi>=3 && gcmp0((GEN)x[mi])) mi--;
! 2027: mi += ex-2;
1.1 noro 2028: pc=cgetg(ly,t_SER); *c=pc;
2029: ps=cgetg(lx,t_SER); *s=ps;
2030: pc[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(x));
2031: pc[2]=un; ps[1]=x[1];
2032: for (i=2; i<ex+2; i++) ps[i]=lcopy((GEN)x[i]);
2033: for (i=3; i<ex2; i++) pc[i]=zero;
2034: for (i=ex2; i<ly; i++)
2035: {
2036: ii=i-ex; av=avma; p1=gzero;
1.2 ! noro 2037: for (j=ex; j<=min(ii-2,mi); j++)
1.1 noro 2038: p1=gadd(p1,gmulgs(gmul((GEN)x[j-ex+2],(GEN)ps[ii-j]),j));
2039: tetpil=avma;
2040: pc[i]=lpile(av,tetpil,gdivgs(p1,2-i));
2041: if (ii<lx)
2042: {
2043: av=avma; p1=gzero;
1.2 ! noro 2044: for (j=ex; j<=min(i-ex2,mi); j++)
1.1 noro 2045: p1=gadd(p1,gmulgs(gmul((GEN)x[j-ex+2],(GEN)pc[i-j]),j));
2046: p1=gdivgs(p1,i-2); tetpil=avma;
2047: ps[i-ex]=lpile(av,tetpil,gadd(p1,(GEN)x[i-ex]));
2048: }
2049: }
2050: return;
2051: /* transc doesn't work for this prototype */
2052: case t_QUAD:
2053: av = avma; p1=gmul(x,realun(prec)); tetpil = avma;
2054: gsincos(p1,s,c,prec);
2055: gptr[0]=s; gptr[1]=c;
2056: gerepilemanysp(av,tetpil,gptr,2);
2057: return;
2058:
2059: case t_POL: case t_RFRAC: case t_RFRACN:
2060: av = avma; p1=tayl(x,gvar(x),precdl); tetpil=avma;
2061: gsincos(p1,s,c,prec);
2062: gptr[0]=s; gptr[1]=c;
2063: gerepilemanysp(av,tetpil,gptr,2);
2064: return;
2065: }
2066: err(typeer,"gsincos");
2067: }
2068:
2069: /********************************************************************/
2070: /** **/
2071: /** FONCTIONS TANGENTE ET COTANGENTE **/
2072: /** **/
2073: /********************************************************************/
2074:
2075: static GEN
2076: mptan(GEN x)
2077: {
1.2 ! noro 2078: gpmem_t av=avma, tetpil;
1.1 noro 2079: GEN s,c;
2080:
2081: mpsincos(x,&s,&c); tetpil=avma;
2082: return gerepile(av,tetpil,divrr(s,c));
2083: }
2084:
2085: GEN
2086: gtan(GEN x, long prec)
2087: {
1.2 ! noro 2088: gpmem_t av, tetpil;
1.1 noro 2089: GEN s,c;
2090:
2091: switch(typ(x))
2092: {
2093: case t_REAL:
2094: return mptan(x);
2095:
2096: case t_SER:
2097: if (gcmp0(x)) return gcopy(x);
2098: if (valp(x)<0) err(negexper,"gtan"); /* fall through */
2099: case t_COMPLEX:
2100: av=avma; gsincos(x,&s,&c,prec); tetpil=avma;
2101: return gerepile(av,tetpil,gdiv(s,c));
2102:
2103: case t_INTMOD: case t_PADIC:
2104: err(typeer,"gtan");
2105: }
2106: return transc(gtan,x,prec);
2107: }
2108:
2109: void
2110: gtanz(GEN x, GEN y)
2111: {
1.2 ! noro 2112: long prec = precision(y);
! 2113: gpmem_t av=avma;
1.1 noro 2114:
2115: if (!prec) err(infprecer,"gtanz");
2116: gaffect(gtan(x,prec),y); avma=av;
2117: }
2118:
2119: static GEN
2120: mpcotan(GEN x)
2121: {
1.2 ! noro 2122: gpmem_t av=avma, tetpil;
1.1 noro 2123: GEN s,c;
2124:
2125: mpsincos(x,&s,&c); tetpil=avma;
2126: return gerepile(av,tetpil,divrr(c,s));
2127: }
2128:
2129: GEN
2130: gcotan(GEN x, long prec)
2131: {
1.2 ! noro 2132: gpmem_t av, tetpil;
1.1 noro 2133: GEN s,c;
2134:
2135: switch(typ(x))
2136: {
2137: case t_REAL:
2138: return mpcotan(x);
2139:
2140: case t_SER:
1.2 ! noro 2141: if (gcmp0(x)) err(talker,"0 argument in cotan");
! 2142: if (valp(x) < 0) err(negexper,"cotan"); /* fall through */
1.1 noro 2143: case t_COMPLEX:
2144: av=avma; gsincos(x,&s,&c,prec); tetpil=avma;
2145: return gerepile(av,tetpil,gdiv(c,s));
2146:
2147: case t_INTMOD: case t_PADIC:
2148: err(typeer,"gcotan");
2149: }
2150: return transc(gcotan,x,prec);
2151: }
2152:
2153: void
2154: gcotanz(GEN x, GEN y)
2155: {
1.2 ! noro 2156: long prec = precision(y);
! 2157: gpmem_t av=avma;
1.1 noro 2158:
2159: if (!prec) err(infprecer,"gcotanz");
2160: gaffect(gcotan(x,prec),y); avma=av;
2161: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to trans1.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath