![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Q.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to Q.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/engine/Q.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
! 2: #include "ca.h"
! 3: #include "base.h"
! 4: #include "inline.h"
! 5:
! 6: void addq(n1,n2,nr)
! 7: Q n1,n2,*nr;
! 8: {
! 9: N nm,dn,nm1,nm2,nm3,dn0,dn1,dn2,g,g1,g0,m;
! 10: int sgn;
! 11: unsigned t,t1;
! 12:
! 13: if ( !n1 )
! 14: *nr = n2;
! 15: else if ( !n2 )
! 16: *nr = n1;
! 17: else if ( INT(n1) )
! 18: if ( INT(n2) ) {
! 19: nm1 = NM(n1); nm2 = NM(n2);
! 20: if ( SGN(n1) == SGN(n2) ) {
! 21: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
! 22: t1 = BD(nm1)[0]; t = t1+BD(nm2)[0];
! 23: if ( t < t1 ) {
! 24: m = NALLOC(2); PL(m) = 2; BD(m)[0] = t; BD(m)[1] = 1;
! 25: } else
! 26: UTON(t,m);
! 27: } else
! 28: addn(NM(n1),NM(n2),&m);
! 29: NTOQ(m,SGN(n1),*nr);
! 30: } else {
! 31: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
! 32: t1 = BD(nm1)[0]; t = t1-BD(nm2)[0];
! 33: if ( !t )
! 34: sgn = 0;
! 35: else if ( t > t1 ) {
! 36: sgn = -1; t = -((int)t); UTON(t,m);
! 37: } else {
! 38: sgn = 1; UTON(t,m);
! 39: }
! 40: } else
! 41: sgn = subn(NM(n1),NM(n2),&m);
! 42: if ( sgn ) {
! 43: if ( SGN(n1) == sgn )
! 44: NTOQ(m,1,*nr);
! 45: else
! 46: NTOQ(m,-1,*nr);
! 47: } else
! 48: *nr = 0;
! 49: }
! 50: } else {
! 51: kmuln(NM(n1),DN(n2),&m);
! 52: if ( SGN(n1) == SGN(n2) ) {
! 53: sgn = SGN(n1); addn(m,NM(n2),&nm);
! 54: } else
! 55: sgn = SGN(n1)*subn(m,NM(n2),&nm);
! 56: if ( sgn ) {
! 57: dn = DN(n2); NDTOQ(nm,dn,sgn,*nr);
! 58: } else
! 59: *nr = 0;
! 60: }
! 61: else if ( INT(n2) ) {
! 62: kmuln(NM(n2),DN(n1),&m);
! 63: if ( SGN(n1) == SGN(n2) ) {
! 64: sgn = SGN(n1); addn(m,NM(n1),&nm);
! 65: } else
! 66: sgn = SGN(n1)*subn(NM(n1),m,&nm);
! 67: if ( sgn ) {
! 68: dn = DN(n1); NDTOQ(nm,dn,sgn,*nr);
! 69: } else
! 70: *nr = 0;
! 71: } else {
! 72: gcdn(DN(n1),DN(n2),&g); divsn(DN(n1),g,&dn1); divsn(DN(n2),g,&dn2);
! 73: kmuln(NM(n1),dn2,&nm1); kmuln(NM(n2),dn1,&nm2);
! 74: if ( SGN(n1) == SGN(n2) ) {
! 75: sgn = SGN(n1); addn(nm1,nm2,&nm3);
! 76: } else
! 77: sgn = SGN(n1)*subn(nm1,nm2,&nm3);
! 78: if ( sgn ) {
! 79: gcdn(nm3,g,&g1); divsn(nm3,g1,&nm); divsn(g,g1,&g0);
! 80: kmuln(dn1,dn2,&dn0); kmuln(g0,dn0,&dn);
! 81: if ( UNIN(dn) )
! 82: NTOQ(nm,sgn,*nr);
! 83: else
! 84: NDTOQ(nm,dn,sgn,*nr);
! 85: } else
! 86: *nr = 0;
! 87: }
! 88: }
! 89:
! 90: void subq(n1,n2,nr)
! 91: Q n1,n2,*nr;
! 92: {
! 93: Q m;
! 94:
! 95: if ( !n1 )
! 96: if ( !n2 )
! 97: *nr = 0;
! 98: else {
! 99: DUPQ(n2,*nr); SGN(*nr) = -SGN(n2);
! 100: }
! 101: else if ( !n2 )
! 102: *nr = n1;
! 103: else if ( n1 == n2 )
! 104: *nr = 0;
! 105: else {
! 106: DUPQ(n2,m); SGN(m) = -SGN(n2); addq(n1,m,nr);
! 107: }
! 108: }
! 109:
! 110: void mulq(n1,n2,nr)
! 111: Q n1,n2,*nr;
! 112: {
! 113: N nm,nm1,nm2,dn,dn1,dn2,g;
! 114: int sgn;
! 115: unsigned u,l;
! 116:
! 117: if ( !n1 || !n2 ) *nr = 0;
! 118: else if ( INT(n1) )
! 119: if ( INT(n2) ) {
! 120: nm1 = NM(n1); nm2 = NM(n2);
! 121: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
! 122: DM(BD(NM(n1))[0],BD(NM(n2))[0],u,l)
! 123: if ( u ) {
! 124: nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = l; BD(nm)[1] = u;
! 125: } else {
! 126: nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = l;
! 127: }
! 128: } else
! 129: kmuln(nm1,nm2,&nm);
! 130: sgn = (SGN(n1)==SGN(n2)?1:-1); NTOQ(nm,sgn,*nr);
! 131: } else {
! 132: gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn);
! 133: kmuln(nm1,NM(n2),&nm); sgn = SGN(n1)*SGN(n2);
! 134: if ( UNIN(dn) )
! 135: NTOQ(nm,sgn,*nr);
! 136: else
! 137: NDTOQ(nm,dn,sgn,*nr);
! 138: }
! 139: else if ( INT(n2) ) {
! 140: gcdn(NM(n2),DN(n1),&g); divsn(NM(n2),g,&nm2); divsn(DN(n1),g,&dn);
! 141: kmuln(nm2,NM(n1),&nm); sgn = SGN(n1)*SGN(n2);
! 142: if ( UNIN(dn) )
! 143: NTOQ(nm,sgn,*nr);
! 144: else
! 145: NDTOQ(nm,dn,sgn,*nr);
! 146: } else {
! 147: gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn2);
! 148: gcdn(DN(n1),NM(n2),&g); divsn(DN(n1),g,&dn1); divsn(NM(n2),g,&nm2);
! 149: kmuln(nm1,nm2,&nm); kmuln(dn1,dn2,&dn); sgn = SGN(n1) * SGN(n2);
! 150: if ( UNIN(dn) )
! 151: NTOQ(nm,sgn,*nr);
! 152: else
! 153: NDTOQ(nm,dn,sgn,*nr);
! 154: }
! 155: }
! 156:
! 157: void divq(n1,n2,nq)
! 158: Q n1,n2,*nq;
! 159: {
! 160: Q m;
! 161:
! 162: if ( !n2 ) {
! 163: error("division by 0");
! 164: *nq = 0;
! 165: return;
! 166: } else if ( !n1 )
! 167: *nq = 0;
! 168: else if ( n1 == n2 )
! 169: *nq = ONE;
! 170: else {
! 171: invq(n2,&m); mulq(n1,m,nq);
! 172: }
! 173: }
! 174:
! 175: void invq(n,nr)
! 176: Q n,*nr;
! 177: {
! 178: N nm,dn;
! 179:
! 180: if ( INT(n) )
! 181: if ( UNIN(NM(n)) )
! 182: if ( SGN(n) > 0 )
! 183: *nr = ONE;
! 184: else {
! 185: nm = ONEN; NTOQ(nm,SGN(n),*nr);
! 186: }
! 187: else {
! 188: nm = ONEN; dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);
! 189: }
! 190: else if ( UNIN(NM(n)) ) {
! 191: nm = DN(n); NTOQ(nm,SGN(n),*nr);
! 192: } else {
! 193: nm = DN(n); dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);
! 194: }
! 195: }
! 196:
! 197: void chsgnq(n,nr)
! 198: Q n,*nr;
! 199: {
! 200: Q t;
! 201:
! 202: if ( !n )
! 203: *nr = 0;
! 204: else {
! 205: DUPQ(n,t); SGN(t) = -SGN(t); *nr = t;
! 206: }
! 207: }
! 208:
! 209: void pwrq(n1,n,nr)
! 210: Q n1,n,*nr;
! 211: {
! 212: N nm,dn;
! 213: int sgn;
! 214:
! 215: if ( !n || UNIQ(n1) )
! 216: *nr = ONE;
! 217: else if ( !n1 )
! 218: *nr = 0;
! 219: else if ( !INT(n) ) {
! 220: error("can't calculate fractional power."); *nr = 0;
! 221: } else if ( MUNIQ(n1) )
! 222: *nr = BD(NM(n))[0]%2 ? n1 : ONE;
! 223: else if ( PL(NM(n)) > 1 ) {
! 224: error("exponent too big."); *nr = 0;
! 225: } else {
! 226: sgn = ((SGN(n1)>0)||EVENN(NM(n))?1:-1);
! 227: pwrn(NM(n1),BD(NM(n))[0],&nm);
! 228: if ( INT(n1) ) {
! 229: if ( UNIN(nm) )
! 230: if ( sgn == 1 )
! 231: *nr = ONE;
! 232: else
! 233: NTOQ(ONEN,sgn,*nr);
! 234: else if ( SGN(n) > 0 )
! 235: NTOQ(nm,sgn,*nr);
! 236: else
! 237: NDTOQ(ONEN,nm,sgn,*nr);
! 238: } else {
! 239: pwrn(DN(n1),BD(NM(n))[0],&dn);
! 240: if ( SGN(n) > 0 )
! 241: NDTOQ(nm,dn,sgn,*nr);
! 242: else if ( UNIN(nm) )
! 243: NTOQ(dn,sgn,*nr);
! 244: else
! 245: NDTOQ(dn,nm,sgn,*nr);
! 246: }
! 247: }
! 248: }
! 249:
! 250: int cmpq(q1,q2)
! 251: Q q1,q2;
! 252: {
! 253: int sgn;
! 254: N t,s;
! 255:
! 256: if ( !q1 )
! 257: if ( !q2 )
! 258: return ( 0 );
! 259: else
! 260: return ( -SGN(q2) );
! 261: else if ( !q2 )
! 262: return ( SGN(q1) );
! 263: else if ( SGN(q1) != SGN(q2) )
! 264: return ( SGN(q1) );
! 265: else if ( INT(q1) )
! 266: if ( INT(q2) ) {
! 267: sgn = cmpn(NM(q1),NM(q2));
! 268: if ( !sgn )
! 269: return ( 0 );
! 270: else
! 271: return ( SGN(q1)==sgn?1:-1 );
! 272: } else {
! 273: kmuln(NM(q1),DN(q2),&t); sgn = cmpn(t,NM(q2));
! 274: return ( SGN(q1) * sgn );
! 275: }
! 276: else if ( INT(q2) ) {
! 277: kmuln(NM(q2),DN(q1),&t); sgn = cmpn(NM(q1),t);
! 278: return ( SGN(q1) * sgn );
! 279: } else {
! 280: kmuln(NM(q1),DN(q2),&t); kmuln(NM(q2),DN(q1),&s); sgn = cmpn(t,s);
! 281: return ( SGN(q1) * sgn );
! 282: }
! 283: }
! 284:
! 285: void remq(n,m,nr)
! 286: Q n,m;
! 287: Q *nr;
! 288: {
! 289: N q,r,s;
! 290:
! 291: if ( !n ) {
! 292: *nr = 0;
! 293: return;
! 294: }
! 295: divn(NM(n),NM(m),&q,&r);
! 296: if ( !r )
! 297: *nr = 0;
! 298: else if ( SGN(n) > 0 )
! 299: NTOQ(r,1,*nr);
! 300: else {
! 301: subn(NM(m),r,&s); NTOQ(s,1,*nr);
! 302: }
! 303: }
! 304:
! 305: void mkbc(n,t)
! 306: int n;
! 307: Q *t;
! 308: {
! 309: int i;
! 310: N c,d;
! 311:
! 312: for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {
! 313: c = NALLOC(1); PL(c) = 1; BD(c)[0] = n-i+1;
! 314: kmuln(NM(t[i-1]),c,&d); divin(d,i,&c); NTOQ(c,1,t[i]);
! 315: }
! 316: for ( ; i <= n; i++ )
! 317: t[i] = t[n-i];
! 318: }
! 319:
! 320: void factorial(n,r)
! 321: int n;
! 322: Q *r;
! 323: {
! 324: Q t,iq,s;
! 325: unsigned int i,m,m0;
! 326: unsigned int c;
! 327:
! 328: if ( !n )
! 329: *r = ONE;
! 330: else if ( n < 0 )
! 331: *r = 0;
! 332: else {
! 333: for ( i = 1, t = ONE; (int)i <= n; ) {
! 334: for ( m0 = m = 1, c = 0; ((int)i <= n) && !c; i++ ) {
! 335: m0 = m; DM(m0,i,c,m)
! 336: }
! 337: if ( c ) {
! 338: m = m0; i--;
! 339: }
! 340: UTOQ(m,iq); mulq(t,iq,&s); t = s;
! 341: }
! 342: *r = t;
! 343: }
! 344: }
! 345:
! 346: void invl(a,mod,ar)
! 347: Q a,mod,*ar;
! 348: {
! 349: Q f1,f2,a1,a2,q,m,s;
! 350: N qn,rn;
! 351:
! 352: a1 = ONE; a2 = 0;
! 353: DUPQ(a,f1); SGN(f1) = 1; f2 = mod;
! 354:
! 355: while ( 1 ) {
! 356: divn(NM(f1),NM(f2),&qn,&rn);
! 357: if ( !qn )
! 358: q = 0;
! 359: else
! 360: NTOQ(qn,1,q);
! 361:
! 362: f1 = f2;
! 363: if ( !rn )
! 364: break;
! 365: else
! 366: NTOQ(rn,1,f2);
! 367:
! 368: mulq(a2,q,&m); subq(a1,m,&s); a1 = a2;
! 369: if ( !s )
! 370: a2 = 0;
! 371: else
! 372: remq(s,mod,&a2);
! 373: }
! 374: if ( SGN(a) < 0 )
! 375: chsgnq(a2,&m);
! 376: else
! 377: m = a2;
! 378:
! 379: if ( SGN(m) >= 0 )
! 380: *ar = m;
! 381: else
! 382: addq(m,mod,ar);
! 383: }
! 384:
! 385: int kara_mag=100;
! 386:
! 387: void kmuln(n1,n2,nr)
! 388: N n1,n2,*nr;
! 389: {
! 390: N n,t,s,m,carry;
! 391: int d,d1,d2,len,i,l;
! 392: int *r,*r0;
! 393:
! 394: if ( !n1 || !n2 ) {
! 395: *nr = 0; return;
! 396: }
! 397: d1 = PL(n1); d2 = PL(n2);
! 398: if ( (d1 < kara_mag) || (d2 < kara_mag) ) {
! 399: muln(n1,n2,nr); return;
! 400: }
! 401: if ( d1 < d2 ) {
! 402: n = n1; n1 = n2; n2 = n;
! 403: d = d1; d1 = d2; d2 = d;
! 404: }
! 405: if ( d2 > (d1+1)/2 ) {
! 406: kmulnmain(n1,n2,nr); return;
! 407: }
! 408: d = (d1/d2)+((d1%d2)!=0);
! 409: len = (d+1)*d2;
! 410: r0 = (int *)ALLOCA(len*sizeof(int));
! 411: bzero((char *)r0,len*sizeof(int));
! 412: for ( carry = 0, i = 0, r = r0; i < d; i++, r += d2 ) {
! 413: extractn(n1,i*d2,d2,&m);
! 414: if ( m ) {
! 415: kmuln(m,n2,&t);
! 416: addn(t,carry,&s);
! 417: copyn(s,d2,r);
! 418: extractn(s,d2,d2,&carry);
! 419: } else {
! 420: copyn(carry,d2,r);
! 421: carry = 0;
! 422: }
! 423: }
! 424: copyn(carry,d2,r);
! 425: for ( l = len - 1; !r0[l]; l-- );
! 426: l++;
! 427: *nr = t = NALLOC(l); PL(t) = l;
! 428: bcopy((char *)r0,(char *)BD(t),l*sizeof(int));
! 429: }
! 430:
! 431: void extractn(n,index,len,nr)
! 432: N n;
! 433: int index,len;
! 434: N *nr;
! 435: {
! 436: unsigned int *m;
! 437: int l;
! 438: N t;
! 439:
! 440: if ( !n ) {
! 441: *nr = 0; return;
! 442: }
! 443: m = BD(n)+index;
! 444: if ( (l = PL(n)-index) >= len ) {
! 445: for ( l = len - 1; (l >= 0) && !m[l]; l-- );
! 446: l++;
! 447: }
! 448: if ( l <= 0 ) {
! 449: *nr = 0; return;
! 450: } else {
! 451: *nr = t = NALLOC(l); PL(t) = l;
! 452: bcopy((char *)m,(char *)BD(t),l*sizeof(int));
! 453: }
! 454: }
! 455:
! 456: void copyn(n,len,p)
! 457: N n;
! 458: int len;
! 459: int *p;
! 460: {
! 461: if ( n )
! 462: bcopy((char *)BD(n),(char *)p,MIN(PL(n),len)*sizeof(int));
! 463: }
! 464:
! 465: void dupn(n,p)
! 466: N n;
! 467: N p;
! 468: {
! 469: if ( n )
! 470: bcopy((char *)n,(char *)p,(1+PL(n))*sizeof(int));
! 471: }
! 472:
! 473: void kmulnmain(n1,n2,nr)
! 474: N n1,n2,*nr;
! 475: {
! 476: int d1,d2,h,sgn,sgn1,sgn2,len;
! 477: N n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
! 478:
! 479: d1 = PL(n1); d2 = PL(n2); h = (d1+1)/2;
! 480: extractn(n1,0,h,&n1lo); extractn(n1,h,d1-h,&n1hi);
! 481: extractn(n2,0,h,&n2lo); extractn(n2,h,d2-h,&n2hi);
! 482: kmuln(n1hi,n2hi,&hi); kmuln(n1lo,n2lo,&lo);
! 483: len = PL(hi)+2*h; t1 = NALLOC(len); PL(t1) = len;
! 484: bzero((char *)BD(t1),len*sizeof(int));
! 485: if ( lo )
! 486: bcopy((char *)BD(lo),(char *)BD(t1),PL(lo)*sizeof(int));
! 487: if ( hi )
! 488: bcopy((char *)BD(hi),(char *)(BD(t1)+2*h),PL(hi)*sizeof(int));
! 489:
! 490: addn(hi,lo,&mid1);
! 491: sgn1 = subn(n1hi,n1lo,&s1); sgn2 = subn(n2lo,n2hi,&s2);
! 492: kmuln(s1,s2,&mid2); sgn = sgn1*sgn2;
! 493: if ( sgn > 0 )
! 494: addn(mid1,mid2,&mid);
! 495: else
! 496: subn(mid1,mid2,&mid);
! 497: if ( mid ) {
! 498: len = PL(mid)+h; t2 = NALLOC(len); PL(t2) = len;
! 499: bzero((char *)BD(t2),len*sizeof(int));
! 500: bcopy((char *)BD(mid),(char *)(BD(t2)+h),PL(mid)*sizeof(int));
! 501: addn(t1,t2,nr);
! 502: } else
! 503: *nr = t1;
! 504: }