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