version 1.4, 2018/11/12 07:59:33 |
version 1.5, 2019/03/13 08:01:05 |
|
|
* DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, |
* DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, |
* PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. |
* PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. |
* |
* |
* $OpenXM: OpenXM_contrib2/asir2018/builtin/dp.c,v 1.3 2018/11/12 04:25:13 noro Exp $ |
* $OpenXM: OpenXM_contrib2/asir2018/builtin/dp.c,v 1.4 2018/11/12 07:59:33 noro Exp $ |
*/ |
*/ |
#include "ca.h" |
#include "ca.h" |
#include "base.h" |
#include "base.h" |
Line 536 void make_reduced(VECT b,int nv) |
|
Line 536 void make_reduced(VECT b,int nv) |
|
|
|
n = b->len; |
n = b->len; |
p = (DL *)BDY(b); |
p = (DL *)BDY(b); |
if ( p[0]->td == 0 ) { |
|
b->len = 1; |
|
return; |
|
} |
|
for ( i = 0; i < n; i++ ) { |
for ( i = 0; i < n; i++ ) { |
pi = p[i]; |
pi = p[i]; |
if ( !pi ) continue; |
if ( !pi ) continue; |
Line 672 P binpoly(P n,int a,int b) |
|
Line 668 P binpoly(P n,int a,int b) |
|
return r; |
return r; |
} |
} |
|
|
void Pdp_monomial_hilbert_poincare(NODE arg,LIST *rp) |
void mhp_to_hf(VL vl,P hp,int n,P *plist,VECT *head,P *hf,Z *den) |
{ |
{ |
LIST g,v; |
P tv,gcd,q,h,hphead,tt,ai,hpoly,nv,bp,w; |
VL vl; |
Z d; |
int m,n,i,hf; |
DCP dc,topdc; |
VECT b,x; |
VECT hfhead; |
NODE t,nd; |
int i,s,qd; |
Z z; |
|
P hp,tv,mt,t1,u,w; |
|
DP a; |
|
DL *p; |
|
P *plist,*r; |
|
Obj val; |
|
|
|
// init_eg(&eg_comp); |
if ( !hp ) { |
if ( get_opt("hf",&val) && val ) hf = 1; |
MKVECT(hfhead,0); *head = hfhead; |
else hf = 0; |
*hf = 0; *den = ONE; |
i_simple = i_all = 0; |
|
g = (LIST)ARG0(arg); v = (LIST)ARG1(arg); |
|
pltovl(v,&vl); |
|
m = length(BDY(g)); MKVECT(b,m); p = (DL *)BDY(b); |
|
for ( t = BDY(g), i = 0; t; t = NEXT(t), i++ ) { |
|
ptod(CO,vl,(P)BDY(t),&a); p[i] = BDY(a)->dl; |
|
} |
|
n = length(BDY(v)); MKVECT(x,n); p = (DL *)BDY(x); |
|
for ( t = BDY(v), i = 0; t; t = NEXT(t), i++ ) { |
|
ptod(CO,vl,(P)BDY(t),&a); p[i] = BDY(a)->dl; |
|
} |
|
/* create (1,1-t,...,(1-t)^n) */ |
|
plist = (P *)MALLOC((n+1)*sizeof(P)); |
|
/* t1 = 1-t */ |
|
makevar("t",&tv); chsgnp(tv,&mt); addp(CO,mt,(P)ONE,&t1); |
|
for ( plist[0] = (P)ONE, i = 1; i <= n; i++ ) |
|
mulp(CO,plist[i-1],t1,&plist[i]); |
|
r = (P *)CALLOC(n+1,sizeof(P)); |
|
make_reduced(b,n); |
|
mhp_rec(b,x,tv,r); |
|
for ( hp = 0, i = 0; i <= n; i++ ) { |
|
mulp(CO,plist[i],r[i],&u); addp(CO,u,hp,&w); hp = w; |
|
} |
|
// printf("all=%d,simple=%d,comp=%.3fsec\n",i_all,i_simple,eg_comp.exectime); |
|
UTOZ(n,z); |
|
if ( !hf ) { |
|
nd = mknode(2,hp,z); |
|
MKLIST(*rp,nd); |
|
} else { |
} else { |
P gcd,q,head,hphead,t,w,ai,hpoly,nv,bp; |
makevar("t",&tv); |
Z den,d; |
|
DCP dc,topdc; |
|
VECT hfhead; |
|
int s,qd,i; |
|
|
|
ezgcdp(CO,hp,plist[n],&gcd); |
ezgcdp(CO,hp,plist[n],&gcd); |
if ( NUM(gcd) ) { |
if ( NUM(gcd) ) { |
s = n; |
s = n; |
Line 734 void Pdp_monomial_hilbert_poincare(NODE arg,LIST *rp) |
|
Line 691 void Pdp_monomial_hilbert_poincare(NODE arg,LIST *rp) |
|
} |
} |
if ( NUM(q) ) qd = 0; |
if ( NUM(q) ) qd = 0; |
else qd = ZTOS(DEG(DC(q))); |
else qd = ZTOS(DEG(DC(q))); |
topdc = 0; |
if ( qd ) { |
for ( i = 0; i < qd; i++ ) { |
topdc = 0; |
NEWDC(dc); NEXT(dc) = topdc; |
for ( i = 0; i < qd; i++ ) { |
ibin(i+s-1,s-1,&COEF(dc)); |
NEWDC(dc); NEXT(dc) = topdc; |
STOZ(i,d); DEG(dc) = d; |
ibin(i+s-1,s-1,&COEF(dc)); |
topdc = dc; |
STOZ(i,d); DEG(dc) = d; |
|
topdc = dc; |
|
} |
|
MKP(VR(tv),topdc,h); |
|
mulp(CO,h,q,&hphead); |
} |
} |
MKP(VR(tv),topdc,head); |
MKVECT(hfhead,qd); |
mulp(CO,head,q,&hphead); |
|
MKVECT(hfhead,qd); |
|
for ( i = 0; i < qd; i++ ) |
for ( i = 0; i < qd; i++ ) |
coefp(hphead,i,(P *)&BDY(hfhead)[i]); |
coefp(hphead,i,(P *)&BDY(hfhead)[i]); |
|
*head = hfhead; |
hpoly = 0; |
hpoly = 0; |
makevar("n",&nv); |
makevar("n",&nv); |
for ( i = 0; i <= qd; i++ ) { |
for ( i = 0; i <= qd; i++ ) { |
coefp(q,i,&ai); |
coefp(q,i,&ai); |
bp = binpoly(nv,s-i-1,s-1); |
bp = binpoly(nv,s-i-1,s-1); |
mulp(CO,ai,bp,&t); |
mulp(CO,ai,bp,&tt); |
addp(CO,hpoly,t,&w); |
addp(CO,hpoly,tt,&w); |
hpoly = w; |
hpoly = w; |
} |
} |
factorialz(s-1,&den); |
*hf = hpoly; |
nd = mknode(5,hp,z,hfhead,hpoly,den); |
factorialz(s-1,den); |
MKLIST(*rp,nd); |
|
} |
} |
} |
} |
|
|
|
/* create (1,1-t,...,(1-t)^n) */ |
|
|
|
P *mhp_prep(int n,P *tv) { |
|
P *plist; |
|
P mt,t1; |
|
int i; |
|
|
|
plist = (P *)MALLOC((n+1)*sizeof(P)); |
|
/* t1 = 1-t */ |
|
makevar("t",tv); chsgnp(*tv,&mt); addp(CO,mt,(P)ONE,&t1); |
|
for ( plist[0] = (P)ONE, i = 1; i <= n; i++ ) |
|
mulp(CO,plist[i-1],t1,&plist[i]); |
|
return plist; |
|
} |
|
|
|
P mhp_ctop(P *r,P *plist,int n) |
|
{ |
|
int i; |
|
P hp,u,w; |
|
|
|
for ( hp = 0, i = 0; i <= n; i++ ) { |
|
mulp(CO,plist[i],r[i],&u); addp(CO,u,hp,&w); hp = w; |
|
} |
|
return hp; |
|
} |
|
|
|
void Pdp_monomial_hilbert_poincare(NODE arg,LIST *rp) |
|
{ |
|
LIST g,v; |
|
VL vl; |
|
int m,n,i; |
|
VECT b,x,hfhead; |
|
NODE t,nd; |
|
Z z,den; |
|
P hp,tv,mt,t1,u,w,hpoly; |
|
DP a; |
|
DL *p; |
|
P *plist,*r; |
|
Obj val; |
|
|
|
i_simple = i_all = 0; |
|
g = (LIST)ARG0(arg); v = (LIST)ARG1(arg); |
|
pltovl(v,&vl); |
|
m = length(BDY(g)); MKVECT(b,m); p = (DL *)BDY(b); |
|
for ( t = BDY(g), i = 0; t; t = NEXT(t), i++ ) { |
|
if ( !BDY(t) ) |
|
p[i] = 0; |
|
else { |
|
ptod(CO,vl,(P)BDY(t),&a); p[i] = BDY(a)->dl; |
|
} |
|
} |
|
n = length(BDY(v)); MKVECT(x,n); p = (DL *)BDY(x); |
|
for ( t = BDY(v), i = 0; t; t = NEXT(t), i++ ) { |
|
ptod(CO,vl,(P)BDY(t),&a); p[i] = BDY(a)->dl; |
|
} |
|
|
|
r = (P *)CALLOC(n+1,sizeof(P)); |
|
plist = mhp_prep(n,&tv); |
|
make_reduced(b,n); |
|
mhp_rec(b,x,tv,r); |
|
hp = mhp_ctop(r,plist,n); |
|
mhp_to_hf(CO,hp,n,plist,&hfhead,&hpoly,&den); |
|
UTOZ(n,z); |
|
nd = mknode(5,hp,z,hfhead,hpoly,den); |
|
MKLIST(*rp,nd); |
|
} |
|
|
#endif |
#endif |
|
|
void Pdp_compute_last_t(NODE arg,LIST *rp) |
void Pdp_compute_last_t(NODE arg,LIST *rp) |