Annotation of OpenXM_contrib2/asir2000/engine/P.c, Revision 1.12
1.2 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.12 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/P.c,v 1.11 2004/08/18 01:10:59 noro Exp $
1.2 noro 49: */
1.1 noro 50: #ifndef FBASE
51: #define FBASE
52: #endif
53:
54: #include "b.h"
55: #include "ca.h"
56:
1.10 noro 57: #include "poly.c"
1.1 noro 58:
59: void ptozp(p,sgn,c,pr)
60: P p;
61: int sgn;
62: Q *c;
63: P *pr;
64: {
1.12 ! noro 65: N nm,dn;
1.1 noro 66:
1.12 ! noro 67: if ( !p ) {
! 68: *c = 0; *pr = 0;
! 69: } else {
! 70: lgp(p,&nm,&dn);
! 71: if ( UNIN(dn) )
! 72: NTOQ(nm,sgn,*c);
! 73: else
! 74: NDTOQ(nm,dn,sgn,*c);
! 75: divcp(p,*c,pr);
! 76: }
1.1 noro 77: }
78:
79: void lgp(p,g,l)
80: P p;
81: N *g,*l;
82: {
1.12 ! noro 83: N g1,g2,l1,l2,l3,l4,tmp;
! 84: DCP dc;
1.1 noro 85:
1.12 ! noro 86: if ( NUM(p) ) {
! 87: *g = NM((Q)p);
! 88: if ( INT((Q)p) )
! 89: *l = ONEN;
! 90: else
! 91: *l = DN((Q)p);
! 92: } else {
! 93: dc = DC(p); lgp(COEF(dc),g,l);
! 94: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
! 95: lgp(COEF(dc),&g1,&l1); gcdn(*g,g1,&g2); *g = g2;
! 96: gcdn(*l,l1,&l2); kmuln(*l,l1,&l3); divn(l3,l2,&l4,&tmp); *l = l4;
! 97: }
! 98: }
1.1 noro 99: }
100:
101: void divcp(f,q,rp)
102: P f;
103: Q q;
104: P *rp;
105: {
1.12 ! noro 106: DCP dc,dcr,dcr0;
1.1 noro 107:
1.12 ! noro 108: if ( !f )
! 109: *rp = 0;
! 110: else if ( NUM(f) )
! 111: DIVNUM(f,q,rp);
! 112: else {
! 113: for ( dc = DC(f), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 114: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
! 115: divcp(COEF(dc),q,&COEF(dcr));
! 116: }
! 117: NEXT(dcr) = 0; MKP(VR(f),dcr0,*rp);
! 118: }
1.1 noro 119: }
120:
121: void diffp(vl,p,v,r)
122: VL vl;
123: P p;
124: V v;
125: P *r;
126: {
1.12 ! noro 127: P t;
! 128: DCP dc,dcr,dcr0;
1.1 noro 129:
1.12 ! noro 130: if ( !p || NUM(p) )
! 131: *r = 0;
! 132: else {
! 133: if ( v == VR(p) ) {
! 134: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 135: if ( !DEG(dc) ) continue;
! 136: MULPQ(COEF(dc),(P)DEG(dc),&t);
! 137: if ( t ) {
! 138: NEXTDC(dcr0,dcr); SUBQ(DEG(dc),ONE,&DEG(dcr));
! 139: COEF(dcr) = t;
! 140: }
! 141: }
! 142: if ( !dcr0 )
! 143: *r = 0;
! 144: else {
! 145: NEXT(dcr) = 0; MKP(v,dcr0,*r);
! 146: }
! 147: } else {
! 148: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 149: diffp(vl,COEF(dc),v,&t);
! 150: if ( t ) {
! 151: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
! 152: }
! 153: }
! 154: if ( !dcr0 )
! 155: *r = 0;
! 156: else {
! 157: NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
! 158: }
! 159: }
! 160: }
1.8 ohara 161: }
162:
163: /* Euler derivation */
164: void ediffp(vl,p,v,r)
165: VL vl;
166: P p;
167: V v;
168: P *r;
169: {
1.12 ! noro 170: P t;
! 171: DCP dc,dcr,dcr0;
1.8 ohara 172:
1.12 ! noro 173: if ( !p || NUM(p) )
! 174: *r = 0;
! 175: else {
! 176: if ( v == VR(p) ) {
! 177: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 178: if ( !DEG(dc) ) continue;
! 179: MULPQ(COEF(dc),(P)DEG(dc),&t);
! 180: if ( t ) {
! 181: NEXTDC(dcr0,dcr);
! 182: DEG(dcr) = DEG(dc);
! 183: COEF(dcr) = t;
! 184: }
! 185: }
! 186: if ( !dcr0 )
! 187: *r = 0;
! 188: else {
! 189: NEXT(dcr) = 0; MKP(v,dcr0,*r);
! 190: }
! 191: } else {
! 192: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 193: ediffp(vl,COEF(dc),v,&t);
! 194: if ( t ) {
! 195: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
! 196: }
! 197: }
! 198: if ( !dcr0 )
! 199: *r = 0;
! 200: else {
! 201: NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
! 202: }
! 203: }
! 204: }
1.1 noro 205: }
206:
207: void coefp(p,d,pr)
208: P p;
209: int d;
210: P *pr;
211: {
1.12 ! noro 212: DCP dc;
! 213: int sgn;
! 214: Q dq;
! 215:
! 216: if ( NUM(p) )
! 217: if ( d == 0 )
! 218: *pr = p;
! 219: else
! 220: *pr = 0;
! 221: else {
! 222: for ( STOQ(d,dq), dc = DC(p); dc; dc = NEXT(dc) )
! 223: if ( (sgn = cmpq(DEG(dc),dq)) > 0 )
! 224: continue;
! 225: else if ( sgn == 0 ) {
! 226: *pr = COEF(dc);
! 227: return;
! 228: } else {
! 229: *pr = 0;
! 230: break;
! 231: }
! 232: *pr = 0;
! 233: }
1.1 noro 234: }
235:
236: int compp(vl,p1,p2)
237: VL vl;
238: P p1,p2;
239: {
1.12 ! noro 240: DCP dc1,dc2;
! 241: V v1,v2;
1.1 noro 242:
1.12 ! noro 243: if ( !p1 )
! 244: return p2 ? -1 : 0;
! 245: else if ( !p2 )
! 246: return 1;
! 247: else if ( NUM(p1) )
! 248: return NUM(p2) ? (*cmpnumt[MAX(NID(p1),NID(p2))])(p1,p2) : -1;
! 249: else if ( NUM(p2) )
! 250: return 1;
! 251: if ( (v1 = VR(p1)) == (v2 = VR(p2)) ) {
! 252: for ( dc1 = DC(p1), dc2 = DC(p2);
! 253: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) )
! 254: switch ( cmpq(DEG(dc1),DEG(dc2)) ) {
! 255: case 1:
! 256: return 1; break;
! 257: case -1:
! 258: return -1; break;
! 259: default:
! 260: switch ( compp(vl,COEF(dc1),COEF(dc2)) ) {
! 261: case 1:
! 262: return 1; break;
! 263: case -1:
! 264: return -1; break;
! 265: default:
! 266: break;
! 267: }
! 268: break;
! 269: }
! 270: return dc1 ? 1 : (dc2 ? -1 : 0);
! 271: } else {
! 272: for ( ; v1 != VR(vl) && v2 != VR(vl); vl = NEXT(vl) );
! 273: return v1 == VR(vl) ? 1 : -1;
! 274: }
1.9 noro 275: }
276:
277: int equalp(vl,p1,p2)
278: VL vl;
279: P p1,p2;
280: {
1.12 ! noro 281: DCP dc1,dc2;
! 282: V v1,v2;
1.9 noro 283:
1.12 ! noro 284: if ( !p1 ) {
! 285: if ( !p2 ) return 1;
! 286: else return 0;
! 287: }
! 288: /* p1 != 0 */
! 289: if ( !p2 ) return 0;
! 290:
! 291: /* p1 != 0, p2 != 0 */
! 292: if ( NUM(p1) ) {
! 293: if ( !NUM(p2) ) return 0;
! 294: /* p1 and p2 are numbers */
! 295: if ( NID((Num)p1) != NID((Num)p2) ) return 0;
! 296: if ( !(*cmpnumt[NID(p1),NID(p1)])(p1,p2) ) return 1;
! 297: return 0;
! 298: }
! 299: if ( VR(p1) != VR(p2) ) return 0;
! 300: for ( dc1 = DC(p1), dc2 = DC(p2);
! 301: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) ) {
! 302: if ( cmpq(DEG(dc1),DEG(dc2)) ) return 0;
! 303: else if ( !equalp(vl,COEF(dc1),COEF(dc2)) ) return 0;
! 304: }
! 305: if ( dc1 || dc2 ) return 0;
! 306: else return 1;
1.1 noro 307: }
308:
309: void csump(vl,p,c)
310: VL vl;
311: P p;
312: Q *c;
313: {
1.12 ! noro 314: DCP dc;
! 315: Q s,s1,s2;
1.1 noro 316:
1.12 ! noro 317: if ( !p || NUM(p) )
! 318: *c = (Q)p;
! 319: else {
! 320: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
! 321: csump(vl,COEF(dc),&s1); addq(s,s1,&s2); s = s2;
! 322: }
! 323: *c = s;
! 324: }
1.1 noro 325: }
326:
327: void degp(v,p,d)
328: V v;
329: P p;
330: Q *d;
331: {
1.12 ! noro 332: DCP dc;
! 333: Q m,m1;
! 334:
! 335: if ( !p || NUM(p) )
! 336: *d = 0;
! 337: else if ( v == VR(p) )
! 338: *d = DEG(DC(p));
! 339: else {
! 340: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
! 341: degp(v,COEF(dc),&m1);
! 342: if ( cmpq(m,m1) < 0 )
! 343: m = m1;
! 344: }
! 345: *d = m;
! 346: }
1.6 noro 347: }
348:
349: void mulpc_trunc(VL vl,P p,P c,VN vn,P *pr);
350: void mulpq_trunc(P p,Q q,VN vn,P *pr);
351: void mulp_trunc(VL vl,P p1,P p2,VN vn,P *pr);
352:
353: void mulp_trunc(VL vl,P p1,P p2,VN vn,P *pr)
354: {
1.12 ! noro 355: register DCP dc,dct,dcr,dcr0;
! 356: V v1,v2;
! 357: P t,s,u;
! 358: int n1,n2,i,d,d1;
! 359:
! 360: if ( !p1 || !p2 ) *pr = 0;
! 361: else if ( NUM(p1) )
! 362: mulpq_trunc(p2,(Q)p1,vn,pr);
! 363: else if ( NUM(p2) )
! 364: mulpq_trunc(p1,(Q)p2,vn,pr);
! 365: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
! 366: for ( ; vn->v && vn->v != v1; vn++ )
! 367: if ( vn->n ) {
! 368: /* p1,p2 do not contain vn->v */
! 369: *pr = 0;
! 370: return;
! 371: }
! 372: if ( !vn->v )
! 373: error("mulp_trunc : invalid vn");
! 374: d = vn->n;
! 375: for ( dc = DC(p2), s = 0; dc; dc = NEXT(dc) ) {
! 376: for ( dcr0 = 0, dct = DC(p1); dct; dct = NEXT(dct) ) {
! 377: d1 = QTOS(DEG(dct))+QTOS(DEG(dc));
! 378: if ( d1 >= d ) {
! 379: mulp_trunc(vl,COEF(dct),COEF(dc),vn+1,&t);
! 380: if ( t ) {
! 381: NEXTDC(dcr0,dcr);
! 382: STOQ(d1,DEG(dcr));
! 383: COEF(dcr) = t;
! 384: }
! 385: }
! 386: }
! 387: if ( dcr0 ) {
! 388: NEXT(dcr) = 0; MKP(v1,dcr0,t);
! 389: addp(vl,s,t,&u); s = u; t = u = 0;
! 390: }
! 391: }
! 392: *pr = s;
! 393: } else {
! 394: while ( v1 != VR(vl) && v2 != VR(vl) )
! 395: vl = NEXT(vl);
! 396: if ( v1 == VR(vl) )
! 397: mulpc_trunc(vl,p1,p2,vn,pr);
! 398: else
! 399: mulpc_trunc(vl,p2,p1,vn,pr);
! 400: }
1.6 noro 401: }
402:
403: void mulpq_trunc(P p,Q q,VN vn,P *pr)
404: {
1.12 ! noro 405: DCP dc,dcr,dcr0;
! 406: P t;
! 407: int i,d;
! 408: V v;
! 409:
! 410: if (!p || !q)
! 411: *pr = 0;
! 412: else if ( NUM(p) ) {
! 413: for ( ; vn->v; vn++ )
! 414: if ( vn->n ) {
! 415: *pr = 0;
! 416: return;
! 417: }
! 418: MULNUM(p,q,pr);
! 419: } else {
! 420: v = VR(p);
! 421: for ( ; vn->v && vn->v != v; vn++ ) {
! 422: if ( vn->n ) {
! 423: /* p does not contain vn->v */
! 424: *pr = 0;
! 425: return;
! 426: }
! 427: }
! 428: if ( !vn->v )
! 429: error("mulpq_trunc : invalid vn");
! 430: d = vn->n;
! 431: for ( dcr0 = 0, dc = DC(p); dc && QTOS(DEG(dc)) >= d; dc = NEXT(dc) ) {
! 432: mulpq_trunc(COEF(dc),q,vn+1,&t);
! 433: if ( t ) {
! 434: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
! 435: }
! 436: }
! 437: if ( dcr0 ) {
! 438: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
! 439: } else
! 440: *pr = 0;
! 441: }
1.6 noro 442: }
443:
444: void mulpc_trunc(VL vl,P p,P c,VN vn,P *pr)
445: {
1.12 ! noro 446: DCP dc,dcr,dcr0;
! 447: P t;
! 448: V v;
! 449: int i,d;
! 450:
! 451: if ( NUM(c) )
! 452: mulpq_trunc(p,(Q)c,vn,pr);
! 453: else {
! 454: v = VR(p);
! 455: for ( ; vn->v && vn->v != v; vn++ )
! 456: if ( vn->n ) {
! 457: /* p,c do not contain vn->v */
! 458: *pr = 0;
! 459: return;
! 460: }
! 461: if ( !vn->v )
! 462: error("mulpc_trunc : invalid vn");
! 463: d = vn->n;
! 464: for ( dcr0 = 0, dc = DC(p); dc && QTOS(DEG(dc)) >= d; dc = NEXT(dc) ) {
! 465: mulp_trunc(vl,COEF(dc),c,vn+1,&t);
! 466: if ( t ) {
! 467: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
! 468: }
! 469: }
! 470: if ( dcr0 ) {
! 471: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
! 472: } else
! 473: *pr = 0;
! 474: }
1.7 noro 475: }
476:
477: void quop_trunc(VL vl,P p1,P p2,VN vn,P *pr)
478: {
1.12 ! noro 479: DCP dc,dcq0,dcq;
! 480: P t,s,m,lc2,qt;
! 481: V v1,v2;
! 482: Q d2;
! 483: VN vn1;
! 484:
! 485: if ( !p1 )
! 486: *pr = 0;
! 487: else if ( NUM(p2) )
! 488: divsp(vl,p1,p2,pr);
! 489: else if ( (v1 = VR(p1)) != (v2 = VR(p2)) ) {
! 490: for ( dcq0 = 0, dc = DC(p1); dc; dc = NEXT(dc) ) {
! 491: NEXTDC(dcq0,dcq);
! 492: DEG(dcq) = DEG(dc);
! 493: quop_trunc(vl,COEF(dc),p2,vn,&COEF(dcq));
! 494: }
! 495: NEXT(dcq) = 0;
! 496: MKP(v1,dcq0,*pr);
! 497: } else {
! 498: d2 = DEG(DC(p2));
! 499: lc2 = COEF(DC(p2));
! 500: t = p1;
! 501: dcq0 = 0;
! 502: /* vn1 = degree list of LC(p2) */
! 503: for ( vn1 = vn; vn1->v != v1; vn1++ );
! 504: vn1++;
! 505: while ( t ) {
! 506: dc = DC(t);
! 507: NEXTDC(dcq0,dcq);
! 508: subq(DEG(dc),d2,&DEG(dcq));
! 509: quop_trunc(vl,COEF(dc),lc2,vn1,&COEF(dcq));
! 510: NEXT(dcq) = 0;
! 511: MKP(v1,dcq,qt);
! 512: mulp_trunc(vl,p2,qt,vn,&m);
! 513: subp(vl,t,m,&s); t = s;
! 514: }
! 515: NEXT(dcq) = 0;
! 516: MKP(v1,dcq0,*pr);
! 517: }
1.1 noro 518: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to P.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine