=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/fctr.c,v retrieving revision 1.1.1.1 retrieving revision 1.24 diff -u -p -r1.1.1.1 -r1.24 --- OpenXM_contrib2/asir2000/builtin/fctr.c 1999/12/03 07:39:07 1.1.1.1 +++ OpenXM_contrib2/asir2000/builtin/fctr.c 2017/02/27 05:14:54 1.24 @@ -1,15 +1,68 @@ -/* $OpenXM: OpenXM/src/asir99/builtin/fctr.c,v 1.1.1.1 1999/11/10 08:12:25 noro Exp $ */ +/* + * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED + * All rights reserved. + * + * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited, + * non-exclusive and royalty-free license to use, copy, modify and + * redistribute, solely for non-commercial and non-profit purposes, the + * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and + * conditions of this Agreement. For the avoidance of doubt, you acquire + * only a limited right to use the SOFTWARE hereunder, and FLL or any + * third party developer retains all rights, including but not limited to + * copyrights, in and to the SOFTWARE. + * + * (1) FLL does not grant you a license in any way for commercial + * purposes. You may use the SOFTWARE only for non-commercial and + * non-profit purposes only, such as academic, research and internal + * business use. + * (2) The SOFTWARE is protected by the Copyright Law of Japan and + * international copyright treaties. If you make copies of the SOFTWARE, + * with or without modification, as permitted hereunder, you shall affix + * to all such copies of the SOFTWARE the above copyright notice. + * (3) An explicit reference to this SOFTWARE and its copyright owner + * shall be made on your publication or presentation in any form of the + * results obtained by use of the SOFTWARE. + * (4) In the event that you modify the SOFTWARE, you shall notify FLL by + * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification + * for such modification or the source code of the modified part of the + * SOFTWARE. + * + * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL + * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND + * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS + * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES' + * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY + * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY. + * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT, + * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY + * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL + * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES + * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES + * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY + * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF + * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART + * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY + * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, + * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. + * + * $OpenXM: OpenXM_contrib2/asir2000/builtin/fctr.c,v 1.23 2004/07/30 07:34:40 noro Exp $ +*/ #include "ca.h" #include "parse.h" void Pfctr(), Pgcd(), Pgcdz(), Plcm(), Psqfr(), Pufctrhint(); -void Pptozp(), Pcont(); +void Pptozp(), Pcont(), Psfcont(); void Pafctr(), Pagcd(); void Pmodsqfr(),Pmodfctr(),Pddd(),Pnewddd(),Pddd_tab(); +void Psfsqfr(),Psffctr(),Psfbfctr(),Psfufctr(),Psfmintdeg(),Psfgcd(); void Pirred_check(), Pnfctr_mod(); +void Pbivariate_hensel_special(); +void sfmintdeg(VL vl,P fx,int dy,int c,P *fr); + struct ftab fctr_tab[] = { - {"fctr",Pfctr,1}, + {"bivariate_hensel_special",Pbivariate_hensel_special,6}, + {"fctr",Pfctr,-2}, {"gcd",Pgcd,-3}, {"gcdz",Pgcdz,2}, {"lcm",Plcm,2}, @@ -17,10 +70,17 @@ struct ftab fctr_tab[] = { {"ufctrhint",Pufctrhint,2}, {"ptozp",Pptozp,1}, {"cont",Pcont,-2}, + {"sfcont",Psfcont,-2}, {"afctr",Pafctr,2}, {"agcd",Pagcd,3}, {"modsqfr",Pmodsqfr,2}, {"modfctr",Pmodfctr,2}, + {"sfsqfr",Psfsqfr,1}, + {"sffctr",Psffctr,1}, + {"sfufctr",Psfufctr,1}, + {"sfbfctr",Psfbfctr,-4}, + {"sfmintdeg",Psfmintdeg,5}, + {"sfgcd",Psfgcd,2}, #if 0 {"ddd",Pddd,2}, {"newddd",Pnewddd,2}, @@ -31,6 +91,40 @@ struct ftab fctr_tab[] = { {0,0,0}, }; +/* bivariate_hensel_special(f(x,y):monic in x,g0(x),h0(y),x,y,d) */ + +void Pbivariate_hensel_special(arg,rp) +NODE arg; +LIST *rp; +{ + DCP dc; + struct oVN vn[2]; + P f,g0,h0,ak,bk,gk,hk; + V vx,vy; + VL nvl; + Q qk,cbd,bb; + int d; + NODE n; + + f = (P)ARG0(arg); + g0 = (P)ARG1(arg); + h0 = (P)ARG2(arg); + vx = VR((P)ARG3(arg)); + vy = VR((P)ARG4(arg)); + d = QTOS((Q)ARG5(arg)); + NEWVL(nvl); nvl->v = vx; + NEWVL(NEXT(nvl)); NEXT(nvl)->v = vy; + NEXT(NEXT(nvl)) = 0; + vn[0].v = vy; vn[0].n = 0; + vn[1].v = 0; vn[1].n = 0; + cbound(nvl,f,&cbd); + addq(cbd,cbd,&bb); + henzq1(g0,h0,bb,&bk,&ak,&qk); + henmv(nvl,vn,f,g0,h0,ak,bk,(P)ONE,(P)ONE,(P)ONE,(P)ONE,qk,d,&gk,&hk); + n = mknode(2,gk,hk); + MKLIST(*rp,n); +} + void Pfctr(arg,rp) NODE arg; LIST *rp; @@ -38,7 +132,12 @@ LIST *rp; DCP dc; asir_assert(ARG0(arg),O_P,"fctr"); - fctrp(CO,(P)ARG0(arg),&dc); + if ( argc(arg) == 1 ) + fctrp(CO,(P)ARG0(arg),&dc); + else { + asir_assert(ARG1(arg),O_P,"fctr"); + fctr_wrt_v_p(CO,(P)ARG0(arg),VR((P)ARG1(arg)),&dc); + } dcptolist(dc,rp); } @@ -58,7 +157,7 @@ P *rp; else if ( !p2 ) *rp = p1; else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) ) - error("gcd : invalid argument"); + gcdprsp(CO,p1,p2,rp); else if ( argc(arg) == 2 ) ezgcdp(CO,p1,p2,rp); else { @@ -198,14 +297,83 @@ P *rp; } } -void Pptozp(arg,rp) +void Psfcont(arg,rp) NODE arg; P *rp; { + DCP dc; + MP mp; + int m; + Obj obj; + P p,p1; + P *l; + V v; + + obj = (Obj)ARG0(arg); + if ( !obj || NUM(obj) ) + *rp = (P)obj; + else if ( OID(obj) == O_P ) { + p = (P)obj; + if ( argc(arg) == 2 ) { + v = VR((P)ARG1(arg)); + change_mvar(CO,p,v,&p1); + if ( VR(p1) != v ) { + *rp = p1; return; + } else + p = p1; + } + for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ ); + l = (P *)ALLOCA(m*sizeof(P)); + for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ ) + l[m] = COEF(dc); + gcdsf(CO,l,m,rp); + } else if ( OID(obj) == O_DP ) { + for ( m = 0, mp = BDY((DP)obj); mp; mp = NEXT(mp), m++ ); + l = (P *)ALLOCA(m*sizeof(P)); + for ( m = 0, mp = BDY((DP)obj); mp; mp = NEXT(mp), m++) + l[m] = mp->c; + gcdsf(CO,l,m,rp); + } +} + +void Pptozp(arg,rp) +NODE arg; +Obj *rp; +{ Q t; + NODE tt,p; + NODE n,n0; + char *key; + P pp; + LIST list; + int get_factor=0; asir_assert(ARG0(arg),O_P,"ptozp"); - ptozp((P)ARG0(arg),1,&t,rp); + + /* analyze the option */ + if ( current_option ) { + for ( tt = current_option; tt; tt = NEXT(tt) ) { + p = BDY((LIST)BDY(tt)); + key = BDY((STRING)BDY(p)); + /* value = (Obj)BDY(NEXT(p)); */ + if ( !strcmp(key,"factor") ) get_factor=1; + else { + error("ptozp: unknown option."); + } + } + } + + ptozp((P)ARG0(arg),1,&t,&pp); + + /* printexpr(NULL,t); */ + /* if the option factor is given, then it returns the answer + in the format [zpoly, num] where num*zpoly is equal to the argument.*/ + if (get_factor) { + n0 = mknode(2,pp,t); + MKLIST(list,n0); + *rp = (Obj)list; + } else + *rp = (Obj)pp; } void Pafctr(arg,rp) @@ -247,33 +415,142 @@ int Divum(); UM *resberle(); +void reduce_sfdc(DCP sfdc, DCP *dc); + void Pmodfctr(arg,rp) NODE arg; LIST *rp; { - DCP dc; - int mod; + DCP dc,dcu; + int mod,i,t; + P p; + Obj u; + VL vl; mod = QTOS((Q)ARG1(arg)); if ( mod < 0 ) error("modfctr : invalid modulus"); - modfctrp(ARG0(arg),mod,NEWDDD,&dc); + p = (P)ARG0(arg); + clctv(CO,p,&vl); + if ( !vl ) { + NEWDC(dc); COEF(dc) = p; DEG(dc) = ONE; NEXT(dc) = 0; + } else if ( !NEXT(vl) ) + modfctrp(ARG0(arg),mod,NEWDDD,&dc); + else { + /* XXX 16384 should be replaced by a macro */ + for ( i = 1, t = mod; t*mod < 16384; t *= mod, i++ ); + current_ff = FF_GFS; + setmod_sf(mod,i); + simp_ff((Obj)p,&u); + mfctrsf(CO,(P)u,&dcu); + reduce_sfdc(dcu,&dc); + } if ( !dc ) { NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0; } dcptolist(dc,rp); } +void Psfgcd(arg,rp) +NODE arg; +LIST *rp; +{ + P ps[2]; + + ps[0] = (P)ARG0(arg); + ps[1] = (P)ARG1(arg); + gcdsf(CO,ps,2,rp); +} + +void Psffctr(arg,rp) +NODE arg; +LIST *rp; +{ + DCP dc; + + mfctrsf(CO,ARG0(arg),&dc); + dcptolist(dc,rp); +} + +void Psfsqfr(arg,rp) +NODE arg; +LIST *rp; +{ + DCP dc; + + sqfrsf(CO,ARG0(arg),&dc); + dcptolist(dc,rp); +} + +void Psfufctr(arg,rp) +NODE arg; +LIST *rp; +{ + DCP dc; + + ufctrsf(ARG0(arg),&dc); + dcptolist(dc,rp); +} + +void Psfbfctr(arg,rp) +NODE arg; +LIST *rp; +{ + V x,y; + DCP dc,dct; + P t; + struct oVL vl1,vl2; + VL vl; + int degbound; + + x = VR((P)ARG1(arg)); + y = VR((P)ARG2(arg)); + vl1.v = x; vl1.next = &vl2; + vl2.v = y; vl2.next = 0; + vl = &vl1; + if ( argc(arg) == 4 ) + degbound = QTOS((Q)ARG3(arg)); + else + degbound = -1; + + sfbfctr((P)ARG0(arg),x,y,degbound,&dc); + for ( dct = dc; dct; dct = NEXT(dct) ) { + reorderp(CO,vl,COEF(dct),&t); COEF(dct) = t; + } + dcptolist(dc,rp); +} + +void Psfmintdeg(arg,rp) +NODE arg; +P *rp; +{ + V x,y; + P r; + struct oVL vl1,vl2; + VL vl; + int dy,c; + + x = VR((P)ARG1(arg)); + y = VR((P)ARG2(arg)); + vl1.v = x; vl1.next = &vl2; + vl2.v = y; vl2.next = 0; + vl = &vl1; + dy = QTOS((Q)ARG3(arg)); + c = QTOS((Q)ARG4(arg)); + sfmintdeg(vl,(P)ARG0(arg),dy,c,&r); + reorderp(CO,vl,r,rp); +} + void Pmodsqfr(arg,rp) NODE arg; LIST *rp; { DCP dc; - if ( !dc ) { + if ( !ARG0(arg) ) { NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0; - } - modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),SQFR,&dc); + } else + modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),SQFR,&dc); dcptolist(dc,rp); } @@ -283,10 +560,10 @@ LIST *rp; { DCP dc; - if ( !dc ) { + if ( !ARG0(arg) ) { NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0; - } - modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),DDD,&dc); + } else + modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),DDD,&dc); dcptolist(dc,rp); } @@ -294,12 +571,12 @@ void Pnewddd(arg,rp) NODE arg; LIST *rp; { - DCP dc; + DCP dc=0; - if ( !dc ) { + if ( !ARG0(arg) ) { NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0; - } - modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),NEWDDD,&dc); + } else + modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),NEWDDD,&dc); dcptolist(dc,rp); } @@ -383,4 +660,61 @@ VECT *rp; for ( i = 0; i < n; i++ ) umtop(v,r[i],(P *)&BDY(result)[i]); *rp = result; +} + +void reduce_sfdc(DCP sfdc,DCP *dcr) +{ + P c,t,s,u,f; + DCP dc0,dc,tdc; + DCP *a; + int i,j,n; + + if ( !current_gfs_ext ) { + /* we simply apply sfptop() */ + for ( dc0 = 0; sfdc; sfdc = NEXT(sfdc) ) { + NEXTDC(dc0,dc); + DEG(dc) = DEG(sfdc); + sfptop(COEF(sfdc),&COEF(dc)); + } + NEXT(dc) = 0; + *dcr = dc0; + return; + } + + if ( NUM(COEF(sfdc)) ) { + sfptop(COEF(sfdc),&c); + sfdc = NEXT(sfdc); + } else + c = (P)ONE; + + for ( n = 0, tdc = sfdc; tdc; tdc = NEXT(tdc), n++ ); + a = (DCP *)ALLOCA(n*sizeof(DCP)); + for ( i = 0, tdc = sfdc; i < n; tdc = NEXT(tdc), i++ ) + a[i] = tdc; + + dc0 = 0; NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = c; + for ( i = 0; i < n; i++ ) { + if ( !a[i] ) + continue; + t = COEF(a[i]); + f = t; + while ( 1 ) { + sf_galois_action(t,ONE,&s); + for ( j = i; j < n; j++ ) + if ( a[j] && !compp(CO,s,COEF(a[j])) ) + break; + if ( j == n ) + error("reduce_sfdc : cannot happen"); + if ( j == i ) { + NEXTDC(dc0,dc); DEG(dc) = DEG(a[i]); + sfptop(f,&COEF(dc)); + break; + } else { + mulp(CO,f,s,&u); f = u; + t = s; + a[j] = 0; + } + } + } + *dcr = dc0; }