Annotation of OpenXM_contrib/pari-2.2/src/basemath/trans1.c, Revision 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