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