Annotation of OpenXM_contrib2/asir2000/builtin/gf.c, Revision 1.16
1.3 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.4 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.3 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.16 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/gf.c,v 1.15 2002/03/15 02:52:09 noro Exp $
1.3 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "parse.h"
52:
53: struct resf_dlist {
1.16 ! noro 54: int ib,id;
1.1 noro 55: };
56:
57: int resf_degtest(int,int *,int,struct resf_dlist *);
58: void uhensel(P,NODE,int,int,NODE *);
1.15 noro 59: void uhensel_incremental(P,NODE,int,int,int,NODE *);
1.1 noro 60: void resf_hensel(int,P,int,P *,ML *);
61: void resf_dtest(P,ML,int,int *,int *,DCP *);
62: void resf_dtest_special(P,ML,int,int *,int *,DCP *);
63: void dtest_special(P,ML,P *);
64: void hensel_special(int,int,P,P *,ML *);
65:
66: void nullspace(UM **,UM,int,int,int *);
67: void nullspace_lm(LM **,int,int *);
68: void nullspace_gf2n(GF2N **,int,int *);
69: void nullspace_gfpn(GFPN **,int,int *);
1.5 noro 70: void nullspace_gfs(GFS **,int,int *);
1.12 noro 71: void nullspace_gfsn(GFSN **,int,int *);
1.1 noro 72: void null_to_sol(int **,int *,int,int,UM *);
73:
74: void showgfmat(UM **,int);
75: void pwr_mod(P,P,V,P,int,N,P *);
76: void rem_mod(P,P,V,P,int,P *);
77:
78: void Pnullspace(),Pgcda_mod(),Pftest(),Presfmain(),Ppwr_mod(),Puhensel();
1.15 noro 79: void Puhensel_incremental();
1.6 noro 80: void Psfuhensel();
1.1 noro 81:
82: void Pnullspace_ff();
83:
84: void Psolve_linear_equation_gf2n();
85: void Plinear_form_to_vect(),Pvect_to_linear_form();
86:
87: void solve_linear_equation_gf2n(GF2N **,int,int,int *);
88: void linear_form_to_array(P,VL,int,Num *);
89: void array_to_linear_form(Num *,VL,int,P *);
1.9 noro 90: void sfuhensel(P,NODE,GFS,int,NODE *);
1.1 noro 91:
92: extern int current_ff;
93:
94: struct ftab gf_tab[] = {
1.16 ! noro 95: {"solve_linear_equation_gf2n",Psolve_linear_equation_gf2n,1},
! 96: {"nullspace",Pnullspace,3},
! 97: {"nullspace_ff",Pnullspace_ff,1},
! 98: /* {"gcda_mod",Pgcda_mod,5}, */
! 99: {"ftest",Pftest,2},
! 100: {"resfmain",Presfmain,4},
! 101: {"pwr_mod",Ppwr_mod,6},
! 102: {"uhensel",Puhensel,4},
! 103: {"uhensel_incremental",Puhensel_incremental,5},
! 104: {"sfuhensel",Psfuhensel,4},
! 105: {0,0,0},
1.1 noro 106: };
107:
108: int resf_degtest(k,in,nb,dlist)
109: int k;
110: int *in;
111: int nb;
112: struct resf_dlist *dlist;
113: {
1.16 ! noro 114: int i,d0,d;
! 115: int dl_i;
! 116: struct resf_dlist *t;
! 117:
! 118: if ( k < nb )
! 119: return 0;
! 120: if ( nb == 1 )
! 121: return 1;
! 122: d0 = 0; d = 0; dl_i = 0;
! 123: for ( i = 0; i < k; i++ ) {
! 124: t = &dlist[in[i]];
! 125: if ( t->ib > dl_i + 1 )
! 126: return 0;
! 127: else if ( t->ib == dl_i )
! 128: d += t->id;
! 129: else if ( !d || (dl_i && d0 != d) )
! 130: return 0;
! 131: else {
! 132: d0 = d; dl_i++; d = t->id;
! 133: }
! 134: }
! 135: if ( dl_i != nb - 1 || d0 != d )
! 136: return 0;
! 137: else
! 138: return 1;
1.1 noro 139: }
140:
141: void Puhensel(arg,rp)
142: NODE arg;
143: LIST *rp;
144: {
1.16 ! noro 145: P f;
! 146: NODE mfl,r;
! 147: int mod,bound;
! 148:
! 149: f = (P)ARG0(arg);
! 150: mfl = BDY((LIST)ARG1(arg));
! 151: mod = QTOS((Q)ARG2(arg));
! 152: bound = QTOS((Q)ARG3(arg));
! 153: uhensel(f,mfl,mod,bound,&r);
! 154: MKLIST(*rp,r);
1.1 noro 155: }
156:
1.15 noro 157: void Puhensel_incremental(arg,rp)
158: NODE arg;
159: LIST *rp;
160: {
1.16 ! noro 161: P f;
! 162: NODE mfl,r;
! 163: int mod,bound,start;
! 164:
! 165: f = (P)ARG0(arg);
! 166: mfl = BDY((LIST)ARG1(arg));
! 167: mod = QTOS((Q)ARG2(arg));
! 168: start = QTOS((Q)ARG3(arg));
! 169: bound = QTOS((Q)ARG4(arg));
! 170: uhensel_incremental(f,mfl,mod,start,bound,&r);
! 171: MKLIST(*rp,r);
1.15 noro 172: }
173:
1.1 noro 174: void uhensel(f,mfl,mod,bound,rp)
175: P f;
176: NODE mfl;
177: int mod,bound;
178: NODE *rp;
179: {
1.16 ! noro 180: ML blist,clist,rlist;
! 181: LUM fl;
! 182: int nf,i;
! 183: P s;
! 184: V v;
! 185: NODE t,top;
! 186:
! 187: nf = length(mfl);
! 188: blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
! 189: for ( i = 0, t = mfl; i < nf; i++, t = NEXT(t) ) {
! 190: blist->c[i] = (pointer)UMALLOC(UDEG((P)BDY(t)));
! 191: ptoum(mod,(P)BDY(t),blist->c[i]);
! 192: }
! 193: gcdgen(f,blist,&clist);
! 194: blist->bound = clist->bound = bound;
! 195: W_LUMALLOC((int)UDEG(f),bound,fl);
! 196: ptolum(mod,bound,f,fl);
! 197: henmain(fl,blist,clist,&rlist);
! 198: v = VR(f);
! 199: for ( i = nf-1, top = 0; i >= 0; i-- ) {
! 200: lumtop(v,mod,bound,rlist->c[i],&s);
! 201: MKNODE(t,s,top); top = t;
! 202: }
! 203: *rp = top;
1.15 noro 204: }
205:
206: void uhensel_incremental(f,mfl,mod,start,bound,rp)
207: P f;
208: NODE mfl;
209: int mod,start,bound;
210: NODE *rp;
211: {
1.16 ! noro 212: ML blist,clist,rlist;
! 213: LUM fl;
! 214: LUM *lblist;
! 215: int nf,i,j,k;
! 216: int **p;
! 217: P s;
! 218: V v;
! 219: NODE t,top;
! 220:
! 221: nf = length(mfl);
! 222: blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
! 223: lblist = (LUM *)MALLOC(nf*sizeof(LUM));
! 224: for ( i = 0, t = mfl; i < nf; i++, t = NEXT(t) ) {
! 225: blist->c[i] = (pointer)UMALLOC(UDEG((P)BDY(t)));
! 226: ptoum(mod,(P)BDY(t),blist->c[i]);
! 227: W_LUMALLOC((int)UDEG((P)BDY(t)),bound,lblist[i]);
! 228: ptolum(mod,start,(P)BDY(t),lblist[i]);
! 229: p = lblist[i]->c;
! 230: for ( j = DEG(lblist[i]); j >= 0; j-- )
! 231: for ( k = start; k < bound; k++ )
! 232: p[j][k] = 0;
! 233: }
! 234: gcdgen(f,blist,&clist);
! 235: clist->bound = bound;
! 236: W_LUMALLOC((int)UDEG(f),bound,fl);
! 237: ptolum(mod,bound,f,fl);
! 238: henmain_incremental(fl,lblist,clist,nf,mod,start,bound);
! 239: v = VR(f);
! 240: for ( i = nf-1, top = 0; i >= 0; i-- ) {
! 241: lumtop_unsigned(v,mod,bound,lblist[i],&s);
! 242: MKNODE(t,s,top); top = t;
! 243: }
! 244: *rp = top;
1.6 noro 245: }
246:
247: void Psfuhensel(arg,rp)
248: NODE arg;
249: LIST *rp;
250: {
1.16 ! noro 251: P f;
! 252: int bound;
! 253: NODE r,mfl;
! 254: GFS ev;
! 255:
! 256: f = (P)ARG0(arg);
! 257: mfl = BDY((LIST)ARG1(arg));
! 258: ev = (GFS)ARG2(arg);
! 259: bound = QTOS((Q)ARG3(arg));
! 260: sfuhensel(f,mfl,ev,bound,&r);
! 261: MKLIST(*rp,r);
1.6 noro 262: }
263:
1.9 noro 264: void sfuhensel(f,mfl,ev,bound,rp)
1.6 noro 265: P f;
1.9 noro 266: NODE mfl;
267: GFS ev;
268: int bound;
1.6 noro 269: NODE *rp;
270: {
1.16 ! noro 271: BM fl;
! 272: BM *r;
! 273: VL vl,nvl;
! 274: int i,fn,dx,dy,d;
! 275: NODE t,top;
! 276: UM fm,hm,q;
! 277: UM *gm;
! 278: V x,y;
! 279: P g,s,u;
! 280:
! 281: clctv(CO,f,&vl);
! 282: if ( !vl || !vl->next || vl->next->next )
! 283: error("sfuhensel : f must be a bivariate poly");
! 284:
! 285: for ( i = 0, t = mfl; t; i++, t = NEXT(t) );
! 286: fn = i;
! 287:
! 288: gm = (UM *)MALLOC(fn*sizeof(UM));
! 289:
! 290: /* XXX : more severe check is necessary */
! 291: x = VR((P)BDY(mfl));
! 292: y = vl->v == x ? vl->next->v : vl->v;
! 293:
! 294: for ( i = 0, t = mfl, d = 0; i < fn; i++, t = NEXT(t) ) {
! 295: gm[i] = (pointer)UMALLOC(getdeg(x,(P)BDY(t)));
! 296: ptosfum((P)BDY(t),gm[i]);
! 297: d += DEG(gm[i]);
! 298: }
! 299:
! 300: /* reorder f if necessary */
! 301: if ( vl->v != x ) {
! 302: reordvar(vl,x,&nvl); reorderp(nvl,vl,f,&g);
! 303: vl = nvl; f = g;
! 304: }
! 305: dx = getdeg(x,f);
! 306: if ( dx != d )
! 307: error("sfuhensel : product of factors has incompatible degree");
! 308:
! 309: dy = getdeg(y,f);
! 310: dy = MAX(dy,bound);
! 311: fl = BMALLOC(dx,dy);
! 312: ptosfbm(dy,f,fl);
! 313: if ( ev ) shiftsfbm(fl,FTOIF(CONT(ev)));
! 314:
! 315: /* fm = fl mod y */
! 316: fm = W_UMALLOC(dx);
! 317: cpyum(COEF(fl)[0],fm);
! 318: hm = W_UMALLOC(dx);
! 319:
! 320: q = W_UMALLOC(dx);
! 321: r = (BM *)MLALLOC(fn*sizeof(BM));
! 322: for ( i = 0; i < fn-1; i++ ) {
! 323: /* fl = gm[i]*hm mod y */
! 324: divsfum(fm,gm[i],hm);
! 325: /* fl is replaced by the cofactor of gk mod y^bound */
! 326: /* r[i] = gk */
! 327: sfhenmain2(fl,gm[i],hm,bound,r+i);
! 328: cpyum(hm,fm);
! 329: }
! 330: /* finally, fl must be the lift of gm[fn-1] */
! 331: r[i] = fl;
! 332:
! 333: for ( i = fn-1, top = 0; i >= 0; i-- ) {
! 334: sfbmtop(r[i],x,y,&s);
! 335: reorderp(CO,vl,s,&u);
! 336: MKNODE(t,u,top); top = t;
! 337: }
! 338: *rp = top;
1.1 noro 339: }
340:
341: void Presfmain(arg,rp)
342: NODE arg;
343: LIST *rp;
344: {
1.16 ! noro 345: P f;
! 346: int mod,n,nb,i,j,k;
! 347: int *nf,*md;
! 348: P *mf;
! 349: NODE mfl,mdl,t,s,u;
! 350: ML list;
! 351: DCP dc;
! 352: int sflag;
! 353:
! 354: f = (P)ARG0(arg);
! 355: mod = QTOS((Q)ARG1(arg));
! 356: mfl = BDY((LIST)ARG2(arg));
! 357: mdl = BDY((LIST)ARG3(arg));
! 358: for ( n = nb = 0, t = mfl; t; nb++, t = NEXT(t) )
! 359: for ( s = BDY((LIST)BDY(t)); s; n++, s = NEXT(s) );
! 360: if ( n == nb ) {
! 361: /* f must be irreducible! */
! 362: NEWDC(dc);
! 363: dc->c = f; dc->d = ONE;
! 364: } else {
! 365: nf = W_ALLOC(nb); md = W_ALLOC(nb); mf = (P *)ALLOCA(n*sizeof(P));
! 366: for ( k = i = 0, t = mfl, u = mdl, sflag = 1; t;
! 367: i++, t = NEXT(t), u = NEXT(u) ) {
! 368: if ( sflag && length(BDY((LIST)BDY(t))) != 2 )
! 369: sflag = 0;
! 370: for ( j = 0, s = BDY((LIST)BDY(t)); s; j++, s = NEXT(s) )
! 371: mf[k++] = (P)BDY(s);
! 372: nf[i] = j; md[i] = QTOS((Q)BDY(u));
! 373: }
! 374: resf_hensel(mod,f,n,mf,&list);
! 375: if ( sflag )
! 376: resf_dtest_special(f,list,nb,nf,md,&dc);
! 377: else
! 378: resf_dtest(f,list,nb,nf,md,&dc);
! 379: }
! 380: dcptolist(dc,rp);
1.1 noro 381: }
382:
383: void resf_hensel(mod,f,nf,mfl,listp)
384: int mod;
385: P f;
386: int nf;
387: P *mfl;
388: ML *listp;
389: {
1.16 ! noro 390: register int i,j;
! 391: int q,n,bound,inv,lc;
! 392: int *p;
! 393: int **pp;
! 394: ML blist,clist,bqlist,cqlist,rlist;
! 395: UM *b;
! 396: LUM fl,tl;
! 397: LUM *l;
! 398: UM w;
! 399:
! 400: w = W_UMALLOC(UDEG(f));
! 401: blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
! 402:
! 403: /* c[0] must have lc(f) */
! 404: blist->c[0] = (pointer)UMALLOC(UDEG(mfl[0]));
! 405: ptoum(mod,mfl[0],w);
! 406: inv = invm(w->c[UDEG(mfl[0])],mod);
! 407: lc = rem(NM((Q)LC(f)),mod);
! 408: if ( SGN((Q)LC(f)) < 0 )
! 409: lc = (mod-lc)%mod;
! 410: lc = dmar(inv,lc,0,mod);
! 411: if ( lc == 1 )
! 412: copyum(w,blist->c[0]);
! 413: else
! 414: mulsum(mod,w,lc,blist->c[0]);
! 415:
! 416: /* c[i] (i=1,...,nf-1) must be monic */
! 417: for ( i = 1; i < nf; i++ ) {
! 418: blist->c[i] = (pointer)UMALLOC(UDEG(mfl[i]));
! 419: ptoum(mod,mfl[i],w);
! 420: inv = invm(w->c[UDEG(mfl[i])],mod);
! 421: if ( inv == 1 )
! 422: copyum(w,blist->c[i]);
! 423: else
! 424: mulsum(mod,w,inv,blist->c[i]);
! 425: }
! 426:
! 427: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
! 428: n = bqlist->n; q = bqlist->mod;
! 429: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
! 430: if ( bound == 1 ) {
! 431: *listp = rlist = MLALLOC(n);
! 432: rlist->n = n; rlist->mod = q; rlist->bound = bound;
! 433: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
! 434: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
! 435: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
! 436: pp[j][0] = p[j];
! 437: }
! 438: } else {
! 439: W_LUMALLOC((int)UDEG(f),bound,fl);
! 440: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
! 441: }
1.1 noro 442: }
443:
444: void resf_dtest(f,list,nb,nfl,mdl,dcp)
445: P f;
446: ML list;
447: int nb;
448: int *nfl,*mdl;
449: DCP *dcp;
450: {
1.16 ! noro 451: int n,np,bound,q;
! 452: int i,j,k;
! 453: int *win;
! 454: P g,factor,cofactor;
! 455: Q csum,csumt;
! 456: DCP dcf,dcf0;
! 457: LUM *c;
! 458: ML wlist;
! 459: struct resf_dlist *dlist;
! 460: int z;
! 461:
! 462: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 463: win = W_ALLOC(np+1);
! 464: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 465: wlist = W_MLALLOC(np); wlist->n = list->n;
! 466: wlist->mod = list->mod; wlist->bound = list->bound;
! 467: c = (LUM *)COEF(wlist);
! 468: bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 469: dlist = (struct resf_dlist *)ALLOCA(np*sizeof(struct resf_dlist));
! 470: for ( i = 0, j = 0; j < nb; j++ )
! 471: for ( k = 0; k < nfl[j]; k++, i++ ) {
! 472: dlist[i].ib = j; dlist[i].id = DEG(c[i])/mdl[j];
! 473: }
! 474: k = nb;
! 475: for ( i = 0; i < nb; i++ )
! 476: win[i] = i+1;
! 477: for ( g = f, dcf = dcf0 = 0, --np, z = 0; ; ) {
1.1 noro 478: #if 0
1.16 ! noro 479: if ( !(z++ % 10000) )
! 480: fputc('.',stderr);
1.1 noro 481: #endif
1.16 ! noro 482: if ( resf_degtest(k,win,nb,dlist) &&
! 483: dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
! 484: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
! 485: g = cofactor;
! 486: ucsump(g,&csum); mulq(csum,(Q)COEF(DC(g)),&csumt); csum = csumt;
! 487: for ( i = 0; i < k - 1; i++ )
! 488: for ( j = win[i] + 1; j < win[i + 1]; j++ ) {
! 489: c[j-i-1] = c[j];
! 490: dlist[j-i-1] = dlist[j];
! 491: }
! 492: for ( j = win[k-1] + 1; j <= np; j++ ) {
! 493: c[j-k] = c[j];
! 494: dlist[j-k] = dlist[j];
! 495: }
! 496: if ( ( np -= k ) < k )
! 497: break;
! 498: if ( np - win[0] + 1 < k )
! 499: if ( ++k > np )
! 500: break;
! 501: else
! 502: for ( i = 0; i < k; i++ )
! 503: win[i] = i + 1;
! 504: else
! 505: for ( i = 1; i < k; i++ )
! 506: win[i] = win[0] + i;
! 507: } else if ( !ncombi(1,np,k,win) )
! 508: if ( k == np )
! 509: break;
! 510: else
! 511: for ( i = 0, ++k; i < k; i++ )
! 512: win[i] = i + 1;
! 513: }
! 514: NEXTDC(dcf0,dcf); COEF(dcf) = g;
! 515: DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;
1.1 noro 516: }
517:
518: void resf_dtest_special(f,list,nb,nfl,mdl,dcp)
519: P f;
520: ML list;
521: int nb;
522: int *nfl,*mdl;
523: DCP *dcp;
524: {
1.16 ! noro 525: int n,np,bound,q;
! 526: int i,j;
! 527: int *win;
! 528: P g,factor,cofactor;
! 529: Q csum,csumt;
! 530: DCP dcf,dcf0;
! 531: LUM *c;
! 532: ML wlist;
! 533: int max;
! 534:
! 535: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 536: win = W_ALLOC(np+1);
! 537: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 538: wlist = W_MLALLOC(np); wlist->n = list->n;
! 539: wlist->mod = list->mod; wlist->bound = list->bound;
! 540: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 541: max = 1<<nb;
! 542: for ( g = f, dcf = dcf0 = 0, i = 0; i < max; i++ ) {
1.1 noro 543: #if 0
1.16 ! noro 544: if ( !(i % 1000) )
! 545: fprintf(stderr,"i=%d\n",i);
1.1 noro 546: #endif
1.16 ! noro 547: for ( j = 0; j < nb; j++ )
! 548: win[j] = (i&(1<<j)) ? 2*j+1 : 2*j;
! 549: if ( dtestmain(g,csum,wlist,nb,win,&factor,&cofactor) ) {
1.1 noro 550: #if 0
1.16 ! noro 551: fprintf(stderr,"success : i=%d\n",i);
1.1 noro 552: #endif
1.16 ! noro 553: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
! 554: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = cofactor;
! 555: NEXT(dcf) = 0;*dcp = dcf0;
! 556: return;
! 557: }
! 558: }
! 559: NEXTDC(dcf0,dcf); COEF(dcf) = g;
! 560: DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;
1.1 noro 561: }
562:
563: #if 0
564: void Pftest(arg,rp)
565: NODE arg;
566: P *rp;
567: {
1.16 ! noro 568: ML list;
! 569: DCP dc;
! 570: P p;
! 571: P *mfl;
! 572:
! 573: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
! 574: hensel_special(4,1,p,mfl,&list);
! 575: dtest_special(p,list,rp);
1.1 noro 576: }
577:
578: void dtest_special(f,list,pr)
579: P f;
580: ML list;
581: P *pr;
582: {
1.16 ! noro 583: int n,np,bound,q;
! 584: int i,j,k;
! 585: int *win;
! 586: P g,factor,cofactor;
! 587: Q csum,csumt;
! 588: DCP dc,dcf,dcf0;
! 589: LUM *c;
! 590: ML wlist;
! 591:
! 592: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 593: win = W_ALLOC(np+1);
! 594: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 595: wlist = W_MLALLOC(np); wlist->n = list->n;
! 596: wlist->mod = list->mod; wlist->bound = list->bound;
! 597: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 598: for ( g = f, i = 130000; i < (1<<20); i++ ) {
1.1 noro 599: #if 0
1.16 ! noro 600: if ( !(i % 1000) )
! 601: fprintf(stderr,"i=%d\n",i);
1.1 noro 602: #endif
1.16 ! noro 603: for ( j = 0; j < 20; j++ )
! 604: win[j] = (i&(1<<j)) ? 2*j+1 : 2*j;
! 605: if ( dtestmain(g,csum,wlist,20,win,&factor,&cofactor) ) {
1.1 noro 606: #if 0
1.16 ! noro 607: fprintf(stderr,"success : i=%d\n",i);
1.1 noro 608: #endif
1.16 ! noro 609: *pr = factor; return;
! 610: }
! 611: }
1.1 noro 612: }
613:
614: void hensel_special(index,count,f,mfl,listp)
615: int index,count;
616: P f;
617: P *mfl;
618: ML *listp;
619: {
1.16 ! noro 620: register int i,j;
! 621: int q,n,t,d,r,u,br,tmp,bound;
! 622: int *c,*p,*m,*w;
! 623: int **pp;
! 624: DCP dc;
! 625: ML blist,clist,bqlist,cqlist,rlist;
! 626: UM *b;
! 627: LUM fl,tl;
! 628: LUM *l;
! 629:
! 630: blist = MLALLOC(40); blist->n = 40; blist->mod = 11;
! 631: for ( i = 0; i < 40; i++ ) {
! 632: blist->c[i] = (pointer)UMALLOC(6);
! 633: ptoum(11,mfl[i],blist->c[i]);
! 634: }
! 635: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
! 636: n = bqlist->n; q = bqlist->mod;
! 637: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
! 638: if ( bound == 1 ) {
! 639: *listp = rlist = MLALLOC(n);
! 640: rlist->n = n; rlist->mod = q; rlist->bound = bound;
! 641: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
! 642: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
! 643: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
! 644: pp[j][0] = p[j];
! 645: }
! 646: } else {
! 647: W_LUMALLOC(UDEG(f),bound,fl);
! 648: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
! 649: }
1.1 noro 650: }
651: #endif
652:
653: #if 0
654: void Pftest(arg,rp)
655: NODE arg;
656: P *rp;
657: {
1.16 ! noro 658: ML list;
! 659: DCP dc;
! 660: P p;
! 661: P *mfl;
! 662:
! 663: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
! 664: hensel_special(2,1,p,mfl,&list);
! 665: dtest_special(p,list,rp);
1.1 noro 666: }
667:
668: void dtest_special(f,list,pr)
669: P f;
670: ML list;
671: P *pr;
672: {
1.16 ! noro 673: int n,np,bound,q;
! 674: int i,j,k,t,b0;
! 675: int *win;
! 676: P g,factor,cofactor;
! 677: Q csum,csumt;
! 678: DCP dc,dcf,dcf0;
! 679: LUM *c;
! 680: ML wlist;
! 681: static int nbits[16] = {0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4};
! 682:
! 683: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 684: win = W_ALLOC(np+1);
! 685: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 686: wlist = W_MLALLOC(np); wlist->n = list->n;
! 687: wlist->mod = list->mod; wlist->bound = list->bound;
! 688: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 689: for ( g = f, i = 0; i < (1<<23); i++ ) {
1.1 noro 690: #if 0
1.16 ! noro 691: if ( !(i % 1000) )
! 692: fprintf(stderr,"i=%d\n",i);
1.1 noro 693: #endif
1.16 ! noro 694: t = i>>20; b0 = nbits[t];
! 695: if ( !b0 )
! 696: continue;
! 697: for ( j = 1; j < 6; j++ ) {
! 698: t = (i>>(20-4*j))&0xf;
! 699: if ( nbits[t] != b0 )
! 700: break;
! 701: }
! 702: if ( j != 6 )
! 703: continue;
! 704: for ( j = k = 0; j < 24; j++ )
! 705: if ( i & (1<<(23-j)) )
! 706: win[k++] = j;
! 707: if ( dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
1.1 noro 708: #if 0
1.16 ! noro 709: fprintf(stderr,"success : i=%d\n",i);
1.1 noro 710: #endif
1.16 ! noro 711: *pr = factor; return;
! 712: }
! 713: }
! 714: *pr = f;
1.1 noro 715: }
716:
717: void hensel_special(index,count,f,mfl,listp)
718: int index,count;
719: P f;
720: P *mfl;
721: ML *listp;
722: {
1.16 ! noro 723: register int i,j;
! 724: int q,n,t,d,r,u,br,tmp,bound;
! 725: int *c,*p,*m,*w;
! 726: int **pp;
! 727: DCP dc;
! 728: ML blist,clist,bqlist,cqlist,rlist;
! 729: UM *b;
! 730: LUM fl,tl;
! 731: LUM *l;
! 732:
! 733: blist = MLALLOC(24); blist->n = 24; blist->mod = 5;
! 734: for ( i = 0; i < 24; i++ ) {
! 735: blist->c[i] = (pointer)UMALLOC(7);
! 736: ptoum(5,mfl[i],blist->c[i]);
! 737: }
! 738: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
! 739: n = bqlist->n; q = bqlist->mod;
! 740: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
! 741: if ( bound == 1 ) {
! 742: *listp = rlist = MLALLOC(n);
! 743: rlist->n = n; rlist->mod = q; rlist->bound = bound;
! 744: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
! 745: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
! 746: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
! 747: pp[j][0] = p[j];
! 748: }
! 749: } else {
! 750: W_LUMALLOC(UDEG(f),bound,fl);
! 751: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
! 752: }
1.1 noro 753: }
754: #endif
755:
756: void Pftest(arg,rp)
757: NODE arg;
758: P *rp;
759: {
1.16 ! noro 760: ML list;
! 761: P p;
! 762: P *mfl;
! 763:
! 764: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
! 765: hensel_special(5,1,p,mfl,&list);
! 766: dtest_special(p,list,rp);
1.1 noro 767: }
768:
769: int nbits(a)
770: int a;
771: {
1.16 ! noro 772: int i,s;
1.1 noro 773:
1.16 ! noro 774: for ( i = 0, s = 0; a && (i < 20); i++, a >>= 1 )
! 775: if ( a & 1 ) s++;
! 776: return s;
1.1 noro 777: }
778:
779: void dtest_special(f,list,pr)
780: P f;
781: ML list;
782: P *pr;
783: {
1.16 ! noro 784: int n,np,bound,q;
! 785: int i,j,k,b0;
! 786: int *win;
! 787: P g,factor,cofactor;
! 788: Q csum,csumt;
! 789: LUM *c;
! 790: ML wlist;
! 791:
! 792: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 793: win = W_ALLOC(np+1);
! 794: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 795: wlist = W_MLALLOC(np); wlist->n = list->n;
! 796: wlist->mod = list->mod; wlist->bound = list->bound;
! 797: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 798: for ( g = f, i = 0; i < (1<<19); i++ ) {
1.1 noro 799: #if 0
1.16 ! noro 800: if ( !(i % 10000) )
! 801: fprintf(stderr,"i=%d\n",i);
1.1 noro 802: #endif
1.16 ! noro 803: b0 = nbits(i>>10);
! 804: if ( !b0 || (nbits(i&0x3ff) != b0) )
! 805: continue;
! 806: for ( j = k = 0; j < 20; j++ )
! 807: if ( i & (1<<(19-j)) )
! 808: win[k++] = j;
! 809: if ( dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
1.1 noro 810: #if 0
1.16 ! noro 811: fprintf(stderr,"success : i=%d\n",i);
1.1 noro 812: #endif
1.16 ! noro 813: *pr = factor; return;
! 814: }
! 815: }
! 816: *pr = f;
1.1 noro 817: }
818:
819: void hensel_special(index,count,f,mfl,listp)
820: int index,count;
821: P f;
822: P *mfl;
823: ML *listp;
824: {
1.16 ! noro 825: register int i,j;
! 826: int q,n,bound;
! 827: int *p;
! 828: int **pp;
! 829: ML blist,clist,bqlist,cqlist,rlist;
! 830: UM *b;
! 831: LUM fl,tl;
! 832: LUM *l;
! 833:
! 834: blist = MLALLOC(20); blist->n = 20; blist->mod = 11;
! 835: for ( i = 0; i < 20; i++ ) {
! 836: blist->c[i] = (pointer)UMALLOC(10);
! 837: ptoum(11,mfl[i],blist->c[i]);
! 838: }
! 839: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
! 840: n = bqlist->n; q = bqlist->mod;
! 841: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
! 842: if ( bound == 1 ) {
! 843: *listp = rlist = MLALLOC(n);
! 844: rlist->n = n; rlist->mod = q; rlist->bound = bound;
! 845: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
! 846: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
! 847: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
! 848: pp[j][0] = p[j];
! 849: }
! 850: } else {
! 851: W_LUMALLOC((int)UDEG(f),bound,fl);
! 852: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
! 853: }
1.1 noro 854: }
855:
856: void Pnullspace(arg,rp)
857: NODE arg;
858: LIST *rp;
859: {
1.16 ! noro 860: int i,j,n,mod;
! 861: MAT mat,r;
! 862: VECT u;
! 863: V v;
! 864: P p,z;
! 865: Q q;
! 866: UM **w;
! 867: UM mp;
! 868: P *t;
! 869: UM *s;
! 870: int *ind;
! 871: NODE n0,n1;
! 872:
! 873: mat = (MAT)ARG0(arg);
! 874: p = (P)ARG1(arg);
! 875: v = VR(p);
! 876: mod = QTOS((Q)ARG2(arg));
! 877: n = mat->row;
! 878: w = (UM **)almat_pointer(n,n);
! 879: for ( i = 0; i < n; i++ )
! 880: for ( j = 0, t = (P *)mat->body[i], s = w[i]; j < n; j++ ) {
! 881: ptomp(mod,t[j],&z);
! 882: s[j] = W_UMALLOC((z&&!NUM(z))?UDEG(z):0);
! 883: mptoum(z,s[j]);
! 884: }
! 885: mp = W_UMALLOC(UDEG(p)); ptoum(mod,p,mp);
! 886: ind = (int *)ALLOCA(n*sizeof(int));
! 887: nullspace(w,mp,mod,n,ind);
! 888: MKMAT(r,n,n);
! 889: for ( i = 0; i < n; i++ )
! 890: for ( j = 0, t = (P *)r->body[i], s = w[i]; j < n; j++ )
! 891: umtop(v,s[j],&t[j]);
! 892: MKVECT(u,n);
! 893: for ( i = 0; i < n; i++ ) {
! 894: STOQ(ind[i],q); u->body[i] = (pointer)q;
! 895: }
! 896: MKNODE(n1,u,0); MKNODE(n0,r,n1); MKLIST(*rp,n0);
1.1 noro 897: }
898:
899: void nullspace(mat,p,mod,n,ind)
900: UM **mat;
901: UM p;
902: int mod,n;
903: int *ind;
904: {
1.16 ! noro 905: int i,j,l,s,d;
! 906: UM q,w,w1,h,inv;
! 907: UM *t,*u;
! 908:
! 909: d = DEG(p); inv = W_UMALLOC(d); q = W_UMALLOC(2*d);
! 910: w = W_UMALLOC(2*d); w1 = W_UMALLOC(2*d); h = W_UMALLOC(d);
! 911: bzero(ind,n*sizeof(int));
! 912: ind[0] = 0;
! 913: for ( i = j = 0; j < n; i++, j++ ) {
! 914: for ( ; j < n; j++ ) {
! 915: for ( l = i; l < n; l++ )
! 916: if ( DEG(mat[l][j])>=0 )
! 917: break;
! 918: if ( l < n ) {
! 919: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 920: } else
! 921: ind[j] = 1;
! 922: }
! 923: if ( j == n )
! 924: break;
! 925: invum(mod,p,mat[i][j],inv);
! 926: for ( s = j, t = mat[i]; s < n; s++ ) {
! 927: mulum(mod,t[s],inv,w);
! 928: DEG(w) = divum(mod,w,p,q);
! 929: cpyum(w,t[s]);
! 930: }
! 931: for ( l = 0; l < n; l++ ) {
! 932: if ( l == i )
! 933: continue;
! 934: u = mat[l]; DEG(w) = -1; subum(mod,w,u[j],h);
! 935: for ( s = j; s < n; s++ ) {
! 936: mulum(mod,h,t[s],w); addum(mod,w,u[s],w1);
! 937: DEG(w1) = divum(mod,w1,p,q); cpyum(w1,u[s]);
! 938: }
! 939: }
! 940: }
1.1 noro 941: }
942:
943: void Pnullspace_ff(arg,rp)
944: NODE arg;
945: LIST *rp;
946: {
1.16 ! noro 947: int i,j,n;
! 948: MAT mat,r;
! 949: VECT u;
! 950: Q q;
! 951: Obj **w;
! 952: Obj *t;
! 953: Obj *s;
! 954: int *ind;
! 955: NODE n0,n1;
! 956:
! 957: mat = (MAT)ARG0(arg);
! 958: n = mat->row;
! 959: w = (Obj **)almat_pointer(n,n);
! 960: for ( i = 0; i < n; i++ )
! 961: for ( j = 0, t = (Obj *)mat->body[i], s = w[i]; j < n; j++ )
! 962: s[j] = t[j];
! 963: ind = (int *)ALLOCA(n*sizeof(int));
! 964: switch ( current_ff ) {
! 965: case FF_GFP:
! 966: nullspace_lm((LM **)w,n,ind); break;
! 967: case FF_GF2N:
! 968: nullspace_gf2n((GF2N **)w,n,ind); break;
! 969: case FF_GFPN:
! 970: nullspace_gfpn((GFPN **)w,n,ind); break;
! 971: case FF_GFS:
! 972: nullspace_gfs((GFS **)w,n,ind); break;
! 973: case FF_GFSN:
! 974: nullspace_gfsn((GFSN **)w,n,ind); break;
! 975: default:
! 976: error("nullspace_ff : current_ff is not set");
! 977: }
! 978: MKMAT(r,n,n);
! 979: for ( i = 0; i < n; i++ )
! 980: for ( j = 0, t = (Obj *)r->body[i], s = w[i]; j < n; j++ )
! 981: t[j] = s[j];
! 982: MKVECT(u,n);
! 983: for ( i = 0; i < n; i++ ) {
! 984: STOQ(ind[i],q); u->body[i] = (pointer)q;
! 985: }
! 986: MKNODE(n1,u,0); MKNODE(n0,r,n1); MKLIST(*rp,n0);
1.1 noro 987: }
988:
989: void nullspace_lm(mat,n,ind)
990: LM **mat;
991: int n;
992: int *ind;
993: {
1.16 ! noro 994: int i,j,l,s;
! 995: Q q,mod;
! 996: N lm;
! 997: LM w,w1,h,inv;
! 998: LM *t,*u;
! 999:
! 1000: getmod_lm(&lm); NTOQ(lm,1,mod);
! 1001:
! 1002: bzero(ind,n*sizeof(int));
! 1003: ind[0] = 0;
! 1004: for ( i = j = 0; j < n; i++, j++ ) {
! 1005: for ( ; j < n; j++ ) {
! 1006: for ( l = i; l < n; l++ ) {
! 1007: simplm(mat[l][j],&w); mat[l][j] = w;
! 1008: if ( mat[l][j] )
! 1009: break;
! 1010: }
! 1011: if ( l < n ) {
! 1012: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 1013: } else
! 1014: ind[j] = 1;
! 1015: }
! 1016: if ( j == n )
! 1017: break;
! 1018: NTOQ(mat[i][j]->body,1,q); invl(q,mod,(Q *)&inv);
! 1019: for ( s = j, t = mat[i]; s < n; s++ ) {
! 1020: mullm(t[s],inv,&w); t[s] = w;
! 1021: }
! 1022: for ( l = 0; l < n; l++ ) {
! 1023: if ( l == i )
! 1024: continue;
! 1025: u = mat[l]; chsgnlm(u[j],&h);
! 1026: for ( s = j; s < n; s++ ) {
! 1027: mullm(h,t[s],&w); addlm(w,u[s],&w1); u[s] = w1;
! 1028: }
! 1029: }
! 1030: }
1.1 noro 1031: }
1032:
1033: void nullspace_gf2n(mat,n,ind)
1034: GF2N **mat;
1035: int n;
1036: int *ind;
1037: {
1.16 ! noro 1038: int i,j,l,s;
! 1039: GF2N w,w1,h,inv;
! 1040: GF2N *t,*u;
! 1041: extern gf2n_lazy;
! 1042:
! 1043: bzero(ind,n*sizeof(int));
! 1044: ind[0] = 0;
! 1045: for ( i = j = 0; j < n; i++, j++ ) {
! 1046: for ( ; j < n; j++ ) {
! 1047: for ( l = i; l < n; l++ ) {
! 1048: simpgf2n(mat[l][j],&w); mat[l][j] = w;
! 1049: if ( mat[l][j] )
! 1050: break;
! 1051: }
! 1052: if ( l < n ) {
! 1053: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 1054: } else
! 1055: ind[j] = 1;
! 1056: }
! 1057: if ( j == n )
! 1058: break;
! 1059: invgf2n(mat[i][j],&inv);
! 1060: for ( s = j, t = mat[i]; s < n; s++ ) {
! 1061: mulgf2n(t[s],inv,&w); t[s] = w;
! 1062: }
! 1063: for ( l = 0; l < n; l++ ) {
! 1064: if ( l == i )
! 1065: continue;
! 1066: u = mat[l]; h = u[j];
! 1067: for ( s = j; s < n; s++ ) {
! 1068: mulgf2n(h,t[s],&w); addgf2n(w,u[s],&w1); u[s] = w1;
! 1069: }
! 1070: }
! 1071: }
1.1 noro 1072: }
1073:
1074: void nullspace_gfpn(mat,n,ind)
1075: GFPN **mat;
1076: int n;
1077: int *ind;
1078: {
1.16 ! noro 1079: int i,j,l,s;
! 1080: GFPN w,w1,h,inv;
! 1081: GFPN *t,*u;
! 1082:
! 1083: bzero(ind,n*sizeof(int));
! 1084: ind[0] = 0;
! 1085: for ( i = j = 0; j < n; i++, j++ ) {
! 1086: for ( ; j < n; j++ ) {
! 1087: for ( l = i; l < n; l++ ) {
! 1088: simpgfpn(mat[l][j],&w); mat[l][j] = w;
! 1089: if ( mat[l][j] )
! 1090: break;
! 1091: }
! 1092: if ( l < n ) {
! 1093: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 1094: } else
! 1095: ind[j] = 1;
! 1096: }
! 1097: if ( j == n )
! 1098: break;
! 1099: divgfpn((GFPN)ONE,(GFPN)mat[i][j],&inv);
! 1100: for ( s = j, t = mat[i]; s < n; s++ ) {
! 1101: mulgfpn(t[s],inv,&w); t[s] = w;
! 1102: }
! 1103: for ( l = 0; l < n; l++ ) {
! 1104: if ( l == i )
! 1105: continue;
! 1106: u = mat[l]; chsgngfpn(u[j],&h);
! 1107: for ( s = j; s < n; s++ ) {
! 1108: mulgfpn(h,t[s],&w); addgfpn(w,u[s],&w1); u[s] = w1;
! 1109: }
! 1110: }
! 1111: }
1.1 noro 1112: }
1.5 noro 1113:
1114: void nullspace_gfs(mat,n,ind)
1115: GFS **mat;
1116: int n;
1117: int *ind;
1118: {
1.16 ! noro 1119: int i,j,l,s;
! 1120: GFS w,w1,h,inv;
! 1121: GFS *t,*u;
! 1122: GFS one;
! 1123:
! 1124: bzero(ind,n*sizeof(int));
! 1125: ind[0] = 0;
! 1126: mqtogfs(ONEM,&one);
! 1127:
! 1128: for ( i = j = 0; j < n; i++, j++ ) {
! 1129: for ( ; j < n; j++ ) {
! 1130: for ( l = i; l < n; l++ )
! 1131: if ( mat[l][j] )
! 1132: break;
! 1133: if ( l < n ) {
! 1134: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 1135: } else
! 1136: ind[j] = 1;
! 1137: }
! 1138: if ( j == n )
! 1139: break;
! 1140: divgfs(one,mat[i][j],&inv);
! 1141: for ( s = j, t = mat[i]; s < n; s++ ) {
! 1142: mulgfs(t[s],inv,&w); t[s] = w;
! 1143: }
! 1144: for ( l = 0; l < n; l++ ) {
! 1145: if ( l == i )
! 1146: continue;
! 1147: u = mat[l];
! 1148: chsgngfs(u[j],&h);
! 1149: for ( s = j; s < n; s++ ) {
! 1150: mulgfs(h,t[s],&w); addgfs(w,u[s],&w1); u[s] = w1;
! 1151: }
! 1152: }
! 1153: }
1.12 noro 1154: }
1155:
1156: void nullspace_gfsn(mat,n,ind)
1157: GFSN **mat;
1158: int n;
1159: int *ind;
1160: {
1.16 ! noro 1161: int i,j,l,s;
! 1162: GFSN w,w1,h,inv;
! 1163: GFSN *t,*u;
! 1164:
! 1165: bzero(ind,n*sizeof(int));
! 1166: ind[0] = 0;
! 1167:
! 1168: for ( i = j = 0; j < n; i++, j++ ) {
! 1169: for ( ; j < n; j++ ) {
! 1170: for ( l = i; l < n; l++ )
! 1171: if ( mat[l][j] )
! 1172: break;
! 1173: if ( l < n ) {
! 1174: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 1175: } else
! 1176: ind[j] = 1;
! 1177: }
! 1178: if ( j == n )
! 1179: break;
! 1180: invgfsn(mat[i][j],&inv);
! 1181: for ( s = j, t = mat[i]; s < n; s++ ) {
! 1182: mulgfsn(t[s],inv,&w); t[s] = w;
! 1183: }
! 1184: for ( l = 0; l < n; l++ ) {
! 1185: if ( l == i )
! 1186: continue;
! 1187: u = mat[l];
! 1188: chsgngfsn(u[j],&h);
! 1189: for ( s = j; s < n; s++ ) {
! 1190: mulgfsn(h,t[s],&w); addgfsn(w,u[s],&w1); u[s] = w1;
! 1191: }
! 1192: }
! 1193: }
1.5 noro 1194: }
1195:
1.1 noro 1196: /* p = a(0)vl[0]+a(1)vl[1]+...+a(m-1)vl[m-1]+a(m) -> array = [a(0) a(1) ... a(m)] */
1197:
1198: void linear_form_to_array(p,vl,m,array)
1199: P p;
1200: VL vl;
1201: int m;
1202: Num *array;
1203: {
1.16 ! noro 1204: int i;
! 1205: DCP dc;
1.1 noro 1206:
1.16 ! noro 1207: bzero((char *)array,(m+1)*sizeof(Num *));
! 1208: for ( i = 0; p && vl; vl = NEXT(vl), i++ ) {
! 1209: if ( ID(p) == O_N )
! 1210: break;
! 1211: else if ( VR(p) == vl->v ) {
! 1212: dc = DC(p);
! 1213: array[i] = (Num)COEF(dc);
! 1214: dc = NEXT(dc);
! 1215: p = dc ? COEF(dc) : 0;
! 1216: }
! 1217: }
! 1218: array[m] = (Num)p;
1.1 noro 1219: }
1220:
1221: void array_to_linear_form(array,vl,m,r)
1222: Num *array;
1223: VL vl;
1224: int m;
1225: P *r;
1226: {
1.16 ! noro 1227: P t;
! 1228: DCP dc0,dc1;
1.1 noro 1229:
1.16 ! noro 1230: if ( !m )
! 1231: *r = (P)array[0];
! 1232: else {
! 1233: array_to_linear_form(array+1,NEXT(vl),m-1,&t);
! 1234: if ( !array[0] )
! 1235: *r = t;
! 1236: else {
! 1237: NEWDC(dc0); DEG(dc0) = ONE; COEF(dc0) = (P)array[0];
! 1238: if ( !t )
! 1239: NEXT(dc0) = 0;
! 1240: else {
! 1241: NEWDC(dc1); DEG(dc1) = 0; COEF(dc1) = t;
! 1242: NEXT(dc1) = 0;
! 1243: NEXT(dc0) = dc1;
! 1244: }
! 1245: MKP(vl->v,dc0,*r);
! 1246: }
! 1247: }
1.1 noro 1248: }
1249:
1250: void Psolve_linear_equation_gf2n(arg,rp)
1251: NODE arg;
1252: LIST *rp;
1253: {
1.16 ! noro 1254: NODE eqs,tn;
! 1255: VL vars,tvl;
! 1256: int i,j,n,m,dim,codim;
! 1257: GF2N **w;
! 1258: int *ind;
! 1259: NODE n0,n1;
! 1260:
! 1261: get_vars(ARG0(arg),&vars);
! 1262: eqs = BDY((LIST)ARG0(arg));
! 1263: for ( n = 0, tn = eqs; tn; tn = NEXT(tn), n++);
! 1264: for ( m = 0, tvl = vars; tvl; tvl = NEXT(tvl), m++);
! 1265: w = (GF2N **)almat_pointer(n,m+1);
! 1266: for ( i = 0, tn = eqs; i < n; i++, tn = NEXT(tn) )
! 1267: linear_form_to_array(BDY(tn),vars,m,(Num *)w[i]);
! 1268: ind = (int *)ALLOCA(m*sizeof(int));
! 1269: solve_linear_equation_gf2n(w,n,m,ind);
! 1270: for ( j = 0, dim = 0; j < m; j++ )
! 1271: if ( ind[j] )
! 1272: dim++;
! 1273: codim = m-dim;
! 1274: for ( i = codim; i < n; i++ )
! 1275: if ( w[i][m] ) {
! 1276: MKLIST(*rp,0); return;
! 1277: }
! 1278: for ( i = 0, n0 = 0; i < codim; i++ ) {
! 1279: NEXTNODE(n0,n1);
! 1280: array_to_linear_form((Num *)w[i],vars,m,(P *)&BDY(n1));
! 1281: }
! 1282: if ( n0 )
! 1283: NEXT(n1) = 0;
! 1284: MKLIST(*rp,n0);
1.1 noro 1285: }
1286:
1287: void solve_linear_equation_gf2n(mat,n,m,ind)
1288: GF2N **mat;
1289: int n;
1290: int *ind;
1291: {
1.16 ! noro 1292: int i,j,l,s;
! 1293: GF2N w,w1,h,inv;
! 1294: GF2N *t,*u;
! 1295: extern gf2n_lazy;
! 1296:
! 1297: bzero(ind,m*sizeof(int));
! 1298: ind[0] = 0;
! 1299: for ( i = j = 0; j < m; i++, j++ ) {
! 1300: for ( ; j < m; j++ ) {
! 1301: for ( l = i; l < n; l++ ) {
! 1302: simpgf2n(mat[l][j],&w); mat[l][j] = w;
! 1303: if ( mat[l][j] )
! 1304: break;
! 1305: }
! 1306: if ( l < n ) {
! 1307: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 1308: } else
! 1309: ind[j] = 1;
! 1310: }
! 1311: if ( j == m )
! 1312: break;
! 1313: invgf2n(mat[i][j],&inv);
! 1314: for ( s = j, t = mat[i]; s <= m; s++ ) {
! 1315: mulgf2n(t[s],inv,&w); t[s] = w;
! 1316: }
! 1317: for ( l = 0; l < n; l++ ) {
! 1318: if ( l == i )
! 1319: continue;
! 1320: u = mat[l]; h = u[j];
! 1321: for ( s = j; s <= m; s++ ) {
! 1322: mulgf2n(h,t[s],&w); addgf2n(w,u[s],&w1); u[s] = w1;
! 1323: }
! 1324: }
! 1325: }
1.1 noro 1326: }
1327:
1328: /*
1329: void null_to_sol(mat,ind,mod,n,r)
1330: int **mat;
1331: int *ind;
1332: int mod,n;
1333: UM *r;
1334: {
1.16 ! noro 1335: int i,j,k,l;
! 1336: int *c;
! 1337: UM w;
! 1338:
! 1339: for ( i = 0, l = 0; i < n; i++ ) {
! 1340: if ( !ind[i] )
! 1341: continue;
! 1342: w = UMALLOC(n);
! 1343: for ( j = k = 0, c = COEF(w); j < n; j++ )
! 1344: if ( ind[j] )
! 1345: c[j] = 0;
! 1346: else
! 1347: c[j] = mat[k++][i];
! 1348: c[i] = mod-1;
! 1349: for ( j = n; j >= 0; j-- )
! 1350: if ( c[j] )
! 1351: break;
! 1352: DEG(w) = j;
! 1353: r[l++] = w;
! 1354: }
1.1 noro 1355: }
1356: */
1357:
1358: void showgfmat(mat,n)
1359: UM **mat;
1360: int n;
1361: {
1.16 ! noro 1362: int i,j,k;
! 1363: int *c;
! 1364: UM p;
! 1365:
! 1366: for ( i = 0; i < n; i++ ) {
! 1367: for ( j = 0; j < n; j++ ) {
! 1368: p = mat[i][j];
! 1369: if ( DEG(p) < 0 )
! 1370: fprintf(asir_out,"0");
! 1371: else
! 1372: for ( p = mat[i][j], k = DEG(p), c = COEF(p); k >= 0; k-- ) {
! 1373: if ( c[k] )
! 1374: fprintf(asir_out,"+%d",c[k]);
! 1375: if ( k > 1 )
! 1376: fprintf(asir_out,"a^%d",k);
! 1377: else if ( k == 1 )
! 1378: fprintf(asir_out,"a",k);
! 1379: }
! 1380: fprintf(asir_out," ");
! 1381: }
! 1382: fprintf(asir_out,"\n");
! 1383: }
1.1 noro 1384: }
1385:
1386: #if 0
1387: void Pgcda_mod(arg,rp)
1388: NODE arg;
1389: P *rp;
1390: {
1.16 ! noro 1391: p1 = (P)ARG0(arg);
! 1392: p2 = (P)ARG1(arg);
! 1393: v = VR((P)ARG2(arg));
! 1394: d = (P)ARG3(arg);
! 1395: m = QTOS((Q)ARG4(arg));
! 1396: reordvar(CO,v,&vl);
! 1397: reorderp(vl,CO,p1,&t); ptomp(m,t,&m1);
! 1398: reorderp(vl,CO,p2,&t); ptomp(m,t,&m2);
! 1399: if ( NUM(m1) || NUM(m2) || VR(m1) != v || VR(m2) != v ) {
! 1400: *rp = ONE; return;
! 1401: }
! 1402: if ( deg(v,m1) >= deg(v,m2) ) {
! 1403: t = m1; m1 = m2; m2 = t;
! 1404: }
! 1405: while ( 1 ) {
! 1406: inva_mod(COEF(DC(m2)),d,m,&inv);
! 1407: NEWDC(dc); NEXT(dc) = 0; MKP(v,dc,h);
! 1408: d0 = deg(v,m1)-deg(v,m2); STOQ(d0,DEG(dc));
! 1409: mulgq(m,d,inv,COEF(DC(m1)),&COEF(dc));
! 1410: mulgp(vl,m,d,m2,h,&t); subgp(vl,m,d,m1,t,&s);
! 1411: }
1.1 noro 1412: }
1413: #endif
1414:
1415: void Ppwr_mod(arg,rp)
1416: NODE arg;
1417: P *rp;
1418: {
1.16 ! noro 1419: P p,a,d,r;
! 1420: int m;
! 1421: Q q;
! 1422: N n;
! 1423:
! 1424: m = QTOS((Q)ARG4(arg)); q = (Q)ARG5(arg); n = q ? NM(q) : 0;
! 1425: ptomp(m,(P)ARG0(arg),&p); ptomp(m,(P)ARG1(arg),&a);
! 1426: ptomp(m,(P)ARG3(arg),&d);
! 1427: pwr_mod(p,a,VR((P)ARG2(arg)),d,m,n,&r);
! 1428: mptop(r,rp);
1.1 noro 1429: }
1430:
1431: void pwr_mod(p,a,v,d,m,n,rp)
1432: P p,a,d;
1433: V v;
1434: int m;
1435: N n;
1436: P *rp;
1437: {
1.16 ! noro 1438: int b;
! 1439: P t,s,r;
! 1440: N n1;
! 1441:
! 1442: if ( !n )
! 1443: *rp = (P)ONEM;
! 1444: else if ( UNIN(n) )
! 1445: *rp = p;
! 1446: else {
! 1447: b = divin(n,2,&n1);
! 1448: pwr_mod(p,a,v,d,m,n1,&t);
! 1449: mulmp(CO,m,t,t,&s); rem_mod(s,a,v,d,m,&r);
! 1450: if ( b ) {
! 1451: mulmp(CO,m,r,p,&t); rem_mod(t,a,v,d,m,rp);
! 1452: } else
! 1453: *rp = r;
! 1454: }
1.1 noro 1455: }
1456:
1457: void rem_mod(p,a,v,d,m,rp)
1458: P p,a,d;
1459: V v;
1460: int m;
1461: P *rp;
1462: {
1.16 ! noro 1463: P tmp,r1,r2;
1.1 noro 1464:
1.16 ! noro 1465: divsrmp(CO,m,p,d,&tmp,&r1);
! 1466: divsrmp(CO,m,r1,a,&tmp,&r2);
! 1467: divsrmp(CO,m,r2,d,&tmp,rp);
1.1 noro 1468: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gf.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / builtin