![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to dist.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 dist.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.8 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.9 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.8 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.21 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.20 2002/01/28 00:54:43 noro Exp $
1.8 noro 49: */
1.1 noro 50: #include "ca.h"
51:
52: #define ORD_REVGRADLEX 0
53: #define ORD_GRADLEX 1
54: #define ORD_LEX 2
55: #define ORD_BREVGRADLEX 3
56: #define ORD_BGRADLEX 4
57: #define ORD_BLEX 5
58: #define ORD_BREVREV 6
59: #define ORD_BGRADREV 7
60: #define ORD_BLEXREV 8
61: #define ORD_ELIM 9
1.12 noro 62: #define ORD_WEYL_ELIM 10
1.13 noro 63: #define ORD_HOMO_WW_DRL 11
1.21 ! noro 64: #define ORD_DRL_ZIGZAG 12
! 65: #define ORD_HOMO_WW_DRL_ZIGZAG 13
! 66:
! 67: int cmpdl_drl_zigzag(), cmpdl_homo_ww_drl_zigzag();
1.1 noro 68:
69: int (*cmpdl)()=cmpdl_revgradlex;
70: int (*primitive_cmpdl[3])() = {cmpdl_revgradlex,cmpdl_gradlex,cmpdl_lex};
71:
1.2 noro 72: int do_weyl;
73:
1.1 noro 74: int dp_nelim,dp_fcoeffs;
75: struct order_spec dp_current_spec;
76: int *dp_dl_work;
77:
1.19 noro 78: int has_fcoef(DP f)
1.1 noro 79: {
80: MP t;
81:
82: if ( !f )
83: return 0;
84: for ( t = BDY(f); t; t = NEXT(t) )
85: if ( has_fcoef_p(t->c) )
86: break;
87: return t ? 1 : 0;
88: }
89:
1.19 noro 90: int has_fcoef_p(P f)
1.1 noro 91: {
92: DCP dc;
93:
94: if ( !f )
95: return 0;
96: else if ( NUM(f) )
1.14 noro 97: return (NID((Num)f) == N_LM
98: || NID((Num)f) == N_GF2N
1.15 noro 99: || NID((Num)f) == N_GFPN
100: || NID((Num)f) == N_GFS) ? 1 : 0;
1.1 noro 101: else {
102: for ( dc = DC(f); dc; dc = NEXT(dc) )
103: if ( has_fcoef_p(COEF(dc)) )
104: return 1;
105: return 0;
106: }
107: }
108:
1.19 noro 109: void initd(struct order_spec *spec)
1.1 noro 110: {
111: switch ( spec->id ) {
112: case 2:
113: cmpdl = cmpdl_matrix;
114: dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
115: break;
116: case 1:
117: cmpdl = cmpdl_order_pair;
118: break;
119: default:
120: switch ( spec->ord.simple ) {
121: case ORD_REVGRADLEX:
122: cmpdl = cmpdl_revgradlex; break;
123: case ORD_GRADLEX:
124: cmpdl = cmpdl_gradlex; break;
125: case ORD_BREVGRADLEX:
126: cmpdl = cmpdl_brevgradlex; break;
127: case ORD_BGRADLEX:
128: cmpdl = cmpdl_bgradlex; break;
129: case ORD_BLEX:
130: cmpdl = cmpdl_blex; break;
131: case ORD_BREVREV:
132: cmpdl = cmpdl_brevrev; break;
133: case ORD_BGRADREV:
134: cmpdl = cmpdl_bgradrev; break;
135: case ORD_BLEXREV:
136: cmpdl = cmpdl_blexrev; break;
137: case ORD_ELIM:
138: cmpdl = cmpdl_elim; break;
1.12 noro 139: case ORD_WEYL_ELIM:
140: cmpdl = cmpdl_weyl_elim; break;
1.13 noro 141: case ORD_HOMO_WW_DRL:
142: cmpdl = cmpdl_homo_ww_drl; break;
1.21 ! noro 143: case ORD_DRL_ZIGZAG:
! 144: cmpdl = cmpdl_drl_zigzag; break;
! 145: case ORD_HOMO_WW_DRL_ZIGZAG:
! 146: cmpdl = cmpdl_homo_ww_drl_zigzag; break;
1.1 noro 147: case ORD_LEX: default:
148: cmpdl = cmpdl_lex; break;
149: }
150: break;
151: }
152: dp_current_spec = *spec;
153: }
154:
1.19 noro 155: void ptod(VL vl,VL dvl,P p,DP *pr)
1.1 noro 156: {
157: int isconst = 0;
1.16 noro 158: int n,i,j,k;
1.1 noro 159: VL tvl;
160: V v;
161: DL d;
162: MP m;
163: DCP dc;
1.16 noro 164: DCP *w;
1.1 noro 165: DP r,s,t,u;
166: P x,c;
167:
168: if ( !p )
169: *pr = 0;
170: else {
171: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
172: if ( NUM(p) ) {
173: NEWDL(d,n);
174: NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
175: } else {
176: for ( i = 0, tvl = dvl, v = VR(p);
177: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
178: if ( !tvl ) {
1.16 noro 179: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
180: w = (DCP *)ALLOCA(k*sizeof(DCP));
181: for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
182: w[j] = dc;
183:
184: for ( j = k-1, s = 0, MKV(v,x); j >= 0; j-- ) {
185: ptod(vl,dvl,COEF(w[j]),&t); pwrp(vl,x,DEG(w[j]),&c);
1.1 noro 186: muldc(vl,t,c,&r); addd(vl,r,s,&t); s = t;
187: }
188: *pr = s;
189: } else {
1.16 noro 190: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
191: w = (DCP *)ALLOCA(k*sizeof(DCP));
192: for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
193: w[j] = dc;
194:
195: for ( j = k-1, s = 0; j >= 0; j-- ) {
196: ptod(vl,dvl,COEF(w[j]),&t);
1.20 noro 197: NEWDL(d,n); d->d[i] = QTOS(DEG(w[j]));
198: d->td = MUL_WEIGHT(d->d[i],i);
1.1 noro 199: NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
1.2 noro 200: comm_muld(vl,t,u,&r); addd(vl,r,s,&t); s = t;
1.1 noro 201: }
202: *pr = s;
203: }
204: }
205: }
1.17 noro 206: #if 0
1.1 noro 207: if ( !dp_fcoeffs && has_fcoef(*pr) )
208: dp_fcoeffs = 1;
1.17 noro 209: #endif
1.1 noro 210: }
211:
1.19 noro 212: void dtop(VL vl,VL dvl,DP p,P *pr)
1.1 noro 213: {
1.16 noro 214: int n,i,j,k;
1.1 noro 215: DL d;
216: MP m;
1.16 noro 217: MP *a;
1.1 noro 218: P r,s,t,u,w;
219: Q q;
220: VL tvl;
221:
222: if ( !p )
223: *pr = 0;
224: else {
1.16 noro 225: for ( k = 0, m = BDY(p); m; m = NEXT(m), k++ );
226: a = (MP *)ALLOCA(k*sizeof(MP));
227: for ( j = 0, m = BDY(p); j < k; m = NEXT(m), j++ )
228: a[j] = m;
229:
230: for ( n = p->nv, j = k-1, s = 0; j >= 0; j-- ) {
231: m = a[j];
1.1 noro 232: t = C(m);
233: if ( NUM(t) && NID((Num)t) == N_M ) {
234: mptop(t,&u); t = u;
235: }
236: for ( i = 0, d = m->dl, tvl = dvl;
237: i < n; tvl = NEXT(tvl), i++ ) {
238: MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
239: mulp(vl,t,u,&w); t = w;
240: }
241: addp(vl,s,t,&u); s = u;
242: }
243: *pr = s;
244: }
245: }
246:
1.19 noro 247: void nodetod(NODE node,DP *dp)
1.1 noro 248: {
249: NODE t;
250: int len,i,td;
251: Q e;
252: DL d;
253: MP m;
254: DP u;
255:
256: for ( t = node, len = 0; t; t = NEXT(t), len++ );
257: NEWDL(d,len);
258: for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
259: e = (Q)BDY(t);
260: if ( !e )
261: d->d[i] = 0;
262: else if ( !NUM(e) || !RATN(e) || !INT(e) )
263: error("nodetod : invalid input");
264: else {
1.20 noro 265: d->d[i] = QTOS((Q)e); td += MUL_WEIGHT(d->d[i],i);
1.1 noro 266: }
267: }
268: d->td = td;
269: NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0;
270: MKDP(len,m,u); u->sugar = td; *dp = u;
271: }
272:
1.19 noro 273: int sugard(MP m)
1.1 noro 274: {
275: int s;
276:
277: for ( s = 0; m; m = NEXT(m) )
278: s = MAX(s,m->dl->td);
279: return s;
280: }
281:
1.19 noro 282: void addd(VL vl,DP p1,DP p2,DP *pr)
1.1 noro 283: {
284: int n;
285: MP m1,m2,mr,mr0;
286: P t;
287:
288: if ( !p1 )
289: *pr = p2;
290: else if ( !p2 )
291: *pr = p1;
292: else {
293: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
294: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
295: case 0:
296: addp(vl,C(m1),C(m2),&t);
297: if ( t ) {
298: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
299: }
300: m1 = NEXT(m1); m2 = NEXT(m2); break;
301: case 1:
302: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
303: m1 = NEXT(m1); break;
304: case -1:
305: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
306: m2 = NEXT(m2); break;
307: }
308: if ( !mr0 )
309: if ( m1 )
310: mr0 = m1;
311: else if ( m2 )
312: mr0 = m2;
313: else {
314: *pr = 0;
315: return;
316: }
317: else if ( m1 )
318: NEXT(mr) = m1;
319: else if ( m2 )
320: NEXT(mr) = m2;
321: else
322: NEXT(mr) = 0;
323: MKDP(NV(p1),mr0,*pr);
324: if ( *pr )
325: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
326: }
327: }
328:
329: /* for F4 symbolic reduction */
330:
1.19 noro 331: void symb_addd(DP p1,DP p2,DP *pr)
1.1 noro 332: {
333: int n;
334: MP m1,m2,mr,mr0;
335:
336: if ( !p1 )
337: *pr = p2;
338: else if ( !p2 )
339: *pr = p1;
340: else {
341: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
342: NEXTMP(mr0,mr); C(mr) = (P)ONE;
343: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
344: case 0:
345: mr->dl = m1->dl;
346: m1 = NEXT(m1); m2 = NEXT(m2); break;
347: case 1:
348: mr->dl = m1->dl;
349: m1 = NEXT(m1); break;
350: case -1:
351: mr->dl = m2->dl;
352: m2 = NEXT(m2); break;
353: }
354: }
355: if ( !mr0 )
356: if ( m1 )
357: mr0 = m1;
358: else if ( m2 )
359: mr0 = m2;
360: else {
361: *pr = 0;
362: return;
363: }
364: else if ( m1 )
365: NEXT(mr) = m1;
366: else if ( m2 )
367: NEXT(mr) = m2;
368: else
369: NEXT(mr) = 0;
370: MKDP(NV(p1),mr0,*pr);
371: if ( *pr )
372: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1.3 noro 373: }
374: }
375:
376: /*
377: * destructive merge of two list
378: *
379: * p1, p2 : list of DL
380: * return : a merged list
381: */
382:
1.19 noro 383: NODE symb_merge(NODE m1,NODE m2,int n)
1.3 noro 384: {
385: NODE top,prev,cur,m,t;
386:
387: if ( !m1 )
388: return m2;
389: else if ( !m2 )
390: return m1;
391: else {
392: switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
393: case 0:
394: top = m1; m = NEXT(m2);
395: break;
396: case 1:
397: top = m1; m = m2;
398: break;
399: case -1:
400: top = m2; m = m1;
401: break;
402: }
403: prev = top; cur = NEXT(top);
404: /* BDY(prev) > BDY(m) always holds */
405: while ( cur && m ) {
406: switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
407: case 0:
408: m = NEXT(m);
409: prev = cur; cur = NEXT(cur);
410: break;
411: case 1:
412: t = NEXT(cur); NEXT(cur) = m; m = t;
413: prev = cur; cur = NEXT(cur);
414: break;
415: case -1:
416: NEXT(prev) = m; m = cur;
417: prev = NEXT(prev); cur = NEXT(prev);
418: break;
1.18 noro 419: }
420: }
421: if ( !cur )
422: NEXT(prev) = m;
423: return top;
424: }
425: }
426:
1.19 noro 427: DLBUCKET symb_merge_bucket(DLBUCKET m1,DLBUCKET m2,int n)
1.18 noro 428: {
429: DLBUCKET top,prev,cur,m,t;
430:
431: if ( !m1 )
432: return m2;
433: else if ( !m2 )
434: return m1;
435: else {
436: if ( m1->td == m2->td ) {
437: top = m1;
438: BDY(top) = symb_merge(BDY(top),BDY(m2),n);
439: m = NEXT(m2);
440: } else if ( m1->td > m2->td ) {
441: top = m1; m = m2;
442: } else {
443: top = m2; m = m1;
444: }
445: prev = top; cur = NEXT(top);
446: /* prev->td > m->td always holds */
447: while ( cur && m ) {
448: if ( cur->td == m->td ) {
449: BDY(cur) = symb_merge(BDY(cur),BDY(m),n);
450: m = NEXT(m);
451: prev = cur; cur = NEXT(cur);
452: } else if ( cur->td > m->td ) {
453: t = NEXT(cur); NEXT(cur) = m; m = t;
454: prev = cur; cur = NEXT(cur);
455: } else {
456: NEXT(prev) = m; m = cur;
457: prev = NEXT(prev); cur = NEXT(prev);
1.3 noro 458: }
459: }
460: if ( !cur )
461: NEXT(prev) = m;
462: return top;
1.1 noro 463: }
464: }
465:
1.19 noro 466: void subd(VL vl,DP p1,DP p2,DP *pr)
1.1 noro 467: {
468: DP t;
469:
470: if ( !p2 )
471: *pr = p1;
472: else {
473: chsgnd(p2,&t); addd(vl,p1,t,pr);
474: }
475: }
476:
1.19 noro 477: void chsgnd(DP p,DP *pr)
1.1 noro 478: {
479: MP m,mr,mr0;
480:
481: if ( !p )
482: *pr = 0;
483: else {
484: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
485: NEXTMP(mr0,mr); chsgnp(C(m),&C(mr)); mr->dl = m->dl;
486: }
487: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
488: if ( *pr )
489: (*pr)->sugar = p->sugar;
490: }
491: }
492:
1.19 noro 493: void muld(VL vl,DP p1,DP p2,DP *pr)
1.1 noro 494: {
1.2 noro 495: if ( ! do_weyl )
496: comm_muld(vl,p1,p2,pr);
497: else
498: weyl_muld(vl,p1,p2,pr);
499: }
500:
1.19 noro 501: void comm_muld(VL vl,DP p1,DP p2,DP *pr)
1.2 noro 502: {
1.1 noro 503: MP m;
504: DP s,t,u;
1.5 noro 505: int i,l,l1;
506: static MP *w;
507: static int wlen;
1.1 noro 508:
509: if ( !p1 || !p2 )
510: *pr = 0;
511: else if ( OID(p1) <= O_P )
512: muldc(vl,p2,(P)p1,pr);
513: else if ( OID(p2) <= O_P )
514: muldc(vl,p1,(P)p2,pr);
515: else {
1.5 noro 516: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
1.4 noro 517: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
1.5 noro 518: if ( l1 < l ) {
519: t = p1; p1 = p2; p2 = t;
520: l = l1;
521: }
522: if ( l > wlen ) {
523: if ( w ) GC_free(w);
524: w = (MP *)MALLOC(l*sizeof(MP));
525: wlen = l;
526: }
1.4 noro 527: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
528: w[i] = m;
529: for ( s = 0, i = l-1; i >= 0; i-- ) {
530: muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
1.1 noro 531: }
1.5 noro 532: bzero(w,l*sizeof(MP));
1.1 noro 533: *pr = s;
534: }
535: }
536:
1.19 noro 537: void muldm(VL vl,DP p,MP m0,DP *pr)
1.1 noro 538: {
539: MP m,mr,mr0;
540: P c;
541: DL d;
542: int n;
543:
544: if ( !p )
545: *pr = 0;
546: else {
547: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
548: m; m = NEXT(m) ) {
549: NEXTMP(mr0,mr);
550: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
551: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
552: else
553: mulp(vl,C(m),c,&C(mr));
554: adddl(n,m->dl,d,&mr->dl);
555: }
556: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
557: if ( *pr )
558: (*pr)->sugar = p->sugar + m0->dl->td;
1.2 noro 559: }
560: }
561:
1.19 noro 562: void weyl_muld(VL vl,DP p1,DP p2,DP *pr)
1.2 noro 563: {
564: MP m;
565: DP s,t,u;
1.4 noro 566: int i,l;
1.5 noro 567: static MP *w;
568: static int wlen;
1.2 noro 569:
570: if ( !p1 || !p2 )
571: *pr = 0;
572: else if ( OID(p1) <= O_P )
573: muldc(vl,p2,(P)p1,pr);
574: else if ( OID(p2) <= O_P )
575: muldc(vl,p1,(P)p2,pr);
576: else {
1.10 noro 577: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
1.5 noro 578: if ( l > wlen ) {
579: if ( w ) GC_free(w);
580: w = (MP *)MALLOC(l*sizeof(MP));
581: wlen = l;
582: }
1.10 noro 583: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
1.4 noro 584: w[i] = m;
585: for ( s = 0, i = l-1; i >= 0; i-- ) {
1.10 noro 586: weyl_muldm(vl,w[i],p2,&t); addd(vl,s,t,&u); s = u;
1.2 noro 587: }
1.5 noro 588: bzero(w,l*sizeof(MP));
1.2 noro 589: *pr = s;
590: }
591: }
592:
1.10 noro 593: /* monomial * polynomial */
594:
1.19 noro 595: void weyl_muldm(VL vl,MP m0,DP p,DP *pr)
1.2 noro 596: {
597: DP r,t,t1;
598: MP m;
1.10 noro 599: DL d0;
600: int n,n2,l,i,j,tlen;
601: static MP *w,*psum;
602: static struct cdl *tab;
1.5 noro 603: static int wlen;
1.10 noro 604: static int rtlen;
1.2 noro 605:
606: if ( !p )
607: *pr = 0;
608: else {
1.4 noro 609: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
1.5 noro 610: if ( l > wlen ) {
611: if ( w ) GC_free(w);
612: w = (MP *)MALLOC(l*sizeof(MP));
613: wlen = l;
614: }
1.4 noro 615: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
616: w[i] = m;
1.10 noro 617:
618: n = NV(p); n2 = n>>1;
619: d0 = m0->dl;
620: for ( i = 0, tlen = 1; i < n2; i++ )
621: tlen *= d0->d[n2+i]+1;
622: if ( tlen > rtlen ) {
623: if ( tab ) GC_free(tab);
624: if ( psum ) GC_free(psum);
625: rtlen = tlen;
626: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
627: psum = (MP *)MALLOC(rtlen*sizeof(MP));
628: }
629: bzero(psum,tlen*sizeof(MP));
630: for ( i = l-1; i >= 0; i-- ) {
631: bzero(tab,tlen*sizeof(struct cdl));
632: weyl_mulmm(vl,m0,w[i],n,tab,tlen);
633: for ( j = 0; j < tlen; j++ ) {
634: if ( tab[j].c ) {
635: NEWMP(m); m->dl = tab[j].d; C(m) = tab[j].c; NEXT(m) = psum[j];
636: psum[j] = m;
637: }
638: }
1.2 noro 639: }
1.10 noro 640: for ( j = tlen-1, r = 0; j >= 0; j-- )
641: if ( psum[j] ) {
642: MKDP(n,psum[j],t); addd(vl,r,t,&t1); r = t1;
643: }
1.2 noro 644: if ( r )
645: r->sugar = p->sugar + m0->dl->td;
646: *pr = r;
647: }
648: }
649:
1.10 noro 650: /* m0 = x0^d0*x1^d1*... * dx0^e0*dx1^e1*... */
651: /* rtab : array of length (e0+1)*(e1+1)*... */
1.2 noro 652:
1.19 noro 653: void weyl_mulmm(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
1.2 noro 654: {
1.19 noro 655: P c,c0,c1;
1.10 noro 656: DL d,d0,d1,dt;
657: int i,j,a,b,k,l,n2,s,min,curlen;
658: struct cdl *p;
659: static Q *ctab;
660: static struct cdl *tab;
1.5 noro 661: static int tablen;
1.10 noro 662: static struct cdl *tmptab;
663: static int tmptablen;
1.2 noro 664:
1.10 noro 665:
666: if ( !m0 || !m1 ) {
667: rtab[0].c = 0;
668: rtab[0].d = 0;
669: return;
670: }
671: c0 = C(m0); c1 = C(m1);
672: mulp(vl,c0,c1,&c);
673: d0 = m0->dl; d1 = m1->dl;
674: n2 = n>>1;
675: curlen = 1;
676: NEWDL(d,n);
677: if ( n & 1 )
678: /* offset of h-degree */
679: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
680: else
681: d->td = 0;
682: rtab[0].c = c;
683: rtab[0].d = d;
684:
685: if ( rtablen > tmptablen ) {
686: if ( tmptab ) GC_free(tmptab);
687: tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
688: tmptablen = rtablen;
689: }
690: for ( i = 0; i < n2; i++ ) {
691: a = d0->d[i]; b = d1->d[n2+i];
692: k = d0->d[n2+i]; l = d1->d[i];
1.20 noro 693:
694: /* degree of xi^a*(Di^k*xi^l)*Di^b */
695: a += l;
696: b += k;
697: s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
698:
1.10 noro 699: if ( !k || !l ) {
700: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
701: if ( p->c ) {
702: dt = p->d;
703: dt->d[i] = a;
704: dt->d[n2+i] = b;
705: dt->td += s;
1.5 noro 706: }
1.10 noro 707: }
708: curlen *= k+1;
709: continue;
710: }
711: if ( k+1 > tablen ) {
712: if ( tab ) GC_free(tab);
713: if ( ctab ) GC_free(ctab);
714: tablen = k+1;
715: tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
716: ctab = (Q *)MALLOC(tablen*sizeof(Q));
717: }
718: /* compute xi^a*(Di^k*xi^l)*Di^b */
719: min = MIN(k,l);
720: mkwc(k,l,ctab);
721: bzero(tab,(k+1)*sizeof(struct cdl));
722: if ( n & 1 )
723: for ( j = 0; j <= min; j++ ) {
724: NEWDL(d,n);
1.20 noro 725: d->d[i] = a-j; d->d[n2+i] = b-j;
1.10 noro 726: d->td = s;
1.20 noro 727: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
1.10 noro 728: tab[j].d = d;
729: tab[j].c = (P)ctab[j];
730: }
731: else
732: for ( j = 0; j <= min; j++ ) {
733: NEWDL(d,n);
1.20 noro 734: d->d[i] = a-j; d->d[n2+i] = b-j;
735: d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
1.10 noro 736: tab[j].d = d;
737: tab[j].c = (P)ctab[j];
738: }
739: bzero(ctab,(min+1)*sizeof(Q));
740: comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
741: curlen *= k+1;
742: bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
743: }
744: }
745:
746: /* direct product of two cdl tables
747: rt[] = [
748: t[0]*t1[0],...,t[n-1]*t1[0],
749: t[0]*t1[1],...,t[n-1]*t1[1],
750: ...
751: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
752: ]
753: */
754:
1.19 noro 755: void comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
1.10 noro 756: {
757: int i,j;
758: struct cdl *p;
759: P c;
760: DL d;
761:
762: bzero(rt,n*n1*sizeof(struct cdl));
763: for ( j = 0, p = rt; j < n1; j++ ) {
764: c = t1[j].c;
765: d = t1[j].d;
766: if ( !c )
767: break;
768: for ( i = 0; i < n; i++, p++ ) {
769: if ( t[i].c ) {
770: mulp(vl,t[i].c,c,&p->c);
771: adddl(nv,t[i].d,d,&p->d);
772: }
1.6 noro 773: }
1.1 noro 774: }
775: }
776:
1.19 noro 777: void muldc(VL vl,DP p,P c,DP *pr)
1.1 noro 778: {
779: MP m,mr,mr0;
780:
781: if ( !p || !c )
782: *pr = 0;
783: else if ( NUM(c) && UNIQ((Q)c) )
784: *pr = p;
785: else if ( NUM(c) && MUNIQ((Q)c) )
786: chsgnd(p,pr);
787: else {
788: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
789: NEXTMP(mr0,mr);
790: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
791: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
792: else
793: mulp(vl,C(m),c,&C(mr));
794: mr->dl = m->dl;
795: }
796: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
797: if ( *pr )
798: (*pr)->sugar = p->sugar;
799: }
800: }
801:
1.19 noro 802: void divsdc(VL vl,DP p,P c,DP *pr)
1.1 noro 803: {
804: MP m,mr,mr0;
805:
806: if ( !c )
807: error("disvsdc : division by 0");
808: else if ( !p )
809: *pr = 0;
810: else {
811: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
812: NEXTMP(mr0,mr); divsp(vl,C(m),c,&C(mr)); mr->dl = m->dl;
813: }
814: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
815: if ( *pr )
816: (*pr)->sugar = p->sugar;
817: }
818: }
819:
1.19 noro 820: void adddl(int n,DL d1,DL d2,DL *dr)
1.1 noro 821: {
822: DL dt;
823: int i;
824:
825: if ( !d1->td )
826: *dr = d2;
827: else if ( !d2->td )
828: *dr = d1;
829: else {
830: *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
831: dt->td = d1->td + d2->td;
832: for ( i = 0; i < n; i++ )
833: dt->d[i] = d1->d[i]+d2->d[i];
834: }
1.11 noro 835: }
836:
837: /* d1 += d2 */
838:
1.19 noro 839: void adddl_destructive(int n,DL d1,DL d2)
1.11 noro 840: {
841: int i;
842:
843: d1->td += d2->td;
844: for ( i = 0; i < n; i++ )
845: d1->d[i] += d2->d[i];
1.1 noro 846: }
847:
1.19 noro 848: int compd(VL vl,DP p1,DP p2)
1.1 noro 849: {
850: int n,t;
851: MP m1,m2;
852:
853: if ( !p1 )
854: return p2 ? -1 : 0;
855: else if ( !p2 )
856: return 1;
857: else {
858: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
859: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
860: if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
861: (t = compp(vl,C(m1),C(m2)) ) )
862: return t;
863: if ( m1 )
864: return 1;
865: else if ( m2 )
866: return -1;
867: else
868: return 0;
869: }
870: }
871:
1.19 noro 872: int cmpdl_lex(int n,DL d1,DL d2)
1.1 noro 873: {
874: int i;
875:
876: for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
877: return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
878: }
879:
1.19 noro 880: int cmpdl_revlex(int n,DL d1,DL d2)
1.1 noro 881: {
882: int i;
883:
884: for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
885: return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
886: }
887:
1.19 noro 888: int cmpdl_gradlex(int n,DL d1,DL d2)
1.1 noro 889: {
890: if ( d1->td > d2->td )
891: return 1;
892: else if ( d1->td < d2->td )
893: return -1;
894: else
895: return cmpdl_lex(n,d1,d2);
896: }
897:
1.19 noro 898: int cmpdl_revgradlex(int n,DL d1,DL d2)
1.1 noro 899: {
1.7 noro 900: register int i;
901: register int *p1,*p2;
902:
1.1 noro 903: if ( d1->td > d2->td )
904: return 1;
905: else if ( d1->td < d2->td )
906: return -1;
1.7 noro 907: else {
908: for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
909: i >= 0 && *p1 == *p2; i--, p1--, p2-- );
910: return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
911: }
1.1 noro 912: }
913:
1.19 noro 914: int cmpdl_blex(int n,DL d1,DL d2)
1.1 noro 915: {
916: int c;
917:
918: if ( c = cmpdl_lex(n-1,d1,d2) )
919: return c;
920: else {
921: c = d1->d[n-1] - d2->d[n-1];
922: return c > 0 ? 1 : c < 0 ? -1 : 0;
923: }
924: }
925:
1.19 noro 926: int cmpdl_bgradlex(int n,DL d1,DL d2)
1.1 noro 927: {
928: int e1,e2,c;
929:
930: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
931: if ( e1 > e2 )
932: return 1;
933: else if ( e1 < e2 )
934: return -1;
935: else {
936: c = cmpdl_lex(n-1,d1,d2);
937: if ( c )
938: return c;
939: else
940: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
941: }
942: }
943:
1.19 noro 944: int cmpdl_brevgradlex(int n,DL d1,DL d2)
1.1 noro 945: {
946: int e1,e2,c;
947:
948: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
949: if ( e1 > e2 )
950: return 1;
951: else if ( e1 < e2 )
952: return -1;
953: else {
954: c = cmpdl_revlex(n-1,d1,d2);
955: if ( c )
956: return c;
957: else
958: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
959: }
960: }
961:
1.19 noro 962: int cmpdl_brevrev(int n,DL d1,DL d2)
1.1 noro 963: {
964: int e1,e2,f1,f2,c,i;
965:
966: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
967: e1 += d1->d[i]; e2 += d2->d[i];
968: }
969: f1 = d1->td - e1; f2 = d2->td - e2;
970: if ( e1 > e2 )
971: return 1;
972: else if ( e1 < e2 )
973: return -1;
974: else {
975: c = cmpdl_revlex(dp_nelim,d1,d2);
976: if ( c )
977: return c;
978: else if ( f1 > f2 )
979: return 1;
980: else if ( f1 < f2 )
981: return -1;
982: else {
983: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
984: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
985: }
986: }
987: }
988:
1.19 noro 989: int cmpdl_bgradrev(int n,DL d1,DL d2)
1.1 noro 990: {
991: int e1,e2,f1,f2,c,i;
992:
993: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
994: e1 += d1->d[i]; e2 += d2->d[i];
995: }
996: f1 = d1->td - e1; f2 = d2->td - e2;
997: if ( e1 > e2 )
998: return 1;
999: else if ( e1 < e2 )
1000: return -1;
1001: else {
1002: c = cmpdl_lex(dp_nelim,d1,d2);
1003: if ( c )
1004: return c;
1005: else if ( f1 > f2 )
1006: return 1;
1007: else if ( f1 < f2 )
1008: return -1;
1009: else {
1010: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1011: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1012: }
1013: }
1014: }
1015:
1.19 noro 1016: int cmpdl_blexrev(int n,DL d1,DL d2)
1.1 noro 1017: {
1018: int e1,e2,f1,f2,c,i;
1019:
1020: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1021: e1 += d1->d[i]; e2 += d2->d[i];
1022: }
1023: f1 = d1->td - e1; f2 = d2->td - e2;
1024: c = cmpdl_lex(dp_nelim,d1,d2);
1025: if ( c )
1026: return c;
1027: else if ( f1 > f2 )
1028: return 1;
1029: else if ( f1 < f2 )
1030: return -1;
1031: else {
1032: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1033: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1034: }
1035: }
1036:
1.19 noro 1037: int cmpdl_elim(int n,DL d1,DL d2)
1.1 noro 1038: {
1039: int e1,e2,i;
1040:
1041: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1042: e1 += d1->d[i]; e2 += d2->d[i];
1043: }
1044: if ( e1 > e2 )
1045: return 1;
1046: else if ( e1 < e2 )
1047: return -1;
1048: else
1049: return cmpdl_revgradlex(n,d1,d2);
1.12 noro 1050: }
1051:
1.19 noro 1052: int cmpdl_weyl_elim(int n,DL d1,DL d2)
1.12 noro 1053: {
1054: int e1,e2,i;
1055:
1056: for ( i = 1, e1 = 0, e2 = 0; i <= dp_nelim; i++ ) {
1057: e1 += d1->d[n-i]; e2 += d2->d[n-i];
1058: }
1059: if ( e1 > e2 )
1060: return 1;
1061: else if ( e1 < e2 )
1062: return -1;
1063: else if ( d1->td > d2->td )
1064: return 1;
1065: else if ( d1->td < d2->td )
1066: return -1;
1067: else return -cmpdl_revlex(n,d1,d2);
1.13 noro 1068: }
1069:
1070: /*
1071: a special ordering
1072: 1. total order
1073: 2. (-w,w) for the first 2*m variables
1074: 3. DRL for the first 2*m variables
1075: */
1076:
1.20 noro 1077: extern int *current_weyl_weight_vector;
1.13 noro 1078:
1.19 noro 1079: int cmpdl_homo_ww_drl(int n,DL d1,DL d2)
1.13 noro 1080: {
1081: int e1,e2,m,i;
1082: int *p1,*p2;
1083:
1084: if ( d1->td > d2->td )
1085: return 1;
1086: else if ( d1->td < d2->td )
1087: return -1;
1088:
1089: m = n>>1;
1.21 ! noro 1090: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
! 1091: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
! 1092: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1.13 noro 1093: }
1094: if ( e1 > e2 )
1095: return 1;
1096: else if ( e1 < e2 )
1097: return -1;
1098:
1099: e1 = d1->td - d1->d[n-1];
1100: e2 = d2->td - d2->d[n-1];
1101: if ( e1 > e2 )
1102: return 1;
1103: else if ( e1 < e2 )
1104: return -1;
1105:
1106: for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
1107: i >= 0 && *p1 == *p2; i--, p1--, p2-- );
1108: return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
1.21 ! noro 1109: }
! 1110:
! 1111: int cmpdl_drl_zigzag(int n,DL d1,DL d2)
! 1112: {
! 1113: int i,t,m;
! 1114: int *p1,*p2;
! 1115:
! 1116: if ( d1->td > d2->td )
! 1117: return 1;
! 1118: else if ( d1->td < d2->td )
! 1119: return -1;
! 1120: else {
! 1121: m = n>>1;
! 1122: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
! 1123: if ( t = p1[m+i] - p2[m+i] ) return t > 0 ? -1 : 1;
! 1124: if ( t = p1[i] - p2[i] ) return t > 0 ? -1 : 1;
! 1125: }
! 1126: return 0;
! 1127: }
! 1128: }
! 1129:
! 1130: int cmpdl_homo_ww_drl_zigzag(int n,DL d1,DL d2)
! 1131: {
! 1132: int e1,e2,m,i,t;
! 1133: int *p1,*p2;
! 1134:
! 1135: if ( d1->td > d2->td )
! 1136: return 1;
! 1137: else if ( d1->td < d2->td )
! 1138: return -1;
! 1139:
! 1140: m = n>>1;
! 1141: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
! 1142: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
! 1143: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
! 1144: }
! 1145: if ( e1 > e2 )
! 1146: return 1;
! 1147: else if ( e1 < e2 )
! 1148: return -1;
! 1149:
! 1150: e1 = d1->td - d1->d[n-1];
! 1151: e2 = d2->td - d2->d[n-1];
! 1152: if ( e1 > e2 )
! 1153: return 1;
! 1154: else if ( e1 < e2 )
! 1155: return -1;
! 1156:
! 1157: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
! 1158: if ( t = p1[m+i] - p2[m+i] ) return t > 0 ? -1 : 1;
! 1159: if ( t = p1[i] - p2[i] ) return t > 0 ? -1 : 1;
! 1160: }
! 1161: return 0;
1.1 noro 1162: }
1163:
1.19 noro 1164: int cmpdl_order_pair(int n,DL d1,DL d2)
1.1 noro 1165: {
1166: int e1,e2,i,j,l;
1167: int *t1,*t2;
1.20 noro 1168: int len,head;
1.1 noro 1169: struct order_pair *pair;
1170:
1171: len = dp_current_spec.ord.block.length;
1172: pair = dp_current_spec.ord.block.order_pair;
1173:
1.20 noro 1174: head = 0;
1.1 noro 1175: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
1176: l = pair[i].length;
1177: switch ( pair[i].order ) {
1178: case 0:
1179: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1.20 noro 1180: e1 += MUL_WEIGHT(t1[j],head+j);
1181: e2 += MUL_WEIGHT(t2[j],head+j);
1.1 noro 1182: }
1183: if ( e1 > e2 )
1184: return 1;
1185: else if ( e1 < e2 )
1186: return -1;
1187: else {
1188: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
1189: if ( j >= 0 )
1190: return t1[j] < t2[j] ? 1 : -1;
1191: }
1192: break;
1193: case 1:
1194: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1.20 noro 1195: e1 += MUL_WEIGHT(t1[j],head+j);
1196: e2 += MUL_WEIGHT(t2[j],head+j);
1.1 noro 1197: }
1198: if ( e1 > e2 )
1199: return 1;
1200: else if ( e1 < e2 )
1201: return -1;
1202: else {
1203: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1204: if ( j < l )
1205: return t1[j] > t2[j] ? 1 : -1;
1206: }
1207: break;
1208: case 2:
1209: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1210: if ( j < l )
1211: return t1[j] > t2[j] ? 1 : -1;
1212: break;
1213: default:
1214: error("cmpdl_order_pair : invalid order"); break;
1215: }
1.20 noro 1216: t1 += l; t2 += l; head += l;
1.1 noro 1217: }
1218: return 0;
1219: }
1220:
1.19 noro 1221: int cmpdl_matrix(int n,DL d1,DL d2)
1.1 noro 1222: {
1223: int *v,*w,*t1,*t2;
1224: int s,i,j,len;
1225: int **matrix;
1226:
1227: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
1228: w[i] = t1[i]-t2[i];
1229: len = dp_current_spec.ord.matrix.row;
1230: matrix = dp_current_spec.ord.matrix.matrix;
1231: for ( j = 0; j < len; j++ ) {
1232: v = matrix[j];
1233: for ( i = 0, s = 0; i < n; i++ )
1234: s += v[i]*w[i];
1235: if ( s > 0 )
1236: return 1;
1237: else if ( s < 0 )
1238: return -1;
1239: }
1240: return 0;
1241: }