Annotation of OpenXM_contrib2/asir2000/builtin/dp.c, Revision 1.7
1.5 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.6 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.5 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.7 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp.c,v 1.6 2000/08/22 05:03:57 noro Exp $
1.5 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "base.h"
52: #include "parse.h"
53:
54: extern int dp_fcoeffs;
55:
56: void Pdp_ord(), Pdp_ptod(), Pdp_dtop();
57: void Pdp_ptozp(), Pdp_ptozp2(), Pdp_red(), Pdp_red2(), Pdp_lcm(), Pdp_redble();
58: void Pdp_sp(), Pdp_hm(), Pdp_ht(), Pdp_hc(), Pdp_rest(), Pdp_td(), Pdp_sugar();
59: void Pdp_cri1(),Pdp_cri2(),Pdp_subd(),Pdp_mod(),Pdp_red_mod(),Pdp_tdiv();
60: void Pdp_prim(),Pdp_red_coef(),Pdp_mag(),Pdp_set_kara(),Pdp_rat();
61: void Pdp_nf(),Pdp_true_nf(),Pdp_nf_ptozp();
62: void Pdp_nf_mod(),Pdp_true_nf_mod();
63: void Pdp_criB(),Pdp_nelim();
64: void Pdp_minp(),Pdp_nf_demand(),Pdp_sp_mod();
65: void Pdp_homo(),Pdp_dehomo();
66: void Pdp_gr_mod_main();
67: void Pdp_gr_main(),Pdp_gr_hm_main(),Pdp_gr_d_main(),Pdp_gr_flags();
68: void Pdp_f4_main(),Pdp_f4_mod_main();
69: void Pdp_gr_print();
70:
71: struct ftab dp_tab[] = {
72: {"dp_ord",Pdp_ord,-1},
73: {"dp_ptod",Pdp_ptod,2},
74: {"dp_dtop",Pdp_dtop,2},
75: {"dp_ptozp",Pdp_ptozp,1},
76: {"dp_ptozp2",Pdp_ptozp2,2},
77: {"dp_prim",Pdp_prim,1},
78: {"dp_redble",Pdp_redble,2},
79: {"dp_subd",Pdp_subd,2},
80: {"dp_red",Pdp_red,3},
81: {"dp_red_mod",Pdp_red_mod,4},
82: {"dp_sp",Pdp_sp,2},
83: {"dp_sp_mod",Pdp_sp_mod,3},
84: {"dp_lcm",Pdp_lcm,2},
85: {"dp_hm",Pdp_hm,1},
86: {"dp_ht",Pdp_ht,1},
87: {"dp_hc",Pdp_hc,1},
88: {"dp_rest",Pdp_rest,1},
89: {"dp_td",Pdp_td,1},
90: {"dp_sugar",Pdp_sugar,1},
91: {"dp_cri1",Pdp_cri1,2},
92: {"dp_cri2",Pdp_cri2,2},
93: {"dp_criB",Pdp_criB,3},
94: {"dp_minp",Pdp_minp,2},
95: {"dp_mod",Pdp_mod,3},
96: {"dp_rat",Pdp_rat,1},
97: {"dp_tdiv",Pdp_tdiv,2},
98: {"dp_red_coef",Pdp_red_coef,2},
99: {"dp_nelim",Pdp_nelim,-1},
100: {"dp_mag",Pdp_mag,1},
101: {"dp_set_kara",Pdp_set_kara,-1},
102: {"dp_nf",Pdp_nf,4},
103: {"dp_true_nf",Pdp_true_nf,4},
104: {"dp_nf_ptozp",Pdp_nf_ptozp,5},
105: {"dp_nf_demand",Pdp_nf_demand,5},
106: {"dp_nf_mod",Pdp_nf_mod,5},
107: {"dp_true_nf_mod",Pdp_true_nf_mod,5},
108: {"dp_homo",Pdp_homo,1},
109: {"dp_dehomo",Pdp_dehomo,1},
110: {"dp_gr_main",Pdp_gr_main,5},
111: /* {"dp_gr_hm_main",Pdp_gr_hm_main,5}, */
112: /* {"dp_gr_d_main",Pdp_gr_d_main,6}, */
113: {"dp_gr_mod_main",Pdp_gr_mod_main,5},
114: {"dp_f4_main",Pdp_f4_main,3},
115: {"dp_f4_mod_main",Pdp_f4_mod_main,4},
116: {"dp_gr_flags",Pdp_gr_flags,-1},
117: {"dp_gr_print",Pdp_gr_print,-1},
118: {0,0,0},
119: };
120:
121: extern int dp_nelim;
122: extern int dp_order_pair_length;
123: extern struct order_pair *dp_order_pair;
124: extern struct order_spec dp_current_spec;
125:
126: void Pdp_ord(arg,rp)
127: NODE arg;
128: Obj *rp;
129: {
130: struct order_spec spec;
131:
132: if ( !arg )
133: *rp = dp_current_spec.obj;
134: else if ( !create_order_spec((Obj)ARG0(arg),&spec) )
135: error("dp_ord : invalid order specification");
136: else {
137: initd(&spec); *rp = spec.obj;
138: }
139: }
140:
141: int create_order_spec(obj,spec)
142: Obj obj;
143: struct order_spec *spec;
144: {
145: int i,j,n,s,row,col;
146: struct order_pair *l;
147: NODE node,t,tn;
148: MAT m;
149: pointer **b;
150: int **w;
151:
152: if ( !obj || NUM(obj) ) {
153: spec->id = 0; spec->obj = obj;
154: spec->ord.simple = QTOS((Q)obj);
155: return 1;
156: } else if ( OID(obj) == O_LIST ) {
157: node = BDY((LIST)obj);
158: for ( n = 0, t = node; t; t = NEXT(t), n++ );
159: l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair));
160: for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) {
161: tn = BDY((LIST)BDY(t)); l[i].order = QTOS((Q)BDY(tn));
162: tn = NEXT(tn); l[i].length = QTOS((Q)BDY(tn));
163: s += l[i].length;
164: }
165: spec->id = 1; spec->obj = obj;
166: spec->ord.block.order_pair = l;
167: spec->ord.block.length = n; spec->nv = s;
168: return 1;
169: } else if ( OID(obj) == O_MAT ) {
170: m = (MAT)obj; row = m->row; col = m->col; b = BDY(m);
171: w = almat(row,col);
172: for ( i = 0; i < row; i++ )
173: for ( j = 0; j < col; j++ )
174: w[i][j] = QTOS((Q)b[i][j]);
175: spec->id = 2; spec->obj = obj;
176: spec->nv = col; spec->ord.matrix.row = row;
177: spec->ord.matrix.matrix = w;
178: return 1;
179: } else
180: return 0;
181: }
182:
183: void homogenize_order(old,n,new)
184: struct order_spec *old,*new;
185: int n;
186: {
187: struct order_pair *l;
188: int length,nv,row,i,j;
189: int **newm,**oldm;
190:
191: switch ( old->id ) {
192: case 0:
193: switch ( old->ord.simple ) {
194: case 0:
195: new->id = 0; new->ord.simple = 0; break;
196: case 1:
197: l = (struct order_pair *)
198: MALLOC_ATOMIC(2*sizeof(struct order_pair));
199: l[0].length = n; l[0].order = 1;
200: l[1].length = 1; l[1].order = 2;
201: new->id = 1;
202: new->ord.block.order_pair = l;
203: new->ord.block.length = 2; new->nv = n+1;
204: break;
205: case 2:
206: new->id = 0; new->ord.simple = 1; break;
207: case 3: case 4: case 5:
208: new->id = 0; new->ord.simple = old->ord.simple+3;
209: dp_nelim = n-1; break;
210: case 6: case 7: case 8: case 9:
211: new->id = 0; new->ord.simple = old->ord.simple; break;
212: default:
213: error("homogenize_order : invalid input");
214: }
215: break;
216: case 1:
217: length = old->ord.block.length;
218: l = (struct order_pair *)
219: MALLOC_ATOMIC((length+1)*sizeof(struct order_pair));
220: bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair));
221: l[length].order = 2; l[length].length = 1;
222: new->id = 1; new->nv = n+1;
223: new->ord.block.order_pair = l;
224: new->ord.block.length = length+1;
225: break;
226: case 2:
227: nv = old->nv; row = old->ord.matrix.row;
228: oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1);
229: for ( i = 0; i <= nv; i++ )
230: newm[0][i] = 1;
231: for ( i = 0; i < row; i++ ) {
232: for ( j = 0; j < nv; j++ )
233: newm[i+1][j] = oldm[i][j];
234: newm[i+1][j] = 0;
235: }
236: new->id = 2; new->nv = nv+1;
237: new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm;
238: break;
239: default:
240: error("homogenize_order : invalid input");
241: }
242: }
243:
244: void Pdp_ptod(arg,rp)
245: NODE arg;
246: DP *rp;
247: {
248: NODE n;
249: VL vl,tvl;
250:
251: asir_assert(ARG0(arg),O_P,"dp_ptod");
252: asir_assert(ARG1(arg),O_LIST,"dp_ptod");
253: for ( vl = 0, n = BDY((LIST)ARG1(arg)); n; n = NEXT(n) ) {
254: if ( !vl ) {
255: NEWVL(vl); tvl = vl;
256: } else {
257: NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
258: }
259: VR(tvl) = VR((P)BDY(n));
260: }
261: if ( vl )
262: NEXT(tvl) = 0;
263: ptod(CO,vl,(P)ARG0(arg),rp);
264: }
265:
266: void Pdp_dtop(arg,rp)
267: NODE arg;
268: P *rp;
269: {
270: NODE n;
271: VL vl,tvl;
272:
273: asir_assert(ARG0(arg),O_DP,"dp_dtop");
274: asir_assert(ARG1(arg),O_LIST,"dp_dtop");
275: for ( vl = 0, n = BDY((LIST)ARG1(arg)); n; n = NEXT(n) ) {
276: if ( !vl ) {
277: NEWVL(vl); tvl = vl;
278: } else {
279: NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
280: }
281: VR(tvl) = VR((P)BDY(n));
282: }
283: if ( vl )
284: NEXT(tvl) = 0;
285: dtop(CO,vl,(DP)ARG0(arg),rp);
286: }
287:
288: extern LIST Dist;
289:
290: void Pdp_ptozp(arg,rp)
291: NODE arg;
292: DP *rp;
293: {
294: asir_assert(ARG0(arg),O_DP,"dp_ptozp");
295: if ( Dist )
296: dp_ptozp_d(BDY(Dist),length(BDY(Dist)),(DP)ARG0(arg),rp);
297: else
298: dp_ptozp((DP)ARG0(arg),rp);
299: }
300:
301: void Pdp_ptozp2(arg,rp)
302: NODE arg;
303: LIST *rp;
304: {
305: DP p0,p1,h,r;
306: NODE n0;
307:
308: p0 = (DP)ARG0(arg); p1 = (DP)ARG1(arg);
309: asir_assert(p0,O_DP,"dp_ptozp2");
310: asir_assert(p1,O_DP,"dp_ptozp2");
311: if ( Dist )
312: dp_ptozp2_d(BDY(Dist),length(BDY(Dist)),p0,p1,&h,&r);
313: else
314: dp_ptozp2(p0,p1,&h,&r);
315: NEWNODE(n0); BDY(n0) = (pointer)h;
316: NEWNODE(NEXT(n0)); BDY(NEXT(n0)) = (pointer)r;
317: NEXT(NEXT(n0)) = 0;
318: MKLIST(*rp,n0);
319: }
320:
321: void Pdp_prim(arg,rp)
322: NODE arg;
323: DP *rp;
324: {
325: DP t;
326:
327: asir_assert(ARG0(arg),O_DP,"dp_prim");
328: dp_prim((DP)ARG0(arg),&t); dp_ptozp(t,rp);
329: }
330:
331: extern int NoGCD;
332:
333: void dp_prim(p,rp)
334: DP p,*rp;
335: {
336: P t,g;
337: DP p1;
338: MP m,mr,mr0;
339: int i,n;
340: P *w;
341: Q *c;
342: Q dvr;
343:
344: if ( !p )
345: *rp = 0;
346: else if ( dp_fcoeffs )
347: *rp = p;
348: else if ( NoGCD )
349: dp_ptozp(p,rp);
350: else {
351: dp_ptozp(p,&p1); p = p1;
352: for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
353: if ( n == 1 ) {
354: m = BDY(p);
355: NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
356: MKDP(p->nv,mr,*rp); (*rp)->sugar = p->sugar;
357: return;
358: }
359: w = (P *)ALLOCA(n*sizeof(P));
360: c = (Q *)ALLOCA(n*sizeof(Q));
361: for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ )
362: if ( NUM(m->c) ) {
363: c[i] = (Q)m->c; w[i] = (P)ONE;
364: } else
365: ptozp(m->c,1,&c[i],&w[i]);
366: qltozl(c,n,&dvr); heu_nezgcdnpz(CO,w,n,&t); mulp(CO,t,(P)dvr,&g);
367: if ( NUM(g) )
368: *rp = p;
369: else {
370: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
371: NEXTMP(mr0,mr); divsp(CO,m->c,g,&mr->c); mr->dl = m->dl;
372: }
373: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
374: }
375: }
376: }
377:
378: void heu_nezgcdnpz(vl,pl,m,pr)
379: VL vl;
380: P *pl,*pr;
381: int m;
382: {
383: int i,r;
384: P gcd,t,s1,s2,u;
385: Q rq;
386:
387: while ( 1 ) {
388: for ( i = 0, s1 = 0; i < m; i++ ) {
389: r = random(); UTOQ(r,rq);
390: mulp(vl,pl[i],(P)rq,&t); addp(vl,s1,t,&u); s1 = u;
391: }
392: for ( i = 0, s2 = 0; i < m; i++ ) {
393: r = random(); UTOQ(r,rq);
394: mulp(vl,pl[i],(P)rq,&t); addp(vl,s2,t,&u); s2 = u;
395: }
396: ezgcdp(vl,s1,s2,&gcd);
397: for ( i = 0; i < m; i++ ) {
398: if ( !divtpz(vl,pl[i],gcd,&t) )
399: break;
400: }
401: if ( i == m )
402: break;
403: }
404: *pr = gcd;
405: }
406:
407: void dp_prim_mod(p,mod,rp)
408: int mod;
409: DP p,*rp;
410: {
411: P t,g;
412: MP m,mr,mr0;
413:
414: if ( !p )
415: *rp = 0;
416: else if ( NoGCD )
417: *rp = p;
418: else {
419: for ( m = BDY(p), g = m->c, m = NEXT(m); m; m = NEXT(m) ) {
420: gcdprsmp(CO,mod,g,m->c,&t); g = t;
421: }
422: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
423: NEXTMP(mr0,mr); divsmp(CO,mod,m->c,g,&mr->c); mr->dl = m->dl;
424: }
425: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
426: }
427: }
428:
429: void Pdp_mod(arg,rp)
430: NODE arg;
431: DP *rp;
432: {
433: DP p;
434: int mod;
435: NODE subst;
436:
437: asir_assert(ARG0(arg),O_DP,"dp_mod");
438: asir_assert(ARG1(arg),O_N,"dp_mod");
439: asir_assert(ARG2(arg),O_LIST,"dp_mod");
440: p = (DP)ARG0(arg); mod = QTOS((Q)ARG1(arg));
441: subst = BDY((LIST)ARG2(arg));
442: dp_mod(p,mod,subst,rp);
443: }
444:
445: void Pdp_rat(arg,rp)
446: NODE arg;
447: DP *rp;
448: {
449: asir_assert(ARG0(arg),O_DP,"dp_rat");
450: dp_rat((DP)ARG0(arg),rp);
451: }
452:
453: void dp_mod(p,mod,subst,rp)
454: DP p;
455: int mod;
456: NODE subst;
457: DP *rp;
458: {
459: MP m,mr,mr0;
460: P t,s,s1;
461: V v;
462: NODE tn;
463:
464: if ( !p )
465: *rp = 0;
466: else {
467: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
468: for ( tn = subst, s = m->c; tn; tn = NEXT(tn) ) {
469: v = VR((P)BDY(tn)); tn = NEXT(tn);
470: substp(CO,s,v,(P)BDY(tn),&s1); s = s1;
471: }
472: ptomp(mod,s,&t);
473: if ( t ) {
474: NEXTMP(mr0,mr); mr->c = t; mr->dl = m->dl;
475: }
476: }
477: if ( mr0 ) {
478: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
479: } else
480: *rp = 0;
481: }
482: }
483:
484: void dp_rat(p,rp)
485: DP p;
486: DP *rp;
487: {
488: MP m,mr,mr0;
489:
490: if ( !p )
491: *rp = 0;
492: else {
493: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
494: NEXTMP(mr0,mr); mptop(m->c,&mr->c); mr->dl = m->dl;
495: }
496: if ( mr0 ) {
497: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
498: } else
499: *rp = 0;
500: }
501: }
502:
503: void Pdp_nf(arg,rp)
504: NODE arg;
505: DP *rp;
506: {
507: NODE b;
508: DP *ps;
509: DP g;
510: int full;
511:
512: asir_assert(ARG0(arg),O_LIST,"dp_nf");
513: asir_assert(ARG1(arg),O_DP,"dp_nf");
514: asir_assert(ARG2(arg),O_VECT,"dp_nf");
515: asir_assert(ARG3(arg),O_N,"dp_nf");
516: if ( !(g = (DP)ARG1(arg)) ) {
517: *rp = 0; return;
518: }
519: b = BDY((LIST)ARG0(arg)); ps = (DP *)BDY((VECT)ARG2(arg));
520: full = (Q)ARG3(arg) ? 1 : 0;
521: dp_nf(b,g,ps,full,rp);
522: }
523:
524: void Pdp_true_nf(arg,rp)
525: NODE arg;
526: LIST *rp;
527: {
528: NODE b,n;
529: DP *ps;
530: DP g;
531: DP nm;
532: P dn;
533: int full;
534:
535: asir_assert(ARG0(arg),O_LIST,"dp_true_nf");
536: asir_assert(ARG1(arg),O_DP,"dp_true_nf");
537: asir_assert(ARG2(arg),O_VECT,"dp_true_nf");
538: asir_assert(ARG3(arg),O_N,"dp_nf");
539: if ( !(g = (DP)ARG1(arg)) ) {
540: nm = 0; dn = (P)ONE;
541: } else {
542: b = BDY((LIST)ARG0(arg)); ps = (DP *)BDY((VECT)ARG2(arg));
543: full = (Q)ARG3(arg) ? 1 : 0;
544: dp_true_nf(b,g,ps,full,&nm,&dn);
545: }
546: NEWNODE(n); BDY(n) = (pointer)nm;
547: NEWNODE(NEXT(n)); BDY(NEXT(n)) = (pointer)dn;
548: NEXT(NEXT(n)) = 0; MKLIST(*rp,n);
549: }
550:
551: void dp_nf(b,g,ps,full,rp)
552: NODE b;
553: DP g;
554: DP *ps;
555: int full;
556: DP *rp;
557: {
1.4 noro 558: DP u,p,d,s,t,dmy1;
1.1 noro 559: P dmy;
560: NODE l;
561: MP m,mr;
562: int i,n;
563: int *wb;
564: int sugar,psugar;
565:
566: if ( !g ) {
567: *rp = 0; return;
568: }
569: for ( n = 0, l = b; l; l = NEXT(l), n++ );
570: wb = (int *)ALLOCA(n*sizeof(int));
571: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
572: wb[i] = QTOS((Q)BDY(l));
573: sugar = g->sugar;
574: for ( d = 0; g; ) {
575: for ( u = 0, i = 0; i < n; i++ ) {
576: if ( dp_redble(g,p = ps[wb[i]]) ) {
1.4 noro 577: dp_red(d,g,p,&t,&u,&dmy,&dmy1);
1.1 noro 578: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
579: sugar = MAX(sugar,psugar);
580: if ( !u ) {
581: if ( d )
582: d->sugar = sugar;
583: *rp = d; return;
584: }
585: d = t;
586: break;
587: }
588: }
589: if ( u )
590: g = u;
591: else if ( !full ) {
592: if ( g ) {
593: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
594: }
595: *rp = g; return;
596: } else {
597: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
598: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
599: addd(CO,d,t,&s); d = s;
600: dp_rest(g,&t); g = t;
601: }
602: }
603: if ( d )
604: d->sugar = sugar;
605: *rp = d;
606: }
607:
608: void dp_true_nf(b,g,ps,full,rp,dnp)
609: NODE b;
610: DP g;
611: DP *ps;
612: int full;
613: DP *rp;
614: P *dnp;
615: {
1.4 noro 616: DP u,p,d,s,t,dmy;
1.1 noro 617: NODE l;
618: MP m,mr;
619: int i,n;
620: int *wb;
621: int sugar,psugar;
622: P dn,tdn,tdn1;
623:
624: dn = (P)ONE;
625: if ( !g ) {
626: *rp = 0; *dnp = dn; return;
627: }
628: for ( n = 0, l = b; l; l = NEXT(l), n++ );
629: wb = (int *)ALLOCA(n*sizeof(int));
630: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
631: wb[i] = QTOS((Q)BDY(l));
632: sugar = g->sugar;
633: for ( d = 0; g; ) {
634: for ( u = 0, i = 0; i < n; i++ ) {
635: if ( dp_redble(g,p = ps[wb[i]]) ) {
1.4 noro 636: dp_red(d,g,p,&t,&u,&tdn,&dmy);
1.1 noro 637: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
638: sugar = MAX(sugar,psugar);
639: if ( !u ) {
640: if ( d )
641: d->sugar = sugar;
642: *rp = d; *dnp = dn; return;
643: } else {
644: d = t;
645: mulp(CO,dn,tdn,&tdn1); dn = tdn1;
646: }
647: break;
648: }
649: }
650: if ( u )
651: g = u;
652: else if ( !full ) {
653: if ( g ) {
654: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
655: }
656: *rp = g; *dnp = dn; return;
657: } else {
658: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
659: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
660: addd(CO,d,t,&s); d = s;
661: dp_rest(g,&t); g = t;
662: }
663: }
664: if ( d )
665: d->sugar = sugar;
666: *rp = d; *dnp = dn;
667: }
668:
669: #define HMAG(p) (p_mag(BDY(p)->c))
670:
671: void Pdp_nf_ptozp(arg,rp)
672: NODE arg;
673: DP *rp;
674: {
675: NODE b;
676: DP g;
677: DP *ps;
678: int full,multiple;
679:
680: asir_assert(ARG0(arg),O_LIST,"dp_nf_ptozp");
681: asir_assert(ARG1(arg),O_DP,"dp_nf_ptozp");
682: asir_assert(ARG2(arg),O_VECT,"dp_nf_ptozp");
683: asir_assert(ARG3(arg),O_N,"dp_nf_ptozp");
684: asir_assert(ARG4(arg),O_N,"dp_nf_ptozp");
685: if ( !(g = (DP)ARG1(arg)) ) {
686: *rp = 0; return;
687: }
688: b = BDY((LIST)ARG0(arg)); ps = (DP *)BDY((VECT)ARG2(arg));
689: full = (Q)ARG3(arg) ? 1 : 0;
690: multiple = QTOS((Q)ARG4(arg));
691: dp_nf_ptozp(b,g,ps,full,multiple,rp);
692: }
693:
694: void dp_nf_ptozp(b,g,ps,full,multiple,rp)
695: NODE b;
696: DP g;
697: DP *ps;
698: int full,multiple;
699: DP *rp;
700: {
1.4 noro 701: DP u,p,d,s,t,dmy1;
1.1 noro 702: P dmy;
703: NODE l;
704: MP m,mr;
705: int i,n;
706: int *wb;
707: int hmag;
708: int sugar,psugar;
709:
710: if ( !g ) {
711: *rp = 0; return;
712: }
713: for ( n = 0, l = b; l; l = NEXT(l), n++ );
714: wb = (int *)ALLOCA(n*sizeof(int));
715: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
716: wb[i] = QTOS((Q)BDY(l));
717: hmag = multiple*HMAG(g);
718: sugar = g->sugar;
719: for ( d = 0; g; ) {
720: for ( u = 0, i = 0; i < n; i++ ) {
721: if ( dp_redble(g,p = ps[wb[i]]) ) {
1.4 noro 722: dp_red(d,g,p,&t,&u,&dmy,&dmy1);
1.1 noro 723: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
724: sugar = MAX(sugar,psugar);
725: if ( !u ) {
726: if ( d )
727: d->sugar = sugar;
728: *rp = d; return;
729: }
730: d = t;
731: break;
732: }
733: }
734: if ( u ) {
735: g = u;
736: if ( d ) {
737: if ( HMAG(d) > hmag ) {
738: dp_ptozp2(d,g,&t,&u); d = t; g = u;
739: hmag = multiple*HMAG(d);
740: }
741: } else {
742: if ( HMAG(g) > hmag ) {
743: dp_ptozp(g,&t); g = t;
744: hmag = multiple*HMAG(g);
745: }
746: }
747: }
748: else if ( !full ) {
749: if ( g ) {
750: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
751: }
752: *rp = g; return;
753: } else {
754: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
755: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
756: addd(CO,d,t,&s); d = s;
757: dp_rest(g,&t); g = t;
758:
759: }
760: }
761: if ( d )
762: d->sugar = sugar;
763: *rp = d;
764: }
765:
766: void Pdp_nf_demand(arg,rp)
767: NODE arg;
768: DP *rp;
769: {
1.4 noro 770: DP g,u,p,d,s,t,dmy1;
1.1 noro 771: P dmy;
772: NODE b,l;
773: DP *hps;
774: MP m,mr;
775: int i,n;
776: int *wb;
777: int full;
778: char *fprefix;
779: int sugar,psugar;
780:
781: asir_assert(ARG0(arg),O_LIST,"dp_nf_demand");
782: asir_assert(ARG1(arg),O_DP,"dp_nf_demand");
783: asir_assert(ARG2(arg),O_N,"dp_nf_demand");
784: asir_assert(ARG3(arg),O_VECT,"dp_nf_demand");
785: asir_assert(ARG4(arg),O_STR,"dp_nf_demand");
786: if ( !(g = (DP)ARG1(arg)) ) {
787: *rp = 0; return;
788: }
789: b = BDY((LIST)ARG0(arg)); full = (Q)ARG2(arg) ? 1 : 0;
790: hps = (DP *)BDY((VECT)ARG3(arg)); fprefix = BDY((STRING)ARG4(arg));
791: for ( n = 0, l = b; l; l = NEXT(l), n++ );
792: wb = (int *)ALLOCA(n*sizeof(int));
793: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
794: wb[i] = QTOS((Q)BDY(l));
795: sugar = g->sugar;
796: for ( d = 0; g; ) {
797: for ( u = 0, i = 0; i < n; i++ ) {
798: if ( dp_redble(g,hps[wb[i]]) ) {
799: FILE *fp;
800: char fname[BUFSIZ];
801:
802: sprintf(fname,"%s%d",fprefix,wb[i]);
803: fprintf(stderr,"loading %s\n",fname);
804: fp = fopen(fname,"r"); skipvl(fp);
805: loadobj(fp,(Obj *)&p); fclose(fp);
1.4 noro 806: dp_red(d,g,p,&t,&u,&dmy,&dmy1);
1.1 noro 807: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
808: sugar = MAX(sugar,psugar);
809: if ( !u ) {
810: if ( d )
811: d->sugar = sugar;
812: *rp = d; return;
813: }
814: d = t;
815: break;
816: }
817: }
818: if ( u )
819: g = u;
820: else if ( !full ) {
821: if ( g ) {
822: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
823: }
824: *rp = g; return;
825: } else {
826: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
827: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
828: addd(CO,d,t,&s); d = s;
829: dp_rest(g,&t); g = t;
830:
831: }
832: }
833: if ( d )
834: d->sugar = sugar;
835: *rp = d;
836: }
837:
838: void Pdp_nf_mod(arg,rp)
839: NODE arg;
840: DP *rp;
841: {
842: NODE b;
843: DP g;
844: DP *ps;
845: int mod,full,ac;
846:
847: ac = argc(arg);
848: asir_assert(ARG0(arg),O_LIST,"dp_nf_mod");
849: asir_assert(ARG1(arg),O_DP,"dp_nf_mod");
850: asir_assert(ARG2(arg),O_VECT,"dp_nf_mod");
851: asir_assert(ARG3(arg),O_N,"dp_nf_mod");
852: asir_assert(ARG4(arg),O_N,"dp_nf_mod");
853: if ( !(g = (DP)ARG1(arg)) ) {
854: *rp = 0; return;
855: }
856: b = BDY((LIST)ARG0(arg)); ps = (DP *)BDY((VECT)ARG2(arg));
857: full = QTOS((Q)ARG3(arg)); mod = QTOS((Q)ARG4(arg));
858: dp_nf_mod_qindex(b,g,ps,mod,full,rp);
859: }
860:
861: void Pdp_true_nf_mod(arg,rp)
862: NODE arg;
863: LIST *rp;
864: {
865: NODE b;
866: DP g,nm;
867: P dn;
868: DP *ps;
869: int mod,full;
870: NODE n;
871:
872: asir_assert(ARG0(arg),O_LIST,"dp_nf_mod");
873: asir_assert(ARG1(arg),O_DP,"dp_nf_mod");
874: asir_assert(ARG2(arg),O_VECT,"dp_nf_mod");
875: asir_assert(ARG3(arg),O_N,"dp_nf_mod");
876: asir_assert(ARG4(arg),O_N,"dp_nf_mod");
877: if ( !(g = (DP)ARG1(arg)) ) {
878: nm = 0; dn = (P)ONEM;
879: } else {
880: b = BDY((LIST)ARG0(arg)); ps = (DP *)BDY((VECT)ARG2(arg));
881: full = QTOS((Q)ARG3(arg)); mod = QTOS((Q)ARG4(arg));
882: dp_true_nf_mod(b,g,ps,mod,full,&nm,&dn);
883: }
884: NEWNODE(n); BDY(n) = (pointer)nm;
885: NEWNODE(NEXT(n)); BDY(NEXT(n)) = (pointer)dn;
886: NEXT(NEXT(n)) = 0; MKLIST(*rp,n);
887: }
888:
889: void dp_nf_mod_qindex(b,g,ps,mod,full,rp)
890: NODE b;
891: DP g;
892: DP *ps;
893: int mod,full;
894: DP *rp;
895: {
896: DP u,p,d,s,t;
897: P dmy;
898: NODE l;
899: MP m,mr;
900: int sugar,psugar;
901:
902: if ( !g ) {
903: *rp = 0; return;
904: }
905: sugar = g->sugar;
906: for ( d = 0; g; ) {
907: for ( u = 0, l = b; l; l = NEXT(l) ) {
908: if ( dp_redble(g,p = ps[QTOS((Q)BDY(l))]) ) {
909: dp_red_mod(d,g,p,mod,&t,&u,&dmy);
910: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
911: sugar = MAX(sugar,psugar);
912: if ( !u ) {
913: if ( d )
914: d->sugar = sugar;
915: *rp = d; return;
916: }
917: d = t;
918: break;
919: }
920: }
921: if ( u )
922: g = u;
923: else if ( !full ) {
924: if ( g ) {
925: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
926: }
927: *rp = g; return;
928: } else {
929: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
930: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
931: addmd(CO,mod,d,t,&s); d = s;
932: dp_rest(g,&t); g = t;
933: }
934: }
935: if ( d )
936: d->sugar = sugar;
937: *rp = d;
938: }
939:
940: void dp_nf_mod(b,g,ps,mod,full,rp)
941: NODE b;
942: DP g;
943: DP *ps;
944: int mod,full;
945: DP *rp;
946: {
947: DP u,p,d,s,t;
948: P dmy;
949: NODE l;
950: MP m,mr;
951: int sugar,psugar;
952:
953: if ( !g ) {
954: *rp = 0; return;
955: }
956: sugar = g->sugar;
957: for ( d = 0; g; ) {
958: for ( u = 0, l = b; l; l = NEXT(l) ) {
959: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
960: dp_red_mod(d,g,p,mod,&t,&u,&dmy);
961: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
962: sugar = MAX(sugar,psugar);
963: if ( !u ) {
964: if ( d )
965: d->sugar = sugar;
966: *rp = d; return;
967: }
968: d = t;
969: break;
970: }
971: }
972: if ( u )
973: g = u;
974: else if ( !full ) {
975: if ( g ) {
976: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
977: }
978: *rp = g; return;
979: } else {
980: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
981: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
982: addmd(CO,mod,d,t,&s); d = s;
983: dp_rest(g,&t); g = t;
984: }
985: }
986: if ( d )
987: d->sugar = sugar;
988: *rp = d;
989: }
990:
991: void dp_true_nf_mod(b,g,ps,mod,full,rp,dnp)
992: NODE b;
993: DP g;
994: DP *ps;
995: int mod,full;
996: DP *rp;
997: P *dnp;
998: {
999: DP u,p,d,s,t;
1000: NODE l;
1001: MP m,mr;
1002: int i,n;
1003: int *wb;
1004: int sugar,psugar;
1005: P dn,tdn,tdn1;
1006:
1007: dn = (P)ONEM;
1008: if ( !g ) {
1009: *rp = 0; *dnp = dn; return;
1010: }
1011: for ( n = 0, l = b; l; l = NEXT(l), n++ );
1012: wb = (int *)ALLOCA(n*sizeof(int));
1013: for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
1014: wb[i] = QTOS((Q)BDY(l));
1015: sugar = g->sugar;
1016: for ( d = 0; g; ) {
1017: for ( u = 0, i = 0; i < n; i++ ) {
1018: if ( dp_redble(g,p = ps[wb[i]]) ) {
1019: dp_red_mod(d,g,p,mod,&t,&u,&tdn);
1020: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
1021: sugar = MAX(sugar,psugar);
1022: if ( !u ) {
1023: if ( d )
1024: d->sugar = sugar;
1025: *rp = d; *dnp = dn; return;
1026: } else {
1027: d = t;
1028: mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1;
1029: }
1030: break;
1031: }
1032: }
1033: if ( u )
1034: g = u;
1035: else if ( !full ) {
1036: if ( g ) {
1037: MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
1038: }
1039: *rp = g; *dnp = dn; return;
1040: } else {
1041: m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
1042: NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
1043: addmd(CO,mod,d,t,&s); d = s;
1044: dp_rest(g,&t); g = t;
1045: }
1046: }
1047: if ( d )
1048: d->sugar = sugar;
1049: *rp = d; *dnp = dn;
1050: }
1051:
1052: void Pdp_tdiv(arg,rp)
1053: NODE arg;
1054: DP *rp;
1055: {
1056: MP m,mr,mr0;
1057: DP p;
1058: Q c;
1059: N d,q,r;
1060: int sgn;
1061:
1062: asir_assert(ARG0(arg),O_DP,"dp_tdiv");
1063: asir_assert(ARG1(arg),O_N,"dp_tdiv");
1064: p = (DP)ARG0(arg); d = NM((Q)ARG1(arg)); sgn = SGN((Q)ARG1(arg));
1065: if ( !p )
1066: *rp = 0;
1067: else {
1068: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1069: divn(NM((Q)m->c),d,&q,&r);
1070: if ( r ) {
1071: *rp = 0; return;
1072: } else {
1073: NEXTMP(mr0,mr); NTOQ(q,SGN((Q)m->c)*sgn,c);
1074: mr->c = (P)c; mr->dl = m->dl;
1075: }
1076: }
1077: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1078: }
1079: }
1080:
1081: void Pdp_red_coef(arg,rp)
1082: NODE arg;
1083: DP *rp;
1084: {
1085: MP m,mr,mr0;
1086: P q,r;
1087: DP p;
1088: P mod;
1089:
1090: p = (DP)ARG0(arg); mod = (P)ARG1(arg);
1091: asir_assert(p,O_DP,"dp_red_coef");
1092: asir_assert(mod,O_P,"dp_red_coef");
1093: if ( !p )
1094: *rp = 0;
1095: else {
1096: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1097: divsrp(CO,m->c,mod,&q,&r);
1098: if ( r ) {
1099: NEXTMP(mr0,mr); mr->c = r; mr->dl = m->dl;
1100: }
1101: }
1102: if ( mr0 ) {
1103: NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
1104: } else
1105: *rp = 0;
1106: }
1107: }
1108:
1109: void qltozl(w,n,dvr)
1110: Q *w,*dvr;
1111: int n;
1112: {
1113: N nm,dn;
1114: N g,l1,l2,l3;
1115: Q c,d;
1116: int i;
1117: struct oVECT v;
1118:
1119: for ( i = 0; i < n; i++ )
1120: if ( w[i] && !INT(w[i]) )
1121: break;
1122: if ( i == n ) {
1123: v.id = O_VECT; v.len = n; v.body = (pointer *)w;
1124: igcdv(&v,dvr); return;
1125: }
1126: c = w[0]; nm = NM(c); dn = INT(c) ? ONEN : DN(c);
1127: for ( i = 1; i < n; i++ ) {
1128: c = w[i]; l1 = INT(c) ? ONEN : DN(c);
1129: gcdn(nm,NM(c),&g); nm = g;
1130: gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn);
1131: }
1132: if ( UNIN(dn) )
1133: NTOQ(nm,1,d);
1134: else
1135: NDTOQ(nm,dn,1,d);
1136: *dvr = d;
1137: }
1138:
1139: int comp_nm(a,b)
1140: Q *a,*b;
1141: {
1142: return cmpn((*a)?NM(*a):0,(*b)?NM(*b):0);
1143: }
1144:
1145: void sortbynm(w,n)
1146: Q *w;
1147: int n;
1148: {
1149: qsort(w,n,sizeof(Q),(int (*)(const void *,const void *))comp_nm);
1150: }
1151:
1152: void Pdp_redble(arg,rp)
1153: NODE arg;
1154: Q *rp;
1155: {
1156: asir_assert(ARG0(arg),O_DP,"dp_redble");
1157: asir_assert(ARG1(arg),O_DP,"dp_redble");
1158: if ( dp_redble((DP)ARG0(arg),(DP)ARG1(arg)) )
1159: *rp = ONE;
1160: else
1161: *rp = 0;
1162: }
1163:
1164: void Pdp_red_mod(arg,rp)
1165: NODE arg;
1166: LIST *rp;
1167: {
1168: DP h,r;
1169: P dmy;
1170: NODE n;
1171:
1172: asir_assert(ARG0(arg),O_DP,"dp_red_mod");
1173: asir_assert(ARG1(arg),O_DP,"dp_red_mod");
1174: asir_assert(ARG2(arg),O_DP,"dp_red_mod");
1175: asir_assert(ARG3(arg),O_N,"dp_red_mod");
1176: dp_red_mod((DP)ARG0(arg),(DP)ARG1(arg),(DP)ARG2(arg),QTOS((Q)ARG3(arg)),
1177: &h,&r,&dmy);
1178: NEWNODE(n); BDY(n) = (pointer)h;
1179: NEWNODE(NEXT(n)); BDY(NEXT(n)) = (pointer)r;
1180: NEXT(NEXT(n)) = 0; MKLIST(*rp,n);
1181: }
1182:
1183: int dp_redble(p1,p2)
1184: DP p1,p2;
1185: {
1186: int i,n;
1187: DL d1,d2;
1188:
1189: d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1190: if ( d1->td < d2->td )
1191: return 0;
1192: else {
1193: for ( i = 0, n = p1->nv; i < n; i++ )
1194: if ( d1->d[i] < d2->d[i] )
1195: return 0;
1196: return 1;
1197: }
1198: }
1199:
1200: void dp_red_mod(p0,p1,p2,mod,head,rest,dnp)
1201: DP p0,p1,p2;
1202: int mod;
1203: DP *head,*rest;
1204: P *dnp;
1205: {
1206: int i,n;
1207: DL d1,d2,d;
1208: MP m;
1209: DP t,s,r,h;
1210: P c1,c2,g,u;
1211:
1212: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1213: NEWDL(d,n); d->td = d1->td - d2->td;
1214: for ( i = 0; i < n; i++ )
1215: d->d[i] = d1->d[i]-d2->d[i];
1216: c1 = (P)BDY(p1)->c; c2 = (P)BDY(p2)->c;
1217: gcdprsmp(CO,mod,c1,c2,&g);
1218: divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u;
1219: if ( NUM(c2) ) {
1220: divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM;
1221: }
1222: NEWMP(m); m->dl = d; chsgnmp(mod,(P)c1,&m->c); NEXT(m) = 0;
1223: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p2,s,&t);
1224: if ( NUM(c2) ) {
1225: addmd(CO,mod,p1,t,&r); h = p0;
1226: } else {
1227: mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h);
1228: }
1229: *head = h; *rest = r; *dnp = c2;
1230: }
1231:
1232: void Pdp_subd(arg,rp)
1233: NODE arg;
1234: DP *rp;
1235: {
1236: DP p1,p2;
1237:
1238: p1 = (DP)ARG0(arg); p2 = (DP)ARG1(arg);
1239: asir_assert(p1,O_DP,"dp_subd");
1240: asir_assert(p2,O_DP,"dp_subd");
1241: dp_subd(p1,p2,rp);
1242: }
1243:
1244: void dp_subd(p1,p2,rp)
1245: DP p1,p2;
1246: DP *rp;
1247: {
1248: int i,n;
1249: DL d1,d2,d;
1250: MP m;
1251: DP s;
1252:
1253: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1254: NEWDL(d,n); d->td = d1->td - d2->td;
1255: for ( i = 0; i < n; i++ )
1256: d->d[i] = d1->d[i]-d2->d[i];
1.2 noro 1257: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
1258: MKDP(n,m,s); s->sugar = d->td;
1259: *rp = s;
1260: }
1261:
1262: void dltod(d,n,rp)
1263: DL d;
1264: int n;
1265: DP *rp;
1266: {
1267: MP m;
1268: DP s;
1269:
1270: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
1271: MKDP(n,m,s); s->sugar = d->td;
1.1 noro 1272: *rp = s;
1273: }
1274:
1275: void Pdp_red(arg,rp)
1276: NODE arg;
1277: LIST *rp;
1278: {
1279: NODE n;
1.4 noro 1280: DP head,rest,dmy1;
1.1 noro 1281: P dmy;
1282:
1283: asir_assert(ARG0(arg),O_DP,"dp_red");
1284: asir_assert(ARG1(arg),O_DP,"dp_red");
1285: asir_assert(ARG2(arg),O_DP,"dp_red");
1.4 noro 1286: dp_red((DP)ARG0(arg),(DP)ARG1(arg),(DP)ARG2(arg),&head,&rest,&dmy,&dmy1);
1.1 noro 1287: NEWNODE(n); BDY(n) = (pointer)head;
1288: NEWNODE(NEXT(n)); BDY(NEXT(n)) = (pointer)rest;
1289: NEXT(NEXT(n)) = 0; MKLIST(*rp,n);
1290: }
1291:
1.4 noro 1292: void dp_red(p0,p1,p2,head,rest,dnp,multp)
1.1 noro 1293: DP p0,p1,p2;
1294: DP *head,*rest;
1295: P *dnp;
1.4 noro 1296: DP *multp;
1.1 noro 1297: {
1298: int i,n;
1299: DL d1,d2,d;
1300: MP m;
1301: DP t,s,r,h;
1302: Q c,c1,c2;
1303: N gn,tn;
1304: P g,a;
1305:
1306: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1307: NEWDL(d,n); d->td = d1->td - d2->td;
1308: for ( i = 0; i < n; i++ )
1309: d->d[i] = d1->d[i]-d2->d[i];
1310: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
1311: if ( dp_fcoeffs ) {
1312: /* do nothing */
1313: } else if ( INT(c1) && INT(c2) ) {
1314: gcdn(NM(c1),NM(c2),&gn);
1315: if ( !UNIN(gn) ) {
1316: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
1317: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
1318: }
1319: } else {
1320: ezgcdpz(CO,(P)c1,(P)c2,&g);
1321: divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
1322: }
1323: NEWMP(m); m->dl = d; chsgnp((P)c1,&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
1.4 noro 1324: *multp = s;
1.3 noro 1325: muld(CO,s,p2,&t); muldc(CO,p1,(P)c2,&s); addd(CO,s,t,&r);
1.1 noro 1326: muldc(CO,p0,(P)c2,&h);
1327: *head = h; *rest = r; *dnp = (P)c2;
1328: }
1329:
1330: void Pdp_sp(arg,rp)
1331: NODE arg;
1332: DP *rp;
1333: {
1334: DP p1,p2;
1335:
1336: p1 = (DP)ARG0(arg); p2 = (DP)ARG1(arg);
1337: asir_assert(p1,O_DP,"dp_sp"); asir_assert(p2,O_DP,"dp_sp");
1338: dp_sp(p1,p2,rp);
1339: }
1340:
1.4 noro 1341: extern int GenTrace;
1342: extern NODE TraceList;
1343:
1.1 noro 1344: void dp_sp(p1,p2,rp)
1345: DP p1,p2;
1346: DP *rp;
1347: {
1348: int i,n,td;
1349: int *w;
1350: DL d1,d2,d;
1351: MP m;
1.4 noro 1352: DP t,s1,s2,u;
1.1 noro 1353: Q c,c1,c2;
1354: N gn,tn;
1355:
1356: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1357: w = (int *)ALLOCA(n*sizeof(int));
1358: for ( i = 0, td = 0; i < n; i++ ) {
1359: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
1360: }
1361:
1362: NEWDL(d,n); d->td = td - d1->td;
1363: for ( i = 0; i < n; i++ )
1364: d->d[i] = w[i] - d1->d[i];
1365: c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c;
1366: if ( INT(c1) && INT(c2) ) {
1367: gcdn(NM(c1),NM(c2),&gn);
1368: if ( !UNIN(gn) ) {
1369: divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c;
1370: divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c;
1371: }
1372: }
1373:
1374: NEWMP(m); m->dl = d; m->c = (P)c2; NEXT(m) = 0;
1.4 noro 1375: MKDP(n,m,s1); s1->sugar = d->td; muld(CO,s1,p1,&t);
1.1 noro 1376:
1377: NEWDL(d,n); d->td = td - d2->td;
1378: for ( i = 0; i < n; i++ )
1379: d->d[i] = w[i] - d2->d[i];
1380: NEWMP(m); m->dl = d; m->c = (P)c1; NEXT(m) = 0;
1.4 noro 1381: MKDP(n,m,s2); s2->sugar = d->td; muld(CO,s2,p2,&u);
1.1 noro 1382:
1383: subd(CO,t,u,rp);
1.4 noro 1384: if ( GenTrace ) {
1385: LIST hist;
1386: NODE node;
1387:
1388: node = mknode(4,ONE,0,s1,ONE);
1389: MKLIST(hist,node);
1390: MKNODE(TraceList,hist,0);
1391:
1392: node = mknode(4,ONE,0,0,ONE);
1393: chsgnd(s2,(DP *)&ARG2(node));
1394: MKLIST(hist,node);
1395: MKNODE(node,hist,TraceList); TraceList = node;
1396: }
1.1 noro 1397: }
1398:
1399: void Pdp_sp_mod(arg,rp)
1400: NODE arg;
1401: DP *rp;
1402: {
1403: DP p1,p2;
1404: int mod;
1405:
1406: p1 = (DP)ARG0(arg); p2 = (DP)ARG1(arg);
1407: asir_assert(p1,O_DP,"dp_sp_mod"); asir_assert(p2,O_DP,"dp_sp_mod");
1408: asir_assert(ARG2(arg),O_N,"dp_sp_mod");
1409: mod = QTOS((Q)ARG2(arg));
1410: dp_sp_mod(p1,p2,mod,rp);
1411: }
1412:
1413: void dp_sp_mod(p1,p2,mod,rp)
1414: DP p1,p2;
1415: int mod;
1416: DP *rp;
1417: {
1418: int i,n,td;
1419: int *w;
1420: DL d1,d2,d;
1421: MP m;
1422: DP t,s,u;
1423:
1424: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1425: w = (int *)ALLOCA(n*sizeof(int));
1426: for ( i = 0, td = 0; i < n; i++ ) {
1427: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
1428: }
1429: NEWDL(d,n); d->td = td - d1->td;
1430: for ( i = 0; i < n; i++ )
1431: d->d[i] = w[i] - d1->d[i];
1432: NEWMP(m); m->dl = d; m->c = (P)BDY(p2)->c; NEXT(m) = 0;
1433: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p1,s,&t);
1434: NEWDL(d,n); d->td = td - d2->td;
1435: for ( i = 0; i < n; i++ )
1436: d->d[i] = w[i] - d2->d[i];
1437: NEWMP(m); m->dl = d; m->c = (P)BDY(p1)->c; NEXT(m) = 0;
1438: MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,p2,s,&u);
1439: submd(CO,mod,t,u,rp);
1440: }
1441:
1442: void Pdp_lcm(arg,rp)
1443: NODE arg;
1444: DP *rp;
1445: {
1446: int i,n,td;
1447: DL d1,d2,d;
1448: MP m;
1449: DP p1,p2;
1450:
1451: p1 = (DP)ARG0(arg); p2 = (DP)ARG1(arg);
1452: asir_assert(p1,O_DP,"dp_lcm"); asir_assert(p2,O_DP,"dp_lcm");
1453: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
1454: NEWDL(d,n);
1455: for ( i = 0, td = 0; i < n; i++ ) {
1456: d->d[i] = MAX(d1->d[i],d2->d[i]); td += d->d[i];
1457: }
1458: d->td = td;
1459: NEWMP(m); m->dl = d; m->c = (P)ONE; NEXT(m) = 0;
1460: MKDP(n,m,*rp); (*rp)->sugar = td; /* XXX */
1461: }
1462:
1463: void Pdp_hm(arg,rp)
1464: NODE arg;
1465: DP *rp;
1466: {
1467: DP p;
1468:
1469: p = (DP)ARG0(arg); asir_assert(p,O_DP,"dp_hm");
1470: dp_hm(p,rp);
1471: }
1472:
1473: void dp_hm(p,rp)
1474: DP p;
1475: DP *rp;
1476: {
1477: MP m,mr;
1478:
1479: if ( !p )
1480: *rp = 0;
1481: else {
1482: m = BDY(p);
1483: NEWMP(mr); mr->dl = m->dl; mr->c = m->c; NEXT(mr) = 0;
1484: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
1485: }
1486: }
1487:
1488: void Pdp_ht(arg,rp)
1489: NODE arg;
1490: DP *rp;
1491: {
1492: DP p;
1493: MP m,mr;
1494:
1495: p = (DP)ARG0(arg); asir_assert(p,O_DP,"dp_ht");
1496: if ( !p )
1497: *rp = 0;
1498: else {
1499: m = BDY(p);
1500: NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0;
1501: MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */
1502: }
1503: }
1504:
1505: void Pdp_hc(arg,rp)
1506: NODE arg;
1507: P *rp;
1508: {
1509: asir_assert(ARG0(arg),O_DP,"dp_hc");
1510: if ( !ARG0(arg) )
1511: *rp = 0;
1512: else
1513: *rp = BDY((DP)ARG0(arg))->c;
1514: }
1515:
1516: void Pdp_rest(arg,rp)
1517: NODE arg;
1518: DP *rp;
1519: {
1520: asir_assert(ARG0(arg),O_DP,"dp_rest");
1521: if ( !ARG0(arg) )
1522: *rp = 0;
1523: else
1524: dp_rest((DP)ARG0(arg),rp);
1525: }
1526:
1527: void dp_rest(p,rp)
1528: DP p,*rp;
1529: {
1530: MP m;
1531:
1532: m = BDY(p);
1533: if ( !NEXT(m) )
1534: *rp = 0;
1535: else {
1536: MKDP(p->nv,NEXT(m),*rp);
1537: if ( *rp )
1538: (*rp)->sugar = p->sugar;
1539: }
1540: }
1541:
1542: void Pdp_td(arg,rp)
1543: NODE arg;
1544: Q *rp;
1545: {
1546: DP p;
1547:
1548: p = (DP)ARG0(arg); asir_assert(p,O_DP,"dp_td");
1549: if ( !p )
1550: *rp = 0;
1551: else
1552: STOQ(BDY(p)->dl->td,*rp);
1553: }
1554:
1555: void Pdp_sugar(arg,rp)
1556: NODE arg;
1557: Q *rp;
1558: {
1559: DP p;
1560:
1561: p = (DP)ARG0(arg); asir_assert(p,O_DP,"dp_sugar");
1562: if ( !p )
1563: *rp = 0;
1564: else
1565: STOQ(p->sugar,*rp);
1566: }
1567:
1568: void Pdp_cri1(arg,rp)
1569: NODE arg;
1570: Q *rp;
1571: {
1572: DP p1,p2;
1573: int *d1,*d2;
1574: int i,n;
1575:
1576: p1 = (DP)ARG0(arg); p2 = (DP)ARG1(arg);
1577: asir_assert(p1,O_DP,"dp_cri1"); asir_assert(p2,O_DP,"dp_cri1");
1578: n = p1->nv; d1 = BDY(p1)->dl->d; d2 = BDY(p2)->dl->d;
1579: for ( i = 0; i < n; i++ )
1580: if ( d1[i] > d2[i] )
1581: break;
1582: *rp = i == n ? ONE : 0;
1583: }
1584:
1585: void Pdp_cri2(arg,rp)
1586: NODE arg;
1587: Q *rp;
1588: {
1589: DP p1,p2;
1590: int *d1,*d2;
1591: int i,n;
1592:
1593: p1 = (DP)ARG0(arg); p2 = (DP)ARG1(arg);
1594: asir_assert(p1,O_DP,"dp_cri2"); asir_assert(p2,O_DP,"dp_cri2");
1595: n = p1->nv; d1 = BDY(p1)->dl->d; d2 = BDY(p2)->dl->d;
1596: for ( i = 0; i < n; i++ )
1597: if ( MIN(d1[i],d2[i]) >= 1 )
1598: break;
1599: *rp = i == n ? ONE : 0;
1600: }
1601:
1602: void Pdp_minp(arg,rp)
1603: NODE arg;
1604: LIST *rp;
1605: {
1606: NODE tn,tn1,d,dd,dd0,p,tp;
1607: LIST l,minp;
1608: DP lcm,tlcm;
1609: int s,ts;
1610:
1611: asir_assert(ARG0(arg),O_LIST,"dp_minp");
1612: d = BDY((LIST)ARG0(arg)); minp = (LIST)BDY(d);
1613: p = BDY(minp); p = NEXT(NEXT(p)); lcm = (DP)BDY(p); p = NEXT(p);
1614: if ( !ARG1(arg) ) {
1615: s = QTOS((Q)BDY(p)); p = NEXT(p);
1616: for ( dd0 = 0, d = NEXT(d); d; d = NEXT(d) ) {
1617: tp = BDY((LIST)BDY(d)); tp = NEXT(NEXT(tp));
1618: tlcm = (DP)BDY(tp); tp = NEXT(tp);
1619: ts = QTOS((Q)BDY(tp)); tp = NEXT(tp);
1620: NEXTNODE(dd0,dd);
1621: if ( ts < s ) {
1622: BDY(dd) = (pointer)minp;
1623: minp = (LIST)BDY(d); lcm = tlcm; s = ts;
1624: } else if ( ts == s ) {
1625: if ( compd(CO,lcm,tlcm) > 0 ) {
1626: BDY(dd) = (pointer)minp;
1627: minp = (LIST)BDY(d); lcm = tlcm; s = ts;
1628: } else
1629: BDY(dd) = BDY(d);
1630: } else
1631: BDY(dd) = BDY(d);
1632: }
1633: } else {
1634: for ( dd0 = 0, d = NEXT(d); d; d = NEXT(d) ) {
1635: tp = BDY((LIST)BDY(d)); tp = NEXT(NEXT(tp));
1636: tlcm = (DP)BDY(tp);
1637: NEXTNODE(dd0,dd);
1638: if ( compd(CO,lcm,tlcm) > 0 ) {
1639: BDY(dd) = (pointer)minp; minp = (LIST)BDY(d); lcm = tlcm;
1640: } else
1641: BDY(dd) = BDY(d);
1642: }
1643: }
1644: if ( dd0 )
1645: NEXT(dd) = 0;
1646: MKLIST(l,dd0); MKNODE(tn,l,0); MKNODE(tn1,minp,tn); MKLIST(*rp,tn1);
1647: }
1648:
1649: void Pdp_criB(arg,rp)
1650: NODE arg;
1651: LIST *rp;
1652: {
1653: NODE d,ij,dd,ddd;
1654: int i,j,s,n;
1655: DP *ps;
1656: DL ts,ti,tj,lij,tdl;
1657:
1658: asir_assert(ARG0(arg),O_LIST,"dp_criB"); d = BDY((LIST)ARG0(arg));
1659: asir_assert(ARG1(arg),O_N,"dp_criB"); s = QTOS((Q)ARG1(arg));
1660: asir_assert(ARG2(arg),O_VECT,"dp_criB"); ps = (DP *)BDY((VECT)ARG2(arg));
1661: if ( !d )
1662: *rp = (LIST)ARG0(arg);
1663: else {
1664: ts = BDY(ps[s])->dl;
1665: n = ps[s]->nv;
1666: NEWDL(tdl,n);
1667: for ( dd = 0; d; d = NEXT(d) ) {
1668: ij = BDY((LIST)BDY(d));
1669: i = QTOS((Q)BDY(ij)); ij = NEXT(ij);
1670: j = QTOS((Q)BDY(ij)); ij = NEXT(ij);
1671: lij = BDY((DP)BDY(ij))->dl;
1672: ti = BDY(ps[i])->dl; tj = BDY(ps[j])->dl;
1673: if ( lij->td != lcm_of_DL(n,lij,ts,tdl)->td
1674: || !dl_equal(n,lij,tdl)
1675: || (lij->td == lcm_of_DL(n,ti,ts,tdl)->td
1676: && dl_equal(n,tdl,lij))
1677: || (lij->td == lcm_of_DL(n,tj,ts,tdl)->td
1678: && dl_equal(n,tdl,lij)) ) {
1679: MKNODE(ddd,BDY(d),dd);
1680: dd = ddd;
1681: }
1682: }
1683: MKLIST(*rp,dd);
1684: }
1685: }
1686:
1687: DL lcm_of_DL(nv,dl1,dl2,dl)
1688: int nv;
1689: DL dl1,dl2;
1690: register DL dl;
1691: {
1692: register int n, *d1, *d2, *d, td;
1693:
1694: if ( !dl ) NEWDL(dl,nv);
1695: d = dl->d, d1 = dl1->d, d2 = dl2->d;
1696: for ( td = 0, n = nv; --n >= 0; d1++, d2++, d++ )
1697: td += (*d = *d1 > *d2 ? *d1 : *d2 );
1698: dl->td = td;
1699: return dl;
1700: }
1701:
1702: int dl_equal(nv,dl1,dl2)
1703: int nv;
1704: DL dl1, dl2;
1705: {
1706: register int *d1, *d2, n;
1707:
1708: if ( dl1->td != dl2->td ) return 0;
1709: for ( d1 = dl1->d, d2 = dl2->d, n = nv; --n >= 0; d1++, d2++ )
1710: if ( *d1 != *d2 ) return 0;
1711: return 1;
1712: }
1713:
1714: void Pdp_nelim(arg,rp)
1715: NODE arg;
1716: Q *rp;
1717: {
1718: if ( arg ) {
1719: asir_assert(ARG0(arg),O_N,"dp_nelim");
1720: dp_nelim = QTOS((Q)ARG0(arg));
1721: }
1722: STOQ(dp_nelim,*rp);
1723: }
1724:
1725: void Pdp_mag(arg,rp)
1726: NODE arg;
1727: Q *rp;
1728: {
1729: DP p;
1730: int s;
1731: MP m;
1732:
1733: p = (DP)ARG0(arg);
1734: asir_assert(p,O_DP,"dp_mag");
1735: if ( !p )
1736: *rp = 0;
1737: else {
1738: for ( s = 0, m = BDY(p); m; m = NEXT(m) )
1739: s += p_mag(m->c);
1740: STOQ(s,*rp);
1741: }
1742: }
1743:
1744: extern int kara_mag;
1745:
1746: void Pdp_set_kara(arg,rp)
1747: NODE arg;
1748: Q *rp;
1749: {
1750: if ( arg ) {
1751: asir_assert(ARG0(arg),O_N,"dp_set_kara");
1752: kara_mag = QTOS((Q)ARG0(arg));
1753: }
1754: STOQ(kara_mag,*rp);
1755: }
1756:
1757: void Pdp_homo(arg,rp)
1758: NODE arg;
1759: DP *rp;
1760: {
1761: asir_assert(ARG0(arg),O_DP,"dp_homo");
1762: dp_homo((DP)ARG0(arg),rp);
1763: }
1764:
1765: void dp_homo(p,rp)
1766: DP p;
1767: DP *rp;
1768: {
1769: MP m,mr,mr0;
1770: int i,n,nv,td;
1771: DL dl,dlh;
1772:
1773: if ( !p )
1774: *rp = 0;
1775: else {
1776: n = p->nv; nv = n + 1;
1777: m = BDY(p); td = sugard(m);
1778: for ( mr0 = 0; m; m = NEXT(m) ) {
1779: NEXTMP(mr0,mr); mr->c = m->c;
1780: dl = m->dl;
1781: mr->dl = dlh = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
1782: dlh->td = td;
1783: for ( i = 0; i < n; i++ )
1784: dlh->d[i] = dl->d[i];
1785: dlh->d[n] = td - dl->td;
1786: }
1787: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1788: }
1789: }
1790:
1791: void Pdp_dehomo(arg,rp)
1792: NODE arg;
1793: DP *rp;
1794: {
1795: asir_assert(ARG0(arg),O_DP,"dp_dehomo");
1796: dp_dehomo((DP)ARG0(arg),rp);
1797: }
1798:
1799: void dp_dehomo(p,rp)
1800: DP p;
1801: DP *rp;
1802: {
1803: MP m,mr,mr0;
1804: int i,n,nv;
1805: DL dl,dlh;
1806:
1807: if ( !p )
1808: *rp = 0;
1809: else {
1810: n = p->nv; nv = n - 1;
1811: m = BDY(p);
1812: for ( mr0 = 0; m; m = NEXT(m) ) {
1813: NEXTMP(mr0,mr); mr->c = m->c;
1814: dlh = m->dl;
1815: mr->dl = dl = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int));
1816: dl->td = dlh->td - dlh->d[nv];
1817: for ( i = 0; i < nv; i++ )
1818: dl->d[i] = dlh->d[i];
1819: }
1820: NEXT(mr) = 0; MKDP(nv,mr0,*rp); (*rp)->sugar = p->sugar;
1821: }
1822: }
1823:
1824: int dp_nt(p)
1825: DP p;
1826: {
1827: int i;
1828: MP m;
1829:
1830: if ( !p )
1831: return 0;
1832: else {
1833: for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ );
1834: return i;
1835: }
1836: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>