version 1.42, 2004/12/13 23:04:16 |
version 1.45, 2005/01/23 14:03:47 |
|
|
* 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.41 2004/12/04 09:39:27 noro Exp $ |
* $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.44 2005/01/12 10:38:07 noro Exp $ |
*/ |
*/ |
#include "ca.h" |
#include "ca.h" |
#include "base.h" |
#include "base.h" |
Line 760 void Psize(NODE arg,LIST *rp) |
|
Line 760 void Psize(NODE arg,LIST *rp) |
|
n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col; |
n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col; |
STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s); |
STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s); |
break; |
break; |
|
case O_IMAT: |
|
n = ((IMAT)ARG0(arg))->row; m = ((IMAT)ARG0(arg))->col; |
|
STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s); |
|
break; |
default: |
default: |
error("size : invalid argument"); break; |
error("size : invalid argument"); break; |
} |
} |
Line 828 void Pinvmat(NODE arg,LIST *rp) |
|
Line 832 void Pinvmat(NODE arg,LIST *rp) |
|
input : a row x col matrix A |
input : a row x col matrix A |
A[I] <-> A[I][0]*x_0+A[I][1]*x_1+... |
A[I] <-> A[I][0]*x_0+A[I][1]*x_1+... |
|
|
output : [B,R,C] |
output : [B,D,R,C] |
B : a rank(A) x col-rank(A) matrix |
B : a rank(A) x col-rank(A) matrix |
|
D : the denominator |
R : a vector of length rank(A) |
R : a vector of length rank(A) |
C : a vector of length col-rank(A) |
C : a vector of length col-rank(A) |
B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+... |
B[I] <-> D*x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+... |
*/ |
*/ |
|
|
void Pgeneric_gauss_elim(NODE arg,LIST *rp) |
void Pgeneric_gauss_elim(NODE arg,LIST *rp) |
Line 1174 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn |
|
Line 1179 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn |
|
int ret; |
int ret; |
struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1; |
struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1; |
int period; |
int period; |
|
int *wx,*ptr; |
|
int wxsize,nsize; |
|
N wn; |
|
Q wq; |
|
|
a0 = (Q **)mat->body; |
a0 = (Q **)mat->body; |
row = mat->row; col = mat->col; |
row = mat->row; col = mat->col; |
Line 1223 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn |
|
Line 1232 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn |
|
init_eg(&eg_mul); init_eg(&eg_inv); |
init_eg(&eg_mul); init_eg(&eg_inv); |
init_eg(&eg_check); init_eg(&eg_intrat); |
init_eg(&eg_check); init_eg(&eg_intrat); |
period = F4_INTRAT_PERIOD; |
period = F4_INTRAT_PERIOD; |
for ( q = ONE, count = 0; ; count++ ) { |
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 > 3 ) |
if ( DP_Print > 3 ) |
fprintf(stderr,"o"); |
fprintf(stderr,"o"); |
/* wc = -b mod md */ |
/* wc = -b mod md */ |
|
get_eg(&tmp0); |
for ( i = 0; i < rank; i++ ) |
for ( i = 0; i < rank; i++ ) |
for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) |
for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) |
if ( u = (Q)bi[j] ) { |
if ( u = (Q)bi[j] ) { |
Line 1236 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn |
|
Line 1250 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn |
|
wi[j] = t; |
wi[j] = t; |
} else |
} else |
wi[j] = 0; |
wi[j] = 0; |
/* wc = A^(-1)wc; wc is normalized */ |
/* wc = A^(-1)wc; wc is not normalized */ |
get_eg(&tmp0); |
solve_by_lu_mod(w,rank,md,wc,ri,0); |
solve_by_lu_mod(w,rank,md,wc,ri); |
/* wx += q*wc */ |
get_eg(&tmp1); |
ptr = wx; |
add_eg(&eg_inv,&tmp0,&tmp1); |
|
/* x = x-q*wc */ |
|
for ( i = 0; i < rank; i++ ) |
for ( i = 0; i < rank; i++ ) |
for ( j = 0, xi = x[i], wi = wc[i]; j < ri; j++ ) { |
for ( j = 0, wi = wc[i]; j < ri; j++ ) { |
STOQ(wi[j],u); mulq(q,u,&s); |
if ( wi[j] ) |
subq(xi[j],s,&u); xi[j] = u; |
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); |
get_eg(&tmp0); |
for ( i = 0; i < rank; i++ ) |
for ( i = 0; i < rank; i++ ) |
for ( j = 0; j < ri; j++ ) { |
for ( j = 0; j < ri; j++ ) { |
Line 1264 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn |
|
Line 1280 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn |
|
add_eg(&eg_mul,&tmp0,&tmp1); |
add_eg(&eg_mul,&tmp0,&tmp1); |
/* q = q*md */ |
/* q = q*md */ |
mulq(q,mdq,&u); q = u; |
mulq(q,mdq,&u); q = u; |
if ( !(count % period) ) { |
if ( count == period ) { |
get_eg(&tmp0); |
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); |
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 ) { |
Line 1288 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn |
|
Line 1317 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn |
|
} |
} |
return rank; |
return rank; |
} |
} |
} else |
} else { |
period *=2; |
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; |
|
} |
} |
} |
} |
} |
} |
} |
Line 1890 int find_lhs_and_lu_mod(unsigned int **a,int row,int c |
|
Line 1925 int find_lhs_and_lu_mod(unsigned int **a,int row,int c |
|
b = a^(-1)b |
b = a^(-1)b |
*/ |
*/ |
|
|
void solve_by_lu_mod(int **a,int n,int md,int **b,int l) |
void solve_by_lu_mod(int **a,int n,int md,int **b,int l,int normalize) |
{ |
{ |
unsigned int *y,*c; |
unsigned int *y,*c; |
int i,j,k; |
int i,j,k; |
Line 1923 void solve_by_lu_mod(int **a,int n,int md,int **b,int |
|
Line 1958 void solve_by_lu_mod(int **a,int n,int md,int **b,int |
|
DMAR(t,a[i][i],0,md,c[i]) |
DMAR(t,a[i][i],0,md,c[i]) |
} |
} |
/* copy c to b[.][k] with normalization */ |
/* copy c to b[.][k] with normalization */ |
for ( i = 0; i < n; i++ ) |
if ( normalize ) |
b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]); |
for ( i = 0; i < n; i++ ) |
|
b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]); |
|
else |
|
for ( i = 0; i < n; i++ ) |
|
b[i][k] = c[i]; |
} |
} |
} |
} |
|
|