OpenXM_contrib2/asir2018/builtin/dp.c - diff

Return to dp.c CVS log

Up to [local] / OpenXM_contrib2 / asir2018 / builtin

Diff for /OpenXM_contrib2/asir2018/builtin/dp.c between version 1.30 and 1.31

version 1.30, 2021/03/10 06:36:20

version 1.31, 2021/12/05 22:41:03

Line 45

* DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,

* PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.

* $OpenXM: OpenXM_contrib2/asir2018/builtin/dp.c,v 1.29 2021/02/28 02:33:15 noro Exp $

* $OpenXM: OpenXM_contrib2/asir2018/builtin/dp.c,v 1.30 2021/03/10 06:36:20 noro Exp $

#include "ca.h"

#include "base.h"

Line 521 void mhp_rec(VECT b,VECT x,P t,P *r)

// get_eg(&eg1); add_eg(&eg_comp,&eg0,&eg1);

}

P mhp_rec_weight(VECT b,VECT x,P t,int *w)

{

int n,i,j,k,l,i2,nv,len,td;

int *d;

Z wj;

P twj,tmp,tmp2,ret,qadd,qcolon;

DL *p,*q;

DL pi,xj,d1;

VECT c;

i_all++;

n = b->len; nv = x->len; p = (DL *)BDY(b);

if ( !n ) {

// I=<0> => HP(t)=1/(1-t^w1)...(1-t^wn) => Q(t)=1

return (P)ONE;

}

if ( n == 1 && p[0]->td == 0 ) {

// I=<1> => HP(t)=0 => Q(t)=0

return 0;

}

for ( i = 0; i < n; i++ ) {

d = p[i]->d;

for ( td = 0, j = 0; j < nv; j++ ) td += d[j];

if (td > 1 ) break;

}

if ( i == n ) {

// I=<xi1,...,xin> => Q(t)=(1-t^wi1)...(1-t^win)

for ( ret = (P)ONE, i = 0; i < n; i++ ) {

d = p[i]->d;

for ( j = 0; j < nv; j++ ) if ( d[j] ) break;

STOZ(w[j],wj); pwrp(CO,t,wj,&tmp);

subp(CO,(P)ONE,tmp,&tmp2); mulp(CO,ret,tmp2,&tmp); ret = tmp;

}

return ret;

}

for ( j = 0, d = p[i]->d; j < nv; j++ )

if ( d[j] ) break;

xj = BDY(x)[j];

MKVECT(c,n); q = (DL *)BDY(c);

for ( i = k = l = 0; i < n; i++ )

if ( p[i]->d[j] ) {

pi = p[i];

NEWDL(d1,nv); d1->td =pi->td - 1;

memcpy(d1->d,pi->d,nv*sizeof(int));

d1->d[j]--;

p[k++] = d1;

} else

q[l++] = p[i];

for ( i = k, i2 = 0; i2 < l; i++, i2++ )

p[i] = q[i2];

/* b=(b[0]/xj,...,b[k-1]/xj,b[k],...b[n-1]) where

b[0],...,b[k-1] are divisible by k */

make_reduced2(b,k,nv);

qcolon = mhp_rec_weight(b,x,t,w);

/* c = (b[0],...,b[l-1],xj) */

q[l] = xj; c->len = l+1;

qadd = mhp_rec_weight(c,x,t,w);

// Q(t)=Qadd+t^wj*Qcolon

STOZ(w[j],wj); pwrp(CO,t,wj,&twj);

mulp(CO,twj,qcolon,&tmp); addp(CO,qadd,tmp,&ret);

return ret;

}

/* (n+a)Cb as a polynomial of n; return (n+a)*...*(n+a-b+1) */

P binpoly(P n,int a,int b)

Line 636 P mhp_ctop(P *r,P *plist,int n)

Line 699 P mhp_ctop(P *r,P *plist,int n)

return hp;

}

LIST dp_monomial_hilbert_poincare(VECT b,VECT x,P *plist)

LIST dp_monomial_hilbert_poincare(VECT b,VECT x)

{

int n;

P *r;

P *r,*plist;

P tv;

P hp,hpoly;

VECT hfhead;

Line 649 LIST dp_monomial_hilbert_poincare(VECT b,VECT x,P *pli

Line 712 LIST dp_monomial_hilbert_poincare(VECT b,VECT x,P *pli

LIST list;

n = x->len;

plist = mhp_prep(n,&tv);

r = (P *)CALLOC(n+1,sizeof(P));

make_reduced(b,n);

makevar("t",&tv);

mhp_rec(b,x,tv,r);

hp = mhp_ctop(r,plist,n);

mhp_to_hf(CO,hp,n,plist,&hfhead,&hpoly);

Line 662 LIST dp_monomial_hilbert_poincare(VECT b,VECT x,P *pli

Line 725 LIST dp_monomial_hilbert_poincare(VECT b,VECT x,P *pli

return list;

}

LIST dp_monomial_hilbert_poincare_weight(VECT b,VECT x,int *w)

{

int n,i;

NODE nd;

LIST list;

P tv,ret;

n = x->len;

make_reduced(b,n);

makevar("t",&tv);

ret = mhp_rec_weight(b,x,tv,w);

nd = mknode(1,ret);

MKLIST(list,nd);

return list;

}

void Pdp_monomial_hilbert_poincare(NODE arg,LIST *rp)

{

LIST g,v;

VL vl;

int m,n,i;

int m,n,i,wlen;

VECT b,x,hfhead,prep;

NODE t,nd;

Z z,den;

P hp,tv,mt,t1,u,w,hpoly;

P hp,tv,mt,t1,u,hpoly;

DP a;

DL *p;

P *plist;

Obj val,ord,weight;

Obj val,ord;

int *w;

struct order_spec *current_spec=0,*spec;

weight = 0;

if ( current_option ) {

if ( peek_option(current_option,"ord",&ord) ) {

current_spec = dp_current_spec;

create_order_spec(0,ord,&spec);

initd(spec);

}

peek_option(current_option,"weight",&weight);

}

i_simple = i_all = 0;

g = (LIST)ARG0(arg); v = (LIST)ARG1(arg);

Line 699 void Pdp_monomial_hilbert_poincare(NODE arg,LIST *rp)

Line 780 void Pdp_monomial_hilbert_poincare(NODE arg,LIST *rp)

for ( t = BDY(v), i = 0; t; t = NEXT(t), i++ ) {

ptod(CO,vl,(P)BDY(t),&a); p[i] = BDY(a)->dl;

}

plist = mhp_prep(n,&tv);

if ( weight ) {

*rp = dp_monomial_hilbert_poincare(b,x,plist);

wlen = length(BDY((LIST)weight));

if ( n != wlen )

error("dp_monomial_hilbert_poincare: inconsistent weight length");

w = (int *)MALLOC(n*sizeof(int));

for ( i = 0, nd = BDY((LIST)weight); i < n; i++, nd = NEXT(nd) )

w[i] = ZTOS((Z)BDY(nd));

} else if ( current_dl_weight_vector )

w = current_dl_weight_vector;

else

w = 0;

if ( w ) {

*rp = dp_monomial_hilbert_poincare_weight(b,x,w);

} else {

*rp = dp_monomial_hilbert_poincare(b,x);

}

if ( current_spec )

initd(current_spec);

}

Line 746 void Pdp_monomial_hilbert_poincare_incremental(NODE ar

Line 841 void Pdp_monomial_hilbert_poincare_incremental(NODE ar

g = BDY((LIST)ARG0(arg)); new = BDY((DP)ARG1(arg))->dl;

data = BDY((LIST)ARG2(arg));

hn = (P)ARG0(data); n = ZTOS((Z)ARG1(data));

plist = (P *)BDY((VECT)ARG4(data));

len = length(g); MKVECT(b,len); p = (DL *)BDY(b);

for ( t = g, i = 0; t; t = NEXT(t), i++ )

p[i] = monomial_colon(BDY((DP)BDY(t))->dl,new,n);

Line 755 void Pdp_monomial_hilbert_poincare_incremental(NODE ar

Line 849 void Pdp_monomial_hilbert_poincare_incremental(NODE ar

NEWDL(dl,n); dl->d[i] = 1; dl->td = 1; BDY(x)[i] = dl;

}

// compute HP(I:new)

list1 = dp_monomial_hilbert_poincare(b,x,plist);

list1 = dp_monomial_hilbert_poincare(b,x);

data1 = BDY((LIST)list1);

hn1 = (P)ARG0(data1);

// HP(I+<new>) = H(I)-t^d*H(I:new), d=tdeg(new)

makevar("t",&tv); UTOZ(new->td,dz);

plist = mhp_prep(n,&tv);

UTOZ(new->td,dz);

pwrp(CO,tv,dz,&td);

mulp(CO,hn1,td,&s);

subp(CO,hn,s,&newhn);

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>