OpenXM_contrib2/asir2018/builtin/array.c - diff

Return to array.c CVS log

Up to [local] / OpenXM_contrib2 / asir2018 / builtin

Diff for /OpenXM_contrib2/asir2018/builtin/array.c between version 1.1 and 1.8

version 1.1, 2018/09/19 05:45:05

version 1.8, 2020/10/06 06:31:19

Line 45

* DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,

* PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.

* $OpenXM$

* $OpenXM: OpenXM_contrib2/asir2018/builtin/array.c,v 1.7 2020/01/09 01:47:40 noro Exp $

#include "ca.h"

#include "base.h"

Line 170 void solve_u(int *,ent **,int,int *,int);

static int *ul,*ll;

static ent **u,**l;

static int modulus;

#if defined(ANDROID)

int getw(FILE *fp)

{

int x;

return (fread((void *)&x, sizeof(x), 1, fp) == 1 ? x : EOF);

}

#endif

void Plusolve_prep(NODE arg,Q *rp)

{

Line 206 void Plusolve_main(NODE arg,VECT *rp)

Line 213 void Plusolve_main(NODE arg,VECT *rp)

v = (VECT)ARG0(arg); len = v->len;

d = (Z *)BDY(v);

rhs = (int *)MALLOC_ATOMIC(len*sizeof(int));

for ( i = 0; i < len; i++ ) rhs[i] = QTOS(d[i]);

for ( i = 0; i < len; i++ ) rhs[i] = ZTOS(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 = (Z *)r->body;

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

STOQ(rhs[i],p[i]);

STOZ(rhs[i],p[i]);

*rp = r;

}

Line 565 void Psepmat_destructive(NODE arg,LIST *rp)

Line 572 void Psepmat_destructive(NODE arg,LIST *rp)

void Psepvect(NODE arg,VECT *rp)

{

sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp);

sepvect((VECT)ARG0(arg),ZTOS((Q)ARG1(arg)),rp);

}

void sepvect(VECT v,int d,VECT *rp)

Line 601 void Pnewvect(NODE arg,VECT *rp)

Line 608 void Pnewvect(NODE arg,VECT *rp)

NODE tn;

asir_assert(ARG0(arg),O_N,"newvect");

len = QTOS((Q)ARG0(arg));

len = ZTOS((Q)ARG0(arg));

if ( len < 0 )

error("newvect : invalid size");

MKVECT(vect,len);

Line 691 void Pnewbytearray(NODE arg,BYTEARRAY *rp)

Line 698 void Pnewbytearray(NODE arg,BYTEARRAY *rp)

case O_N:

if ( !RATN(ARG0(arg)) )

error("newbytearray : invalid argument");

len = QTOS((Q)ARG0(arg));

len = ZTOS((Q)ARG0(arg));

if ( len < 0 )

error("newbytearray : invalid size");

MKBYTEARRAY(array,len);

Line 701 void Pnewbytearray(NODE arg,BYTEARRAY *rp)

Line 708 void Pnewbytearray(NODE arg,BYTEARRAY *rp)

}

} else if ( ac == 2 ) {

asir_assert(ARG0(arg),O_N,"newbytearray");

len = QTOS((Q)ARG0(arg));

len = ZTOS((Q)ARG0(arg));

if ( len < 0 )

error("newbytearray : invalid size");

MKBYTEARRAY(array,len);

Line 715 void Pnewbytearray(NODE arg,BYTEARRAY *rp)

Line 722 void Pnewbytearray(NODE arg,BYTEARRAY *rp)

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));

vb[i] = (unsigned char)ZTOS((Q)BDY(tn));

}

break;

case O_STR:

Line 746 void Pmemoryplot_to_coord(NODE arg,LIST *rp)

Line 753 void Pmemoryplot_to_coord(NODE arg,LIST *rp)

asir_assert(ARG0(arg),O_LIST,"memoryplot_to_coord");

arg = BDY((LIST)ARG0(arg));

len = QTOS((Q)ARG0(arg));

len = ZTOS((Q)ARG0(arg));

blen = (len+7)/8;

y = QTOS((Q)ARG1(arg));

y = ZTOS((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);

STOZ(i,iq); STOZ(j,jq);

n = mknode(2,iq,jq);

MKLIST(l,n);

BDY(r) = l;

Line 775 void Pnewmat(NODE arg,MAT *rp)

Line 782 void Pnewmat(NODE arg,MAT *rp)

asir_assert(ARG0(arg),O_N,"newmat");

asir_assert(ARG1(arg),O_N,"newmat");

row = QTOS((Q)ARG0(arg)); col = QTOS((Q)ARG1(arg));

row = ZTOS((Q)ARG0(arg)); col = ZTOS((Q)ARG1(arg));

if ( row < 0 || col < 0 )

error("newmat : invalid size");

MKMAT(m,row,col);

Line 967 void Premainder(NODE arg,Obj *rp)

Line 974 void Premainder(NODE arg,Obj *rp)

cmp(md,(P)a,(P *)rp); break;

case O_VECT:

/* strage spec */

smd = QTOS(md);

smd = ZTOS(md);

v = (VECT)a; n = v->len; vb = v->body;

MKVECT(w,n); wb = w->body;

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

Line 1015 void Psremainder(NODE arg,Obj *rp)

Line 1022 void Psremainder(NODE arg,Obj *rp)

case O_P:

cmp(md,(P)a,(P *)rp); break;

case O_VECT:

smd = QTOS(md);

smd = ZTOS(md);

v = (VECT)a; n = v->len; vb = v->body;

MKVECT(w,n); wb = w->body;

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

Line 1049 void Psize(NODE arg,LIST *rp)

Line 1056 void Psize(NODE arg,LIST *rp)

switch (OID(ARG0(arg))) {

case O_VECT:

n = ((VECT)ARG0(arg))->len;

STOQ(n,q); MKNODE(t,q,0);

STOZ(n,q); MKNODE(t,q,0);

break;

case O_MAT:

n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col;

STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);

STOZ(m,q); MKNODE(s,q,0); STOZ(n,q); MKNODE(t,q,s);

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);

STOZ(m,q); MKNODE(s,q,0); STOZ(n,q); MKNODE(t,q,s);

break;

default:

error("size : invalid argument"); break;

Line 1080 void Pdet(NODE arg,P *rp)

Line 1087 void Pdet(NODE arg,P *rp)

else if ( argc(arg) == 1 )

detp(CO,(P **)BDY(m),m->row,rp);

else {

n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);

n = m->row; mod = ZTOS((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++ )

Line 1109 void Pinvmat(NODE arg,LIST *rp)

Line 1116 void Pinvmat(NODE arg,LIST *rp)

nd = mknode(2,r,dn);

MKLIST(*rp,nd);

} else {

n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);

n = m->row; mod = ZTOS((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++ )

Line 1169 void Pgeneric_gauss_elim(NODE arg,LIST *rp)

Line 1176 void Pgeneric_gauss_elim(NODE arg,LIST *rp)

MKVECT(rind,rank);

MKVECT(cind,t);

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

STOQ(ri[i],q);

STOZ(ri[i],q);

BDY(rind)[i] = (pointer)q;

}

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

STOQ(ci[i],q);

STOZ(ci[i],q);

BDY(cind)[i] = (pointer)q;

}

n0 = mknode(4,nm,dn,rind,cind);

Line 1195 void Pindep_rows_mod(NODE arg,VECT *rp)

Line 1202 void Pindep_rows_mod(NODE arg,VECT *rp)

asir_assert(ARG0(arg),O_MAT,"indep_rows_mod");

asir_assert(ARG1(arg),O_N,"indep_rows_mod");

m = (MAT)ARG0(arg); md = QTOS((Z)ARG1(arg));

m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));

row = m->row; col = m->col; tmat = (Z **)m->body;

wmat = (int **)almat(row,col);

Line 1211 void Pindep_rows_mod(NODE arg,VECT *rp)

Line 1218 void Pindep_rows_mod(NODE arg,VECT *rp)

MKVECT(rind,rank);

rib = (Z *)rind->body;

for ( j = 0; j < rank; j++ ) {

STOQ(rowstat[j],rib[j]);

STOZ(rowstat[j],rib[j]);

}

*rp = rind;

}

Line 1229 void Pindep_rows_mod(NODE arg,VECT *rp)

Line 1236 void Pindep_rows_mod(NODE arg,VECT *rp)

B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...

#if SIZEOF_LONG==8

void Pgeneric_gauss_elim_mod64(NODE arg,LIST *rp)

{

NODE n0;

MAT m,mat;

VECT rind,cind,rnum;

Z **tmat;

mp_limb_t **wmat,**row0;

Z *rib,*cib,*rnb;

int *colstat;

mp_limb_t *p;

Z q;

mp_limb_t md;

int i,j,k,l,row,col,t,rank;

asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod64");

asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod64");

m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));

row = m->row; col = m->col; tmat = (Z **)m->body;

wmat = (mp_limb_t **)almat64(row,col);

row0 = (mp_limb_t **)ALLOCA(row*sizeof(mp_limb_t *));

for ( i = 0; i < row; i++ ) row0[i] = wmat[i];

colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));

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

for ( j = 0; j < col; j++ )

wmat[i][j] = remqi64((Q)tmat[i][j],md);

rank = generic_gauss_elim_mod64(wmat,row,col,md,colstat);

MKVECT(rnum,rank);

rnb = (Z *)rnum->body;

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

for ( j = 0, p = wmat[i]; j < row; j++ )

if ( p == row0[j] )

STOZ(j,rnb[i]);

MKMAT(mat,rank,col-rank);

tmat = (Z **)mat->body;

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

for ( j = k = 0; j < col; j++ )

if ( !colstat[j] ) {

UTOZ(wmat[i][j],tmat[i][k]); k++;

}

MKVECT(rind,rank);

MKVECT(cind,col-rank);

rib = (Z *)rind->body; cib = (Z *)cind->body;

for ( j = k = l = 0; j < col; j++ )

if ( colstat[j] ) {

STOZ(j,rib[k]); k++;

} else {

STOZ(j,cib[l]); l++;

}

n0 = mknode(4,mat,rind,cind,rnum);

MKLIST(*rp,n0);

}

#endif

void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)

{

NODE n0;

Line 1239 void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)

Line 1305 void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)

Z *rib,*cib,*rnb;

int *colstat,*p;

Z q;

long mdl;

int md,i,j,k,l,row,col,t,rank;

asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");

asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");

m = (MAT)ARG0(arg); md = QTOS((Z)ARG1(arg));

#if SIZEOF_LONG==8

mdl = ZTOS((Z)ARG1(arg));

if ( mdl >= ((mp_limb_t)1)<<32 ) {

Pgeneric_gauss_elim_mod64(arg,rp);

return;

}

#endif

m = (MAT)ARG0(arg);

md = ZTOS((Z)ARG1(arg));

row = m->row; col = m->col; tmat = (Z **)m->body;

wmat = (int **)almat(row,col);

Line 1261 void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)

Line 1336 void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)

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

for ( j = 0, p = wmat[i]; j < row; j++ )

if ( p == row0[j] )

STOQ(j,rnb[i]);

STOZ(j,rnb[i]);

MKMAT(mat,rank,col-rank);

tmat = (Z **)mat->body;

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

for ( j = k = 0; j < col; j++ )

if ( !colstat[j] ) {

UTOQ(wmat[i][j],tmat[i][k]); k++;

UTOZ(wmat[i][j],tmat[i][k]); k++;

}

MKVECT(rind,rank);

Line 1276 void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)

Line 1351 void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)

rib = (Z *)rind->body; cib = (Z *)cind->body;

for ( j = k = l = 0; j < col; j++ )

if ( colstat[j] ) {

STOQ(j,rib[k]); k++;

STOZ(j,rib[k]); k++;

} else {

STOQ(j,cib[l]); l++;

STOZ(j,cib[l]); l++;

}

n0 = mknode(4,mat,rind,cind,rnum);

MKLIST(*rp,n0);

}

void Pleqm(NODE arg,VECT *rp)

{

MAT m;

Line 1296 void Pleqm(NODE arg,VECT *rp)

Line 1372 void Pleqm(NODE arg,VECT *rp)

asir_assert(ARG0(arg),O_MAT,"leqm");

asir_assert(ARG1(arg),O_N,"leqm");

m = (MAT)ARG0(arg); md = QTOS((Z)ARG1(arg));

m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));

row = m->row; col = m->col; mat = m->body;

wmat = (int **)almat(row,col);

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

Line 1311 void Pleqm(NODE arg,VECT *rp)

Line 1387 void Pleqm(NODE arg,VECT *rp)

n = col - 1;

MKVECT(vect,n);

for ( i = 0, v = (Z *)vect->body; i < n; i++ ) {

t = (md-wmat[i][n])%md; STOQ(t,v[i]);

t = (md-wmat[i][n])%md; STOZ(t,v[i]);

}

*rp = vect;

}

Line 1354 int gauss_elim_mod(int **mat,int row,int col,int md)

Line 1430 int gauss_elim_mod(int **mat,int row,int col,int md)

return -1;

}

struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,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,eg_back,eg_fore;

struct oEGT eg_conv;

#if 0

Line 1467 void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)

Line 1543 void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)

int f4_nocheck;

#define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;

#define ONE_STEP1 if ( ( zzz = *s ) != 0 ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;

void reduce_reducers_mod(int **mat,int row,int col,int md)

{

Line 1483 void reduce_reducers_mod(int **mat,int row,int col,int

Line 1559 void reduce_reducers_mod(int **mat,int row,int col,int

ind[i] = j;

for ( l = i-1; l >= 0; l-- ) {

/* reduce mat[i] by mat[l] */

if ( hc = t[ind[l]] ) {

if ( ( hc = t[ind[l]] ) != 0 ) {

/* mat[i] = mat[i]-hc*mat[l] */

j = ind[l];

s = mat[l]+j;

Line 1509 void reduce_reducers_mod(int **mat,int row,int col,int

Line 1585 void reduce_reducers_mod(int **mat,int row,int col,int

ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1

}

for ( ; k > 0; k-- ) {

if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;

if ( ( zzz = *s ) != 0 ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;

}

Line 1544 void pre_reduce_mod(int **mat,int row,int col,int nred

Line 1620 void pre_reduce_mod(int **mat,int row,int col,int nred

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]] ) {

if ( ( hc = t[ind[l]] ) != 0 ) {

/* 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++ )

Line 1558 void pre_reduce_mod(int **mat,int row,int col,int nred

Line 1634 void pre_reduce_mod(int **mat,int row,int col,int nred

t = mat[i];

for ( l = nred-1; l >= 0; l-- ) {

/* reduce mat[i] by mat[l] */

if ( hc = t[ind[l]] ) {

if ( ( hc = t[ind[l]] ) != 0 ) {

/* 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++ )

Line 1582 void reduce_sp_by_red_mod(int *sp,int **redmat,int *in

Line 1658 void reduce_sp_by_red_mod(int *sp,int **redmat,int *in

/* reduce the spolys by redmat */

for ( i = nred-1; i >= 0; i-- ) {

/* reduce sp by redmat[i] */

if ( hc = sp[ind[i]] ) {

if ( ( hc = sp[ind[i]] ) != 0 ) {

/* 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++;

if ( ( zzz = *s ) != 0 ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;

}

Line 1635 void red_by_vect(int m,unsigned int *p,unsigned int *r

Line 1711 void red_by_vect(int m,unsigned int *p,unsigned int *r

}

#if defined(__GNUC__) && SIZEOF_LONG==8

#if SIZEOF_LONG==8

/* 64bit vector += UNIT vector(normalized) */

/* mp_limb_t vector += U32 vector(normalized)*U32 */

void red_by_vect64(int m, U64 *p,unsigned int *c,U64 *r,unsigned int hc,int len)

void red_by_vect64(int m, mp_limb_t *p,unsigned int *c,mp_limb_t *r,unsigned int hc,int len)

{

U64 t;

mp_limb_t t;

/* (p[0],c[0]) is normalized */

*p++ = 0; *c++ = 0; r++; len--;

Line 1651 void red_by_vect64(int m, U64 *p,unsigned int *c,U64 *

Line 1727 void red_by_vect64(int m, U64 *p,unsigned int *c,U64 *

*p = t;

}

/* mp_limb_t vector = (mp_limb_t vector+mp_limb_t vector*mp_limb_t)%mp_limb_t */

void red_by_vect64mod(mp_limb_t m, mp_limb_t *p,mp_limb_t *r,mp_limb_t hc,int len)

{

*p++ = 0; r++; len--;

for ( ; len; len--, r++, p++ )

if ( *r )

*p = muladdmod64(*r,hc,*p,m);

}

int generic_gauss_elim_mod64(mp_limb_t **mat,int row,int col,mp_limb_t md,int *colstat)

{

int i,j,k,l,rank;

mp_limb_t inv,a;

mp_limb_t *t,*pivot,*pk;

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 = invmod64(pivot[j],md);

for ( k = j, pk = pivot+k; k < col; k++, pk++ )

if ( *pk )

*pk = mulmod64(*pk,inv,md);

for ( i = rank+1; i < row; i++ ) {

t = mat[i];

if ( ( a = t[j] ) != 0 )

red_by_vect64mod(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];

if ( ( a = t[j] ) != 0 )

red_by_vect64mod(md,t+j,pivot+j,md-a,col-j);

}

l--;

}

return rank;

}

int find_lhs_and_lu_mod64(mp_limb_t **a,int row,int col,

mp_limb_t md,int **rinfo,int **cinfo)

{

int i,j,k,d;

int *rp,*cp;

mp_limb_t *t,*pivot;

mp_limb_t 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 = invmod64(pivot[k],md);

for ( i = d+1; i < row; i++ ) {

t = a[i];

if ( (m = t[k]) != 0 ) {

t[k] = mulmod64(inv,m,md);

for ( j = k+1, m = md - t[k]; j < col; j++ )

if ( pivot[j] ) {

t[j] = muladdmod64(m,pivot[j],t[j],md);

}

d++;

}

return d;

}

int lu_mod64(mp_limb_t **a,int n,mp_limb_t md,int **rinfo)

{

int i,j,k;

int *rp;

mp_limb_t *t,*pivot;

mp_limb_t 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 = invmod64(pivot[k],md);

for ( i = k+1; i < n; i++ ) {

t = a[i];

if ( (m = t[k]) != 0 ) {

t[k] = mulmod64(inv,m,md);

for ( j = k+1, m = md - t[k]; j < n; j++ )

if ( pivot[j] ) {

t[j] = muladdmod64(m,pivot[j],t[j],md);

}

return 1;

}

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_mod64(mp_limb_t **a,int n,mp_limb_t md,mp_limb_signed_t **b,int l,int normalize)

{

mp_limb_t *y,*c;

int i,j,k;

mp_limb_t t,m,m2;

y = (mp_limb_t *)MALLOC_ATOMIC(n*sizeof(mp_limb_t));

c = (mp_limb_t *)MALLOC_ATOMIC(n*sizeof(mp_limb_t));

m2 = md/2;

for ( k = 0; k < l; k++ ) {

/* copy b[.][k] to c */

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

c[i] = 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];

t = muladdmod64(m,y[j],t,md);

}

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];

t = muladdmod64(m,c[j],t,md);

}

/* a[i][i] = 1/U[i][i] */

c[i] = mulmod64(t,a[i][i],md);

}

/* copy c to b[.][k] with normalization */

if ( normalize )

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

b[i][k] = (mp_limb_signed_t)(c[i]>m2 ? c[i]-md : c[i]);

else

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

b[i][k] = (mp_limb_signed_t)c[i];

}

#endif

void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)

Line 1692 void reduce_sp_by_red_mod_compress (int *sp,CDP *redma

Line 1944 void reduce_sp_by_red_mod_compress (int *sp,CDP *redma

for ( i = nred-1; i >= 0; i-- ) {

/* reduce sp by redmat[i] */

usp[ind[i]] %= md;

if ( hc = usp[ind[i]] ) {

if ( ( hc = usp[ind[i]] ) != 0 ) {

/* sp = sp-hc*redmat[i] */

hc = md-hc;

ri = redmat[i];

Line 1738 int generic_gauss_elim_mod(int **mat0,int row,int col,

Line 1990 int generic_gauss_elim_mod(int **mat0,int row,int col,

}

for ( i = rank+1; i < row; i++ ) {

t = mat[i];

if ( a = t[j] )

if ( ( a = t[j] ) != 0 )

red_by_vect(md,t+j,pivot+j,md-a,col-j);

}

rank++;

Line 1749 int generic_gauss_elim_mod(int **mat0,int row,int col,

Line 2001 int generic_gauss_elim_mod(int **mat0,int row,int col,

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

t = mat[i];

t[j] %= md;

if ( a = t[j] )

if ( ( a = t[j] ) != 0 )

red_by_vect(md,t+j,pivot+j,md-a,col-j);

}

l--;

Line 1798 int generic_gauss_elim_mod2(int **mat0,int row,int col

Line 2050 int generic_gauss_elim_mod2(int **mat0,int row,int col

}

for ( i = rank+1; i < row; i++ ) {

t = mat[i];

if ( a = t[j] )

if ( ( a = t[j] ) != 0 )

red_by_vect(md,t+j,pivot+j,md-a,col-j);

}

rank++;

Line 1809 int generic_gauss_elim_mod2(int **mat0,int row,int col

Line 2061 int generic_gauss_elim_mod2(int **mat0,int row,int col

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

t = mat[i];

t[j] %= md;

if ( a = t[j] )

if ( ( a = t[j] ) != 0 )

red_by_vect(md,t+j,pivot+j,md-a,col-j);

}

l--;

Line 1854 int indep_rows_mod(int **mat0,int row,int col,int md,i

Line 2106 int indep_rows_mod(int **mat0,int row,int col,int md,i

}

for ( i = rank+1; i < row; i++ ) {

t = mat[i];

if ( a = t[j] )

if ( ( a = t[j] ) != 0 )

red_by_vect(md,t+j,pivot+j,md-a,col-j);

}

rank++;

Line 1888 int generic_gauss_elim_sf(int **mat0,int row,int col,i

Line 2140 int generic_gauss_elim_sf(int **mat0,int row,int col,i

*pk = _mulsf(*pk,inv);

for ( i = rank+1; i < row; i++ ) {

t = mat[i];

if ( a = t[j] )

if ( ( a = t[j] ) != 0 )

red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);

}

rank++;

Line 1898 int generic_gauss_elim_sf(int **mat0,int row,int col,i

Line 2150 int generic_gauss_elim_sf(int **mat0,int row,int col,i

pivot = mat[l];

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

t = mat[i];

if ( a = t[j] )

if ( ( a = t[j] ) != 0 )

red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);

}

l--;

Line 1934 int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)

Line 2186 int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)

pivot[k] = inv = invm(pivot[k],md);

for ( i = k+1; i < row; i++ ) {

t = a[i];

if ( m = t[k] ) {

if ( ( m = t[k] ) != 0 ) {

DMAR(inv,m,0,md,t[k])

for ( j = k+1, m = md - t[k]; j < col; j++ )

if ( pivot[j] ) {

Line 1990 int find_lhs_and_lu_mod(unsigned int **a,int row,int c

Line 2242 int find_lhs_and_lu_mod(unsigned int **a,int row,int c

pivot[k] = inv = invm(pivot[k],md);

for ( i = d+1; i < row; i++ ) {

t = a[i];

if ( m = t[k] ) {

if ( ( m = t[k] ) != 0 ) {

DMAR(inv,m,0,md,t[k])

for ( j = k+1, m = md - t[k]; j < col; j++ )

if ( pivot[j] ) {

Line 2025 int lu_mod(unsigned int **a,int n,unsigned int md,int

Line 2277 int lu_mod(unsigned int **a,int n,unsigned int md,int

inv = invm(pivot[k],md);

for ( i = k+1; i < n; i++ ) {

t = a[i];

if ( m = t[k] ) {

if ( ( m = t[k] ) != 0 ) {

DMAR(inv,m,0,md,t[k])

for ( j = k+1, m = md - t[k]; j < n; j++ )

if ( pivot[j] ) {

Line 2054 void solve_by_lu_mod(int **a,int n,int md,int **b,int

Line 2306 void solve_by_lu_mod(int **a,int n,int md,int **b,int

int i,j,k;

unsigned int t,m,m2;

y = (int *)MALLOC_ATOMIC(n*sizeof(int));

y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));

c = (int *)MALLOC_ATOMIC(n*sizeof(int));

c = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));

m2 = md>>1;

for ( k = 0; k < l; k++ ) {

/* copy b[.][k] to c */

Line 2102 void Pleqm1(NODE arg,VECT *rp)

Line 2354 void Pleqm1(NODE arg,VECT *rp)

asir_assert(ARG0(arg),O_MAT,"leqm1");

asir_assert(ARG1(arg),O_N,"leqm1");

m = (MAT)ARG0(arg); md = QTOS((Z)ARG1(arg));

m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));

row = m->row; col = m->col; mat = m->body;

wmat = (int **)almat(row,col);

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

Line 2117 void Pleqm1(NODE arg,VECT *rp)

Line 2369 void Pleqm1(NODE arg,VECT *rp)

n = col - 1;

MKVECT(vect,n);

for ( i = 0, v = (Z *)vect->body; i < n; i++ ) {

t = (md-wmat[i][n])%md; STOQ(t,v[i]);

t = (md-wmat[i][n])%md; STOZ(t,v[i]);

}

*rp = vect;

}

Line 2173 void Pgeninvm(NODE arg,LIST *rp)

Line 2425 void Pgeninvm(NODE arg,LIST *rp)

asir_assert(ARG0(arg),O_MAT,"leqm1");

asir_assert(ARG1(arg),O_N,"leqm1");

m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));

m = (MAT)ARG0(arg); md = ZTOS((Q)ARG1(arg));

row = m->row; col = m->col; mat = m->body;

wmat = (unsigned int **)almat(row,col+row);

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

Line 2188 void Pgeninvm(NODE arg,LIST *rp)

Line 2440 void Pgeninvm(NODE arg,LIST *rp)

MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);

for ( i = 0, tmat = (Z **)mat1->body; i < col; i++ )

for ( j = 0; j < row; j++ )

UTOQ(wmat[i][j+col],tmat[i][j]);

UTOZ(wmat[i][j+col],tmat[i][j]);

for ( tmat = (Z **)mat2->body; i < row; i++ )

for ( j = 0; j < row; j++ )

UTOQ(wmat[i][j+col],tmat[i-col][j]);

UTOZ(wmat[i][j+col],tmat[i-col][j]);

MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);

}

Line 2215 int gauss_elim_geninv_mod(unsigned int **mat,int row,i

Line 2467 int gauss_elim_geninv_mod(unsigned int **mat,int row,i

pivot[k] = dmar(pivot[k],inv,0,md);

for ( i = j+1; i < row; i++ ) {

t = mat[i];

if ( a = t[j] )

if ( ( a = t[j] ) != 0 )

for ( k = j, a = md - a; k < m; k++ )

t[k] = dmar(pivot[k],a,t[k],md);

}

Line 2224 int gauss_elim_geninv_mod(unsigned int **mat,int row,i

Line 2476 int gauss_elim_geninv_mod(unsigned int **mat,int row,i

pivot = mat[j];

for ( i = j-1; i >= 0; i-- ) {

t = mat[i];

if ( a = t[j] )

if ( ( a = t[j] ) != 0 )

for ( k = j, a = md - a; k < m; k++ )

t[k] = dmar(pivot[k],a,t[k],md);

}

Line 2244 void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)

Line 2496 void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)

lu = (GFMMAT)ARG0(arg);

perm = (Z *)BDY((VECT)ARG1(arg));

rhs = (Z *)BDY((VECT)ARG2(arg));

md = (unsigned int)QTOS((Z)ARG3(arg));

md = (unsigned int)ZTOS((Z)ARG3(arg));

n = lu->col;

b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));

sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));

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

b[i] = QTOS(rhs[QTOS(perm[i])]);

b[i] = ZTOS(rhs[ZTOS(perm[i])]);

solve_by_lu_gfmmat(lu,md,b,sol);

MKVECT(r,n);

for ( i = 0, v = (Z *)r->body; i < n; i++ )

UTOQ(sol[i],v[i]);

UTOZ(sol[i],v[i]);

*rp = r;

}

Line 2308 void Plu_mat(NODE arg,LIST *rp)

Line 2560 void Plu_mat(NODE arg,LIST *rp)

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]);

STOZ(iperm[i],v[i]);

n0 = mknode(3,lu,dn,perm);

MKLIST(*rp,n0);

}

Line 2327 void Plu_gfmmat(NODE arg,LIST *rp)

Line 2579 void Plu_gfmmat(NODE arg,LIST *rp)

asir_assert(ARG0(arg),O_MAT,"lu_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)ZTOS((Q)ARG1(arg));

mat_to_gfmmat(m,md,&mm);

row = m->row;

col = m->col;

Line 2338 void Plu_gfmmat(NODE arg,LIST *rp)

Line 2590 void Plu_gfmmat(NODE arg,LIST *rp)

else {

MKVECT(perm,row);

for ( i = 0, v = (Z *)perm->body; i < row; i++ )

STOQ(iperm[i],v[i]);

STOZ(iperm[i],v[i]);

n0 = mknode(2,mm,perm);

}

MKLIST(*rp,n0);

Line 2351 void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)

Line 2603 void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)

asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");

asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");

m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));

m = (MAT)ARG0(arg); md = (unsigned int)ZTOS((Q)ARG1(arg));

mat_to_gfmmat(m,md,rp);

}

Line 2390 void Pgeninvm_swap(NODE arg,LIST *rp)

Line 2642 void Pgeninvm_swap(NODE arg,LIST *rp)

asir_assert(ARG0(arg),O_MAT,"geninvm_swap");

asir_assert(ARG1(arg),O_N,"geninvm_swap");

m = (MAT)ARG0(arg); md = QTOS((Z)ARG1(arg));

m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));

row = m->row; col = m->col; mat = m->body;

wmat = (unsigned int **)almat(row,col+row);

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

Line 2406 void Pgeninvm_swap(NODE arg,LIST *rp)

Line 2658 void Pgeninvm_swap(NODE arg,LIST *rp)

MKMAT(mat1,col,col);

for ( i = 0, tmat = (Z **)mat1->body; i < col; i++ )

for ( j = 0; j < col; j++ )

UTOQ(invmat[i][j],tmat[i][j]);

UTOZ(invmat[i][j],tmat[i][j]);

MKVECT(vect1,row);

for ( i = 0, tvect = (Z *)vect1->body; i < row; i++ )

STOQ(index[i],tvect[i]);

STOZ(index[i],tvect[i]);

MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);

}

Line 2442 int gauss_elim_geninv_mod_swap(unsigned int **mat,int

Line 2694 int gauss_elim_geninv_mod_swap(unsigned int **mat,int

pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);

for ( i = j+1; i < row; i++ ) {

t = mat[i];

if ( a = t[j] )

if ( ( a = t[j] ) != 0 )

for ( k = j, a = md - a; k < m; k++ )

if ( pivot[k] )

t[k] = dmar(pivot[k],a,t[k],md);

Line 2452 int gauss_elim_geninv_mod_swap(unsigned int **mat,int

Line 2704 int gauss_elim_geninv_mod_swap(unsigned int **mat,int

pivot = mat[j];

for ( i = j-1; i >= 0; i-- ) {

t = mat[i];

if ( a = t[j] )

if ( ( a = t[j] ) != 0 )

for ( k = j, a = md - a; k < m; k++ )

if ( pivot[k] )

t[k] = dmar(pivot[k],a,t[k],md);

Line 2488 void Pgeninv_sf_swap(NODE arg,LIST *rp)

Line 2740 void Pgeninv_sf_swap(NODE arg,LIST *rp)

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] )

if ( ( q = (GFS)mat[i][j] ) != 0 )

wmat[i][j] = FTOIF(CONT(q));

wmat[i][col+i] = _onesf();

}

Line 2499 void Pgeninv_sf_swap(NODE arg,LIST *rp)

Line 2751 void Pgeninv_sf_swap(NODE arg,LIST *rp)

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] ) {

if ( ( t = invmat[i][j] ) != 0 ) {

MKGFS(IFTOF(t),tmat[i][j]);

}

MKVECT(vect1,row);

for ( i = 0, tvect = (Z *)vect1->body; i < row; i++ )

STOQ(index[i],tvect[i]);

STOZ(index[i],tvect[i]);

MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);

}

Line 2536 int gauss_elim_geninv_sf_swap(int **mat,int row,int co

Line 2788 int gauss_elim_geninv_sf_swap(int **mat,int row,int co

pivot[k] = _mulsf(pivot[k],inv);

for ( i = j+1; i < row; i++ ) {

t = mat[i];

if ( a = t[j] )

if ( ( a = t[j] ) != 0 )

for ( k = j, a = _chsgnsf(a); k < m; k++ )

if ( pivot[k] ) {

u = _mulsf(pivot[k],a);

Line 2548 int gauss_elim_geninv_sf_swap(int **mat,int row,int co

Line 2800 int gauss_elim_geninv_sf_swap(int **mat,int row,int co

pivot = mat[j];

for ( i = j-1; i >= 0; i-- ) {

t = mat[i];

if ( a = t[j] )

if ( ( a = t[j] ) != 0 )

for ( k = j, a = _chsgnsf(a); k < m; k++ )

if ( pivot[k] ) {

u = _mulsf(pivot[k],a);

Line 2596 void Pnbpoly_up2(NODE arg,GF2N *rp)

Line 2848 void Pnbpoly_up2(NODE arg,GF2N *rp)

int m,type,ret;

UP2 r;

m = QTOS((Z)ARG0(arg));

m = ZTOS((Z)ARG0(arg));

type = QTOS((Z)ARG1(arg));

type = ZTOS((Z)ARG1(arg));

ret = generate_ONB_polynomial(&r,m,type);

if ( ret == 0 )

MKGF2N(r,*rp);

Line 2611 void Px962_irredpoly_up2(NODE arg,GF2N *rp)

Line 2863 void Px962_irredpoly_up2(NODE arg,GF2N *rp)

GF2N prev;

UP2 r;

m = QTOS((Q)ARG0(arg));

m = ZTOS((Q)ARG0(arg));

prev = (GF2N)ARG1(arg);

if ( !prev ) {

w = (m>>5)+1; NEWUP2(r,w); r->w = 0;

Line 2636 void Pirredpoly_up2(NODE arg,GF2N *rp)

Line 2888 void Pirredpoly_up2(NODE arg,GF2N *rp)

GF2N prev;

UP2 r;

m = QTOS((Q)ARG0(arg));

m = ZTOS((Q)ARG0(arg));

prev = (GF2N)ARG1(arg);

if ( !prev ) {

w = (m>>5)+1; NEWUP2(r,w); r->w = 0;

Line 2665 void Pmat_swap_row_destructive(NODE arg, MAT *m)

Line 2917 void Pmat_swap_row_destructive(NODE arg, MAT *m)

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));

i1 = ZTOS((Q)ARG1(arg));

i2 = QTOS((Q)ARG2(arg));

i2 = ZTOS((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];

Line 2686 void Pmat_swap_col_destructive(NODE arg, MAT *m)

Line 2938 void Pmat_swap_col_destructive(NODE arg, MAT *m)

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));

j1 = ZTOS((Q)ARG1(arg));

j2 = QTOS((Q)ARG2(arg));

j2 = ZTOS((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;

Line 2721 int generate_ONB_polynomial(UP2 *rp,int m,int type)

Line 2973 int generate_ONB_polynomial(UP2 *rp,int m,int type)

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

f->b[i] = 0xffffffff;

/* mask the top word if necessary */

if ( r = (m+1)&31 )

if ( ( r = (m+1)&31 ) != 0 )

f->b[w-1] &= (1<<r)-1;

return 0;

break;

Line 2999 void Pnd_det(NODE arg,P *rp)

Line 3251 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);

nd_det(ZTOS((Q)ARG1(arg)),ARG0(arg),rp);

}

void Pmat_col(NODE arg,VECT *rp)

Line 3011 void Pmat_col(NODE arg,VECT *rp)

Line 3263 void Pmat_col(NODE arg,VECT *rp)

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));

j = ZTOS((Q)ARG1(arg));

if ( j < 0 || j >= mat->col) {

error("mat_col : Out of range");

}

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