OpenXM_contrib2/asir2000/engine/nd.c - diff

Return to nd.c CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / engine

Diff for /OpenXM_contrib2/asir2000/engine/nd.c between version 1.145 and 1.154

version 1.145, 2006/08/26 05:38:06

version 1.154, 2008/02/28 01:25:31

Line 1

/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.144 2006/08/26 03:09:55 noro Exp $ */

/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.153 2007/04/30 02:11:12 noro Exp $ */

#include "nd.h"

Line 8 int nd_dcomp;

NM _nm_free_list;

ND _nd_free_list;

ND_pairs _ndp_free_list;

NODE nd_hcf;

static NODE nd_subst;

static VL nd_vc;

static int nd_ntrans;

static int nd_nalg;

#if 0

Line 47 static int nd_demand;

Line 50 static int nd_demand;

NumberField get_numberfield();

UINT *nd_det_compute_bound(NDV **dm,int n,int j);

void nd_det_reconstruct(NDV **dm,int n,int j,NDV d);

void nd_heu_nezgcdnpz(VL vl,P *pl,int m,int full,P *pr);

int nd_monic(int m,ND *p);

NDV plain_vect_to_ndv_q(Q *mat,int col,UINT *s0vect);

extern int Denominator,DP_Multiple;

void nd_free_private_storage()

{

_nm_free_list = 0;

Line 207 void ndl_homogenize(UINT *d,UINT *r,int obpe,EPOS oepo

Line 213 void ndl_homogenize(UINT *d,UINT *r,int obpe,EPOS oepo

{

int w,i,e,n,omask0;

omask0 = (1<<obpe)-1;

omask0 = obpe==32?0xffffffff:((1<<obpe)-1);

n = nd_nvar-1;

ndl_zero(r);

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

Line 976 ND nd_add(int mod,ND p1,ND p2)

Line 982 ND nd_add(int mod,ND p1,ND p2)

int nm_destructive_add_q(NM *m1,NM *m2,NM *mr0,NM *mr)

{

NM s;

Q t;

P t;

int can;

addq(CQ(*m1),CQ(*m2),&t);

addp(nd_vc,CP(*m1),CP(*m2),&t);

s = *m1; *m1 = NEXT(*m1);

if ( t ) {

can = 1; NEXTNM2(*mr0,*mr,s); CQ(*mr) = (t);

can = 1; NEXTNM2(*mr0,*mr,s); CP(*mr) = (t);

} else {

can = 2; FREENM(s);

}

Line 995 ND nd_add_q(ND p1,ND p2)

Line 1001 ND nd_add_q(ND p1,ND p2)

int n,c,can;

ND r;

NM m1,m2,mr0,mr,s;

Q t;

P t;

if ( !p1 ) return p2;

else if ( !p2 ) return p1;

Line 1008 ND nd_add_q(ND p1,ND p2)

Line 1014 ND nd_add_q(ND p1,ND p2)

#if defined(__x86_64__)

can += nm_destructive_add_q(&m1,&m2,&mr0,&mr);

#else

addq(CQ(m1),CQ(m2),&t);

addp(nd_vc,CP(m1),CP(m2),&t);

s = m1; m1 = NEXT(m1);

if ( t ) {

can++; NEXTNM2(mr0,mr,s); CQ(mr) = (t);

can++; NEXTNM2(mr0,mr,s); CP(mr) = (t);

} else {

can += 2; FREENM(s);

}

Line 1088 ND nd_add_sf(ND p1,ND p2)

Line 1094 ND nd_add_sf(ND p1,ND p2)

}

ND nd_reduce2(int mod,ND d,ND g,NDV p,NM mul,NDC dn)

{

int c,c1,c2;

Q cg,cred,gcd,tq;

P cgp,credp,gcdp;

Obj tr,tr1;

if ( mod == -1 )

CM(mul) = _mulsf(_invsf(HCM(p)),_chsgnsf(HCM(g)));

else if ( mod ) {

c1 = invm(HCM(p),mod); c2 = mod-HCM(g);

DMAR(c1,c2,0,mod,c); CM(mul) = c;

} else if ( nd_vc ) {

ezgcdpz(nd_vc,HCP(g),HCP(p),&gcdp);

divsp(nd_vc,HCP(g),gcdp,&cgp); divsp(nd_vc,HCP(p),gcdp,&credp);

chsgnp(cgp,&CP(mul));

nd_mul_c_q(d,credp); nd_mul_c_q(g,credp);

if ( dn ) {

mulr(nd_vc,(Obj)dn->r,(Obj)credp,&tr);

reductr(nd_vc,tr,&tr1); dn->r = (R)tr1;

}

} else {

igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred);

chsgnq(cg,&CQ(mul));

nd_mul_c_q(d,(P)cred); nd_mul_c_q(g,(P)cred);

if ( dn ) {

mulq(dn->z,cred,&tq); dn->z = tq;

}

return nd_add(mod,g,ndv_mul_nm(mod,mul,p));

}

/* ret=1 : success, ret=0 : overflow */

int nd_nf(int mod,ND g,NDV *ps,int full,NDC dn,ND *rp)

int nd_nf(int mod,ND d,ND g,NDV *ps,int full,NDC dn,ND *rp)

{

ND d;

NM m,mrd,tail;

NM mul;

int n,sugar,psugar,sugar0,stat,index;

Line 1100 int nd_nf(int mod,ND g,NDV *ps,int full,NDC dn,ND *rp)

Line 1137 int nd_nf(int mod,ND g,NDV *ps,int full,NDC dn,ND *rp)

NDV p,red;

Q cg,cred,gcd,tq,qq;

double hmag;

P tp,tp1;

Obj tr,tr1;

if ( dn ) {

if ( mod )

dn->m = 1;

else if ( nd_vc )

dn->r = (R)ONE;

else

dn->z = ONE;

}

if ( !g ) {

*rp = 0;

*rp = d;

return 1;

}

if ( !mod ) hmag = ((double)p_mag((P)HCQ(g)))*nd_scale;

if ( !mod ) hmag = ((double)p_mag(HCP(g)))*nd_scale;

sugar0 = sugar = SG(g);

n = NV(g);

mul = (NM)ALLOCA(sizeof(struct oNM)+(nd_wpd-1)*sizeof(UINT));

for ( d = 0; g; ) {

if ( d )

for ( tail = BDY(d); NEXT(tail); tail = NEXT(tail) );

for ( ; g; ) {

index = ndl_find_reducer(HDL(g));

if ( index >= 0 ) {

h = nd_psh[index];

Line 1125 int nd_nf(int mod,ND g,NDV *ps,int full,NDC dn,ND *rp)

Line 1168 int nd_nf(int mod,ND g,NDV *ps,int full,NDC dn,ND *rp)

nd_free(g); nd_free(d);

return 0;

}

if ( nd_demand )

p = nd_demand ? ndv_load(index) : ps[index];

p = ndv_load(index);

g = nd_reduce2(mod,d,g,p,mul,dn);

else

p = ps[index];

if ( mod == -1 )

CM(mul) = _mulsf(_invsf(HCM(p)),_chsgnsf(HCM(g)));

else if ( mod ) {

c1 = invm(HCM(p),mod); c2 = mod-HCM(g);

DMAR(c1,c2,0,mod,c); CM(mul) = c;

} else {

igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred);

chsgnq(cg,&CQ(mul));

nd_mul_c_q(d,cred); nd_mul_c_q(g,cred);

if ( dn ) {

mulq(dn->z,cred,&tq); dn->z = tq;

}

g = nd_add(mod,g,ndv_mul_nm(mod,mul,p));

sugar = MAX(sugar,SG(p)+TD(DL(mul)));

if ( !mod && hmag && g && ((double)(p_mag((P)HCQ(g))) > hmag) ) {

if ( !mod && g && ((double)(p_mag(HCP(g))) > hmag) ) {

tq = HCQ(g);

tp = HCP(g);

nd_removecont2(d,g);

if ( dn ) {

divq(tq,HCQ(g),&qq); divq(dn->z,qq,&tq); dn->z = tq;

if ( nd_vc ) {

divsp(nd_vc,tp,HCP(g),&tp1);

divr(nd_vc,(Obj)dn->r,(Obj)tp1,&tr);

reductr(nd_vc,(Obj)tr,&tr1); dn->r = (R)tr1;

} else {

divq((Q)tp,HCQ(g),&qq);

divq(dn->z,qq,&tq); dn->z = tq;

}

hmag = ((double)p_mag((P)HCQ(g)))*nd_scale;

hmag = ((double)p_mag(HCP(g)))*nd_scale;

}

} else if ( !full ) {

*rp = g;

Line 1231 int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *rp

Line 1265 int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *rp

} else {

igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred);

chsgnq(cg,&CQ(mul));

nd_mul_c_q(d,cred);

nd_mul_c_q(d,(P)cred);

mulq_pbucket(bucket,cred);

g = bucket->body[hindex];

gmag = (double)p_mag((P)HCQ(g));

Line 1286 int ndv_check_candidate(NODE input,int obpe,int oadv,E

Line 1320 int ndv_check_candidate(NODE input,int obpe,int oadv,E

NODE t,s;

union oNDC dn;

ndv_setup(0,0,cand,0);

ndv_setup(0,0,cand,0,1);

n = length(cand);

/* membercheck : list is a subset of Id(cand) ? */

Line 1297 again:

Line 1331 again:

else

r = (NDV)BDY(t);

d = ndvtond(0,r);

stat = nd_nf(0,d,nd_ps,0,0,&nf);

stat = nd_nf(0,0,d,nd_ps,0,0,&nf);

if ( !stat ) {

nd_reconstruct(0,0);

goto again;

Line 1404 void mulq_pbucket(PGeoBucket g,Q c)

Line 1438 void mulq_pbucket(PGeoBucket g,Q c)

int k;

for ( k = 0; k <= g->m; k++ )

nd_mul_c_q(g->body[k],c);

nd_mul_c_q(g->body[k],(P)c);

}

NM remove_head_pbucket_symbolic(PGeoBucket g)

Line 1539 ND normalize_pbucket(int mod,PGeoBucket g)

Line 1573 ND normalize_pbucket(int mod,PGeoBucket g)

return r;

}

#if 0

void register_hcf(NDV p)

{

DCP dc,t;

P hc,h;

int c;

NODE l,l1,prev;

hc = p->body->c.p;

if ( !nd_vc || NUM(hc) ) return;

fctrp(nd_vc,hc,&dc);

for ( t = dc; t; t = NEXT(t) ) {

h = t->c;

if ( NUM(h) ) continue;

for ( prev = 0, l = nd_hcf; l; prev = l, l = NEXT(l) ) {

c = compp(nd_vc,h,(P)BDY(l));

if ( c >= 0 ) break;

}

if ( !l || c > 0 ) {

MKNODE(l1,h,l);

if ( !prev )

nd_hcf = l1;

else

NEXT(prev) = l1;

}

#else

void register_hcf(NDV p)

{

DCP dc,t;

P hc,h,q;

Q dmy;

int c;

NODE l,l1,prev;

hc = p->body->c.p;

if ( NUM(hc) ) return;

ptozp(hc,1,&dmy,&h);

#if 1

for ( l = nd_hcf; l; l = NEXT(l) ) {

while ( 1 ) {

if ( divtpz(nd_vc,h,(P)BDY(l),&q) ) h = q;

else break;

}

if ( NUM(h) ) return;

#endif

for ( prev = 0, l = nd_hcf; l; prev = l, l = NEXT(l) ) {

c = compp(nd_vc,h,(P)BDY(l));

if ( c >= 0 ) break;

}

if ( !l || c > 0 ) {

MKNODE(l1,h,l);

if ( !prev )

nd_hcf = l1;

else

NEXT(prev) = l1;

}

#endif

int do_diagonalize(int sugar,int m)

{

int i,nh,stat;

Line 1546 int do_diagonalize(int sugar,int m)

Line 1642 int do_diagonalize(int sugar,int m)

ND h,nf,s,head;

NDV nfv;

Q q,num,den;

union oNDC dn;

P nm,nmp,dn,mnp,dnp;

for ( i = nd_psn-1; i >= 0 && SG(nd_psh[i]) == sugar; i-- ) {

if ( nd_demand )

Line 1555 int do_diagonalize(int sugar,int m)

Line 1651 int do_diagonalize(int sugar,int m)

nfv = nd_ps[i];

s = ndvtond(m,nfv);

s = nd_separate_head(s,&head);

stat = nd_nf(m,s,nd_ps,1,&dn,&nf);

stat = nd_nf(m,head,s,nd_ps,1,0,&nf);

if ( !stat ) return 0;

if ( !m ) {

NTOQ(NM(dn.z),SGN(dn.z),num);

mulq(HCQ(head),num,&q); HCQ(head) = q;

if ( DN(dn.z) ) {

NTOQ(DN(dn.z),1,den);

nd_mul_c_q(nf,den);

}

nf = nd_add(m,head,nf);

ndv_free(nfv);

nd_removecont(m,nf);

nfv = ndtondv(m,nf);

nd_free(nf);

nd_bound[i] = ndv_compute_bound(nfv);

if ( !m ) register_hcf(nfv);

if ( nd_demand ) {

ndv_save(nfv,i);

ndv_free(nfv);

Line 1623 again:

Line 1711 again:

goto again;

}

#if USE_GEOBUCKET

stat = m?nd_nf_pbucket(m,h,nd_ps,!Top,&nf):nd_nf(m,h,nd_ps,!Top,0,&nf);

stat = m?nd_nf_pbucket(m,h,nd_ps,!Top,&nf):nd_nf(m,0,h,nd_ps,!Top,0,&nf);

#else

stat = nd_nf(m,h,nd_ps,!Top,0,&nf);

stat = nd_nf(m,0,h,nd_ps,!Top,0,&nf);

#endif

if ( !stat ) {

NEXT(l) = d; d = l;

Line 1664 again:

Line 1752 again:

else

for ( t = g; t; t = NEXT(t) )

BDY(t) = (pointer)nd_ps[(int)BDY(t)];

if ( !checkonly && DP_Print ) { printf("nd_gb done.\n"); fflush(stdout); }

return g;

}

Line 1674 int do_diagonalize_trace(int sugar,int m)

Line 1763 int do_diagonalize_trace(int sugar,int m)

ND h,nf,nfq,s,head;

NDV nfv,nfqv;

Q q,den,num;

union oNDC dn;

for ( i = nd_psn-1; i >= 0 && SG(nd_psh[i]) == sugar; i-- ) {

/* for nd_ps */

Line 1694 int do_diagonalize_trace(int sugar,int m)

Line 1782 int do_diagonalize_trace(int sugar,int m)

nfv = nd_ps_trace[i];

s = ndvtond(0,nfv);

s = nd_separate_head(s,&head);

stat = nd_nf(0,s,nd_ps_trace,1,&dn,&nf);

stat = nd_nf(0,head,s,nd_ps_trace,1,0,&nf);

if ( !stat ) return 0;

NTOQ(NM(dn.z),SGN(dn.z),num);

mulq(HCQ(head),num,&q); HCQ(head) = q;

if ( DN(dn.z) ) {

NTOQ(DN(dn.z),1,den);

nd_mul_c_q(nf,den);

}

nf = nd_add(0,head,nf);

ndv_free(nfv);

nd_removecont(0,nf);

nfv = ndtondv(0,nf);

nd_free(nf);

nd_bound[i] = ndv_compute_bound(nfv);

register_hcf(nfv);

if ( nd_demand ) {

ndv_save(nfv,i);

ndv_free(nfv);

Line 1720 int do_diagonalize_trace(int sugar,int m)

Line 1802 int do_diagonalize_trace(int sugar,int m)

static struct oEGT eg_invdalg;

struct oEGT eg_le;

void nd_subst_vector(VL vl,P p,NODE subst,P *r)

{

NODE tn;

P p1;

for ( tn = subst; tn; tn = NEXT(NEXT(tn)) ) {

substp(vl,p,BDY(tn),BDY(NEXT(tn)),&p1); p = p1;

}

*r = p;

}

NODE nd_gb_trace(int m,int ishomo)

{

int i,nh,sugar,stat;

Line 1729 NODE nd_gb_trace(int m,int ishomo)

Line 1822 NODE nd_gb_trace(int m,int ishomo)

ND h,nf,nfq,s,head;

NDV nfv,nfqv;

Q q,den,num;

P hc;

union oNDC dn;

struct oEGT eg_monic,egm0,egm1;

int diag_count = 0;

Line 1748 again:

Line 1842 again:

if ( SG(l) != sugar ) {

#if 1

if ( ishomo ) {

if ( DP_Print > 2 ) fprintf(asir_out,"|");

stat = do_diagonalize_trace(sugar,m);

if ( DP_Print > 2 ) fprintf(asir_out,"|");

diag_count = 0;

if ( !stat ) {

NEXT(l) = d; d = l;

Line 1769 again:

Line 1865 again:

#if USE_GEOBUCKET

stat = nd_nf_pbucket(m,h,nd_ps,!Top,&nf);

#else

stat = nd_nf(m,h,nd_ps,!Top,0,&nf);

stat = nd_nf(m,0,h,nd_ps,!Top,0,&nf);

#endif

if ( !stat ) {

NEXT(l) = d; d = l;

Line 1782 again:

Line 1878 again:

} else

nfq = 0;

if ( !nfq ) {

if ( !nd_sp(0,1,l,&h) || !nd_nf(0,h,nd_ps_trace,!Top,0,&nfq) ) {

if ( !nd_sp(0,1,l,&h) || !nd_nf(0,0,h,nd_ps_trace,!Top,0,&nfq) ) {

NEXT(l) = d; d = l;

d = nd_reconstruct(1,d);

goto again;

Line 1790 again:

Line 1886 again:

}

if ( nfq ) {

/* m|HC(nfq) => failure */

if ( !rem(NM(HCQ(nfq)),m) ) return 0;

if ( nd_vc ) {

nd_subst_vector(nd_vc,HCP(nfq),nd_subst,&hc); q = (Q)hc;

} else

q = HCQ(nfq);

if ( !rem(NM(q),m) ) return 0;

if ( DP_Print ) { printf("+"); fflush(stdout); }

if ( nd_nalg ) {

Line 1807 again:

Line 1907 again:

nh = ndv_newps(0,nfv,nfqv);

if ( ishomo && ++diag_count == diag_period ) {

diag_count = 0;

if ( DP_Print > 2 ) fprintf(asir_out,"|");

stat = do_diagonalize_trace(sugar,m);

if ( DP_Print > 2 ) fprintf(asir_out,"|");

if ( !stat ) {

NEXT(l) = d; d = l;

d = nd_reconstruct(1,d);

Line 1835 again:

Line 1937 again:

print_eg("invdalg",&eg_invdalg);

print_eg("le",&eg_le);

}

if ( DP_Print ) { printf("nd_gb_trace done.\n"); fflush(stdout); }

return g;

}

Line 1865 NODE ndv_reduceall(int m,NODE f)

Line 1968 NODE ndv_reduceall(int m,NODE f)

(int (*)(const void *,const void *))ndv_compare);

for ( t = f, i = 0; t; i++, t = NEXT(t) ) BDY(t) = (pointer)w[i];

#endif

ndv_setup(m,0,f,0);

ndv_setup(m,0,f,0,1);

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

g = ndvtond(m,nd_ps[i]);

g = nd_separate_head(g,&head);

stat = nd_nf(m,g,nd_ps,1,&dn,&nf);

stat = nd_nf(m,head,g,nd_ps,1,0,&nf);

if ( !stat )

nd_reconstruct(0,0);

else {

if ( DP_Print ) { printf("."); fflush(stdout); }

if ( !m ) {

NTOQ(NM(dn.z),SGN(dn.z),num);

mulq(HCQ(head),num,&q); HCQ(head) = q;

if ( DN(dn.z) ) {

NTOQ(DN(dn.z),1,den);

nd_mul_c_q(nf,den);

}

nf = nd_add(m,head,nf);

ndv_free(nd_ps[i]);

nd_removecont(m,nf);

nd_ps[i] = ndtondv(m,nf); nd_free(nf);

Line 2187 int ndv_newps(int m,NDV a,NDV aq)

Line 2281 int ndv_newps(int m,NDV a,NDV aq)

nd_ps[nd_psn] = a;

if ( aq ) {

nd_ps_trace[nd_psn] = aq;

register_hcf(aq);

nd_bound[nd_psn] = ndv_compute_bound(aq);

SG(r) = SG(aq); ndl_copy(HDL(aq),DL(r));

} else {

if ( !m ) register_hcf(a);

nd_bound[nd_psn] = ndv_compute_bound(a);

SG(r) = SG(a); ndl_copy(HDL(a),DL(r));

}

Line 2205 int ndv_newps(int m,NDV a,NDV aq)

Line 2301 int ndv_newps(int m,NDV a,NDV aq)

return nd_psn++;

}

void ndv_setup(int mod,int trace,NODE f,int dont_sort)

void ndv_setup(int mod,int trace,NODE f,int dont_sort,int dont_removecont)

{

int i,j,td,len,max;

NODE s,s0,f0;

Line 2233 void ndv_setup(int mod,int trace,NODE f,int dont_sort)

Line 2329 void ndv_setup(int mod,int trace,NODE f,int dont_sort)

nd_ps_trace = (NDV *)MALLOC(nd_pslen*sizeof(NDV));

nd_psh = (RHist *)MALLOC(nd_pslen*sizeof(RHist));

nd_bound = (UINT **)MALLOC(nd_pslen*sizeof(UINT *));

nd_hcf = 0;

if ( trace && nd_vc )

makesubst(nd_vc,&nd_subst);

else

nd_subst = 0;

if ( !nd_red )

nd_red = (RHist *)MALLOC(REDTAB_LEN*sizeof(RHist));

for ( i = 0; i < REDTAB_LEN; i++ ) nd_red[i] = 0;

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

if ( trace ) {

a = nd_ps_trace[i] = ndv_dup(0,w[i]);

ndv_removecont(0,a);

if ( !dont_removecont) ndv_removecont(0,a);

register_hcf(a);

am = nd_ps[i] = ndv_dup(mod,a);

ndv_mod(mod,am);

ndv_removecont(mod,am);

} else {

a = nd_ps[i] = ndv_dup(mod,w[i]);

ndv_removecont(mod,a);

if ( mod || !dont_removecont ) ndv_removecont(mod,a);

if ( !mod ) register_hcf(a);

}

NEWRHist(r); SG(r) = HTD(a); ndl_copy(HDL(a),DL(r));

nd_bound[i] = ndv_compute_bound(a);

Line 2386 void nd_gr(LIST f,LIST v,int m,int f4,struct order_spe

Line 2490 void nd_gr(LIST f,LIST v,int m,int f4,struct order_spe

if ( !m && Demand ) nd_demand = 1;

else nd_demand = 0;

if ( DP_Multiple )

nd_scale = ((double)DP_Multiple)/(double)(Denominator?Denominator:1);

#if 0

ndv_alloc = 0;

#endif

get_vars((Obj)f,&fv); pltovl(v,&vv);

get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc);

for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ );

switch ( ord->id ) {

case 1:

Line 2431 void nd_gr(LIST f,LIST v,int m,int f4,struct order_spe

Line 2537 void nd_gr(LIST f,LIST v,int m,int f4,struct order_spe

if ( b ) { NEXTNODE(fd0,fd); BDY(fd) = (pointer)b; }

}

if ( fd0 ) NEXT(fd) = 0;

ndv_setup(m,0,fd0,0);

ndv_setup(m,0,fd0,0,0);

x = f4?nd_f4(m):nd_gb(m,ishomo,0);

nd_demand = 0;

x = ndv_reducebase(x);

Line 2464 void nd_gr_postproc(LIST f,LIST v,int m,struct order_s

Line 2570 void nd_gr_postproc(LIST f,LIST v,int m,struct order_s

NumberField nf;

struct order_spec *ord1;

get_vars((Obj)f,&fv); pltovl(v,&vv);

get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc);

for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ );

switch ( ord->id ) {

case 1:

Line 2506 void nd_gr_postproc(LIST f,LIST v,int m,struct order_s

Line 2612 void nd_gr_postproc(LIST f,LIST v,int m,struct order_s

if ( b ) { NEXTNODE(fd0,fd); BDY(fd) = (pointer)b; }

}

if ( fd0 ) NEXT(fd) = 0;

ndv_setup(m,0,fd0,0);

ndv_setup(m,0,fd0,0,1);

for ( x = 0, i = 0; i < nd_psn; i++ )

x = update_base(x,i);

if ( do_check ) {

Line 2550 void nd_gr_trace(LIST f,LIST v,int trace,int homo,int

Line 2656 void nd_gr_trace(LIST f,LIST v,int trace,int homo,int

struct order_spec *ord1;

struct oEGT eg_check,eg0,eg1;

get_vars((Obj)f,&fv); pltovl(v,&vv);

if ( DP_Multiple )

nd_scale = ((double)DP_Multiple)/(double)(Denominator?Denominator:1);

get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc);

for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ );

switch ( ord->id ) {

case 1:

Line 2621 void nd_gr_trace(LIST f,LIST v,int trace,int homo,int

Line 2730 void nd_gr_trace(LIST f,LIST v,int trace,int homo,int

while ( 1 ) {

if ( Demand )

nd_demand = 1;

ndv_setup(m,1,fd0,0);

ndv_setup(m,1,fd0,0,0);

cand = f4?nd_f4_trace(m):nd_gb_trace(m,ishomo || homo);

if ( !cand ) {

/* failure */

Line 2800 void nd_removecont(int mod,ND p)

Line 2909 void nd_removecont(int mod,ND p)

v.len = n;

v.body = (pointer *)w;

for ( m = BDY(p), i = 0; i < n; m = NEXT(m), i++ ) w[i] = CQ(m);

removecont_array(w,n);

removecont_array((P *)w,n,1);

for ( m = BDY(p), i = 0; i < n; m = NEXT(m), i++ ) CQ(m) = w[i];

}

Line 2814 void nd_removecont2(ND p1,ND p2)

Line 2923 void nd_removecont2(ND p1,ND p2)

struct oVECT v;

N q,r;

if ( !p1 ) {

nd_removecont(0,p2); return;

} else if ( !p2 ) {

nd_removecont(0,p1); return;

}

n1 = nd_length(p1);

n2 = nd_length(p2);

n = n1+n2;

w = (Q *)ALLOCA(n*sizeof(Q));

v.len = n;

v.body = (pointer *)w;

for ( m = BDY(p1), i = 0; i < n1; m = NEXT(m), i++ ) w[i] = CQ(m);

i = 0;

for ( m = BDY(p2); i < n; m = NEXT(m), i++ ) w[i] = CQ(m);

if ( p1 )

removecont_array(w,n);

for ( m = BDY(p1); i < n1; m = NEXT(m), i++ ) w[i] = CQ(m);

for ( m = BDY(p1), i = 0; i < n1; m = NEXT(m), i++ ) CQ(m) = w[i];

if ( p2 )

for ( m = BDY(p2); i < n; m = NEXT(m), i++ ) CQ(m) = w[i];

for ( m = BDY(p2); i < n; m = NEXT(m), i++ ) w[i] = CQ(m);

removecont_array((P *)w,n,1);

i = 0;

if ( p1 )

for ( m = BDY(p1); i < n1; m = NEXT(m), i++ ) CQ(m) = w[i];

if ( p2 )

for ( m = BDY(p2); i < n; m = NEXT(m), i++ ) CQ(m) = w[i];

}

void ndv_removecont(int mod,NDV p)

{

int i,len;

int i,len,all_p;

Q *w;

Q *c;

P *w;

Q dvr,t;

P g,cont,tp;

NMV m;

if ( mod == -1 )

Line 2845 void ndv_removecont(int mod,NDV p)

Line 2957 void ndv_removecont(int mod,NDV p)

ndv_mul_c(mod,p,invm(HCM(p),mod));

else {

len = p->len;

w = (Q *)ALLOCA(len*sizeof(Q));

w = (P *)ALLOCA(len*sizeof(P));

for ( m = BDY(p), i = 0; i < len; NMV_ADV(m), i++ ) w[i] = CQ(m);

c = (Q *)ALLOCA(len*sizeof(Q));

sortbynm(w,len);

for ( m = BDY(p), all_p = 1, i = 0; i < len; NMV_ADV(m), i++ ) {

qltozl(w,len,&dvr);

ptozp(CP(m),1,&c[i],&w[i]);

for ( m = BDY(p), i = 0; i < len; NMV_ADV(m), i++ ) {

all_p = all_p && !NUM(w[i]);

divq(CQ(m),dvr,&t); CQ(m) = t;

}

if ( all_p ) {

qltozl(c,len,&dvr); nd_heu_nezgcdnpz(nd_vc,w,len,1,&g);

mulp(nd_vc,(P)dvr,g,&cont);

for ( m = BDY(p), i = 0; i < len; NMV_ADV(m), i++ ) {

divsp(nd_vc,CP(m),cont,&tp); CP(m) = tp;

}

} else {

sortbynm((Q *)c,len);

qltozl((Q *)c,len,&dvr);

for ( m = BDY(p), i = 0; i < len; NMV_ADV(m), i++ ) {

divsp(nd_vc,CP(m),(P)dvr,&tp); CP(m) = tp;

}

Line 2901 void ndv_dehomogenize(NDV p,struct order_spec *ord)

Line 3025 void ndv_dehomogenize(NDV p,struct order_spec *ord)

NV(p)--;

}

void removecont_array(Q *c,int n)

void nd_heu_nezgcdnpz(VL vl,P *pl,int m,int full,P *pr)

{

int i;

P *tpl,*tpl1;

NODE l;

P h,gcd,t;

tpl = (P *)ALLOCA(m*sizeof(P));

tpl1 = (P *)ALLOCA(m*sizeof(P));

bcopy(pl,tpl,m*sizeof(P));

gcd = (P)ONE;

for ( l = nd_hcf; l; l = NEXT(l) ) {

h = (P)BDY(l);

while ( 1 ) {

for ( i = 0; i < m; i++ )

if ( !divtpz(vl,tpl[i],h,&tpl1[i]) )

break;

if ( i == m ) {

bcopy(tpl1,tpl,m*sizeof(P));

mulp(vl,gcd,h,&t); gcd = t;

} else

break;

}

if ( DP_Print > 2 ){fprintf(asir_out,"[%d]",nmonop(gcd)); fflush(asir_out);}

if ( full ) {

heu_nezgcdnpz(vl,tpl,m,&t);

mulp(vl,gcd,t,pr);

} else

*pr = gcd;

}

void removecont_array(P *p,int n,int full)

{

int all_p,all_q,i;

Q *c;

P *w;

P t,s;

for ( all_q = 1, i = 0; i < n; i++ )

all_q = all_q && NUM(p[i]);

if ( all_q ) {

removecont_array_q((Q *)p,n);

} else {

c = (Q *)ALLOCA(n*sizeof(Q));

w = (P *)ALLOCA(n*sizeof(P));

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

ptozp(p[i],1,&c[i],&w[i]);

}

removecont_array_q(c,n);

nd_heu_nezgcdnpz(nd_vc,w,n,full,&t);

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

divsp(nd_vc,w[i],t,&s); mulp(nd_vc,s,(P)c[i],&p[i]);

}

void removecont_array_q(Q *c,int n)

{

struct oVECT v;

Q d0,d1,a,u,u1,gcd;

int i,j;

Line 2952 void nd_mul_c(int mod,ND p,int mul)

Line 3133 void nd_mul_c(int mod,ND p,int mul)

}

void nd_mul_c_q(ND p,Q mul)

void nd_mul_c_q(ND p,P mul)

{

NM m;

Q c;

P c;

if ( !p ) return;

if ( UNIQ(mul) ) return;

for ( m = BDY(p); m; m = NEXT(m) ) {

mulq(CQ(m),mul,&c); CQ(m) = c;

mulp(nd_vc,CP(m),mul,&c); CP(m) = c;

}

Line 3250 ND nd_copy(ND p)

Line 3431 ND nd_copy(ND p)

int nd_sp(int mod,int trace,ND_pairs p,ND *rp)

{

NM m;

NM m1,m2;

NDV p1,p2;

ND t1,t2;

UINT *lcm;

P gp,tp;

Q g,t;

int td;

if ( !mod && nd_demand ) {

Line 3266 int nd_sp(int mod,int trace,ND_pairs p,ND *rp)

Line 3449 int nd_sp(int mod,int trace,ND_pairs p,ND *rp)

}

lcm = LCM(p);

NEWNM(m);

NEWNM(m1); ndl_sub(lcm,HDL(p1),DL(m1));

CQ(m) = HCQ(p2);

if ( ndl_check_bound2(p->i1,DL(m1)) ) {

ndl_sub(lcm,HDL(p1),DL(m));

FREENM(m1); return 0;

if ( ndl_check_bound2(p->i1,DL(m)) )

return 0;

t1 = ndv_mul_nm(mod,m,p1);

if ( mod == -1 ) CM(m) = _chsgnsf(HCM(p1));

else if ( mod ) CM(m) = mod-HCM(p1);

else chsgnq(HCQ(p1),&CQ(m));

ndl_sub(lcm,HDL(p2),DL(m));

if ( ndl_check_bound2(p->i2,DL(m)) ) {

nd_free(t1);

return 0;

}

t2 = ndv_mul_nm(mod,m,p2);

NEWNM(m2); ndl_sub(lcm,HDL(p2),DL(m2));

if ( ndl_check_bound2(p->i2,DL(m2)) ) {

FREENM(m1); FREENM(m2); return 0;

}

if ( mod == -1 ) {

CM(m1) = HCM(p2); CM(m2) = _chsgnsf(HCM(p1));

} else if ( mod ) {

CM(m1) = HCM(p2); CM(m2) = mod-HCM(p1);

} else if ( nd_vc ) {

ezgcdpz(nd_vc,HCP(p1),HCP(p2),&gp);

divsp(nd_vc,HCP(p2),gp,&CP(m1));

divsp(nd_vc,HCP(p1),gp,&tp); chsgnp(tp,&CP(m2));

} else {

igcd_cofactor(HCQ(p1),HCQ(p2),&g,&t,&CQ(m1)); chsgnq(t,&CQ(m2));

}

t1 = ndv_mul_nm(mod,m1,p1); t2 = ndv_mul_nm(mod,m2,p2);

*rp = nd_add(mod,t1,t2);

FREENM(m);

FREENM(m1); FREENM(m2);

return 1;

}

Line 3484 ND ndv_mul_nm(int mod,NM m0,NDV p)

Line 3673 ND ndv_mul_nm(int mod,NM m0,NDV p)

NMV m;

UINT *d,*dt,*dm;

int c,n,td,i,c1,c2,len;

Q q;

P q;

ND r;

if ( !p ) return 0;

Line 3516 ND ndv_mul_nm(int mod,NM m0,NDV p)

Line 3705 ND ndv_mul_nm(int mod,NM m0,NDV p)

ndl_add(DL(m),d,DL(mr));

}

} else {

q = CQ(m0);

q = CP(m0);

for ( i = 0; i < len; i++, NMV_ADV(m) ) {

NEXTNM(mr0,mr);

mulq(CQ(m),q,&CQ(mr));

mulp(nd_vc,CP(m),q,&CP(mr));

ndl_add(DL(m),d,DL(mr));

}

Line 3662 void ndv_mod(int mod,NDV p)

Line 3851 void ndv_mod(int mod,NDV p)

NMV t,d;

int r,s,u;

int i,len,dlen;

P cp;

Q c;

Obj gfs;

if ( !p ) return;

Line 3678 void ndv_mod(int mod,NDV p)

Line 3869 void ndv_mod(int mod,NDV p)

}

else

for ( t = d = BDY(p), i = 0; i < len; i++, NMV_ADV(t) ) {

r = rem(NM(CQ(t)),mod);

if ( nd_vc ) {

nd_subst_vector(nd_vc,CP(t),nd_subst,&cp);

c = (Q)cp;

} else

c = CQ(t);

r = rem(NM(c),mod);

if ( r ) {

if ( SGN(CQ(t)) < 0 )

if ( SGN(c) < 0 )

r = mod-r;

if ( DN(CQ(t)) ) {

if ( DN(c) ) {

s = rem(DN(CQ(t)),mod);

s = rem(DN(c),mod);

if ( !s )

error("ndv_mod : division by 0");

s = invm(s,mod);

Line 4043 void nd_nf_p(P f,LIST g,LIST v,int m,struct order_spec

Line 4239 void nd_nf_p(P f,LIST g,LIST v,int m,struct order_spec

VL vv,tv;

int stat,nvar,max,e;

union oNDC dn;

Q cont;

P pp;

if ( !f ) {

*rp = 0;

Line 4068 void nd_nf_p(P f,LIST g,LIST v,int m,struct order_spec

Line 4266 void nd_nf_p(P f,LIST g,LIST v,int m,struct order_spec

/* conversion to ndv */

for ( in0 = 0, t = BDY(g); t; t = NEXT(t) ) {

NEXTNODE(in0,in);

BDY(in) = (pointer)ptondv(CO,vv,(P)BDY(t));

ptozp((P)BDY(t),1,&cont,&pp);

BDY(in) = (pointer)ptondv(CO,vv,pp);

if ( m ) ndv_mod(m,(NDV)BDY(in));

}

NEXTNODE(in0,in);

Line 4076 void nd_nf_p(P f,LIST g,LIST v,int m,struct order_spec

Line 4275 void nd_nf_p(P f,LIST g,LIST v,int m,struct order_spec

if ( m ) ndv_mod(m,(NDV)BDY(in));

NEXT(in) = 0;

/* dont sort */

/* dont sort, dont removecont */

ndv_setup(m,0,in0,1);

ndv_setup(m,0,in0,1,1);

nd_psn--;

nd_scale=2;

while ( 1 ) {

nd = (pointer)ndvtond(m,nd_ps[nd_psn]);

stat = nd_nf(m,nd,nd_ps,1,0,&nf);

stat = nd_nf(m,0,nd,nd_ps,1,0,&nf);

if ( !stat ) {

nd_psn++;

nd_reconstruct(0,0);

Line 4275 int ndv_reduce_vect_q(Q *svect,int trace,int col,IndAr

Line 4474 int ndv_reduce_vect_q(Q *svect,int trace,int col,IndAr

if ( j == col ) break;

if ( hmag && ((double)p_mag((P)svect[j]) > hmag) ) {

nz = compress_array(svect,cvect,col);

removecont_array(cvect,nz);

removecont_array((P *)cvect,nz,1);

expand_array(svect,cvect,nz);

hmag = ((double)p_mag((P)svect[j]))*nd_scale;

}

nz = compress_array(svect,cvect,col);

removecont_array(cvect,nz);

removecont_array((P *)cvect,nz,1);

expand_array(svect,cvect,nz);

if ( DP_Print ) {

fprintf(asir_out,"-"); fflush(asir_out);

Line 5607 void nd_det(int mod,MAT f,P *rp)

Line 5806 void nd_det(int mod,MAT f,P *rp)

sgn = -sgn;

}

bound = nd_det_compute_bound(dm,n,j);

if ( ndl_check_bound(bound,bound) )

for ( k = 0; k < nd_nvar; k++ )

if ( bound[k]*2 > nd_mask0 ) break;

if ( k < nd_nvar )

nd_det_reconstruct(dm,n,j,d);

for ( i = j+1, mj = dm[j], mjj = mj[j]; i < n; i++ ) {

Line 5737 void nd_det_reconstruct(NDV **dm,int n,int j,NDV d)

Line 5938 void nd_det_reconstruct(NDV **dm,int n,int j,NDV d)

#endif

}

/* returns a UINT array containing degree bounds */

UINT *nd_det_compute_bound(NDV **dm,int n,int j)

{

UINT *d0,*d1,*d,*t,*r;

int k,l;

int k,l,i;

d0 = (UINT *)ALLOCA(nd_wpd*sizeof(UINT));

d0 = (UINT *)MALLOC(nd_nvar*sizeof(UINT));

d1 = (UINT *)ALLOCA(nd_wpd*sizeof(UINT));

for ( k = 0; k < nd_nvar; k++ ) d0[k] = 0;

for ( k = 0; k < nd_wpd; k++ ) d0[k] = 0;

for ( k = j; k < n; k++ )

for ( l = j; l < n; l++ )

if ( dm[k][l] ) {

d = ndv_compute_bound(dm[k][l]);

ndl_lcm(d,d0,d1);

for ( i = 0; i < nd_nvar; i++ )

t = d1; d1 = d0; d0 = t;

d0[i] = MAX(d0[i],d[i]);

}

r = (UINT *)ALLOCA(nd_wpd*sizeof(UINT));

return d0;

for ( k = 0; k < nd_wpd; k++ ) r[k] = d0[k];

return r;

}

DL nd_separate_d(UINT *d,UINT *trans)

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