| version 1.2, 2000/03/14 05:25:43 |
version 1.51, 2006/03/16 10:08:20 |
|
|
| /* $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.1.1.1 1999/12/03 07:39:07 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.50 2006/01/05 00:21:20 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> |
| |
#include <unistd.h> |
| |
|
| |
#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 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 *); |
|
| int generic_gauss_elim(MAT ,MAT *,Q *,int **,int **); |
|
| |
|
| int gauss_elim_mod(int **,int,int,int); |
|
| 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 Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(); |
| |
void Pinvmat(); |
| |
void Pnewbytearray(),Pmemoryplot_to_coord(); |
| |
|
| |
void Pgeneric_gauss_elim(); |
| void Pgeneric_gauss_elim_mod(); |
void Pgeneric_gauss_elim_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 38 void Px962_irredpoly_up2(); |
|
| Line 87 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(); |
| |
|
| struct ftab array_tab[] = { |
struct ftab array_tab[] = { |
| {"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}, |
| {"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 66 struct ftab array_tab[] = { |
|
| Line 136 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}, |
| {0,0,0}, |
{0,0,0}, |
| }; |
}; |
| |
|
| int comp_obj(a,b) |
int comp_obj(Obj *a,Obj *b) |
| Obj *a,*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; |
| |
|
| int generic_comp_obj(a,b) |
int generic_comp_obj(Obj *a,Obj *b) |
| Obj *a,*b; |
|
| { |
{ |
| Q r; |
Q r; |
| |
|
|
|
| } |
} |
| |
|
| |
|
| 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 ) |
| |
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); |
MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n); |
| 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; |
|
|
| 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 ) { |
|
|
| *rp = vect; |
*rp = vect; |
| } |
} |
| |
|
| void Pnewmat(arg,rp) |
void Pvect(NODE arg,VECT *rp) { |
| NODE arg; |
int len,i,r; |
| 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 ) { |
| |
/* ARG0(arg) must be a filename */ |
| |
asir_assert(ARG0(arg),O_STR,"newbytearray"); |
| |
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); |
| |
} 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; |
|
|
| MKLIST(*rp,n); |
MKLIST(*rp,n); |
| } |
} |
| |
|
| void Premainder(arg,rp) |
void Pltov(NODE arg,VECT *rp) |
| NODE arg; |
|
| Obj *rp; |
|
| { |
{ |
| |
NODE n; |
| |
VECT v; |
| |
int len,i; |
| |
|
| |
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,j,k,l,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); |
| |
} |
| |
|
| |
/* |
| |
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 1073 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 1087 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 706 int **rindp,**cindp; |
|
| Line 1117 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); |
|
|
| else |
else |
| cind[l++] = j; |
cind[l++] = j; |
| get_eg(&tmp0); |
get_eg(&tmp0); |
| if ( gensolve_check(mat,*nm,*dn,rind,cind) ) |
if ( gensolve_check(mat,*nm,*dn,rind,cind) ) { |
| |
get_eg(&tmp1); |
| |
add_eg(&eg_gschk,&tmp0,&tmp1); |
| |
add_eg(&eg_gschk_split,&tmp0,&tmp1); |
| |
if ( DP_Print ) { |
| |
print_eg("Mod",&eg_mod_split); |
| |
print_eg("Elim",&eg_elim_split); |
| |
print_eg("ChRem",&eg_chrem_split); |
| |
print_eg("IntRat",&eg_intrat_split); |
| |
print_eg("Check",&eg_gschk_split); |
| |
fflush(asir_out); |
| |
} |
| |
return rank; |
| |
} |
| |
} |
| |
} |
| |
} |
| |
} |
| |
|
| |
int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp) |
| |
{ |
| |
MAT bmat,xmat; |
| |
Q **a0,**a,**b,**x,**nm; |
| |
Q *ai,*bi,*xi; |
| |
int row,col; |
| |
int **w; |
| |
int *wi; |
| |
int **wc; |
| |
Q mdq,q,s,u; |
| |
N tn; |
| |
int ind,md,i,j,k,l,li,ri,rank; |
| |
unsigned int t; |
| |
int *cinfo,*rinfo; |
| |
int *rind,*cind; |
| |
int count; |
| |
int ret; |
| |
struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1; |
| |
int period; |
| |
int *wx,*ptr; |
| |
int wxsize,nsize; |
| |
N wn; |
| |
Q wq; |
| |
|
| |
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); |
| |
} |
| |
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] ) { |
| |
/* 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 { |
| |
/* 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 */ |
| |
/* ri = col - rank */ |
| |
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((ri)*sizeof(int)); |
| |
|
| |
init_eg(&eg_mul); init_eg(&eg_inv); |
| |
init_eg(&eg_check); init_eg(&eg_intrat); |
| |
period = F4_INTRAT_PERIOD; |
| |
nsize = period; |
| |
wxsize = rank*ri*nsize; |
| |
wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int)); |
| |
for ( i = 0; i < wxsize; i++ ) wx[i] = 0; |
| |
for ( q = ONE, count = 0; ; ) { |
| |
if ( DP_Print ) |
| |
fprintf(stderr,"o"); |
| |
/* wc = -b mod md */ |
| |
get_eg(&tmp0); |
| |
for ( i = 0; i < rank; i++ ) |
| |
for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) |
| |
if ( u = (Q)bi[j] ) { |
| |
t = rem(NM(u),md); |
| |
if ( t && SGN(u) > 0 ) |
| |
t = (md - t) % md; |
| |
wi[j] = t; |
| |
} else |
| |
wi[j] = 0; |
| |
/* wc = A^(-1)wc; wc is not normalized */ |
| |
solve_by_lu_mod(w,rank,md,wc,ri,0); |
| |
/* wx += q*wc */ |
| |
ptr = wx; |
| |
for ( i = 0; i < rank; i++ ) |
| |
for ( j = 0, wi = wc[i]; j < ri; j++ ) { |
| |
if ( wi[j] ) |
| |
muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr); |
| |
ptr += nsize; |
| |
} |
| |
count++; |
| get_eg(&tmp1); |
get_eg(&tmp1); |
| add_eg(&eg_gschk,&tmp0,&tmp1); |
add_eg(&eg_inv,&tmp0,&tmp1); |
| add_eg(&eg_gschk_split,&tmp0,&tmp1); |
get_eg(&tmp0); |
| if ( Print ) { |
for ( i = 0; i < rank; i++ ) |
| print_eg("Mod",&eg_mod_split); |
for ( j = 0; j < ri; j++ ) { |
| print_eg("Elim",&eg_elim_split); |
inner_product_mat_int_mod(a,wc,rank,i,j,&u); |
| print_eg("ChRem",&eg_chrem_split); |
addq(b[i][j],u,&s); |
| print_eg("IntRat",&eg_intrat_split); |
if ( s ) { |
| print_eg("Check",&eg_gschk_split); |
t = divin(NM(s),md,&tn); |
| fflush(asir_out); |
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; |
| |
} |
| } |
} |
| return rank; |
|
| } |
} |
| |
} |
| |
} |
| |
|
| |
int generic_gauss_elim_hensel_dalg(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp) |
| |
{ |
| |
MAT bmat,xmat; |
| |
Q **a0,**a,**b,**x,**nm; |
| |
Q *ai,*bi,*xi; |
| |
int row,col; |
| |
int **w; |
| |
int *wi; |
| |
int **wc; |
| |
Q mdq,q,s,u; |
| |
N tn; |
| |
int ind,md,i,j,k,l,li,ri,rank; |
| |
unsigned int t; |
| |
int *cinfo,*rinfo; |
| |
int *rind,*cind; |
| |
int count; |
| |
int ret; |
| |
struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1; |
| |
int period; |
| |
int *wx,*ptr; |
| |
int wxsize,nsize; |
| |
N wn; |
| |
Q wq; |
| |
NumberField nf; |
| |
DP *mb; |
| |
DP m; |
| |
int col1; |
| |
|
| |
nf = get_numberfield(); |
| |
mb = nf->mb; |
| |
a0 = (Q **)mat->body; |
| |
row = mat->row; col = mat->col; |
| |
w = (int **)almat(row,col); |
| |
for ( ind = 0; ; ind++ ) { |
| |
md = get_lprime(ind); |
| |
STOQ(md,mdq); |
| |
for ( i = 0; i < row; i++ ) |
| |
for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ ) |
| |
if ( q = (Q)ai[j] ) { |
| |
t = rem(NM(q),md); |
| |
if ( t && SGN(q) < 0 ) |
| |
t = (md - t) % md; |
| |
wi[j] = t; |
| |
} else |
| |
wi[j] = 0; |
| |
|
| |
if ( DP_Print ) { |
| |
fprintf(asir_out,"LU decomposition.."); fflush(asir_out); |
| } |
} |
| |
rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo); |
| |
if ( DP_Print ) { |
| |
fprintf(asir_out,"done.\n"); fflush(asir_out); |
| |
} |
| |
for ( i = 0; i < col-1; i++ ) { |
| |
if ( !cinfo[i] ) { |
| |
m = mb[i]; |
| |
for ( j = i+1; j < col-1; j++ ) |
| |
if ( dp_redble(mb[j],m) ) |
| |
cinfo[j] = -1; |
| |
} |
| |
} |
| |
a = (Q **)almat_pointer(rank,rank); /* lhs mat */ |
| |
MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */ |
| |
for ( j = li = ri = 0; j < col; j++ ) |
| |
if ( cinfo[j] > 0 ) { |
| |
/* the column is in lhs */ |
| |
for ( i = 0; i < rank; i++ ) { |
| |
w[i][li] = w[i][j]; |
| |
a[i][li] = a0[rinfo[i]][j]; |
| |
} |
| |
li++; |
| |
} else if ( !cinfo[j] ) { |
| |
/* the column is in rhs */ |
| |
for ( i = 0; i < rank; i++ ) |
| |
b[i][ri] = a0[rinfo[i]][j]; |
| |
ri++; |
| |
} |
| |
|
| |
/* solve Ax+B=0; A: rank x rank, B: rank x ri */ |
| |
MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body; |
| |
MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body; |
| |
/* use the right part of w as work area */ |
| |
wc = (int **)almat(rank,ri); |
| |
for ( i = 0; i < rank; i++ ) |
| |
wc[i] = w[i]+rank; |
| |
*rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); |
| |
*cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int)); |
| |
init_eg(&eg_mul); init_eg(&eg_inv); |
| |
init_eg(&eg_check); init_eg(&eg_intrat); |
| |
period = F4_INTRAT_PERIOD; |
| |
nsize = period; |
| |
wxsize = rank*ri*nsize; |
| |
wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int)); |
| |
for ( i = 0; i < wxsize; i++ ) wx[i] = 0; |
| |
for ( q = ONE, count = 0; ; ) { |
| |
if ( DP_Print ) |
| |
fprintf(stderr,"o"); |
| |
/* wc = -b mod md */ |
| |
get_eg(&tmp0); |
| |
for ( i = 0; i < rank; i++ ) |
| |
for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) |
| |
if ( u = (Q)bi[j] ) { |
| |
t = rem(NM(u),md); |
| |
if ( t && SGN(u) > 0 ) |
| |
t = (md - t) % md; |
| |
wi[j] = t; |
| |
} else |
| |
wi[j] = 0; |
| |
/* wc = A^(-1)wc; wc is not normalized */ |
| |
solve_by_lu_mod(w,rank,md,wc,ri,0); |
| |
/* wx += q*wc */ |
| |
ptr = wx; |
| |
for ( i = 0; i < rank; i++ ) |
| |
for ( j = 0, wi = wc[i]; j < ri; j++ ) { |
| |
if ( wi[j] ) |
| |
muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr); |
| |
ptr += nsize; |
| |
} |
| |
count++; |
| |
get_eg(&tmp1); |
| |
add_eg(&eg_inv,&tmp0,&tmp1); |
| |
get_eg(&tmp0); |
| |
for ( i = 0; i < rank; i++ ) |
| |
for ( j = 0; j < ri; j++ ) { |
| |
inner_product_mat_int_mod(a,wc,rank,i,j,&u); |
| |
addq(b[i][j],u,&s); |
| |
if ( s ) { |
| |
t = divin(NM(s),md,&tn); |
| |
if ( t ) |
| |
error("generic_gauss_elim_hensel:incosistent"); |
| |
NTOQ(tn,SGN(s),b[i][j]); |
| |
} else |
| |
b[i][j] = 0; |
| |
} |
| |
get_eg(&tmp1); |
| |
add_eg(&eg_mul,&tmp0,&tmp1); |
| |
/* q = q*md */ |
| |
mulq(q,mdq,&u); q = u; |
| |
if ( count == period ) { |
| |
get_eg(&tmp0); |
| |
ptr = wx; |
| |
for ( i = 0; i < rank; i++ ) |
| |
for ( j = 0, xi = x[i]; j < ri; |
| |
j++, ptr += nsize ) { |
| |
for ( k = nsize-1; k >= 0 && !ptr[k]; k-- ); |
| |
if ( k >= 0 ) { |
| |
wn = NALLOC(k+1); |
| |
PL(wn) = k+1; |
| |
for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l]; |
| |
NTOQ(wn,1,wq); |
| |
subq(xi[j],wq,&u); xi[j] = u; |
| |
} |
| |
} |
| |
ret = intmtoratm_q(xmat,NM(q),*nmmat,dn); |
| |
get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1); |
| |
if ( ret ) { |
| |
for ( j = k = l = 0; j < col; j++ ) |
| |
if ( cinfo[j] > 0 ) |
| |
rind[k++] = j; |
| |
else if ( !cinfo[j] ) |
| |
cind[l++] = j; |
| |
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; |
| |
} |
| |
} |
| |
} |
| } |
} |
| } |
} |
| |
|
| 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 880 int *rind,*cind; |
|
| Line 1663 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; |
|
|
| N u,unm,udn; |
N u,unm,udn; |
| int sgn,ret; |
int sgn,ret; |
| |
|
| |
if ( UNIN(md) ) |
| |
return 0; |
| row = mat->row; col = mat->col; |
row = mat->row; col = mat->col; |
| bshiftn(md,1,&t); |
bshiftn(md,1,&t); |
| isqrt(t,&s); |
isqrt(t,&s); |
|
|
| return 1; |
return 1; |
| } |
} |
| |
|
| int generic_gauss_elim_mod(mat,row,col,md,colstat) |
/* mat->body = Q ** */ |
| int **mat; |
|
| int row,col,md; |
int intmtoratm_q(MAT mat,N md,MAT nm,Q *dn) |
| int *colstat; |
|
| { |
{ |
| |
N t,s,b; |
| |
Q dn0,dn1,nm1,q; |
| |
int i,j,k,l,row,col; |
| |
Q **rmat; |
| |
Q **tmat; |
| |
Q *tmi; |
| |
Q *nmk; |
| |
N u,unm,udn; |
| |
int sgn,ret; |
| |
|
| |
if ( UNIN(md) ) |
| |
return 0; |
| |
row = mat->row; col = mat->col; |
| |
bshiftn(md,1,&t); |
| |
isqrt(t,&s); |
| |
bshiftn(s,64,&b); |
| |
if ( !b ) |
| |
b = ONEN; |
| |
dn0 = ONE; |
| |
tmat = (Q **)mat->body; |
| |
rmat = (Q **)nm->body; |
| |
for ( i = 0; i < row; i++ ) |
| |
for ( j = 0, tmi = tmat[i]; j < col; j++ ) |
| |
if ( tmi[j] ) { |
| |
muln(NM(tmi[j]),NM(dn0),&s); |
| |
remn(s,md,&u); |
| |
ret = inttorat(u,md,b,&sgn,&unm,&udn); |
| |
if ( !ret ) |
| |
return 0; |
| |
else { |
| |
if ( SGN(tmi[j])<0 ) |
| |
sgn = -sgn; |
| |
NTOQ(unm,sgn,nm1); |
| |
NTOQ(udn,1,dn1); |
| |
if ( !UNIQ(dn1) ) { |
| |
for ( k = 0; k < i; k++ ) |
| |
for ( l = 0, nmk = rmat[k]; l < col; l++ ) { |
| |
mulq(nmk[l],dn1,&q); nmk[l] = q; |
| |
} |
| |
for ( l = 0, nmk = rmat[i]; l < j; l++ ) { |
| |
mulq(nmk[l],dn1,&q); nmk[l] = q; |
| |
} |
| |
} |
| |
rmat[i][j] = nm1; |
| |
mulq(dn0,dn1,&q); dn0 = q; |
| |
} |
| |
} |
| |
*dn = dn0; |
| |
return 1; |
| |
} |
| |
|
| |
#define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++; |
| |
|
| |
void reduce_reducers_mod(int **mat,int row,int col,int md) |
| |
{ |
| |
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) |
| |
{ |
| |
register unsigned int up,lo; |
| |
unsigned int 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; |
| |
} |
| |
} |
| |
|
| |
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 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_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; |
| } |
} |
| } |
} |
| } |
} |
|
|
| return 1; |
return 1; |
| } |
} |
| |
|
| void Pleqm1(arg,rp) |
/* |
| NODE arg; |
Input |
| VECT *rp; |
a: a row x col matrix |
| |
md : a modulus |
| |
|
| |
Output: |
| |
return : d = the rank of mat |
| |
a[0..(d-1)][0..(d-1)] : LU decomposition (a[i][i] = 1/U[i][i]) |
| |
rinfo: array of length row |
| |
cinfo: array of length col |
| |
i-th row in new a <-> rinfo[i]-th row in old a |
| |
cinfo[j]=1 <=> j-th column is contained in the LU decomp. |
| |
*/ |
| |
|
| |
int find_lhs_and_lu_mod(unsigned int **a,int row,int col, |
| |
unsigned int md,int **rinfo,int **cinfo) |
| { |
{ |
| |
int i,j,k,d; |
| |
int *rp,*cp; |
| |
unsigned int *t,*pivot; |
| |
unsigned int inv,m; |
| |
|
| |
*rinfo = rp = (int *)MALLOC_ATOMIC(row*sizeof(int)); |
| |
*cinfo = cp = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
| |
for ( i = 0; i < row; i++ ) |
| |
rp[i] = i; |
| |
for ( k = 0, d = 0; k < col; k++ ) { |
| |
for ( i = d; i < row && !a[i][k]; i++ ); |
| |
if ( i == row ) { |
| |
cp[k] = 0; |
| |
continue; |
| |
} else |
| |
cp[k] = 1; |
| |
if ( i != d ) { |
| |
j = rp[i]; rp[i] = rp[d]; rp[d] = j; |
| |
t = a[i]; a[i] = a[d]; a[d] = t; |
| |
} |
| |
pivot = a[d]; |
| |
pivot[k] = inv = invm(pivot[k],md); |
| |
for ( i = d+1; i < row; i++ ) { |
| |
t = a[i]; |
| |
if ( m = t[k] ) { |
| |
DMAR(inv,m,0,md,t[k]) |
| |
for ( j = k+1, m = md - t[k]; j < col; j++ ) |
| |
if ( pivot[j] ) { |
| |
unsigned int tj; |
| |
DMAR(m,pivot[j],t[j],md,tj) |
| |
t[j] = tj; |
| |
} |
| |
} |
| |
} |
| |
d++; |
| |
} |
| |
return d; |
| |
} |
| |
|
| |
/* |
| |
Input |
| |
a : n x n matrix; a result of LU-decomposition |
| |
md : modulus |
| |
b : n x l matrix |
| |
Output |
| |
b = a^(-1)b |
| |
*/ |
| |
|
| |
void solve_by_lu_mod(int **a,int n,int md,int **b,int l,int normalize) |
| |
{ |
| |
unsigned int *y,*c; |
| |
int i,j,k; |
| |
unsigned int t,m,m2; |
| |
|
| |
y = (int *)MALLOC_ATOMIC(n*sizeof(int)); |
| |
c = (int *)MALLOC_ATOMIC(n*sizeof(int)); |
| |
m2 = md>>1; |
| |
for ( k = 0; k < l; k++ ) { |
| |
/* copy b[.][k] to c */ |
| |
for ( i = 0; i < n; i++ ) |
| |
c[i] = (unsigned int)b[i][k]; |
| |
/* solve Ly=c */ |
| |
for ( i = 0; i < n; i++ ) { |
| |
for ( t = c[i], j = 0; j < i; j++ ) |
| |
if ( a[i][j] ) { |
| |
m = md - a[i][j]; |
| |
DMAR(m,y[j],t,md,t) |
| |
} |
| |
y[i] = t; |
| |
} |
| |
/* solve Uc=y */ |
| |
for ( i = n-1; i >= 0; i-- ) { |
| |
for ( t = y[i], j =i+1; j < n; j++ ) |
| |
if ( a[i][j] ) { |
| |
m = md - a[i][j]; |
| |
DMAR(m,c[j],t,md,t) |
| |
} |
| |
/* a[i][i] = 1/U[i][i] */ |
| |
DMAR(t,a[i][i],0,md,c[i]) |
| |
} |
| |
/* copy c to b[.][k] with normalization */ |
| |
if ( normalize ) |
| |
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(NODE arg,VECT *rp) |
| |
{ |
| MAT m; |
MAT m; |
| VECT vect; |
VECT vect; |
| pointer **mat; |
pointer **mat; |
|
|
| } |
} |
| } |
} |
| |
|
| 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 1131 int row,col,md; |
|
| Line 2326 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 1213 int row,col,md; |
|
| Line 2404 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 1276 unsigned int *x; |
|
| Line 2462 unsigned int *x; |
|
| } |
} |
| } |
} |
| |
|
| void Plu_gfmmat(arg,rp) |
void Plu_gfmmat(NODE arg,LIST *rp) |
| NODE arg; |
|
| LIST *rp; |
|
| { |
{ |
| MAT m; |
MAT m; |
| GFMMAT mm; |
GFMMAT mm; |
|
|
| 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; |
|
|
| 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; |
|
|
| NTOQ(sum,sgn,*r); |
NTOQ(sum,sgn,*r); |
| } |
} |
| |
|
| void Pmul_mat_vect_int(arg,rp) |
/* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */ |
| NODE arg; |
|
| VECT *rp; |
void inner_product_mat_int_mod(Q **a,int **b,int n,int k,int l,Q *r) |
| { |
{ |
| |
int la,lb,i; |
| |
int sgn,sgn1; |
| |
N wm,wma,sum,t; |
| |
Q aki; |
| |
int bil,bilsgn; |
| |
struct oN tn; |
| |
|
| |
for ( la = 0, i = 0; i < n; i++ ) { |
| |
if ( aki = a[k][i] ) |
| |
if ( DN(aki) ) |
| |
error("inner_product_int : invalid argument"); |
| |
else |
| |
la = MAX(PL(NM(aki)),la); |
| |
} |
| |
lb = 1; |
| |
sgn = 0; |
| |
sum= NALLOC(la+lb+2); |
| |
bzero((char *)sum,(la+lb+3)*sizeof(unsigned int)); |
| |
wm = NALLOC(la+lb+2); |
| |
wma = NALLOC(la+lb+2); |
| |
for ( i = 0; i < n; i++ ) { |
| |
if ( !(aki = a[k][i]) || !(bil = b[i][l]) ) |
| |
continue; |
| |
tn.p = 1; |
| |
if ( bil > 0 ) { |
| |
tn.b[0] = bil; bilsgn = 1; |
| |
} else { |
| |
tn.b[0] = -bil; bilsgn = -1; |
| |
} |
| |
_muln(NM(aki),&tn,wm); |
| |
sgn1 = SGN(aki)*bilsgn; |
| |
if ( !sgn ) { |
| |
sgn = sgn1; |
| |
t = wm; wm = sum; sum = t; |
| |
} else if ( sgn == sgn1 ) { |
| |
_addn(sum,wm,wma); |
| |
if ( !PL(wma) ) |
| |
sgn = 0; |
| |
t = wma; wma = sum; sum = t; |
| |
} else { |
| |
/* sgn*sum+sgn1*wm = sgn*(sum-wm) */ |
| |
sgn *= _subn(sum,wm,wma); |
| |
t = wma; wma = sum; sum = t; |
| |
} |
| |
} |
| |
GC_free(wm); |
| |
GC_free(wma); |
| |
if ( !sgn ) { |
| |
GC_free(sum); |
| |
*r = 0; |
| |
} else |
| |
NTOQ(sum,sgn,*r); |
| |
} |
| |
|
| |
void Pmul_mat_vect_int(NODE arg,VECT *rp) |
| |
{ |
| MAT mat; |
MAT mat; |
| VECT vect,r; |
VECT vect,r; |
| int row,col,i; |
int row,col,i; |
|
|
| 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 |
|
|
| } |
} |
| /* exhausted */ |
/* exhausted */ |
| return 1; |
return 1; |
| |
} |
| |
|
| |
void printqmat(Q **mat,int row,int col) |
| |
{ |
| |
int i,j; |
| |
|
| |
for ( i = 0; i < row; i++ ) { |
| |
for ( j = 0; j < col; j++ ) { |
| |
printnum((Num)mat[i][j]); printf(" "); |
| |
} |
| |
printf("\n"); |
| |
} |
| |
} |
| |
|
| |
void printimat(int **mat,int row,int col) |
| |
{ |
| |
int i,j; |
| |
|
| |
for ( i = 0; i < row; i++ ) { |
| |
for ( j = 0; j < col; j++ ) { |
| |
printf("%d ",mat[i][j]); |
| |
} |
| |
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); |
| } |
} |
![[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