Annotation of OpenXM_contrib2/asir2018/engine/P.c, Revision 1.1
1.1 ! 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
! 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
! 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: *
! 48: * $OpenXM$
! 49: */
! 50: #ifndef FBASE
! 51: #define FBASE
! 52: #endif
! 53:
! 54: #include "b.h"
! 55: #include "ca.h"
! 56:
! 57: #include "poly.c"
! 58:
! 59: void divcp(P f,Q q,P *rp)
! 60: {
! 61: DCP dc,dcr,dcr0;
! 62:
! 63: if ( !f )
! 64: *rp = 0;
! 65: else if ( NUM(f) )
! 66: DIVNUM(f,q,rp);
! 67: else {
! 68: for ( dc = DC(f), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 69: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
! 70: divcp(COEF(dc),q,&COEF(dcr));
! 71: }
! 72: NEXT(dcr) = 0; MKP(VR(f),dcr0,*rp);
! 73: }
! 74: }
! 75:
! 76: void diffp(VL vl,P p,V v,P *r)
! 77: {
! 78: P t;
! 79: DCP dc,dcr,dcr0;
! 80:
! 81: if ( !p || NUM(p) )
! 82: *r = 0;
! 83: else {
! 84: if ( v == VR(p) ) {
! 85: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 86: if ( !DEG(dc) ) continue;
! 87: MULPQ(COEF(dc),(P)DEG(dc),&t);
! 88: if ( t ) {
! 89: NEXTDC(dcr0,dcr); subz(DEG(dc),ONE,&DEG(dcr));
! 90: COEF(dcr) = t;
! 91: }
! 92: }
! 93: if ( !dcr0 )
! 94: *r = 0;
! 95: else {
! 96: NEXT(dcr) = 0; MKP(v,dcr0,*r);
! 97: }
! 98: } else {
! 99: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 100: diffp(vl,COEF(dc),v,&t);
! 101: if ( t ) {
! 102: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
! 103: }
! 104: }
! 105: if ( !dcr0 )
! 106: *r = 0;
! 107: else {
! 108: NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
! 109: }
! 110: }
! 111: }
! 112: }
! 113:
! 114: /* Euler derivation */
! 115: void ediffp(VL vl,P p,V v,P *r)
! 116: {
! 117: P t;
! 118: DCP dc,dcr,dcr0;
! 119:
! 120: if ( !p || NUM(p) )
! 121: *r = 0;
! 122: else {
! 123: if ( v == VR(p) ) {
! 124: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 125: if ( !DEG(dc) ) continue;
! 126: MULPQ(COEF(dc),(P)DEG(dc),&t);
! 127: if ( t ) {
! 128: NEXTDC(dcr0,dcr);
! 129: DEG(dcr) = DEG(dc);
! 130: COEF(dcr) = t;
! 131: }
! 132: }
! 133: if ( !dcr0 )
! 134: *r = 0;
! 135: else {
! 136: NEXT(dcr) = 0; MKP(v,dcr0,*r);
! 137: }
! 138: } else {
! 139: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 140: ediffp(vl,COEF(dc),v,&t);
! 141: if ( t ) {
! 142: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
! 143: }
! 144: }
! 145: if ( !dcr0 )
! 146: *r = 0;
! 147: else {
! 148: NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
! 149: }
! 150: }
! 151: }
! 152: }
! 153:
! 154: void coefp(P p,int d,P *pr)
! 155: {
! 156: DCP dc;
! 157: int sgn;
! 158: Z dq;
! 159:
! 160: if ( NUM(p) )
! 161: if ( d == 0 )
! 162: *pr = p;
! 163: else
! 164: *pr = 0;
! 165: else {
! 166: for ( STOQ(d,dq), dc = DC(p); dc; dc = NEXT(dc) )
! 167: if ( (sgn = cmpz(DEG(dc),dq)) > 0 )
! 168: continue;
! 169: else if ( sgn == 0 ) {
! 170: *pr = COEF(dc);
! 171: return;
! 172: } else {
! 173: *pr = 0;
! 174: break;
! 175: }
! 176: *pr = 0;
! 177: }
! 178: }
! 179:
! 180: int compp(VL vl,P p1,P p2)
! 181: {
! 182: DCP dc1,dc2;
! 183: V v1,v2;
! 184:
! 185: if ( !p1 )
! 186: return p2 ? -1 : 0;
! 187: else if ( !p2 )
! 188: return 1;
! 189: else if ( NUM(p1) )
! 190: return NUM(p2) ? (*cmpnumt[MAX(NID(p1),NID(p2))])(p1,p2) : -1;
! 191: else if ( NUM(p2) )
! 192: return 1;
! 193: if ( (v1 = VR(p1)) == (v2 = VR(p2)) ) {
! 194: for ( dc1 = DC(p1), dc2 = DC(p2);
! 195: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) )
! 196: switch ( cmpz(DEG(dc1),DEG(dc2)) ) {
! 197: case 1:
! 198: return 1; break;
! 199: case -1:
! 200: return -1; break;
! 201: default:
! 202: switch ( compp(vl,COEF(dc1),COEF(dc2)) ) {
! 203: case 1:
! 204: return 1; break;
! 205: case -1:
! 206: return -1; break;
! 207: default:
! 208: break;
! 209: }
! 210: break;
! 211: }
! 212: return dc1 ? 1 : (dc2 ? -1 : 0);
! 213: } else {
! 214: for ( ; v1 != VR(vl) && v2 != VR(vl); vl = NEXT(vl) );
! 215: return v1 == VR(vl) ? 1 : -1;
! 216: }
! 217: }
! 218:
! 219: int equalp(VL vl,P p1,P p2)
! 220: {
! 221: DCP dc1,dc2;
! 222: V v1,v2;
! 223:
! 224: if ( !p1 ) {
! 225: if ( !p2 ) return 1;
! 226: else return 0;
! 227: }
! 228: /* p1 != 0 */
! 229: if ( !p2 ) return 0;
! 230:
! 231: /* p1 != 0, p2 != 0 */
! 232: if ( NUM(p1) ) {
! 233: if ( !NUM(p2) ) return 0;
! 234: /* p1 and p2 are numbers */
! 235: if ( NID((Num)p1) != NID((Num)p2) ) return 0;
! 236: if ( !(*cmpnumt[NID(p1),NID(p1)])(p1,p2) ) return 1;
! 237: return 0;
! 238: }
! 239: if ( VR(p1) != VR(p2) ) return 0;
! 240: for ( dc1 = DC(p1), dc2 = DC(p2);
! 241: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) ) {
! 242: if ( cmpz(DEG(dc1),DEG(dc2)) ) return 0;
! 243: else if ( !equalp(vl,COEF(dc1),COEF(dc2)) ) return 0;
! 244: }
! 245: if ( dc1 || dc2 ) return 0;
! 246: else return 1;
! 247: }
! 248:
! 249: void csump(VL vl,P p,Q *c)
! 250: {
! 251: DCP dc;
! 252: Q s,s1,s2;
! 253:
! 254: if ( !p || NUM(p) )
! 255: *c = (Q)p;
! 256: else {
! 257: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
! 258: csump(vl,COEF(dc),&s1); addq(s,s1,&s2); s = s2;
! 259: }
! 260: *c = s;
! 261: }
! 262: }
! 263:
! 264: void degp(V v,P p,Z *d)
! 265: {
! 266: DCP dc;
! 267: Z m,m1;
! 268:
! 269: if ( !p || NUM(p) )
! 270: *d = 0;
! 271: else if ( v == VR(p) )
! 272: *d = DEG(DC(p));
! 273: else {
! 274: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
! 275: degp(v,COEF(dc),&m1);
! 276: if ( cmpz(m,m1) < 0 )
! 277: m = m1;
! 278: }
! 279: *d = m;
! 280: }
! 281: }
! 282:
! 283: void mulpc_trunc(VL vl,P p,P c,VN vn,P *pr);
! 284: void mulpq_trunc(P p,Q q,VN vn,P *pr);
! 285: void mulp_trunc(VL vl,P p1,P p2,VN vn,P *pr);
! 286:
! 287: void mulp_trunc(VL vl,P p1,P p2,VN vn,P *pr)
! 288: {
! 289: register DCP dc,dct,dcr,dcr0;
! 290: V v1,v2;
! 291: P t,s,u;
! 292: int n1,n2,i,d,d1;
! 293:
! 294: if ( !p1 || !p2 ) *pr = 0;
! 295: else if ( NUM(p1) )
! 296: mulpq_trunc(p2,(Q)p1,vn,pr);
! 297: else if ( NUM(p2) )
! 298: mulpq_trunc(p1,(Q)p2,vn,pr);
! 299: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
! 300: for ( ; vn->v && vn->v != v1; vn++ )
! 301: if ( vn->n ) {
! 302: /* p1,p2 do not contain vn->v */
! 303: *pr = 0;
! 304: return;
! 305: }
! 306: if ( !vn->v )
! 307: error("mulp_trunc : invalid vn");
! 308: d = vn->n;
! 309: for ( dc = DC(p2), s = 0; dc; dc = NEXT(dc) ) {
! 310: for ( dcr0 = 0, dct = DC(p1); dct; dct = NEXT(dct) ) {
! 311: d1 = QTOS(DEG(dct))+QTOS(DEG(dc));
! 312: if ( d1 >= d ) {
! 313: mulp_trunc(vl,COEF(dct),COEF(dc),vn+1,&t);
! 314: if ( t ) {
! 315: NEXTDC(dcr0,dcr);
! 316: STOQ(d1,DEG(dcr));
! 317: COEF(dcr) = t;
! 318: }
! 319: }
! 320: }
! 321: if ( dcr0 ) {
! 322: NEXT(dcr) = 0; MKP(v1,dcr0,t);
! 323: addp(vl,s,t,&u); s = u; t = u = 0;
! 324: }
! 325: }
! 326: *pr = s;
! 327: } else {
! 328: while ( v1 != VR(vl) && v2 != VR(vl) )
! 329: vl = NEXT(vl);
! 330: if ( v1 == VR(vl) )
! 331: mulpc_trunc(vl,p1,p2,vn,pr);
! 332: else
! 333: mulpc_trunc(vl,p2,p1,vn,pr);
! 334: }
! 335: }
! 336:
! 337: void mulpq_trunc(P p,Q q,VN vn,P *pr)
! 338: {
! 339: DCP dc,dcr,dcr0;
! 340: P t;
! 341: int i,d;
! 342: V v;
! 343:
! 344: if (!p || !q)
! 345: *pr = 0;
! 346: else if ( NUM(p) ) {
! 347: for ( ; vn->v; vn++ )
! 348: if ( vn->n ) {
! 349: *pr = 0;
! 350: return;
! 351: }
! 352: MULNUM(p,q,pr);
! 353: } else {
! 354: v = VR(p);
! 355: for ( ; vn->v && vn->v != v; vn++ ) {
! 356: if ( vn->n ) {
! 357: /* p does not contain vn->v */
! 358: *pr = 0;
! 359: return;
! 360: }
! 361: }
! 362: if ( !vn->v )
! 363: error("mulpq_trunc : invalid vn");
! 364: d = vn->n;
! 365: for ( dcr0 = 0, dc = DC(p); dc && QTOS(DEG(dc)) >= d; dc = NEXT(dc) ) {
! 366: mulpq_trunc(COEF(dc),q,vn+1,&t);
! 367: if ( t ) {
! 368: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
! 369: }
! 370: }
! 371: if ( dcr0 ) {
! 372: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
! 373: } else
! 374: *pr = 0;
! 375: }
! 376: }
! 377:
! 378: void mulpc_trunc(VL vl,P p,P c,VN vn,P *pr)
! 379: {
! 380: DCP dc,dcr,dcr0;
! 381: P t;
! 382: V v;
! 383: int i,d;
! 384:
! 385: if ( NUM(c) )
! 386: mulpq_trunc(p,(Q)c,vn,pr);
! 387: else {
! 388: v = VR(p);
! 389: for ( ; vn->v && vn->v != v; vn++ )
! 390: if ( vn->n ) {
! 391: /* p,c do not contain vn->v */
! 392: *pr = 0;
! 393: return;
! 394: }
! 395: if ( !vn->v )
! 396: error("mulpc_trunc : invalid vn");
! 397: d = vn->n;
! 398: for ( dcr0 = 0, dc = DC(p); dc && QTOS(DEG(dc)) >= d; dc = NEXT(dc) ) {
! 399: mulp_trunc(vl,COEF(dc),c,vn+1,&t);
! 400: if ( t ) {
! 401: NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
! 402: }
! 403: }
! 404: if ( dcr0 ) {
! 405: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
! 406: } else
! 407: *pr = 0;
! 408: }
! 409: }
! 410:
! 411: void quop_trunc(VL vl,P p1,P p2,VN vn,P *pr)
! 412: {
! 413: DCP dc,dcq0,dcq;
! 414: P t,s,m,lc2,qt;
! 415: V v1,v2;
! 416: Z d2;
! 417: VN vn1;
! 418:
! 419: if ( !p1 )
! 420: *pr = 0;
! 421: else if ( NUM(p2) )
! 422: divsp(vl,p1,p2,pr);
! 423: else if ( (v1 = VR(p1)) != (v2 = VR(p2)) ) {
! 424: for ( dcq0 = 0, dc = DC(p1); dc; dc = NEXT(dc) ) {
! 425: NEXTDC(dcq0,dcq);
! 426: DEG(dcq) = DEG(dc);
! 427: quop_trunc(vl,COEF(dc),p2,vn,&COEF(dcq));
! 428: }
! 429: NEXT(dcq) = 0;
! 430: MKP(v1,dcq0,*pr);
! 431: } else {
! 432: d2 = DEG(DC(p2));
! 433: lc2 = COEF(DC(p2));
! 434: t = p1;
! 435: dcq0 = 0;
! 436: /* vn1 = degree list of LC(p2) */
! 437: for ( vn1 = vn; vn1->v != v1; vn1++ );
! 438: vn1++;
! 439: while ( t ) {
! 440: dc = DC(t);
! 441: NEXTDC(dcq0,dcq);
! 442: subz(DEG(dc),d2,&DEG(dcq));
! 443: quop_trunc(vl,COEF(dc),lc2,vn1,&COEF(dcq));
! 444: NEXT(dcq) = 0;
! 445: MKP(v1,dcq,qt);
! 446: mulp_trunc(vl,p2,qt,vn,&m);
! 447: subp(vl,t,m,&s); t = s;
! 448: }
! 449: NEXT(dcq) = 0;
! 450: MKP(v1,dcq0,*pr);
! 451: }
! 452: }
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 / asir2018 / engine