Return to array.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / builtin |
version 1.48, 2005/11/27 05:37:53 | version 1.50, 2006/01/05 00:21:20 | ||
---|---|---|---|
|
|
||
* 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/asir2000/builtin/array.c,v 1.47 2005/11/27 00:07:05 noro Exp $ | * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.49 2005/12/21 23:18:15 noro Exp $ | ||
*/ | */ | ||
#include "ca.h" | #include "ca.h" | ||
#include "base.h" | #include "base.h" | ||
|
|
||
void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(); | void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(); | ||
void Pinvmat(); | void Pinvmat(); | ||
void Pnewbytearray(); | void Pnewbytearray(),Pmemoryplot_to_coord(); | ||
void Pgeneric_gauss_elim(); | void Pgeneric_gauss_elim(); | ||
void Pgeneric_gauss_elim_mod(); | void Pgeneric_gauss_elim_mod(); | ||
|
|
||
{"matr",Pmat,-99999999}, | {"matr",Pmat,-99999999}, | ||
{"matc",Pmatc,-99999999}, | {"matc",Pmatc,-99999999}, | ||
{"newbytearray",Pnewbytearray,-2}, | {"newbytearray",Pnewbytearray,-2}, | ||
{"memoryplot_to_coord",Pmemoryplot_to_coord,1}, | |||
{"sepmat_destructive",Psepmat_destructive,2}, | {"sepmat_destructive",Psepmat_destructive,2}, | ||
{"sepvect",Psepvect,2}, | {"sepvect",Psepvect,2}, | ||
{"qsort",Pqsort,-2}, | {"qsort",Pqsort,-2}, | ||
|
|
||
*rp = array; | *rp = array; | ||
} | } | ||
#define MEMORY_GETPOINT(a,len,x,y) (((a)[(len)*(y)+((x)>>3)])&(1<<((x)&7))) | |||
void Pmemoryplot_to_coord(NODE arg,LIST *rp) | |||
{ | |||
int len,blen,y,i,j; | |||
char *a; | |||
NODE r0,r,n; | |||
LIST l; | |||
BYTEARRAY ba; | |||
Q iq,jq; | |||
asir_assert(ARG0(arg),O_LIST,"memoryplot_to_coord"); | |||
arg = BDY((LIST)ARG0(arg)); | |||
len = QTOS((Q)ARG0(arg)); | |||
blen = (len+7)/8; | |||
y = QTOS((Q)ARG1(arg)); | |||
ba = (BYTEARRAY)ARG2(arg); a = ba->body; | |||
r0 = 0; | |||
for ( j = 0; j < y; j++ ) | |||
for ( i = 0; i < len; i++ ) | |||
if ( MEMORY_GETPOINT(a,blen,i,j) ) { | |||
NEXTNODE(r0,r); | |||
STOQ(i,iq); STOQ(j,jq); | |||
n = mknode(2,iq,jq); | |||
MKLIST(l,n); | |||
BDY(r) = l; | |||
} | |||
if ( r0 ) NEXT(r) = 0; | |||
MKLIST(*rp,r0); | |||
} | |||
void Pnewmat(NODE arg,MAT *rp) | void Pnewmat(NODE arg,MAT *rp) | ||
{ | { | ||
int row,col; | int row,col; | ||
|
|
||
ret = intmtoratm_q(xmat,NM(q),*nmmat,dn); | ret = intmtoratm_q(xmat,NM(q),*nmmat,dn); | ||
get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1); | get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1); | ||
if ( ret ) { | if ( ret ) { | ||
rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); | |||
cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int)); | |||
for ( j = k = l = 0; j < col; j++ ) | for ( j = k = l = 0; j < col; j++ ) | ||
if ( cinfo[j] ) | if ( cinfo[j] ) | ||
rind[k++] = j; | rind[k++] = j; | ||
else | else | ||
cind[l++] = j; | |||
get_eg(&tmp0); | |||
ret = gensolve_check(mat,*nmmat,*dn,rind,cind); | |||
get_eg(&tmp1); add_eg(&eg_check,&tmp0,&tmp1); | |||
if ( ret ) { | |||
if ( DP_Print > 3 ) { | |||
fprintf(stderr,"\n"); | |||
print_eg("INV",&eg_inv); | |||
print_eg("MUL",&eg_mul); | |||
print_eg("INTRAT",&eg_intrat); | |||
print_eg("CHECK",&eg_check); | |||
fflush(asir_out); | |||
} | |||
*rindp = rind; | |||
*cindp = cind; | |||
for ( j = k = 0; j < col; j++ ) | |||
if ( !cinfo[j] ) | |||
cind[k++] = j; | |||
return rank; | |||
} | |||
} else { | |||
period = period*3/2; | |||
count = 0; | |||
nsize += period; | |||
wxsize += rank*ri*nsize; | |||
wx = (int *)REALLOC(wx,wxsize*sizeof(int)); | |||
for ( i = 0; i < wxsize; i++ ) wx[i] = 0; | |||
} | |||
} | |||
} | |||
} | |||
} | |||
int generic_gauss_elim_hensel_dalg(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp) | |||
{ | |||
MAT bmat,xmat; | |||
Q **a0,**a,**b,**x,**nm; | |||
Q *ai,*bi,*xi; | |||
int row,col; | |||
int **w; | |||
int *wi; | |||
int **wc; | |||
Q mdq,q,s,u; | |||
N tn; | |||
int ind,md,i,j,k,l,li,ri,rank; | |||
unsigned int t; | |||
int *cinfo,*rinfo; | |||
int *rind,*cind; | |||
int count; | |||
int ret; | |||
struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1; | |||
int period; | |||
int *wx,*ptr; | |||
int wxsize,nsize; | |||
N wn; | |||
Q wq; | |||
NumberField nf; | |||
DP *mb; | |||
DP m; | |||
int col1; | |||
nf = get_numberfield(); | |||
mb = nf->mb; | |||
a0 = (Q **)mat->body; | |||
row = mat->row; col = mat->col; | |||
w = (int **)almat(row,col); | |||
for ( ind = 0; ; ind++ ) { | |||
md = get_lprime(ind); | |||
STOQ(md,mdq); | |||
for ( i = 0; i < row; i++ ) | |||
for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ ) | |||
if ( q = (Q)ai[j] ) { | |||
t = rem(NM(q),md); | |||
if ( t && SGN(q) < 0 ) | |||
t = (md - t) % md; | |||
wi[j] = t; | |||
} else | |||
wi[j] = 0; | |||
if ( DP_Print ) { | |||
fprintf(asir_out,"LU decomposition.."); fflush(asir_out); | |||
} | |||
rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo); | |||
if ( DP_Print ) { | |||
fprintf(asir_out,"done.\n"); fflush(asir_out); | |||
} | |||
for ( i = 0; i < col-1; i++ ) { | |||
if ( !cinfo[i] ) { | |||
m = mb[i]; | |||
for ( j = i+1; j < col-1; j++ ) | |||
if ( dp_redble(mb[j],m) ) | |||
cinfo[j] = -1; | |||
} | |||
} | |||
a = (Q **)almat_pointer(rank,rank); /* lhs mat */ | |||
MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */ | |||
for ( j = li = ri = 0; j < col; j++ ) | |||
if ( cinfo[j] > 0 ) { | |||
/* the column is in lhs */ | |||
for ( i = 0; i < rank; i++ ) { | |||
w[i][li] = w[i][j]; | |||
a[i][li] = a0[rinfo[i]][j]; | |||
} | |||
li++; | |||
} else if ( !cinfo[j] ) { | |||
/* the column is in rhs */ | |||
for ( i = 0; i < rank; i++ ) | |||
b[i][ri] = a0[rinfo[i]][j]; | |||
ri++; | |||
} | |||
/* solve Ax+B=0; A: rank x rank, B: rank x ri */ | |||
MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body; | |||
MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body; | |||
/* use the right part of w as work area */ | |||
wc = (int **)almat(rank,ri); | |||
for ( i = 0; i < rank; i++ ) | |||
wc[i] = w[i]+rank; | |||
*rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); | |||
*cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int)); | |||
init_eg(&eg_mul); init_eg(&eg_inv); | |||
init_eg(&eg_check); init_eg(&eg_intrat); | |||
period = F4_INTRAT_PERIOD; | |||
nsize = period; | |||
wxsize = rank*ri*nsize; | |||
wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int)); | |||
for ( i = 0; i < wxsize; i++ ) wx[i] = 0; | |||
for ( q = ONE, count = 0; ; ) { | |||
if ( DP_Print ) | |||
fprintf(stderr,"o"); | |||
/* wc = -b mod md */ | |||
get_eg(&tmp0); | |||
for ( i = 0; i < rank; i++ ) | |||
for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) | |||
if ( u = (Q)bi[j] ) { | |||
t = rem(NM(u),md); | |||
if ( t && SGN(u) > 0 ) | |||
t = (md - t) % md; | |||
wi[j] = t; | |||
} else | |||
wi[j] = 0; | |||
/* wc = A^(-1)wc; wc is not normalized */ | |||
solve_by_lu_mod(w,rank,md,wc,ri,0); | |||
/* wx += q*wc */ | |||
ptr = wx; | |||
for ( i = 0; i < rank; i++ ) | |||
for ( j = 0, wi = wc[i]; j < ri; j++ ) { | |||
if ( wi[j] ) | |||
muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr); | |||
ptr += nsize; | |||
} | |||
count++; | |||
get_eg(&tmp1); | |||
add_eg(&eg_inv,&tmp0,&tmp1); | |||
get_eg(&tmp0); | |||
for ( i = 0; i < rank; i++ ) | |||
for ( j = 0; j < ri; j++ ) { | |||
inner_product_mat_int_mod(a,wc,rank,i,j,&u); | |||
addq(b[i][j],u,&s); | |||
if ( s ) { | |||
t = divin(NM(s),md,&tn); | |||
if ( t ) | |||
error("generic_gauss_elim_hensel:incosistent"); | |||
NTOQ(tn,SGN(s),b[i][j]); | |||
} else | |||
b[i][j] = 0; | |||
} | |||
get_eg(&tmp1); | |||
add_eg(&eg_mul,&tmp0,&tmp1); | |||
/* q = q*md */ | |||
mulq(q,mdq,&u); q = u; | |||
if ( count == period ) { | |||
get_eg(&tmp0); | |||
ptr = wx; | |||
for ( i = 0; i < rank; i++ ) | |||
for ( j = 0, xi = x[i]; j < ri; | |||
j++, ptr += nsize ) { | |||
for ( k = nsize-1; k >= 0 && !ptr[k]; k-- ); | |||
if ( k >= 0 ) { | |||
wn = NALLOC(k+1); | |||
PL(wn) = k+1; | |||
for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l]; | |||
NTOQ(wn,1,wq); | |||
subq(xi[j],wq,&u); xi[j] = u; | |||
} | |||
} | |||
ret = intmtoratm_q(xmat,NM(q),*nmmat,dn); | |||
get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1); | |||
if ( ret ) { | |||
for ( j = k = l = 0; j < col; j++ ) | |||
if ( cinfo[j] > 0 ) | |||
rind[k++] = j; | |||
else if ( !cinfo[j] ) | |||
cind[l++] = j; | cind[l++] = j; | ||
get_eg(&tmp0); | get_eg(&tmp0); | ||
ret = gensolve_check(mat,*nmmat,*dn,rind,cind); | ret = gensolve_check(mat,*nmmat,*dn,rind,cind); |