| version 1.3, 2000/04/20 02:20:15 |
version 1.75, 2017/09/17 02:34:02 |
|
|
| /* $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.2 2000/03/14 05:25:43 noro Exp $ */ |
/* |
| |
* Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED |
| |
* All rights reserved. |
| |
* |
| |
* FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited, |
| |
* non-exclusive and royalty-free license to use, copy, modify and |
| |
* redistribute, solely for non-commercial and non-profit purposes, the |
| |
* computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and |
| |
* conditions of this Agreement. For the avoidance of doubt, you acquire |
| |
* only a limited right to use the SOFTWARE hereunder, and FLL or any |
| |
* third party developer retains all rights, including but not limited to |
| |
* copyrights, in and to the SOFTWARE. |
| |
* |
| |
* (1) FLL does not grant you a license in any way for commercial |
| |
* purposes. You may use the SOFTWARE only for non-commercial and |
| |
* non-profit purposes only, such as academic, research and internal |
| |
* business use. |
| |
* (2) The SOFTWARE is protected by the Copyright Law of Japan and |
| |
* international copyright treaties. If you make copies of the SOFTWARE, |
| |
* with or without modification, as permitted hereunder, you shall affix |
| |
* to all such copies of the SOFTWARE the above copyright notice. |
| |
* (3) An explicit reference to this SOFTWARE and its copyright owner |
| |
* shall be made on your publication or presentation in any form of the |
| |
* results obtained by use of the SOFTWARE. |
| |
* (4) In the event that you modify the SOFTWARE, you shall notify FLL by |
| |
* e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification |
| |
* for such modification or the source code of the modified part of the |
| |
* SOFTWARE. |
| |
* |
| |
* THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL |
| |
* MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND |
| |
* EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS |
| |
* FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES' |
| |
* RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY |
| |
* MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY. |
| |
* UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT, |
| |
* OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY |
| |
* DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL |
| |
* DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES |
| |
* ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES |
| |
* FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY |
| |
* DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF |
| |
* SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART |
| |
* OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY |
| |
* DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, |
| |
* PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. |
| |
* |
| |
* $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.74 2017/09/15 01:52:51 noro Exp $ |
| |
*/ |
| #include "ca.h" |
#include "ca.h" |
| #include "base.h" |
#include "base.h" |
| #include "parse.h" |
#include "parse.h" |
| #include "inline.h" |
#include "inline.h" |
| /* |
|
| |
#include <sys/types.h> |
| |
#include <sys/stat.h> |
| |
#if !defined(_MSC_VER) |
| |
#include <unistd.h> |
| |
#endif |
| |
|
| |
#define F4_INTRAT_PERIOD 8 |
| |
|
| |
#if 0 |
| #undef DMAR |
#undef DMAR |
| #define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d); |
#define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d); |
| */ |
#endif |
| |
|
| extern int Print; /* XXX */ |
extern int DP_Print; /* XXX */ |
| |
|
| void inner_product_mat_int_mod(Q **,int **,int,int,int,Q *); |
|
| void solve_by_lu_mod(int **,int,int,int **,int); |
|
| void solve_by_lu_gfmmat(GFMMAT,unsigned int,unsigned int *,unsigned int *); |
|
| int lu_gfmmat(GFMMAT,unsigned int,int *); |
|
| void mat_to_gfmmat(MAT,unsigned int,GFMMAT *); |
|
| |
|
| int generic_gauss_elim_mod(int **,int,int,int,int *); |
void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(), Ptriangleq(); |
| int generic_gauss_elim(MAT ,MAT *,Q *,int **,int **); |
void Pinvmat(); |
| |
void Pnewbytearray(),Pmemoryplot_to_coord(); |
| |
|
| int gauss_elim_mod(int **,int,int,int); |
void Pgeneric_gauss_elim(); |
| int gauss_elim_mod1(int **,int,int,int); |
|
| int gauss_elim_geninv_mod(unsigned int **,int,int,int); |
|
| int gauss_elim_geninv_mod_swap(unsigned int **,int,int,unsigned int,unsigned int ***,int **); |
|
| void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(); |
|
| |
|
| void Pgeneric_gauss_elim_mod(); |
void Pgeneric_gauss_elim_mod(); |
| |
|
| |
void Pindep_rows_mod(); |
| |
|
| void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat(); |
void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat(); |
| void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol(); |
void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol(), Pltov(); |
| |
void Pgeninv_sf_swap(); |
| void sepvect(); |
void sepvect(); |
| void Pmulmat_gf2n(); |
void Pmulmat_gf2n(); |
| void Pbconvmat_gf2n(); |
void Pbconvmat_gf2n(); |
| Line 40 void Px962_irredpoly_up2(); |
|
| Line 91 void Px962_irredpoly_up2(); |
|
| void Pirredpoly_up2(); |
void Pirredpoly_up2(); |
| void Pnbpoly_up2(); |
void Pnbpoly_up2(); |
| void Pqsort(); |
void Pqsort(); |
| |
void Pexponent_vector(); |
| |
void Pmat_swap_row_destructive(); |
| |
void Pmat_swap_col_destructive(); |
| |
void Pvect(); |
| |
void Pmat(); |
| |
void Pmatc(); |
| |
void Pnd_det(); |
| |
void Plu_mat(); |
| |
void Pmat_col(); |
| |
void Plusolve_prep(); |
| |
void Plusolve_main(); |
| |
|
| struct ftab array_tab[] = { |
struct ftab array_tab[] = { |
| |
{"lu_mat",Plu_mat,1}, |
| {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4}, |
{"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4}, |
| {"lu_gfmmat",Plu_gfmmat,2}, |
{"lu_gfmmat",Plu_gfmmat,2}, |
| {"mat_to_gfmmat",Pmat_to_gfmmat,2}, |
{"mat_to_gfmmat",Pmat_to_gfmmat,2}, |
| |
{"generic_gauss_elim",Pgeneric_gauss_elim,1}, |
| {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2}, |
{"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2}, |
| |
{"indep_rows_mod",Pindep_rows_mod,2}, |
| {"newvect",Pnewvect,-2}, |
{"newvect",Pnewvect,-2}, |
| |
{"vect",Pvect,-99999999}, |
| |
{"vector",Pnewvect,-2}, |
| |
{"exponent_vector",Pexponent_vector,-99999999}, |
| {"newmat",Pnewmat,-3}, |
{"newmat",Pnewmat,-3}, |
| |
{"matrix",Pnewmat,-3}, |
| |
{"mat",Pmat,-99999999}, |
| |
{"matr",Pmat,-99999999}, |
| |
{"matc",Pmatc,-99999999}, |
| |
{"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}, |
| {"vtol",Pvtol,1}, |
{"vtol",Pvtol,1}, |
| |
{"ltov",Pltov,1}, |
| {"size",Psize,1}, |
{"size",Psize,1}, |
| {"det",Pdet,-2}, |
{"det",Pdet,-2}, |
| |
{"nd_det",Pnd_det,-2}, |
| |
{"invmat",Pinvmat,-2}, |
| {"leqm",Pleqm,2}, |
{"leqm",Pleqm,2}, |
| {"leqm1",Pleqm1,2}, |
{"leqm1",Pleqm1,2}, |
| {"geninvm",Pgeninvm,2}, |
{"geninvm",Pgeninvm,2}, |
| {"geninvm_swap",Pgeninvm_swap,2}, |
{"geninvm_swap",Pgeninvm_swap,2}, |
| |
{"geninv_sf_swap",Pgeninv_sf_swap,1}, |
| {"remainder",Premainder,2}, |
{"remainder",Premainder,2}, |
| {"sremainder",Psremainder,2}, |
{"sremainder",Psremainder,2}, |
| {"mulmat_gf2n",Pmulmat_gf2n,1}, |
{"mulmat_gf2n",Pmulmat_gf2n,1}, |
| Line 68 struct ftab array_tab[] = { |
|
| Line 146 struct ftab array_tab[] = { |
|
| {"x962_irredpoly_up2",Px962_irredpoly_up2,2}, |
{"x962_irredpoly_up2",Px962_irredpoly_up2,2}, |
| {"irredpoly_up2",Pirredpoly_up2,2}, |
{"irredpoly_up2",Pirredpoly_up2,2}, |
| {"nbpoly_up2",Pnbpoly_up2,2}, |
{"nbpoly_up2",Pnbpoly_up2,2}, |
| |
{"mat_swap_row_destructive",Pmat_swap_row_destructive,3}, |
| |
{"mat_swap_col_destructive",Pmat_swap_col_destructive,3}, |
| |
{"mat_col",Pmat_col,2}, |
| |
{"lusolve_prep",Plusolve_prep,1}, |
| |
{"lusolve_main",Plusolve_main,1}, |
| |
{"triangleq",Ptriangleq,1}, |
| {0,0,0}, |
{0,0,0}, |
| }; |
}; |
| |
|
| int comp_obj(a,b) |
typedef struct _ent { int j; unsigned int e; } ent; |
| Obj *a,*b; |
|
| |
ent *get_row(FILE *,int *l); |
| |
void put_row(FILE *out,int l,ent *a); |
| |
void lu_elim(int *l,ent **a,int k,int i,int mul,int mod); |
| |
void lu_append(int *,ent **,int *,int,int,int); |
| |
void solve_l(int *,ent **,int,int *,int); |
| |
void solve_u(int *,ent **,int,int *,int); |
| |
|
| |
|
| |
static int *ul,*ll; |
| |
static ent **u,**l; |
| |
static int modulus; |
| |
|
| |
void Plusolve_prep(NODE arg,Q *rp) |
| { |
{ |
| |
char *fname; |
| |
FILE *in; |
| |
int len,i,rank; |
| |
int *rhs; |
| |
|
| |
fname = BDY((STRING)ARG0(arg)); |
| |
in = fopen(fname,"r"); |
| |
modulus = getw(in); |
| |
len = getw(in); |
| |
ul = (int *)MALLOC_ATOMIC(len*sizeof(int)); |
| |
u = (ent **)MALLOC(len*sizeof(ent *)); |
| |
ll = (int *)MALLOC_ATOMIC(len*sizeof(int)); |
| |
l = (ent **)MALLOC(len*sizeof(ent *)); |
| |
for ( i = 0; i < len; i++ ) { |
| |
u[i] = get_row(in,&ul[i]); |
| |
} |
| |
for ( i = 0; i < len; i++ ) { |
| |
l[i] = get_row(in,&ll[i]); |
| |
} |
| |
fclose(in); |
| |
*rp = ONE; |
| |
} |
| |
|
| |
void Plusolve_main(NODE arg,VECT *rp) |
| |
{ |
| |
Q *d,*p; |
| |
VECT v,r; |
| |
int len,i; |
| |
int *rhs; |
| |
|
| |
v = (VECT)ARG0(arg); len = v->len; |
| |
d = (Q *)BDY(v); |
| |
rhs = (int *)MALLOC_ATOMIC(len*sizeof(int)); |
| |
for ( i = 0; i < len; i++ ) rhs[i] = QTOS(d[i]); |
| |
solve_l(ll,l,len,rhs,modulus); |
| |
solve_u(ul,u,len,rhs,modulus); |
| |
NEWVECT(r); r->len = len; |
| |
r->body = (pointer *)MALLOC(len*sizeof(pointer)); |
| |
p = (Q *)r->body; |
| |
for ( i = 0; i < len; i++ ) |
| |
STOQ(rhs[i],p[i]); |
| |
*rp = r; |
| |
} |
| |
|
| |
ent *get_row(FILE *in,int *l) |
| |
{ |
| |
int len,i; |
| |
ent *a; |
| |
|
| |
*l = len = getw(in); |
| |
a = (ent *)MALLOC_ATOMIC(len*sizeof(ent)); |
| |
for ( i = 0; i < len; i++ ) { |
| |
a[i].j = getw(in); |
| |
a[i].e = getw(in); |
| |
} |
| |
return a; |
| |
} |
| |
|
| |
void lu_gauss(int *ul,ent **u,int *ll,ent **l,int n,int mod) |
| |
{ |
| |
int i,j,k,s,mul; |
| |
unsigned int inv; |
| |
int *ll2; |
| |
|
| |
ll2 = (int *)MALLOC_ATOMIC(n*sizeof(int)); |
| |
for ( i = 0; i < n; i++ ) ll2[i] = 0; |
| |
for ( i = 0; i < n; i++ ) { |
| |
fprintf(stderr,"i=%d\n",i); |
| |
inv = invm(u[i][0].e,mod); |
| |
for ( k = i+1; k < n; k++ ) |
| |
if ( u[k][0].j == n-i ) { |
| |
s = u[k][0].e; |
| |
DMAR(s,inv,0,mod,mul); |
| |
lu_elim(ul,u,k,i,mul,mod); |
| |
lu_append(ll,l,ll2,k,i,mul); |
| |
} |
| |
} |
| |
} |
| |
|
| |
#define INITLEN 10 |
| |
|
| |
void lu_append(int *l,ent **a,int *l2,int k,int i,int mul) |
| |
{ |
| |
int len; |
| |
ent *p; |
| |
|
| |
len = l[k]; |
| |
if ( !len ) { |
| |
a[k] = p = (ent *)MALLOC_ATOMIC(INITLEN*sizeof(ent)); |
| |
p[0].j = i; p[0].e = mul; |
| |
l[k] = 1; l2[k] = INITLEN; |
| |
} else { |
| |
if ( l2[k] == l[k] ) { |
| |
l2[k] *= 2; |
| |
a[k] = REALLOC(a[k],l2[k]*sizeof(ent)); |
| |
} |
| |
p =a[k]; |
| |
p[l[k]].j = i; p[l[k]].e = mul; |
| |
l[k]++; |
| |
} |
| |
} |
| |
|
| |
/* a[k] = a[k]-mul*a[i] */ |
| |
|
| |
void lu_elim(int *l,ent **a,int k,int i,int mul,int mod) |
| |
{ |
| |
ent *ak,*ai,*w; |
| |
int lk,li,j,m,p,q,r,s,t,j0; |
| |
|
| |
ak = a[k]; ai = a[i]; lk = l[k]; li = l[i]; |
| |
w = (ent *)alloca((lk+li)*sizeof(ent)); |
| |
p = 0; q = 0; j = 0; |
| |
mul = mod-mul; |
| |
while ( p < lk && q < li ) { |
| |
if ( ak[p].j > ai[q].j ) { |
| |
w[j] = ak[p]; j++; p++; |
| |
} else if ( ak[p].j < ai[q].j ) { |
| |
w[j].j = ai[q].j; |
| |
t = ai[q].e; |
| |
DMAR(t,mul,0,mod,r); |
| |
w[j].e = r; |
| |
j++; q++; |
| |
} else { |
| |
t = ai[q].e; s = ak[p].e; |
| |
DMAR(t,mul,s,mod,r); |
| |
if ( r ) { |
| |
w[j].j = ai[q].j; w[j].e = r; j++; |
| |
} |
| |
p++; q++; |
| |
} |
| |
} |
| |
if ( q == li ) |
| |
while ( p < lk ) { |
| |
w[j] = ak[p]; j++; p++; |
| |
} |
| |
else if ( p == lk ) |
| |
while ( q < li ) { |
| |
w[j].j = ai[q].j; |
| |
t = ai[q].e; |
| |
DMAR(t,mul,0,mod,r); |
| |
w[j].e = r; |
| |
j++; q++; |
| |
} |
| |
if ( j <= lk ) { |
| |
for ( m = 0; m < j; m++ ) ak[m] = w[m]; |
| |
} else { |
| |
a[k] = ak = (ent *)MALLOC_ATOMIC(j*sizeof(ent)); |
| |
for ( m = 0; m < j; m++ ) ak[m] = w[m]; |
| |
} |
| |
l[k] = j; |
| |
} |
| |
|
| |
void solve_l(int *ll,ent **l,int n,int *rhs,int mod) |
| |
{ |
| |
int j,k,s,len; |
| |
ent *p; |
| |
|
| |
for ( j = 0; j < n; j++ ) { |
| |
len = ll[j]; p = l[j]; |
| |
for ( k = 0, s = 0; k < len; k++ ) |
| |
s = dmar(p[k].e,rhs[p[k].j],s,mod); |
| |
rhs[j] -= s; |
| |
if ( rhs[j] < 0 ) rhs[j] += mod; |
| |
} |
| |
} |
| |
|
| |
void solve_u(int *ul,ent **u,int n,int *rhs,int mod) |
| |
{ |
| |
int j,k,s,len,inv; |
| |
ent *p; |
| |
|
| |
for ( j = n-1; j >= 0; j-- ) { |
| |
len = ul[j]; p = u[j]; |
| |
for ( k = 1, s = 0; k < len; k++ ) |
| |
s = dmar(p[k].e,rhs[p[k].j],s,mod); |
| |
rhs[j] -= s; |
| |
if ( rhs[j] < 0 ) rhs[j] += mod; |
| |
inv = invm((unsigned int)p[0].e,mod); |
| |
rhs[j] = dmar(rhs[j],inv,0,mod); |
| |
} |
| |
} |
| |
|
| |
int comp_obj(Obj *a,Obj *b) |
| |
{ |
| return arf_comp(CO,*a,*b); |
return arf_comp(CO,*a,*b); |
| } |
} |
| |
|
| static FUNC generic_comp_obj_func; |
static FUNC generic_comp_obj_func; |
| static NODE generic_comp_obj_arg; |
static NODE generic_comp_obj_arg; |
| |
static NODE generic_comp_obj_option; |
| |
|
| int generic_comp_obj(a,b) |
int generic_comp_obj(Obj *a,Obj *b) |
| Obj *a,*b; |
|
| { |
{ |
| Q r; |
Q r; |
| |
|
| BDY(generic_comp_obj_arg)=(pointer)(*a); |
BDY(generic_comp_obj_arg)=(pointer)(*a); |
| BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b); |
BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b); |
| r = (Q)bevalf(generic_comp_obj_func,generic_comp_obj_arg); |
r = (Q)bevalf_with_opts(generic_comp_obj_func,generic_comp_obj_arg,generic_comp_obj_option); |
| if ( !r ) |
if ( !r ) |
| return 0; |
return 0; |
| else |
else |
|
|
| } |
} |
| |
|
| |
|
| void Pqsort(arg,rp) |
void Pqsort(NODE arg,LIST *rp) |
| NODE arg; |
|
| VECT *rp; |
|
| { |
{ |
| VECT vect; |
VECT vect; |
| char buf[BUFSIZ]; |
NODE n,n1; |
| char *fname; |
|
| NODE n; |
|
| P p; |
P p; |
| V v; |
V v; |
| |
FUNC func; |
| |
int len,i; |
| |
pointer *a; |
| |
Obj t; |
| |
|
| asir_assert(ARG0(arg),O_VECT,"qsort"); |
t = ARG0(arg); |
| vect = (VECT)ARG0(arg); |
if (OID(t) == O_LIST) { |
| |
n = (NODE)BDY((LIST)t); |
| |
len = length(n); |
| |
MKVECT(vect,len); |
| |
for ( i = 0; i < len; i++, n = NEXT(n) ) { |
| |
BDY(vect)[i] = BDY(n); |
| |
} |
| |
|
| |
}else if (OID(t) != O_VECT) { |
| |
error("qsort : invalid argument"); |
| |
}else { |
| |
vect = (VECT)t; |
| |
} |
| if ( argc(arg) == 1 ) |
if ( argc(arg) == 1 ) |
| qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj); |
qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj); |
| else { |
else { |
|
|
| if ( !p || OID(p)!=2 ) |
if ( !p || OID(p)!=2 ) |
| error("qsort : invalid argument"); |
error("qsort : invalid argument"); |
| v = VR(p); |
v = VR(p); |
| if ( (int)v->attr != V_SR ) |
gen_searchf(NAME(v),&func); |
| error("qsort : no such function"); |
if ( !func ) { |
| generic_comp_obj_func = (FUNC)v->priv; |
if ( (int)v->attr != V_SR ) |
| MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n); |
error("qsort : no such function"); |
| |
func = (FUNC)v->priv; |
| |
} |
| |
generic_comp_obj_func = func; |
| |
MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n); |
| |
generic_comp_obj_option = current_option; |
| qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj); |
qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj); |
| } |
} |
| *rp = vect; |
if (OID(t) == O_LIST) { |
| |
a = BDY(vect); |
| |
for ( i = len - 1, n = 0; i >= 0; i-- ) { |
| |
MKNODE(n1,a[i],n); n = n1; |
| |
} |
| |
MKLIST(*rp,n); |
| |
}else { |
| |
*rp = (LIST)vect; |
| |
} |
| } |
} |
| |
|
| void PNBmul_gf2n(arg,rp) |
void PNBmul_gf2n(NODE arg,GF2N *rp) |
| NODE arg; |
|
| GF2N *rp; |
|
| { |
{ |
| GF2N a,b; |
GF2N a,b; |
| GF2MAT mat; |
GF2MAT mat; |
|
|
| } |
} |
| } |
} |
| |
|
| void Pmul_vect_mat_gf2n(arg,rp) |
void Pmul_vect_mat_gf2n(NODE arg,GF2N *rp) |
| NODE arg; |
|
| GF2N *rp; |
|
| { |
{ |
| GF2N a; |
GF2N a; |
| GF2MAT mat; |
GF2MAT mat; |
|
|
| } |
} |
| } |
} |
| |
|
| void Pbconvmat_gf2n(arg,rp) |
void Pbconvmat_gf2n(NODE arg,LIST *rp) |
| NODE arg; |
|
| LIST *rp; |
|
| { |
{ |
| P p0,p1; |
P p0,p1; |
| int to; |
int to; |
|
|
| MKLIST(*rp,n0); |
MKLIST(*rp,n0); |
| } |
} |
| |
|
| void Pmulmat_gf2n(arg,rp) |
void Pmulmat_gf2n(NODE arg,GF2MAT *rp) |
| NODE arg; |
|
| GF2MAT *rp; |
|
| { |
{ |
| GF2MAT m; |
GF2MAT m; |
| |
|
|
|
| *rp = m; |
*rp = m; |
| } |
} |
| |
|
| void Psepmat_destructive(arg,rp) |
void Psepmat_destructive(NODE arg,LIST *rp) |
| NODE arg; |
|
| LIST *rp; |
|
| { |
{ |
| MAT mat,mat1; |
MAT mat,mat1; |
| int i,j,row,col; |
int i,j,row,col; |
|
|
| sgn = SGN(ent); |
sgn = SGN(ent); |
| divn(nm,mod,&quo,&rem); |
divn(nm,mod,&quo,&rem); |
| /* if ( quo != nm && rem != nm ) */ |
/* if ( quo != nm && rem != nm ) */ |
| /* GC_free(nm); */ |
/* GCFREE(nm); */ |
| /* GC_free(ent); */ |
/* GCFREE(ent); */ |
| NTOQ(rem,sgn,a[i][j]); NTOQ(quo,sgn,a1[i][j]); |
NTOQ(rem,sgn,a[i][j]); NTOQ(quo,sgn,a1[i][j]); |
| } |
} |
| MKNODE(n1,mat1,0); MKNODE(n0,mat,n1); |
MKNODE(n1,mat1,0); MKNODE(n0,mat,n1); |
| MKLIST(*rp,n0); |
MKLIST(*rp,n0); |
| } |
} |
| |
|
| void Psepvect(arg,rp) |
void Psepvect(NODE arg,VECT *rp) |
| NODE arg; |
|
| VECT *rp; |
|
| { |
{ |
| sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp); |
sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp); |
| } |
} |
| |
|
| void sepvect(v,d,rp) |
void sepvect(VECT v,int d,VECT *rp) |
| VECT v; |
|
| int d; |
|
| VECT *rp; |
|
| { |
{ |
| int i,j,k,n,q,q1,r; |
int i,j,k,n,q,q1,r; |
| pointer *pv,*pw,*pu; |
pointer *pv,*pw,*pu; |
|
|
| } |
} |
| } |
} |
| |
|
| void Pnewvect(arg,rp) |
void Pnewvect(NODE arg,VECT *rp) |
| NODE arg; |
|
| VECT *rp; |
|
| { |
{ |
| int len,i,r; |
int len,i,r; |
| VECT vect; |
VECT vect; |
|
|
| |
|
| asir_assert(ARG0(arg),O_N,"newvect"); |
asir_assert(ARG0(arg),O_N,"newvect"); |
| len = QTOS((Q)ARG0(arg)); |
len = QTOS((Q)ARG0(arg)); |
| if ( len <= 0 ) |
if ( len < 0 ) |
| error("newvect : invalid size"); |
error("newvect : invalid size"); |
| MKVECT(vect,len); |
MKVECT(vect,len); |
| if ( argc(arg) == 2 ) { |
if ( argc(arg) == 2 ) { |
| list = (LIST)ARG1(arg); |
list = (LIST)ARG1(arg); |
| asir_assert(list,O_LIST,"newvect"); |
asir_assert(list,O_LIST,"newvect"); |
| |
#if 0 |
| for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ); |
for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ); |
| if ( r > len ) { |
if ( r > len ) { |
| *rp = vect; |
*rp = vect; |
| return; |
return; |
| } |
} |
| |
#endif |
| for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) ) |
for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) ) |
| vb[i] = (pointer)BDY(tn); |
vb[i] = (pointer)BDY(tn); |
| } |
} |
| *rp = vect; |
*rp = vect; |
| } |
} |
| |
|
| void Pnewmat(arg,rp) |
void Pvect(NODE arg,VECT *rp) { |
| NODE arg; |
int len,i; |
| MAT *rp; |
VECT vect; |
| |
pointer *vb; |
| |
NODE tn; |
| |
|
| |
if ( !arg ) { |
| |
*rp =0; |
| |
return; |
| |
} |
| |
|
| |
for (len = 0, tn = arg; tn; tn = NEXT(tn), len++); |
| |
if ( len == 1 ) { |
| |
if ( ARG0(arg) != 0 ) { |
| |
switch ( OID(ARG0(arg)) ) { |
| |
case O_VECT: |
| |
*rp = ARG0(arg); |
| |
return; |
| |
case O_LIST: |
| |
for ( len = 0, tn = ARG0(arg); tn; tn = NEXT(tn), len++ ); |
| |
MKVECT(vect,len-1); |
| |
for ( i = 0, tn = BDY((LIST)ARG0(arg)), vb =BDY(vect); |
| |
tn; i++, tn = NEXT(tn) ) |
| |
vb[i] = (pointer)BDY(tn); |
| |
*rp=vect; |
| |
return; |
| |
} |
| |
} |
| |
} |
| |
MKVECT(vect,len); |
| |
for ( i = 0, tn = arg, vb = BDY(vect); tn; i++, tn = NEXT(tn) ) |
| |
vb[i] = (pointer)BDY(tn); |
| |
*rp = vect; |
| |
} |
| |
|
| |
void Pexponent_vector(NODE arg,DP *rp) |
| { |
{ |
| |
nodetod(arg,rp); |
| |
} |
| |
|
| |
void Pnewbytearray(NODE arg,BYTEARRAY *rp) |
| |
{ |
| |
int len,i,r; |
| |
BYTEARRAY array; |
| |
unsigned char *vb; |
| |
char *str; |
| |
LIST list; |
| |
NODE tn; |
| |
int ac; |
| |
struct stat sbuf; |
| |
char *fname; |
| |
FILE *fp; |
| |
|
| |
ac = argc(arg); |
| |
if ( ac == 1 ) { |
| |
if ( !OID((Obj)ARG0(arg)) ) error("newbytearray : invalid argument"); |
| |
switch ( OID((Obj)ARG0(arg)) ) { |
| |
case O_STR: |
| |
fname = BDY((STRING)ARG0(arg)); |
| |
fp = fopen(fname,"rb"); |
| |
if ( !fp ) error("newbytearray : fopen failed"); |
| |
if ( stat(fname,&sbuf) < 0 ) |
| |
error("newbytearray : stat failed"); |
| |
len = sbuf.st_size; |
| |
MKBYTEARRAY(array,len); |
| |
fread(BDY(array),len,sizeof(char),fp); |
| |
break; |
| |
case O_N: |
| |
if ( !RATN(ARG0(arg)) ) |
| |
error("newbytearray : invalid argument"); |
| |
len = QTOS((Q)ARG0(arg)); |
| |
if ( len < 0 ) |
| |
error("newbytearray : invalid size"); |
| |
MKBYTEARRAY(array,len); |
| |
break; |
| |
default: |
| |
error("newbytearray : invalid argument"); |
| |
} |
| |
} else if ( ac == 2 ) { |
| |
asir_assert(ARG0(arg),O_N,"newbytearray"); |
| |
len = QTOS((Q)ARG0(arg)); |
| |
if ( len < 0 ) |
| |
error("newbytearray : invalid size"); |
| |
MKBYTEARRAY(array,len); |
| |
if ( !ARG1(arg) ) |
| |
error("newbytearray : invalid initialization"); |
| |
switch ( OID((Obj)ARG1(arg)) ) { |
| |
case O_LIST: |
| |
list = (LIST)ARG1(arg); |
| |
asir_assert(list,O_LIST,"newbytearray"); |
| |
for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ); |
| |
if ( r <= len ) { |
| |
for ( i = 0, tn = BDY(list), vb = BDY(array); tn; |
| |
i++, tn = NEXT(tn) ) |
| |
vb[i] = (unsigned char)QTOS((Q)BDY(tn)); |
| |
} |
| |
break; |
| |
case O_STR: |
| |
str = BDY((STRING)ARG1(arg)); |
| |
r = strlen(str); |
| |
if ( r <= len ) |
| |
bcopy(str,BDY(array),r); |
| |
break; |
| |
default: |
| |
if ( !ARG1(arg) ) |
| |
error("newbytearray : invalid initialization"); |
| |
} |
| |
} else |
| |
error("newbytearray : invalid argument"); |
| |
*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; |
| |
unsigned 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) |
| |
{ |
| int row,col; |
int row,col; |
| int i,j,r,c; |
int i,j,r,c; |
| NODE tn,sn; |
NODE tn,sn; |
|
|
| asir_assert(ARG0(arg),O_N,"newmat"); |
asir_assert(ARG0(arg),O_N,"newmat"); |
| asir_assert(ARG1(arg),O_N,"newmat"); |
asir_assert(ARG1(arg),O_N,"newmat"); |
| row = QTOS((Q)ARG0(arg)); col = QTOS((Q)ARG1(arg)); |
row = QTOS((Q)ARG0(arg)); col = QTOS((Q)ARG1(arg)); |
| if ( row <= 0 || col <= 0 ) |
if ( row < 0 || col < 0 ) |
| error("newmat : invalid size"); |
error("newmat : invalid size"); |
| MKMAT(m,row,col); |
MKMAT(m,row,col); |
| if ( argc(arg) == 3 ) { |
if ( argc(arg) == 3 ) { |
|
|
| *rp = m; |
*rp = m; |
| } |
} |
| |
|
| void Pvtol(arg,rp) |
void Pmat(NODE arg, MAT *rp) |
| NODE arg; |
|
| LIST *rp; |
|
| { |
{ |
| |
int row,col; |
| |
int i; |
| |
MAT m; |
| |
pointer **mb; |
| |
pointer *ent; |
| |
NODE tn, sn; |
| |
VECT v; |
| |
|
| |
if ( !arg ) { |
| |
*rp =0; |
| |
return; |
| |
} |
| |
|
| |
for (row = 0, tn = arg; tn; tn = NEXT(tn), row++); |
| |
if ( row == 1 ) { |
| |
if ( OID(ARG0(arg)) == O_MAT ) { |
| |
*rp=ARG0(arg); |
| |
return; |
| |
} else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) { |
| |
error("mat : invalid argument"); |
| |
} |
| |
} |
| |
if ( OID(ARG0(arg)) == O_VECT ) { |
| |
v = ARG0(arg); |
| |
col = v->len; |
| |
} else if ( OID(ARG0(arg)) == O_LIST ) { |
| |
for (col = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), col++); |
| |
} else { |
| |
error("mat : invalid argument"); |
| |
} |
| |
|
| |
MKMAT(m,row,col); |
| |
for (row = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), row++) { |
| |
if ( BDY(tn) == 0 ) { |
| |
error("mat : invalid argument"); |
| |
} else if ( OID(BDY(tn)) == O_VECT ) { |
| |
v = tn->body; |
| |
ent = BDY(v); |
| |
for (i = 0; i < v->len; i++ ) mb[row][i] = (Obj)ent[i]; |
| |
} else if ( OID(BDY(tn)) == O_LIST ) { |
| |
for (col = 0, sn = BDY((LIST)BDY(tn)); sn; col++, sn = NEXT(sn) ) |
| |
mb[row][col] = (pointer)BDY(sn); |
| |
} else { |
| |
error("mat : invalid argument"); |
| |
} |
| |
} |
| |
*rp = m; |
| |
} |
| |
|
| |
void Pmatc(NODE arg, MAT *rp) |
| |
{ |
| |
int row,col; |
| |
int i; |
| |
MAT m; |
| |
pointer **mb; |
| |
pointer *ent; |
| |
NODE tn, sn; |
| |
VECT v; |
| |
|
| |
if ( !arg ) { |
| |
*rp =0; |
| |
return; |
| |
} |
| |
|
| |
for (col = 0, tn = arg; tn; tn = NEXT(tn), col++); |
| |
if ( col == 1 ) { |
| |
if ( OID(ARG0(arg)) == O_MAT ) { |
| |
*rp=ARG0(arg); |
| |
return; |
| |
} else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) { |
| |
error("matc : invalid argument"); |
| |
} |
| |
} |
| |
if ( OID(ARG0(arg)) == O_VECT ) { |
| |
v = ARG0(arg); |
| |
row = v->len; |
| |
} else if ( OID(ARG0(arg)) == O_LIST ) { |
| |
for (row = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), row++); |
| |
} else { |
| |
error("matc : invalid argument"); |
| |
} |
| |
|
| |
MKMAT(m,row,col); |
| |
for (col = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), col++) { |
| |
if ( BDY(tn) == 0 ) { |
| |
error("matc : invalid argument"); |
| |
} else if ( OID(BDY(tn)) == O_VECT ) { |
| |
v = tn->body; |
| |
ent = BDY(v); |
| |
for (i = 0; i < v->len; i++ ) mb[i][col] = (Obj)ent[i]; |
| |
} else if ( OID(BDY(tn)) == O_LIST ) { |
| |
for (row = 0, sn = BDY((LIST)BDY(tn)); sn; row++, sn = NEXT(sn) ) |
| |
mb[row][col] = (pointer)BDY(sn); |
| |
} else { |
| |
error("matc : invalid argument"); |
| |
} |
| |
} |
| |
*rp = m; |
| |
} |
| |
|
| |
void Pvtol(NODE arg,LIST *rp) |
| |
{ |
| NODE n,n1; |
NODE n,n1; |
| VECT v; |
VECT v; |
| pointer *a; |
pointer *a; |
| int len,i; |
int len,i; |
| |
|
| |
if ( OID(ARG0(arg)) == O_LIST ) { |
| |
*rp = ARG0(arg); |
| |
return; |
| |
} |
| asir_assert(ARG0(arg),O_VECT,"vtol"); |
asir_assert(ARG0(arg),O_VECT,"vtol"); |
| v = (VECT)ARG0(arg); len = v->len; a = BDY(v); |
v = (VECT)ARG0(arg); len = v->len; a = BDY(v); |
| for ( i = len - 1, n = 0; i >= 0; i-- ) { |
for ( i = len - 1, n = 0; i >= 0; i-- ) { |
|
|
| MKLIST(*rp,n); |
MKLIST(*rp,n); |
| } |
} |
| |
|
| void Premainder(arg,rp) |
void Pltov(NODE arg,VECT *rp) |
| NODE arg; |
|
| Obj *rp; |
|
| { |
{ |
| |
NODE n; |
| |
VECT v,v0; |
| |
int len,i; |
| |
|
| |
if ( OID(ARG0(arg)) == O_VECT ) { |
| |
v0 = (VECT)ARG0(arg); len = v0->len; |
| |
MKVECT(v,len); |
| |
for ( i = 0; i < len; i++ ) { |
| |
BDY(v)[i] = BDY(v0)[i]; |
| |
} |
| |
*rp = v; |
| |
return; |
| |
} |
| |
asir_assert(ARG0(arg),O_LIST,"ltov"); |
| |
n = (NODE)BDY((LIST)ARG0(arg)); |
| |
len = length(n); |
| |
MKVECT(v,len); |
| |
for ( i = 0; i < len; i++, n = NEXT(n) ) |
| |
BDY(v)[i] = BDY(n); |
| |
*rp = v; |
| |
} |
| |
|
| |
void Premainder(NODE arg,Obj *rp) |
| |
{ |
| Obj a; |
Obj a; |
| VECT v,w; |
VECT v,w; |
| MAT m,l; |
MAT m,l; |
|
|
| } |
} |
| } |
} |
| |
|
| void Psremainder(arg,rp) |
void Psremainder(NODE arg,Obj *rp) |
| NODE arg; |
|
| Obj *rp; |
|
| { |
{ |
| Obj a; |
Obj a; |
| VECT v,w; |
VECT v,w; |
|
|
| } |
} |
| } |
} |
| |
|
| void Psize(arg,rp) |
void Psize(NODE arg,LIST *rp) |
| NODE arg; |
|
| LIST *rp; |
|
| { |
{ |
| |
|
| int n,m; |
int n,m; |
|
|
| 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; |
| } |
} |
|
|
| MKLIST(*rp,t); |
MKLIST(*rp,t); |
| } |
} |
| |
|
| void Pdet(arg,rp) |
void Pdet(NODE arg,P *rp) |
| NODE arg; |
|
| P *rp; |
|
| { |
{ |
| MAT m; |
MAT m; |
| int n,i,j,mod; |
int n,i,j,mod; |
|
|
| } |
} |
| } |
} |
| |
|
| |
void Pinvmat(NODE arg,LIST *rp) |
| |
{ |
| |
MAT m,r; |
| |
int n,i,j,mod; |
| |
P dn; |
| |
P **mat,**imat,**w; |
| |
NODE nd; |
| |
|
| |
m = (MAT)ARG0(arg); |
| |
asir_assert(m,O_MAT,"invmat"); |
| |
if ( m->row != m->col ) |
| |
error("invmat : non-square matrix"); |
| |
else if ( argc(arg) == 1 ) { |
| |
n = m->row; |
| |
invmatp(CO,(P **)BDY(m),n,&imat,&dn); |
| |
NEWMAT(r); r->row = n; r->col = n; r->body = (pointer **)imat; |
| |
nd = mknode(2,r,dn); |
| |
MKLIST(*rp,nd); |
| |
} else { |
| |
n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m); |
| |
w = (P **)almat_pointer(n,n); |
| |
for ( i = 0; i < n; i++ ) |
| |
for ( j = 0; j < n; j++ ) |
| |
ptomp(mod,mat[i][j],&w[i][j]); |
| |
#if 0 |
| |
detmp(CO,mod,w,n,&d); |
| |
mptop(d,rp); |
| |
#else |
| |
error("not implemented yet"); |
| |
#endif |
| |
} |
| |
} |
| |
|
| /* |
/* |
| 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,D,R,C] |
| |
B : a rank(A) x col-rank(A) matrix |
| |
D : the denominator |
| |
R : a vector of length rank(A) |
| |
C : a vector of length col-rank(A) |
| |
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) |
| |
{ |
| |
NODE n0,opt,p; |
| |
MAT m,nm; |
| |
int *ri,*ci; |
| |
VECT rind,cind; |
| |
Q dn,q; |
| |
int i,row,col,t,rank; |
| |
int is_hensel = 0; |
| |
char *key; |
| |
Obj value; |
| |
|
| |
if ( current_option ) { |
| |
for ( opt = current_option; opt; opt = NEXT(opt) ) { |
| |
p = BDY((LIST)BDY(opt)); |
| |
key = BDY((STRING)BDY(p)); |
| |
value = (Obj)BDY(NEXT(p)); |
| |
if ( !strcmp(key,"hensel") && value ) { |
| |
is_hensel = value ? 1 : 0; |
| |
break; |
| |
} |
| |
} |
| |
} |
| |
asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim"); |
| |
m = (MAT)ARG0(arg); |
| |
row = m->row; col = m->col; |
| |
if ( is_hensel ) |
| |
rank = generic_gauss_elim_hensel(m,&nm,&dn,&ri,&ci); |
| |
else |
| |
rank = generic_gauss_elim(m,&nm,&dn,&ri,&ci); |
| |
t = col-rank; |
| |
MKVECT(rind,rank); |
| |
MKVECT(cind,t); |
| |
for ( i = 0; i < rank; i++ ) { |
| |
STOQ(ri[i],q); |
| |
BDY(rind)[i] = (pointer)q; |
| |
} |
| |
for ( i = 0; i < t; i++ ) { |
| |
STOQ(ci[i],q); |
| |
BDY(cind)[i] = (pointer)q; |
| |
} |
| |
n0 = mknode(4,nm,dn,rind,cind); |
| |
MKLIST(*rp,n0); |
| |
} |
| |
|
| |
void Pindep_rows_mod(NODE arg,VECT *rp) |
| |
{ |
| |
MAT m,mat; |
| |
VECT rind; |
| |
Q **tmat; |
| |
int **wmat,**row0; |
| |
Q *rib; |
| |
int *rowstat,*p; |
| |
Q q; |
| |
int md,i,j,k,l,row,col,t,rank; |
| |
|
| |
asir_assert(ARG0(arg),O_MAT,"indep_rows_mod"); |
| |
asir_assert(ARG1(arg),O_N,"indep_rows_mod"); |
| |
m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg)); |
| |
row = m->row; col = m->col; tmat = (Q **)m->body; |
| |
wmat = (int **)almat(row,col); |
| |
|
| |
row0 = (int **)ALLOCA(row*sizeof(int *)); |
| |
for ( i = 0; i < row; i++ ) row0[i] = wmat[i]; |
| |
|
| |
rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int)); |
| |
for ( i = 0; i < row; i++ ) |
| |
for ( j = 0; j < col; j++ ) |
| |
if ( q = (Q)tmat[i][j] ) { |
| |
t = rem(NM(q),md); |
| |
if ( t && SGN(q) < 0 ) |
| |
t = (md - t) % md; |
| |
wmat[i][j] = t; |
| |
} else |
| |
wmat[i][j] = 0; |
| |
rank = indep_rows_mod(wmat,row,col,md,rowstat); |
| |
|
| |
MKVECT(rind,rank); |
| |
rib = (Q *)rind->body; |
| |
for ( j = 0; j < rank; j++ ) { |
| |
STOQ(rowstat[j],rib[j]); |
| |
} |
| |
*rp = rind; |
| |
} |
| |
|
| |
/* |
| |
input : a row x col matrix A |
| |
A[I] <-> A[I][0]*x_0+A[I][1]*x_1+... |
| |
|
| output : [B,R,C] |
output : [B,R,C] |
| B : a rank(A) x col-rank(A) matrix |
B : a rank(A) x col-rank(A) matrix |
| 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) |
| |
RN : a vector of length rank(A) indicating useful rows |
| |
|
| B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+... |
B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+... |
| */ |
*/ |
| |
|
| void Pgeneric_gauss_elim_mod(arg,rp) |
void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp) |
| NODE arg; |
|
| LIST *rp; |
|
| { |
{ |
| NODE n0; |
NODE n0; |
| MAT m,mat; |
MAT m,mat; |
| VECT rind,cind; |
VECT rind,cind,rnum; |
| Q **tmat; |
Q **tmat; |
| int **wmat; |
int **wmat,**row0; |
| Q *rib,*cib; |
Q *rib,*cib,*rnb; |
| int *colstat; |
int *colstat,*p; |
| Q q; |
Q q; |
| int md,i,j,k,l,row,col,t,n,rank; |
int md,i,j,k,l,row,col,t,rank; |
| |
|
| asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod"); |
asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod"); |
| asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod"); |
asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod"); |
| m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg)); |
m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg)); |
| row = m->row; col = m->col; tmat = (Q **)m->body; |
row = m->row; col = m->col; tmat = (Q **)m->body; |
| wmat = (int **)almat(row,col); |
wmat = (int **)almat(row,col); |
| |
|
| |
row0 = (int **)ALLOCA(row*sizeof(int *)); |
| |
for ( i = 0; i < row; i++ ) row0[i] = wmat[i]; |
| |
|
| colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
| for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) |
| for ( j = 0; j < col; j++ ) |
for ( j = 0; j < col; j++ ) |
|
|
| wmat[i][j] = 0; |
wmat[i][j] = 0; |
| rank = generic_gauss_elim_mod(wmat,row,col,md,colstat); |
rank = generic_gauss_elim_mod(wmat,row,col,md,colstat); |
| |
|
| |
MKVECT(rnum,rank); |
| |
rnb = (Q *)rnum->body; |
| |
for ( i = 0; i < rank; i++ ) |
| |
for ( j = 0, p = wmat[i]; j < row; j++ ) |
| |
if ( p == row0[j] ) |
| |
STOQ(j,rnb[i]); |
| |
|
| MKMAT(mat,rank,col-rank); |
MKMAT(mat,rank,col-rank); |
| tmat = (Q **)mat->body; |
tmat = (Q **)mat->body; |
| for ( i = 0; i < rank; i++ ) |
for ( i = 0; i < rank; i++ ) |
|
|
| } else { |
} else { |
| STOQ(j,cib[l]); l++; |
STOQ(j,cib[l]); l++; |
| } |
} |
| n0 = mknode(3,mat,rind,cind); |
n0 = mknode(4,mat,rind,cind,rnum); |
| MKLIST(*rp,n0); |
MKLIST(*rp,n0); |
| } |
} |
| |
|
| void Pleqm(arg,rp) |
void Pleqm(NODE arg,VECT *rp) |
| NODE arg; |
|
| VECT *rp; |
|
| { |
{ |
| MAT m; |
MAT m; |
| VECT vect; |
VECT vect; |
|
|
| } |
} |
| } |
} |
| |
|
| int gauss_elim_mod(mat,row,col,md) |
int gauss_elim_mod(int **mat,int row,int col,int md) |
| int **mat; |
|
| int row,col,md; |
|
| { |
{ |
| int i,j,k,inv,a,n; |
int i,j,k,inv,a,n; |
| int *t,*pivot; |
int *t,*pivot; |
|
| Line 1356 int row,col,md; |
|
| t = mat[i]; |
t = mat[i]; |
| if ( i != j && (a = t[j]) ) |
if ( i != j && (a = t[j]) ) |
| for ( k = j, a = md - a; k <= n; k++ ) { |
for ( k = j, a = md - a; k <= n; k++ ) { |
| |
unsigned int tk; |
| /* t[k] = dmar(pivot[k],a,t[k],md); */ |
/* t[k] = dmar(pivot[k],a,t[k],md); */ |
| DMAR(pivot[k],a,t[k],md,t[k]) |
DMAR(pivot[k],a,t[k],md,tk) |
| |
t[k] = tk; |
| } |
} |
| } |
} |
| } |
} |
|
| Line 1370 int row,col,md; |
|
| return -1; |
return -1; |
| } |
} |
| |
|
| struct oEGT eg_mod,eg_elim,eg_chrem,eg_gschk,eg_intrat,eg_symb; |
struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb; |
| |
struct oEGT eg_conv; |
| |
|
| int generic_gauss_elim(mat,nm,dn,rindp,cindp) |
int generic_gauss_elim(MAT mat,MAT *nm,Q *dn,int **rindp,int **cindp) |
| MAT mat; |
|
| MAT *nm; |
|
| Q *dn; |
|
| int **rindp,**cindp; |
|
| { |
{ |
| int **wmat; |
int **wmat; |
| Q **bmat; |
Q **bmat; |
| Line 708 int **rindp,**cindp; |
|
| Line 1400 int **rindp,**cindp; |
|
| colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
| wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
| for ( ind = 0; ; ind++ ) { |
for ( ind = 0; ; ind++ ) { |
| if ( Print ) { |
if ( DP_Print ) { |
| fprintf(asir_out,"."); fflush(asir_out); |
fprintf(asir_out,"."); fflush(asir_out); |
| } |
} |
| md = lprime[ind]; |
md = get_lprime(ind); |
| get_eg(&tmp0); |
get_eg(&tmp0); |
| for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) |
| for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ ) |
for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ ) |
|
|
| } |
} |
| } else { |
} else { |
| if ( rank < rank0 ) { |
if ( rank < rank0 ) { |
| if ( Print ) { |
if ( DP_Print ) { |
| fprintf(asir_out,"lower rank matrix; continuing...\n"); |
fprintf(asir_out,"lower rank matrix; continuing...\n"); |
| fflush(asir_out); |
fflush(asir_out); |
| } |
} |
| continue; |
continue; |
| } else if ( rank > rank0 ) { |
} else if ( rank > rank0 ) { |
| if ( Print ) { |
if ( DP_Print ) { |
| fprintf(asir_out,"higher rank matrix; resetting...\n"); |
fprintf(asir_out,"higher rank matrix; resetting...\n"); |
| fflush(asir_out); |
fflush(asir_out); |
| } |
} |
|
|
| } else { |
} else { |
| for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ ); |
for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ ); |
| if ( j < col ) { |
if ( j < col ) { |
| if ( Print ) { |
if ( DP_Print ) { |
| fprintf(asir_out,"inconsitent colstat; resetting...\n"); |
fprintf(asir_out,"inconsitent colstat; resetting...\n"); |
| fflush(asir_out); |
fflush(asir_out); |
| } |
} |
|
|
| add_eg(&eg_chrem_split,&tmp0,&tmp1); |
add_eg(&eg_chrem_split,&tmp0,&tmp1); |
| |
|
| get_eg(&tmp0); |
get_eg(&tmp0); |
| ret = intmtoratm(crmat,m1,*nm,dn); |
if ( ind % F4_INTRAT_PERIOD ) |
| |
ret = 0; |
| |
else |
| |
ret = intmtoratm(crmat,m1,*nm,dn); |
| get_eg(&tmp1); |
get_eg(&tmp1); |
| add_eg(&eg_intrat,&tmp0,&tmp1); |
add_eg(&eg_intrat,&tmp0,&tmp1); |
| add_eg(&eg_intrat_split,&tmp0,&tmp1); |
add_eg(&eg_intrat_split,&tmp0,&tmp1); |
|
|
| get_eg(&tmp1); |
get_eg(&tmp1); |
| add_eg(&eg_gschk,&tmp0,&tmp1); |
add_eg(&eg_gschk,&tmp0,&tmp1); |
| add_eg(&eg_gschk_split,&tmp0,&tmp1); |
add_eg(&eg_gschk_split,&tmp0,&tmp1); |
| if ( Print ) { |
if ( DP_Print ) { |
| print_eg("Mod",&eg_mod_split); |
print_eg("Mod",&eg_mod_split); |
| print_eg("Elim",&eg_elim_split); |
print_eg("Elim",&eg_elim_split); |
| print_eg("ChRem",&eg_chrem_split); |
print_eg("ChRem",&eg_chrem_split); |
|
|
| } |
} |
| } |
} |
| |
|
| int generic_gauss_elim_hensel(mat,nmmat,dn,rindp,cindp) |
void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm); |
| MAT mat; |
|
| MAT *nmmat; |
/* XXX broken */ |
| Q *dn; |
void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm) |
| int **rindp,**cindp; |
|
| { |
{ |
| |
Q **a0,**b; |
| |
Q *aiq; |
| |
N **a; |
| |
N *ai; |
| |
Q q,q1,dn2,a1,q0,bik; |
| |
MAT m; |
| |
unsigned int md; |
| |
int n,ind,i,j,rank,t,inv,t1,ret,min,k; |
| |
int **w; |
| |
int *wi,*rinfo0,*rinfo; |
| |
N m1,m2,m3,u,s; |
| |
|
| |
a0 = (Q **)mat->body; |
| |
n = mat->row; |
| |
if ( n != mat->col ) |
| |
error("lu_dec_cr : non-square matrix"); |
| |
w = (int **)almat(n,n); |
| |
MKMAT(m,n,n); |
| |
a = (N **)m->body; |
| |
UTON(1,m1); |
| |
rinfo0 = 0; |
| |
ind = 0; |
| |
while ( 1 ) { |
| |
md = get_lprime(ind); |
| |
/* mat mod md */ |
| |
for ( i = 0; i < n; i++ ) |
| |
for ( j = 0, aiq = a0[i], wi = w[i]; j < n; j++ ) |
| |
if ( q = aiq[j] ) { |
| |
t = rem(NM(q),md); |
| |
if ( t && SGN(q) < 0 ) |
| |
t = (md - t) % md; |
| |
wi[j] = t; |
| |
} else |
| |
wi[j] = 0; |
| |
|
| |
if ( !lu_mod((unsigned int **)w,n,md,&rinfo) ) continue; |
| |
printf("."); fflush(stdout); |
| |
if ( !rinfo0 ) |
| |
*perm = rinfo0 = rinfo; |
| |
else { |
| |
for ( i = 0; i < n; i++ ) |
| |
if ( rinfo[i] != rinfo0[i] ) break; |
| |
if ( i < n ) continue; |
| |
} |
| |
if ( UNIN(m1) ) { |
| |
for ( i = 0; i < n; i++ ) |
| |
for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ ) { |
| |
UTON(wi[j],u); ai[j] = u; |
| |
} |
| |
UTON(md,m1); |
| |
} else { |
| |
inv = invm(rem(m1,md),md); |
| |
UTON(md,m2); muln(m1,m2,&m3); |
| |
for ( i = 0; i < n; i++ ) |
| |
for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ ) |
| |
if ( ai[i] ) { |
| |
/* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */ |
| |
t = rem(ai[j],md); |
| |
if ( wi[j] >= t ) |
| |
t = wi[j]-t; |
| |
else |
| |
t = md-(t-wi[j]); |
| |
DMAR(t,inv,0,md,t1) |
| |
UTON(t1,u); |
| |
muln(m1,u,&s); |
| |
addn(ai[j],s,&u); ai[j] = u; |
| |
} else if ( wi[j] ) { |
| |
/* f3 = m1*(m1 mod m2)^(-1)*f2 */ |
| |
DMAR(wi[j],inv,0,md,t) |
| |
UTON(t,u); |
| |
muln(m1,u,&s); ai[j] = s; |
| |
} |
| |
m1 = m3; |
| |
} |
| |
if ( (++ind%8) == 0 ) { |
| |
ret = intmtoratm(m,m1,lu,dn); |
| |
if ( ret ) { |
| |
b = (Q **)lu->body; |
| |
mulq(*dn,*dn,&dn2); |
| |
for ( i = 0; i < n; i++ ) { |
| |
for ( j = 0; j < n; j++ ) { |
| |
q = 0; |
| |
min = MIN(i,j); |
| |
for ( k = 0; k <= min; k++ ) { |
| |
bik = k==i ? *dn : b[i][k]; |
| |
mulq(bik,b[k][j],&q0); |
| |
addq(q,q0,&q1); q = q1; |
| |
} |
| |
mulq(a0[rinfo0[i]][j],dn2,&q1); |
| |
if ( cmpq(q,q1) ) break; |
| |
} |
| |
if ( j < n ) break; |
| |
} |
| |
if ( i == n ) |
| |
return; |
| |
} |
| |
} |
| |
} |
| |
} |
| |
|
| |
void nmat(N **m,int n) |
| |
{ |
| |
int i,j; |
| |
|
| |
for ( i = 0; i < n; i++ ) { |
| |
for ( j = 0; j < n; j++ ) { |
| |
printn(m[i][j]); printf(" "); |
| |
} |
| |
printf("\n"); |
| |
} |
| |
} |
| |
|
| |
int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp) |
| |
{ |
| MAT bmat,xmat; |
MAT bmat,xmat; |
| Q **a0,**a,**b,**x,**nm; |
Q **a0,**a,**b,**x,**nm; |
| Q *ai,*bi,*xi; |
Q *ai,*bi,*xi; |
| Line 850 int **rindp,**cindp; |
|
| Line 1658 int **rindp,**cindp; |
|
| int *cinfo,*rinfo; |
int *cinfo,*rinfo; |
| int *rind,*cind; |
int *rind,*cind; |
| int count; |
int count; |
| struct oEGT eg_mul,eg_inv,tmp0,tmp1; |
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; |
| |
|
| a0 = (Q **)mat->body; |
a0 = (Q **)mat->body; |
| row = mat->row; col = mat->col; |
row = mat->row; col = mat->col; |
| w = (int **)almat(row,col); |
w = (int **)almat(row,col); |
| for ( ind = 0; ; ind++ ) { |
for ( ind = 0; ; ind++ ) { |
| md = lprime[ind]; |
md = get_lprime(ind); |
| STOQ(md,mdq); |
STOQ(md,mdq); |
| for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) |
| for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ ) |
for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ ) |
| Line 868 int **rindp,**cindp; |
|
| Line 1682 int **rindp,**cindp; |
|
| } else |
} else |
| wi[j] = 0; |
wi[j] = 0; |
| |
|
| rank = find_lhs_and_lu_mod(w,row,col,md,&rinfo,&cinfo); |
if ( DP_Print > 3 ) { |
| |
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 > 3 ) { |
| |
fprintf(asir_out,"done.\n"); fflush(asir_out); |
| |
} |
| a = (Q **)almat_pointer(rank,rank); /* lhs mat */ |
a = (Q **)almat_pointer(rank,rank); /* lhs mat */ |
| MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */ |
MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */ |
| for ( j = li = ri = 0; j < col; j++ ) |
for ( j = li = ri = 0; j < col; j++ ) |
| Line 898 int **rindp,**cindp; |
|
| Line 1718 int **rindp,**cindp; |
|
| *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int)); |
*cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int)); |
| |
|
| init_eg(&eg_mul); init_eg(&eg_inv); |
init_eg(&eg_mul); init_eg(&eg_inv); |
| for ( q = ONE, count = 0; ; count++ ) { |
init_eg(&eg_check); init_eg(&eg_intrat); |
| fprintf(stderr,"."); |
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 > 3 ) |
| |
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 910 int **rindp,**cindp; |
|
| Line 1738 int **rindp,**cindp; |
|
| 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 */ |
| |
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); |
get_eg(&tmp1); |
| add_eg(&eg_inv,&tmp0,&tmp1); |
add_eg(&eg_inv,&tmp0,&tmp1); |
| /* x = x-q*wc */ |
get_eg(&tmp0); |
| 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; j < ri; j++ ) { |
| STOQ(wi[j],u); mulq(q,u,&s); |
inner_product_mat_int_mod(a,wc,rank,i,j,&u); |
| subq(xi[j],s,&u); xi[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 ) { |
| |
rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); |
| |
cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int)); |
| |
for ( j = k = l = 0; j < col; j++ ) |
| |
if ( cinfo[j] ) |
| |
rind[k++] = j; |
| |
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,DP *mb,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 m; |
| |
int col1; |
| |
|
| |
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); |
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++ ) |
| |
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++ ) { |
for ( j = 0; j < ri; j++ ) { |
| inner_product_mat_int_mod(a,wc,rank,i,j,&u); |
inner_product_mat_int_mod(a,wc,rank,i,j,&u); |
| addq(b[i][j],u,&s); |
addq(b[i][j],u,&s); |
| Line 938 int **rindp,**cindp; |
|
| Line 1960 int **rindp,**cindp; |
|
| 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 % 2) && intmtoratm_q(xmat,NM(q),*nmmat,dn) ) { |
if ( count == period ) { |
| for ( j = k = l = 0; j < col; j++ ) |
get_eg(&tmp0); |
| if ( cinfo[j] ) |
ptr = wx; |
| rind[k++] = j; |
for ( i = 0; i < rank; i++ ) |
| else |
for ( j = 0, xi = x[i]; j < ri; |
| cind[l++] = j; |
j++, ptr += nsize ) { |
| if ( gensolve_check(mat,*nmmat,*dn,rind,cind) ) { |
for ( k = nsize-1; k >= 0 && !ptr[k]; k-- ); |
| fprintf(stderr,"\n"); |
if ( k >= 0 ) { |
| print_eg("INV",&eg_inv); |
wn = NALLOC(k+1); |
| print_eg("MUL",&eg_mul); |
PL(wn) = k+1; |
| fflush(asir_out); |
for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l]; |
| return rank; |
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; |
| |
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); |
| |
} |
| |
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; |
| } |
} |
| } |
} |
| } |
} |
| Line 958 int **rindp,**cindp; |
|
| Line 2012 int **rindp,**cindp; |
|
| |
|
| int f4_nocheck; |
int f4_nocheck; |
| |
|
| int gensolve_check(mat,nm,dn,rind,cind) |
int gensolve_check(MAT mat,MAT nm,Q dn,int *rind,int *cind) |
| MAT mat,nm; |
|
| Q dn; |
|
| int *rind,*cind; |
|
| { |
{ |
| int row,col,rank,clen,i,j,k,l; |
int row,col,rank,clen,i,j,k,l; |
| Q s,t,u; |
Q s,t; |
| Q *w; |
Q *w; |
| Q *mati,*nmk; |
Q *mati,*nmk; |
| |
|
| Line 1009 int *rind,*cind; |
|
| Line 2060 int *rind,*cind; |
|
| |
|
| /* assuming 0 < c < m */ |
/* assuming 0 < c < m */ |
| |
|
| int inttorat(c,m,b,sgnp,nmp,dnp) |
int inttorat(N c,N m,N b,int *sgnp,N *nmp,N *dnp) |
| N c,m,b; |
|
| int *sgnp; |
|
| N *nmp,*dnp; |
|
| { |
{ |
| Q qq,t,u1,v1,r1,nm; |
Q qq,t,u1,v1,r1; |
| N q,r,u2,v2,r2; |
N q,u2,v2,r2; |
| |
|
| u1 = 0; v1 = ONE; u2 = m; v2 = c; |
u1 = 0; v1 = ONE; u2 = m; v2 = c; |
| while ( cmpn(v2,b) >= 0 ) { |
while ( cmpn(v2,b) >= 0 ) { |
|
|
| |
|
| /* mat->body = N ** */ |
/* mat->body = N ** */ |
| |
|
| int intmtoratm(mat,md,nm,dn) |
int intmtoratm(MAT mat,N md,MAT nm,Q *dn) |
| MAT mat; |
|
| N md; |
|
| MAT nm; |
|
| Q *dn; |
|
| { |
{ |
| N t,s,b; |
N t,s,b; |
| Q bound,dn0,dn1,nm1,q,tq; |
Q dn0,dn1,nm1,q; |
| int i,j,k,l,row,col; |
int i,j,k,l,row,col; |
| Q **rmat; |
Q **rmat; |
| N **tmat; |
N **tmat; |
|
|
| |
|
| /* mat->body = Q ** */ |
/* mat->body = Q ** */ |
| |
|
| int intmtoratm_q(mat,md,nm,dn) |
int intmtoratm_q(MAT mat,N md,MAT nm,Q *dn) |
| MAT mat; |
|
| N md; |
|
| MAT nm; |
|
| Q *dn; |
|
| { |
{ |
| N t,s,b; |
N t,s,b; |
| Q bound,dn0,dn1,nm1,q,tq; |
Q dn0,dn1,nm1,q; |
| int i,j,k,l,row,col; |
int i,j,k,l,row,col; |
| Q **rmat; |
Q **rmat; |
| Q **tmat; |
Q **tmat; |
|
|
| return 1; |
return 1; |
| } |
} |
| |
|
| int generic_gauss_elim_mod(mat,row,col,md,colstat) |
#define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++; |
| int **mat; |
|
| int row,col,md; |
void reduce_reducers_mod(int **mat,int row,int col,int md) |
| int *colstat; |
|
| { |
{ |
| |
int i,j,k,l,hc,zzz; |
| |
int *t,*s,*tj,*ind; |
| |
|
| |
/* reduce the reducers */ |
| |
ind = (int *)ALLOCA(row*sizeof(int)); |
| |
for ( i = 0; i < row; i++ ) { |
| |
t = mat[i]; |
| |
for ( j = 0; j < col && !t[j]; j++ ); |
| |
/* register the position of the head term */ |
| |
ind[i] = j; |
| |
for ( l = i-1; l >= 0; l-- ) { |
| |
/* reduce mat[i] by mat[l] */ |
| |
if ( hc = t[ind[l]] ) { |
| |
/* mat[i] = mat[i]-hc*mat[l] */ |
| |
j = ind[l]; |
| |
s = mat[l]+j; |
| |
tj = t+j; |
| |
hc = md-hc; |
| |
k = col-j; |
| |
for ( ; k >= 64; k -= 64 ) { |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
| |
} |
| |
for ( ; k > 0; k-- ) { |
| |
if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++; |
| |
} |
| |
} |
| |
} |
| |
} |
| |
} |
| |
|
| |
/* |
| |
mat[i] : reducers (i=0,...,nred-1) |
| |
spolys (i=nred,...,row-1) |
| |
mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order |
| |
1. reduce the reducers |
| |
2. reduce spolys by the reduced reducers |
| |
*/ |
| |
|
| |
void pre_reduce_mod(int **mat,int row,int col,int nred,int md) |
| |
{ |
| |
int i,j,k,l,hc,inv; |
| |
int *t,*s,*tk,*ind; |
| |
|
| |
#if 1 |
| |
/* reduce the reducers */ |
| |
ind = (int *)ALLOCA(row*sizeof(int)); |
| |
for ( i = 0; i < nred; i++ ) { |
| |
/* make mat[i] monic and mat[i] by mat[0],...,mat[i-1] */ |
| |
t = mat[i]; |
| |
for ( j = 0; j < col && !t[j]; j++ ); |
| |
/* register the position of the head term */ |
| |
ind[i] = j; |
| |
inv = invm(t[j],md); |
| |
for ( k = j; k < col; k++ ) |
| |
if ( t[k] ) |
| |
DMAR(t[k],inv,0,md,t[k]) |
| |
for ( l = i-1; l >= 0; l-- ) { |
| |
/* reduce mat[i] by mat[l] */ |
| |
if ( hc = t[ind[l]] ) { |
| |
/* mat[i] = mat[i]-hc*mat[l] */ |
| |
for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k; |
| |
k < col; k++, tk++, s++ ) |
| |
if ( *s ) |
| |
DMAR(*s,hc,*tk,md,*tk) |
| |
} |
| |
} |
| |
} |
| |
/* reduce the spolys */ |
| |
for ( i = nred; i < row; i++ ) { |
| |
t = mat[i]; |
| |
for ( l = nred-1; l >= 0; l-- ) { |
| |
/* reduce mat[i] by mat[l] */ |
| |
if ( hc = t[ind[l]] ) { |
| |
/* mat[i] = mat[i]-hc*mat[l] */ |
| |
for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k; |
| |
k < col; k++, tk++, s++ ) |
| |
if ( *s ) |
| |
DMAR(*s,hc,*tk,md,*tk) |
| |
} |
| |
} |
| |
} |
| |
#endif |
| |
} |
| |
/* |
| |
mat[i] : reducers (i=0,...,nred-1) |
| |
mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order |
| |
*/ |
| |
|
| |
void reduce_sp_by_red_mod(int *sp,int **redmat,int *ind,int nred,int col,int md) |
| |
{ |
| |
int i,j,k,hc,zzz; |
| |
int *s,*tj; |
| |
|
| |
/* reduce the spolys by redmat */ |
| |
for ( i = nred-1; i >= 0; i-- ) { |
| |
/* reduce sp by redmat[i] */ |
| |
if ( hc = sp[ind[i]] ) { |
| |
/* sp = sp-hc*redmat[i] */ |
| |
j = ind[i]; |
| |
hc = md-hc; |
| |
s = redmat[i]+j; |
| |
tj = sp+j; |
| |
for ( k = col-j; k > 0; k-- ) { |
| |
if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++; |
| |
} |
| |
} |
| |
} |
| |
} |
| |
|
| |
/* |
| |
mat[i] : compressed reducers (i=0,...,nred-1) |
| |
mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order |
| |
*/ |
| |
|
| |
void red_by_compress(int m,unsigned int *p,unsigned int *r, |
| |
unsigned int *ri,unsigned int hc,int len) |
| |
{ |
| |
unsigned int up,lo; |
| |
unsigned int dmy; |
| |
unsigned int *pj; |
| |
|
| |
p[*ri] = 0; r++; ri++; |
| |
for ( len--; len; len--, r++, ri++ ) { |
| |
pj = p+ *ri; |
| |
DMA(*r,hc,*pj,up,lo); |
| |
if ( up ) { |
| |
DSAB(m,up,lo,dmy,*pj); |
| |
} else |
| |
*pj = lo; |
| |
} |
| |
} |
| |
|
| |
/* p -= hc*r */ |
| |
|
| |
void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len) |
| |
{ |
| |
unsigned int up,lo,dmy; |
| |
|
| |
*p++ = 0; r++; len--; |
| |
for ( ; len; len--, r++, p++ ) |
| |
if ( *r ) { |
| |
DMA(*r,hc,*p,up,lo); |
| |
if ( up ) { |
| |
DSAB(m,up,lo,dmy,*p); |
| |
} else |
| |
*p = lo; |
| |
} |
| |
} |
| |
|
| |
#if defined(__GNUC__) && SIZEOF_LONG==8 |
| |
/* 64bit vector += UNIT vector(normalized) */ |
| |
|
| |
void red_by_vect64(int m, U64 *p,unsigned int *c,U64 *r,unsigned int hc,int len) |
| |
{ |
| |
U64 t; |
| |
|
| |
/* (p[0],c[0]) is normalized */ |
| |
*p++ = 0; *c++ = 0; r++; len--; |
| |
for ( ; len; len--, r++, p++, c++ ) |
| |
if ( *r ) { |
| |
t = (*p)+(*r)*hc; |
| |
if ( t < *p ) (*c)++; |
| |
*p = t; |
| |
} |
| |
} |
| |
#endif |
| |
|
| |
void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len) |
| |
{ |
| |
*p++ = 0; r++; len--; |
| |
for ( ; len; len--, r++, p++ ) |
| |
if ( *r ) |
| |
*p = _addsf(_mulsf(*r,hc),*p); |
| |
} |
| |
|
| |
extern GZ current_mod_lf; |
| |
extern int current_mod_lf_size; |
| |
|
| |
void red_by_vect_lf(mpz_t *p,mpz_t *r,mpz_t hc,int len) |
| |
{ |
| |
mpz_set_ui(*p++,0); r++; len--; |
| |
for ( ; len; len--, r++, p++ ) { |
| |
mpz_addmul(*p,*r,hc); |
| |
#if 0 |
| |
if ( mpz_size(*p) > current_mod_lf_size ) |
| |
mpz_mod(*p,*p,BDY(current_mod_lf)); |
| |
#endif |
| |
} |
| |
} |
| |
|
| |
|
| |
extern unsigned int **psca; |
| |
|
| |
void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind, |
| |
int nred,int col,int md) |
| |
{ |
| |
int i,len; |
| |
CDP ri; |
| |
unsigned int hc; |
| |
unsigned int *usp; |
| |
|
| |
usp = (unsigned int *)sp; |
| |
/* reduce the spolys by redmat */ |
| |
for ( i = nred-1; i >= 0; i-- ) { |
| |
/* reduce sp by redmat[i] */ |
| |
usp[ind[i]] %= md; |
| |
if ( hc = usp[ind[i]] ) { |
| |
/* sp = sp-hc*redmat[i] */ |
| |
hc = md-hc; |
| |
ri = redmat[i]; |
| |
len = ri->len; |
| |
red_by_compress(md,usp,psca[ri->psindex],ri->body,hc,len); |
| |
} |
| |
} |
| |
for ( i = 0; i < col; i++ ) |
| |
if ( usp[i] >= (unsigned int)md ) |
| |
usp[i] %= md; |
| |
} |
| |
|
| |
#define ONE_STEP2 if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++; |
| |
|
| |
int generic_gauss_elim_mod(int **mat0,int row,int col,int md,int *colstat) |
| |
{ |
| int i,j,k,l,inv,a,rank; |
int i,j,k,l,inv,a,rank; |
| int *t,*pivot; |
unsigned int *t,*pivot,*pk; |
| |
unsigned int **mat; |
| |
|
| |
mat = (unsigned int **)mat0; |
| for ( rank = 0, j = 0; j < col; j++ ) { |
for ( rank = 0, j = 0; j < col; j++ ) { |
| for ( i = rank; i < row && !mat[i][j]; i++ ); |
for ( i = rank; i < row; i++ ) |
| |
mat[i][j] %= md; |
| |
for ( i = rank; i < row; i++ ) |
| |
if ( mat[i][j] ) |
| |
break; |
| if ( i == row ) { |
if ( i == row ) { |
| colstat[j] = 0; |
colstat[j] = 0; |
| continue; |
continue; |
|
|
| } |
} |
| pivot = mat[rank]; |
pivot = mat[rank]; |
| inv = invm(pivot[j],md); |
inv = invm(pivot[j],md); |
| for ( k = j; k < col; k++ ) |
for ( k = j, pk = pivot+k; k < col; k++, pk++ ) |
| if ( pivot[k] ) { |
if ( *pk ) { |
| DMAR(pivot[k],inv,0,md,pivot[k]) |
if ( *pk >= (unsigned int)md ) |
| |
*pk %= md; |
| |
DMAR(*pk,inv,0,md,*pk) |
| } |
} |
| for ( i = rank+1; i < row; i++ ) { |
for ( i = rank+1; i < row; i++ ) { |
| t = mat[i]; |
t = mat[i]; |
| if ( a = t[j] ) |
if ( a = t[j] ) |
| for ( k = j, a = md - a; k < col; k++ ) |
red_by_vect(md,t+j,pivot+j,md-a,col-j); |
| if ( pivot[k] ) { |
|
| DMAR(pivot[k],a,t[k],md,t[k]) |
|
| } |
|
| } |
} |
| rank++; |
rank++; |
| } |
} |
|
|
| pivot = mat[l]; |
pivot = mat[l]; |
| for ( i = 0; i < l; i++ ) { |
for ( i = 0; i < l; i++ ) { |
| t = mat[i]; |
t = mat[i]; |
| |
t[j] %= md; |
| if ( a = t[j] ) |
if ( a = t[j] ) |
| for ( k = j, a = md-a; k < col; k++ ) |
red_by_vect(md,t+j,pivot+j,md-a,col-j); |
| if ( pivot[k] ) { |
|
| DMAR(pivot[k],a,t[k],md,t[k]) |
|
| } |
|
| } |
} |
| l--; |
l--; |
| } |
} |
| |
for ( j = 0, l = 0; l < rank; j++ ) |
| |
if ( colstat[j] ) { |
| |
t = mat[l]; |
| |
for ( k = j; k < col; k++ ) |
| |
if ( t[k] >= (unsigned int)md ) |
| |
t[k] %= md; |
| |
l++; |
| |
} |
| return rank; |
return rank; |
| } |
} |
| |
|
| |
int generic_gauss_elim_mod2(int **mat0,int row,int col,int md,int *colstat,int *rowstat) |
| |
{ |
| |
int i,j,k,l,inv,a,rank; |
| |
unsigned int *t,*pivot,*pk; |
| |
unsigned int **mat; |
| |
|
| |
for ( i = 0; i < row; i++ ) rowstat[i] = i; |
| |
mat = (unsigned int **)mat0; |
| |
for ( rank = 0, j = 0; j < col; j++ ) { |
| |
for ( i = rank; i < row; i++ ) |
| |
mat[i][j] %= md; |
| |
for ( i = rank; i < row; i++ ) |
| |
if ( mat[i][j] ) |
| |
break; |
| |
if ( i == row ) { |
| |
colstat[j] = 0; |
| |
continue; |
| |
} else |
| |
colstat[j] = 1; |
| |
if ( i != rank ) { |
| |
t = mat[i]; mat[i] = mat[rank]; mat[rank] = t; |
| |
k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k; |
| |
} |
| |
pivot = mat[rank]; |
| |
inv = invm(pivot[j],md); |
| |
for ( k = j, pk = pivot+k; k < col; k++, pk++ ) |
| |
if ( *pk ) { |
| |
if ( *pk >= (unsigned int)md ) |
| |
*pk %= md; |
| |
DMAR(*pk,inv,0,md,*pk) |
| |
} |
| |
for ( i = rank+1; i < row; i++ ) { |
| |
t = mat[i]; |
| |
if ( a = t[j] ) |
| |
red_by_vect(md,t+j,pivot+j,md-a,col-j); |
| |
} |
| |
rank++; |
| |
} |
| |
for ( j = col-1, l = rank-1; j >= 0; j-- ) |
| |
if ( colstat[j] ) { |
| |
pivot = mat[l]; |
| |
for ( i = 0; i < l; i++ ) { |
| |
t = mat[i]; |
| |
t[j] %= md; |
| |
if ( a = t[j] ) |
| |
red_by_vect(md,t+j,pivot+j,md-a,col-j); |
| |
} |
| |
l--; |
| |
} |
| |
for ( j = 0, l = 0; l < rank; j++ ) |
| |
if ( colstat[j] ) { |
| |
t = mat[l]; |
| |
for ( k = j; k < col; k++ ) |
| |
if ( t[k] >= (unsigned int)md ) |
| |
t[k] %= md; |
| |
l++; |
| |
} |
| |
return rank; |
| |
} |
| |
|
| |
int indep_rows_mod(int **mat0,int row,int col,int md,int *rowstat) |
| |
{ |
| |
int i,j,k,l,inv,a,rank; |
| |
unsigned int *t,*pivot,*pk; |
| |
unsigned int **mat; |
| |
|
| |
for ( i = 0; i < row; i++ ) rowstat[i] = i; |
| |
mat = (unsigned int **)mat0; |
| |
for ( rank = 0, j = 0; j < col; j++ ) { |
| |
for ( i = rank; i < row; i++ ) |
| |
mat[i][j] %= md; |
| |
for ( i = rank; i < row; i++ ) |
| |
if ( mat[i][j] ) |
| |
break; |
| |
if ( i == row ) continue; |
| |
if ( i != rank ) { |
| |
t = mat[i]; mat[i] = mat[rank]; mat[rank] = t; |
| |
k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k; |
| |
} |
| |
pivot = mat[rank]; |
| |
inv = invm(pivot[j],md); |
| |
for ( k = j, pk = pivot+k; k < col; k++, pk++ ) |
| |
if ( *pk ) { |
| |
if ( *pk >= (unsigned int)md ) |
| |
*pk %= md; |
| |
DMAR(*pk,inv,0,md,*pk) |
| |
} |
| |
for ( i = rank+1; i < row; i++ ) { |
| |
t = mat[i]; |
| |
if ( a = t[j] ) |
| |
red_by_vect(md,t+j,pivot+j,md-a,col-j); |
| |
} |
| |
rank++; |
| |
} |
| |
return rank; |
| |
} |
| |
|
| |
int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat) |
| |
{ |
| |
int i,j,k,l,inv,a,rank; |
| |
unsigned int *t,*pivot,*pk; |
| |
unsigned int **mat; |
| |
|
| |
mat = (unsigned int **)mat0; |
| |
for ( rank = 0, j = 0; j < col; j++ ) { |
| |
for ( i = rank; i < row; i++ ) |
| |
if ( mat[i][j] ) |
| |
break; |
| |
if ( i == row ) { |
| |
colstat[j] = 0; |
| |
continue; |
| |
} else |
| |
colstat[j] = 1; |
| |
if ( i != rank ) { |
| |
t = mat[i]; mat[i] = mat[rank]; mat[rank] = t; |
| |
} |
| |
pivot = mat[rank]; |
| |
inv = _invsf(pivot[j]); |
| |
for ( k = j, pk = pivot+k; k < col; k++, pk++ ) |
| |
if ( *pk ) |
| |
*pk = _mulsf(*pk,inv); |
| |
for ( i = rank+1; i < row; i++ ) { |
| |
t = mat[i]; |
| |
if ( a = t[j] ) |
| |
red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j); |
| |
} |
| |
rank++; |
| |
} |
| |
for ( j = col-1, l = rank-1; j >= 0; j-- ) |
| |
if ( colstat[j] ) { |
| |
pivot = mat[l]; |
| |
for ( i = 0; i < l; i++ ) { |
| |
t = mat[i]; |
| |
if ( a = t[j] ) |
| |
red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j); |
| |
} |
| |
l--; |
| |
} |
| |
return rank; |
| |
} |
| |
|
| /* LU decomposition; a[i][i] = 1/U[i][i] */ |
/* LU decomposition; a[i][i] = 1/U[i][i] */ |
| |
|
| int lu_gfmmat(mat,md,perm) |
int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm) |
| GFMMAT mat; |
|
| unsigned int md; |
|
| int *perm; |
|
| { |
{ |
| int row,col; |
int row,col; |
| int i,j,k,l; |
int i,j,k; |
| unsigned int *t,*pivot; |
unsigned int *t,*pivot; |
| unsigned int **a; |
unsigned int **a; |
| unsigned int inv,m; |
unsigned int inv,m; |
|
|
| DMAR(inv,m,0,md,t[k]) |
DMAR(inv,m,0,md,t[k]) |
| for ( j = k+1, m = md - t[k]; j < col; j++ ) |
for ( j = k+1, m = md - t[k]; j < col; j++ ) |
| if ( pivot[j] ) { |
if ( pivot[j] ) { |
| DMAR(m,pivot[j],t[j],md,t[j]) |
unsigned int tj; |
| |
|
| |
DMAR(m,pivot[j],t[j],md,tj) |
| |
t[j] = tj; |
| } |
} |
| } |
} |
| } |
} |
|
|
| cinfo[j]=1 <=> j-th column is contained in the LU decomp. |
cinfo[j]=1 <=> j-th column is contained in the LU decomp. |
| */ |
*/ |
| |
|
| int find_lhs_and_lu_mod(a,row,col,md,rinfo,cinfo) |
int find_lhs_and_lu_mod(unsigned int **a,int row,int col, |
| unsigned int **a; |
unsigned int md,int **rinfo,int **cinfo) |
| unsigned int md; |
|
| int **rinfo,**cinfo; |
|
| { |
{ |
| int i,j,k,l,d; |
int i,j,k,d; |
| int *rp,*cp; |
int *rp,*cp; |
| unsigned int *t,*pivot; |
unsigned int *t,*pivot; |
| unsigned int inv,m; |
unsigned int inv,m; |
| Line 1288 int **rinfo,**cinfo; |
|
| Line 2715 int **rinfo,**cinfo; |
|
| DMAR(inv,m,0,md,t[k]) |
DMAR(inv,m,0,md,t[k]) |
| for ( j = k+1, m = md - t[k]; j < col; j++ ) |
for ( j = k+1, m = md - t[k]; j < col; j++ ) |
| if ( pivot[j] ) { |
if ( pivot[j] ) { |
| DMAR(m,pivot[j],t[j],md,t[j]) |
unsigned int tj; |
| |
DMAR(m,pivot[j],t[j],md,tj) |
| |
t[j] = tj; |
| } |
} |
| } |
} |
| } |
} |
| Line 1297 int **rinfo,**cinfo; |
|
| Line 2726 int **rinfo,**cinfo; |
|
| return d; |
return d; |
| } |
} |
| |
|
| |
int lu_mod(unsigned int **a,int n,unsigned int md,int **rinfo) |
| |
{ |
| |
int i,j,k; |
| |
int *rp; |
| |
unsigned int *t,*pivot; |
| |
unsigned int inv,m; |
| |
|
| |
*rinfo = rp = (int *)MALLOC_ATOMIC(n*sizeof(int)); |
| |
for ( i = 0; i < n; i++ ) rp[i] = i; |
| |
for ( k = 0; k < n; k++ ) { |
| |
for ( i = k; i < n && !a[i][k]; i++ ); |
| |
if ( i == n ) return 0; |
| |
if ( i != k ) { |
| |
j = rp[i]; rp[i] = rp[k]; rp[k] = j; |
| |
t = a[i]; a[i] = a[k]; a[k] = t; |
| |
} |
| |
pivot = a[k]; |
| |
inv = invm(pivot[k],md); |
| |
for ( i = k+1; i < n; i++ ) { |
| |
t = a[i]; |
| |
if ( m = t[k] ) { |
| |
DMAR(inv,m,0,md,t[k]) |
| |
for ( j = k+1, m = md - t[k]; j < n; j++ ) |
| |
if ( pivot[j] ) { |
| |
unsigned int tj; |
| |
DMAR(m,pivot[j],t[j],md,tj) |
| |
t[j] = tj; |
| |
} |
| |
} |
| |
} |
| |
} |
| |
return 1; |
| |
} |
| |
|
| /* |
/* |
| Input |
Input |
| a : n x n matrix; a result of LU-decomposition |
a : n x n matrix; a result of LU-decomposition |
| Line 1306 int **rinfo,**cinfo; |
|
| Line 2769 int **rinfo,**cinfo; |
|
| b = a^(-1)b |
b = a^(-1)b |
| */ |
*/ |
| |
|
| void solve_by_lu_mod(a,n,md,b,l) |
void solve_by_lu_mod(int **a,int n,int md,int **b,int l,int normalize) |
| int **a; |
|
| int n; |
|
| int md; |
|
| int **b; |
|
| int l; |
|
| { |
{ |
| unsigned int *y,*c; |
unsigned int *y,*c; |
| int i,j,k; |
int i,j,k; |
|
|
| 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]; |
| } |
} |
| } |
} |
| |
|
| void Pleqm1(arg,rp) |
void Pleqm1(NODE arg,VECT *rp) |
| NODE arg; |
|
| VECT *rp; |
|
| { |
{ |
| MAT m; |
MAT m; |
| VECT vect; |
VECT vect; |
|
|
| } |
} |
| } |
} |
| |
|
| gauss_elim_mod1(mat,row,col,md) |
int gauss_elim_mod1(int **mat,int row,int col,int md) |
| int **mat; |
|
| int row,col,md; |
|
| { |
{ |
| int i,j,k,inv,a,n; |
int i,j,k,inv,a,n; |
| int *t,*pivot; |
int *t,*pivot; |
| Line 1429 int row,col,md; |
|
| Line 2887 int row,col,md; |
|
| return -1; |
return -1; |
| } |
} |
| |
|
| void Pgeninvm(arg,rp) |
void Pgeninvm(NODE arg,LIST *rp) |
| NODE arg; |
|
| LIST *rp; |
|
| { |
{ |
| MAT m; |
MAT m; |
| pointer **mat; |
pointer **mat; |
|
|
| MKMAT(mat1,col,row); MKMAT(mat2,row-col,row); |
MKMAT(mat1,col,row); MKMAT(mat2,row-col,row); |
| for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ ) |
for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ ) |
| for ( j = 0; j < row; j++ ) |
for ( j = 0; j < row; j++ ) |
| STOQ(wmat[i][j+col],tmat[i][j]); |
UTOQ(wmat[i][j+col],tmat[i][j]); |
| for ( tmat = (Q **)mat2->body; i < row; i++ ) |
for ( tmat = (Q **)mat2->body; i < row; i++ ) |
| for ( j = 0; j < row; j++ ) |
for ( j = 0; j < row; j++ ) |
| STOQ(wmat[i][j+col],tmat[i-col][j]); |
UTOQ(wmat[i][j+col],tmat[i-col][j]); |
| MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1); |
MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1); |
| } |
} |
| } |
} |
| |
|
| int gauss_elim_geninv_mod(mat,row,col,md) |
int gauss_elim_geninv_mod(unsigned int **mat,int row,int col,int md) |
| unsigned int **mat; |
|
| int row,col,md; |
|
| { |
{ |
| int i,j,k,inv,a,n,m; |
int i,j,k,inv,a,n,m; |
| unsigned int *t,*pivot; |
unsigned int *t,*pivot; |
| Line 1511 int row,col,md; |
|
| Line 2965 int row,col,md; |
|
| return 0; |
return 0; |
| } |
} |
| |
|
| void Psolve_by_lu_gfmmat(arg,rp) |
void Psolve_by_lu_gfmmat(NODE arg,VECT *rp) |
| NODE arg; |
|
| VECT *rp; |
|
| { |
{ |
| GFMMAT lu; |
GFMMAT lu; |
| Q *perm,*rhs,*v; |
Q *perm,*rhs,*v; |
|
|
| solve_by_lu_gfmmat(lu,md,b,sol); |
solve_by_lu_gfmmat(lu,md,b,sol); |
| MKVECT(r,n); |
MKVECT(r,n); |
| for ( i = 0, v = (Q *)r->body; i < n; i++ ) |
for ( i = 0, v = (Q *)r->body; i < n; i++ ) |
| STOQ(sol[i],v[i]); |
UTOQ(sol[i],v[i]); |
| *rp = r; |
*rp = r; |
| } |
} |
| |
|
| void solve_by_lu_gfmmat(lu,md,b,x) |
void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md, |
| GFMMAT lu; |
unsigned int *b,unsigned int *x) |
| unsigned int md; |
|
| unsigned int *b; |
|
| unsigned int *x; |
|
| { |
{ |
| int n; |
int n; |
| unsigned int **a; |
unsigned int **a; |
| Line 1574 unsigned int *x; |
|
| Line 3023 unsigned int *x; |
|
| } |
} |
| } |
} |
| |
|
| void Plu_gfmmat(arg,rp) |
void Plu_mat(NODE arg,LIST *rp) |
| NODE arg; |
|
| LIST *rp; |
|
| { |
{ |
| |
MAT m,lu; |
| |
Q dn; |
| |
Q *v; |
| |
int n,i; |
| |
int *iperm; |
| |
VECT perm; |
| |
NODE n0; |
| |
|
| |
asir_assert(ARG0(arg),O_MAT,"lu_mat"); |
| |
m = (MAT)ARG0(arg); |
| |
n = m->row; |
| |
MKMAT(lu,n,n); |
| |
lu_dec_cr(m,lu,&dn,&iperm); |
| |
MKVECT(perm,n); |
| |
for ( i = 0, v = (Q *)perm->body; i < n; i++ ) |
| |
STOQ(iperm[i],v[i]); |
| |
n0 = mknode(3,lu,dn,perm); |
| |
MKLIST(*rp,n0); |
| |
} |
| |
|
| |
void Plu_gfmmat(NODE arg,LIST *rp) |
| |
{ |
| MAT m; |
MAT m; |
| GFMMAT mm; |
GFMMAT mm; |
| unsigned int md; |
unsigned int md; |
|
|
| VECT perm; |
VECT perm; |
| NODE n0; |
NODE n0; |
| |
|
| asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat"); |
asir_assert(ARG0(arg),O_MAT,"lu_gfmmat"); |
| asir_assert(ARG1(arg),O_N,"mat_to_gfmmat"); |
asir_assert(ARG1(arg),O_N,"lu_gfmmat"); |
| m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg)); |
m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg)); |
| mat_to_gfmmat(m,md,&mm); |
mat_to_gfmmat(m,md,&mm); |
| row = m->row; |
row = m->row; |
|
|
| MKLIST(*rp,n0); |
MKLIST(*rp,n0); |
| } |
} |
| |
|
| void Pmat_to_gfmmat(arg,rp) |
void Pmat_to_gfmmat(NODE arg,GFMMAT *rp) |
| NODE arg; |
|
| GFMMAT *rp; |
|
| { |
{ |
| MAT m; |
MAT m; |
| unsigned int md; |
unsigned int md; |
|
|
| mat_to_gfmmat(m,md,rp); |
mat_to_gfmmat(m,md,rp); |
| } |
} |
| |
|
| void mat_to_gfmmat(m,md,rp) |
void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp) |
| MAT m; |
|
| unsigned int md; |
|
| GFMMAT *rp; |
|
| { |
{ |
| unsigned int **wmat; |
unsigned int **wmat; |
| unsigned int t; |
unsigned int t; |
|
|
| TOGFMMAT(row,col,wmat,*rp); |
TOGFMMAT(row,col,wmat,*rp); |
| } |
} |
| |
|
| void Pgeninvm_swap(arg,rp) |
void Pgeninvm_swap(NODE arg,LIST *rp) |
| NODE arg; |
|
| LIST *rp; |
|
| { |
{ |
| MAT m; |
MAT m; |
| pointer **mat; |
pointer **mat; |
|
|
| } |
} |
| } |
} |
| |
|
| gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp) |
int gauss_elim_geninv_mod_swap(unsigned int **mat,int row,int col,unsigned int md, |
| unsigned int **mat; |
unsigned int ***invmatp,int **indexp) |
| int row,col; |
|
| unsigned int md; |
|
| unsigned int ***invmatp; |
|
| int **indexp; |
|
| { |
{ |
| int i,j,k,inv,a,n,m; |
int i,j,k,inv,a,n,m; |
| unsigned int *t,*pivot,*s; |
unsigned int *t,*pivot,*s; |
|
|
| return 0; |
return 0; |
| } |
} |
| |
|
| |
void Pgeninv_sf_swap(NODE arg,LIST *rp) |
| |
{ |
| |
MAT m; |
| |
GFS **mat,**tmat; |
| |
Q *tvect; |
| |
GFS q; |
| |
int **wmat,**invmat; |
| |
int *index; |
| |
unsigned int t; |
| |
int i,j,row,col,status; |
| |
MAT mat1; |
| |
VECT vect1; |
| |
NODE node1,node2; |
| |
|
| |
asir_assert(ARG0(arg),O_MAT,"geninv_sf_swap"); |
| |
m = (MAT)ARG0(arg); |
| |
row = m->row; col = m->col; mat = (GFS **)m->body; |
| |
wmat = (int **)almat(row,col+row); |
| |
for ( i = 0; i < row; i++ ) { |
| |
bzero((char *)wmat[i],(col+row)*sizeof(int)); |
| |
for ( j = 0; j < col; j++ ) |
| |
if ( q = (GFS)mat[i][j] ) |
| |
wmat[i][j] = FTOIF(CONT(q)); |
| |
wmat[i][col+i] = _onesf(); |
| |
} |
| |
status = gauss_elim_geninv_sf_swap(wmat,row,col,&invmat,&index); |
| |
if ( status > 0 ) |
| |
*rp = 0; |
| |
else { |
| |
MKMAT(mat1,col,col); |
| |
for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ ) |
| |
for ( j = 0; j < col; j++ ) |
| |
if ( t = invmat[i][j] ) { |
| |
MKGFS(IFTOF(t),tmat[i][j]); |
| |
} |
| |
MKVECT(vect1,row); |
| |
for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ ) |
| |
STOQ(index[i],tvect[i]); |
| |
MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1); |
| |
} |
| |
} |
| |
|
| |
int gauss_elim_geninv_sf_swap(int **mat,int row,int col, |
| |
int ***invmatp,int **indexp) |
| |
{ |
| |
int i,j,k,inv,a,n,m,u; |
| |
int *t,*pivot,*s; |
| |
int *index; |
| |
int **invmat; |
| |
|
| |
n = col; m = row+col; |
| |
*indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int)); |
| |
for ( i = 0; i < row; i++ ) |
| |
index[i] = i; |
| |
for ( j = 0; j < n; j++ ) { |
| |
for ( i = j; i < row && !mat[i][j]; i++ ); |
| |
if ( i == row ) { |
| |
*indexp = 0; *invmatp = 0; return 1; |
| |
} |
| |
if ( i != j ) { |
| |
t = mat[i]; mat[i] = mat[j]; mat[j] = t; |
| |
k = index[i]; index[i] = index[j]; index[j] = k; |
| |
} |
| |
pivot = mat[j]; |
| |
inv = _invsf(pivot[j]); |
| |
for ( k = j; k < m; k++ ) |
| |
if ( pivot[k] ) |
| |
pivot[k] = _mulsf(pivot[k],inv); |
| |
for ( i = j+1; i < row; i++ ) { |
| |
t = mat[i]; |
| |
if ( a = t[j] ) |
| |
for ( k = j, a = _chsgnsf(a); k < m; k++ ) |
| |
if ( pivot[k] ) { |
| |
u = _mulsf(pivot[k],a); |
| |
t[k] = _addsf(u,t[k]); |
| |
} |
| |
} |
| |
} |
| |
for ( j = n-1; j >= 0; j-- ) { |
| |
pivot = mat[j]; |
| |
for ( i = j-1; i >= 0; i-- ) { |
| |
t = mat[i]; |
| |
if ( a = t[j] ) |
| |
for ( k = j, a = _chsgnsf(a); k < m; k++ ) |
| |
if ( pivot[k] ) { |
| |
u = _mulsf(pivot[k],a); |
| |
t[k] = _addsf(u,t[k]); |
| |
} |
| |
} |
| |
} |
| |
*invmatp = invmat = (int **)almat(col,col); |
| |
for ( i = 0; i < col; i++ ) |
| |
for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ ) |
| |
s[j] = t[col+index[j]]; |
| |
return 0; |
| |
} |
| |
|
| void _addn(N,N,N); |
void _addn(N,N,N); |
| int _subn(N,N,N); |
int _subn(N,N,N); |
| void _muln(N,N,N); |
void _muln(N,N,N); |
| |
|
| void inner_product_int(a,b,n,r) |
void inner_product_int(Q *a,Q *b,int n,Q *r) |
| Q *a,*b; |
|
| int n; |
|
| Q *r; |
|
| { |
{ |
| int la,lb,i; |
int la,lb,i; |
| int sgn,sgn1; |
int sgn,sgn1; |
|
|
| t = wma; wma = sum; sum = t; |
t = wma; wma = sum; sum = t; |
| } |
} |
| } |
} |
| GC_free(wm); |
GCFREE(wm); |
| GC_free(wma); |
GCFREE(wma); |
| if ( !sgn ) { |
if ( !sgn ) { |
| GC_free(sum); |
GCFREE(sum); |
| *r = 0; |
*r = 0; |
| } else |
} else |
| NTOQ(sum,sgn,*r); |
NTOQ(sum,sgn,*r); |
|
|
| |
|
| /* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */ |
/* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */ |
| |
|
| void inner_product_mat_int_mod(a,b,n,k,l,r) |
void inner_product_mat_int_mod(Q **a,int **b,int n,int k,int l,Q *r) |
| Q **a; |
|
| int **b; |
|
| int n,k,l; |
|
| Q *r; |
|
| { |
{ |
| int la,lb,i; |
int la,lb,i; |
| int sgn,sgn1; |
int sgn,sgn1; |
|
|
| t = wma; wma = sum; sum = t; |
t = wma; wma = sum; sum = t; |
| } |
} |
| } |
} |
| GC_free(wm); |
GCFREE(wm); |
| GC_free(wma); |
GCFREE(wma); |
| if ( !sgn ) { |
if ( !sgn ) { |
| GC_free(sum); |
GCFREE(sum); |
| *r = 0; |
*r = 0; |
| } else |
} else |
| NTOQ(sum,sgn,*r); |
NTOQ(sum,sgn,*r); |
| } |
} |
| |
|
| void Pmul_mat_vect_int(arg,rp) |
void Pmul_mat_vect_int(NODE arg,VECT *rp) |
| NODE arg; |
|
| VECT *rp; |
|
| { |
{ |
| MAT mat; |
MAT mat; |
| VECT vect,r; |
VECT vect,r; |
|
|
| row = mat->row; |
row = mat->row; |
| col = mat->col; |
col = mat->col; |
| MKVECT(r,row); |
MKVECT(r,row); |
| for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) { |
| inner_product_int(mat->body[i],vect->body,col,&r->body[i]); |
inner_product_int((Q *)mat->body[i],(Q *)vect->body,col,(Q *)&r->body[i]); |
| |
} |
| *rp = r; |
*rp = r; |
| } |
} |
| |
|
| void Pnbpoly_up2(arg,rp) |
void Pnbpoly_up2(NODE arg,GF2N *rp) |
| NODE arg; |
|
| GF2N *rp; |
|
| { |
{ |
| int m,type,ret; |
int m,type,ret; |
| UP2 r; |
UP2 r; |
|
|
| *rp = 0; |
*rp = 0; |
| } |
} |
| |
|
| void Px962_irredpoly_up2(arg,rp) |
void Px962_irredpoly_up2(NODE arg,GF2N *rp) |
| NODE arg; |
|
| GF2N *rp; |
|
| { |
{ |
| int m,type,ret,w; |
int m,ret,w; |
| GF2N prev; |
GF2N prev; |
| UP2 r; |
UP2 r; |
| |
|
|
|
| bzero((char *)r->b,w*sizeof(unsigned int)); |
bzero((char *)r->b,w*sizeof(unsigned int)); |
| } |
} |
| } |
} |
| ret = _generate_irreducible_polynomial(r,m,type); |
ret = _generate_irreducible_polynomial(r,m); |
| if ( ret == 0 ) |
if ( ret == 0 ) |
| MKGF2N(r,*rp); |
MKGF2N(r,*rp); |
| else |
else |
| *rp = 0; |
*rp = 0; |
| } |
} |
| |
|
| void Pirredpoly_up2(arg,rp) |
void Pirredpoly_up2(NODE arg,GF2N *rp) |
| NODE arg; |
|
| GF2N *rp; |
|
| { |
{ |
| int m,type,ret,w; |
int m,ret,w; |
| GF2N prev; |
GF2N prev; |
| UP2 r; |
UP2 r; |
| |
|
|
|
| bzero((char *)r->b,w*sizeof(unsigned int)); |
bzero((char *)r->b,w*sizeof(unsigned int)); |
| } |
} |
| } |
} |
| ret = _generate_good_irreducible_polynomial(r,m,type); |
ret = _generate_good_irreducible_polynomial(r,m); |
| if ( ret == 0 ) |
if ( ret == 0 ) |
| MKGF2N(r,*rp); |
MKGF2N(r,*rp); |
| else |
else |
| *rp = 0; |
*rp = 0; |
| } |
} |
| |
|
| |
void Pmat_swap_row_destructive(NODE arg, MAT *m) |
| |
{ |
| |
int i1,i2; |
| |
pointer *t; |
| |
MAT mat; |
| |
|
| |
asir_assert(ARG0(arg),O_MAT,"mat_swap_row_destructive"); |
| |
asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive"); |
| |
asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive"); |
| |
mat = (MAT)ARG0(arg); |
| |
i1 = QTOS((Q)ARG1(arg)); |
| |
i2 = QTOS((Q)ARG2(arg)); |
| |
if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row ) |
| |
error("mat_swap_row_destructive : Out of range"); |
| |
t = mat->body[i1]; |
| |
mat->body[i1] = mat->body[i2]; |
| |
mat->body[i2] = t; |
| |
*m = mat; |
| |
} |
| |
|
| |
void Pmat_swap_col_destructive(NODE arg, MAT *m) |
| |
{ |
| |
int j1,j2,i,n; |
| |
pointer *mi; |
| |
pointer t; |
| |
MAT mat; |
| |
|
| |
asir_assert(ARG0(arg),O_MAT,"mat_swap_col_destructive"); |
| |
asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive"); |
| |
asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive"); |
| |
mat = (MAT)ARG0(arg); |
| |
j1 = QTOS((Q)ARG1(arg)); |
| |
j2 = QTOS((Q)ARG2(arg)); |
| |
if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col ) |
| |
error("mat_swap_col_destructive : Out of range"); |
| |
n = mat->row; |
| |
for ( i = 0; i < n; i++ ) { |
| |
mi = mat->body[i]; |
| |
t = mi[j1]; mi[j1] = mi[j2]; mi[j2] = t; |
| |
} |
| |
*m = mat; |
| |
} |
| /* |
/* |
| * f = type 'type' normal polynomial of degree m if exists |
* f = type 'type' normal polynomial of degree m if exists |
| * IEEE P1363 A.7.2 |
* IEEE P1363 A.7.2 |
|
|
| return 1; |
return 1; |
| } |
} |
| |
|
| printqmat(mat,row,col) |
void printqmat(Q **mat,int row,int col) |
| Q **mat; |
|
| int row,col; |
|
| { |
{ |
| int i,j; |
int i,j; |
| |
|
| for ( i = 0; i < row; i++ ) { |
for ( i = 0; i < row; i++ ) { |
| for ( j = 0; j < col; j++ ) { |
for ( j = 0; j < col; j++ ) { |
| printnum(mat[i][j]); printf(" "); |
printnum((Num)mat[i][j]); printf(" "); |
| } |
} |
| printf("\n"); |
printf("\n"); |
| } |
} |
| } |
} |
| |
|
| printimat(mat,row,col) |
void printimat(int **mat,int row,int col) |
| int **mat; |
|
| int row,col; |
|
| { |
{ |
| int i,j; |
int i,j; |
| |
|
|
|
| } |
} |
| printf("\n"); |
printf("\n"); |
| } |
} |
| |
} |
| |
|
| |
void Pnd_det(NODE arg,P *rp) |
| |
{ |
| |
if ( argc(arg) == 1 ) |
| |
nd_det(0,ARG0(arg),rp); |
| |
else |
| |
nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp); |
| |
} |
| |
|
| |
void Pmat_col(NODE arg,VECT *rp) |
| |
{ |
| |
int i,j,n; |
| |
MAT mat; |
| |
VECT vect; |
| |
|
| |
asir_assert(ARG0(arg),O_MAT,"mat_col"); |
| |
asir_assert(ARG1(arg),O_N,"mat_col"); |
| |
mat = (MAT)ARG0(arg); |
| |
j = QTOS((Q)ARG1(arg)); |
| |
if ( j < 0 || j >= mat->col) { |
| |
error("mat_col : Out of range"); |
| |
} |
| |
n = mat->row; |
| |
MKVECT(vect,n); |
| |
for(i=0; i<n; i++) { |
| |
BDY(vect)[i] = BDY(mat)[i][j]; |
| |
} |
| |
*rp = vect; |
| |
} |
| |
|
| |
NODE triangleq(NODE e) |
| |
{ |
| |
int n,i,k; |
| |
V v; |
| |
VL vl; |
| |
P *p; |
| |
NODE r,r1; |
| |
|
| |
n = length(e); |
| |
p = (P *)MALLOC(n*sizeof(P)); |
| |
for ( i = 0; i < n; i++, e = NEXT(e) ) p[i] = (P)BDY(e); |
| |
i = 0; |
| |
while ( 1 ) { |
| |
for ( ; i < n && !p[i]; i++ ); |
| |
if ( i == n ) break; |
| |
if ( OID(p[i]) == O_N ) return 0; |
| |
v = p[i]->v; |
| |
for ( k = i+1; k < n; k++ ) |
| |
if ( p[k] ) { |
| |
if ( OID(p[k]) == O_N ) return 0; |
| |
if ( p[k]->v == v ) p[k] = 0; |
| |
} |
| |
i++; |
| |
} |
| |
for ( r = 0, i = 0; i < n; i++ ) { |
| |
if ( p[i] ) { |
| |
MKNODE(r1,p[i],r); r = r1; |
| |
} |
| |
} |
| |
return r; |
| |
} |
| |
|
| |
void Ptriangleq(NODE arg,LIST *rp) |
| |
{ |
| |
NODE ret; |
| |
|
| |
asir_assert(ARG0(arg),O_LIST,"sparseleq"); |
| |
ret = triangleq(BDY((LIST)ARG0(arg))); |
| |
MKLIST(*rp,ret); |
| } |
} |
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to array.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / builtin