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