/**************************************************************************** ** ** This file is part of GAP, a system for computational discrete algebra. ** ** Copyright of GAP belongs to its developers, whose names are too numerous ** to list here. Please refer to the COPYRIGHT file for details. ** ** SPDX-License-Identifier: GPL-2.0-or-later ** ** This file contains the functions of the package with the operations for ** generic lists.
*/
/**************************************************************************** ** *F EqListList(<listL>,<listR>) . . . . . . . . . test if two lists are equal ** ** 'EqListList' returns 'true' if the two lists <listL> and <listR> are ** equal and 'false' otherwise. This is a generic function, which works for ** all types of lists.
*/ Int EqListList (
Obj listL,
Obj listR )
{ Int lenL; // length of the left operand Int lenR; // length of the right operand
Obj elmL; // element of the left operand
Obj elmR; // element of the right operand Int i; // loop variable
// get the lengths of the lists and compare them
lenL = LEN_LIST( listL );
lenR = LEN_LIST( listR ); if ( lenL != lenR ) { return 0;
}
// loop over the elements and compare them for ( i = 1; i <= lenL; i++ ) {
elmL = ELMV0_LIST( listL, i );
elmR = ELMV0_LIST( listR, i ); if ( elmL == 0 && elmR != 0 ) { return 0;
} elseif ( elmR == 0 && elmL != 0 ) { return 0;
} elseif ( ! EQ( elmL, elmR ) ) { return 0;
}
}
// no differences found, the lists are equal return 1;
}
/**************************************************************************** ** *F LtListList(<listL>,<listR>) . . . . . . . . . test if two lists are equal ** ** 'LtListList' returns 'true' if the list <listL> is less than the list ** <listR> and 'false' otherwise. This is a generic function, which works ** for all types of lists.
*/ Int LtListList (
Obj listL,
Obj listR )
{ Int lenL; // length of the left operand Int lenR; // length of the right operand
Obj elmL; // element of the left operand
Obj elmR; // element of the right operand Int i; // loop variable
// get the lengths of the lists and compare them
lenL = LEN_LIST( listL );
lenR = LEN_LIST( listR );
// loop over the elements and compare them for ( i = 1; i <= lenL && i <= lenR; i++ ) {
elmL = ELMV0_LIST( listL, i );
elmR = ELMV0_LIST( listR, i ); if ( elmL == 0 && elmR != 0 ) { return 1;
} elseif ( elmR == 0 && elmL != 0 ) { return 0;
} elseif ( ! EQ( elmL, elmR ) ) { return LT( elmL, elmR );
}
}
// reached the end of at least one list return (lenL < lenR);
}
/**************************************************************************** ** *F InList(<objL>,<listR>) . . . . . . . . . . . . membership test for lists ** ** 'InList' returns a nonzero value if <objL> is a member of the list ** <listR>, and zero otherwise.
*/ staticInt InList(Obj objL, Obj listR)
{ return Fail != POS_LIST(listR, objL, INTOBJ_INT(0));
}
/**************************************************************************** ** *F SumSclList(<listL>,<listR>) . . . . . . . . . sum of a scalar and a list *F SumListScl(<listL>,<listR>) . . . . . . . . . sum of a list and a scalar *F SumListList<listL>,<listR>) . . . . . . . . . . . . . sum of two lists ** ** 'SumSclList' is a generic function for the first kind of sum, that of a ** scalar and a list. ** ** 'SumListScl' is a generic function for the second kind of sum, that of a ** list and a scalar. ** ** 'SumListList' is a generic function for the third kind of sum, that of ** two lists.
*/
Obj SumSclList (
Obj listL,
Obj listR )
{
Obj listS; // sum, result
Obj elmS; // one element of sum list
Obj elmR; // one element of right operand Int len; // length Int i; // loop variable
// make the result list
len = LEN_LIST( listR );
listS = NEW_PLIST_WITH_MUTABILITY( IS_MUTABLE_OBJ(listL) || IS_MUTABLE_OBJ(listR),
T_PLIST, len );
SET_LEN_PLIST( listS, len );
// loop over the entries and add for ( i = 1; i <= len; i++ ) {
elmR = ELMV0_LIST( listR, i ); if (elmR)
{
elmS = SUM( listL, elmR );
SET_ELM_PLIST( listS, i, elmS );
CHANGED_BAG( listS );
}
}
return listS;
}
Obj SumListScl (
Obj listL,
Obj listR )
{
Obj listS; // sum, result
Obj elmS; // one element of sum list
Obj elmL; // one element of left operand Int len; // length Int i; // loop variable
// make the result list
len = LEN_LIST( listL );
listS = NEW_PLIST_WITH_MUTABILITY( IS_MUTABLE_OBJ(listR) || IS_MUTABLE_OBJ(listL),
T_PLIST, len );
SET_LEN_PLIST( listS, len );
// loop over the entries and add for ( i = 1; i <= len; i++ ) {
elmL = ELMV0_LIST( listL, i ); if (elmL)
{
elmS = SUM( elmL, listR );
SET_ELM_PLIST( listS, i, elmS );
CHANGED_BAG( listS );
}
}
return listS;
}
Obj SumListList (
Obj listL,
Obj listR )
{
Obj listS; // sum, result
Obj elmS; // one element of the sum
Obj elmL; // one element of the left list
Obj elmR; // one element of the right list Int lenL,lenR, lenS;// lengths Int i; // loop variable
UInt mutS;
/**************************************************************************** ** *F ZeroList(<list>) . . . . . . . . . . . . . . . . . . . . zero of a list *F ZeroListDefault(<list>) . . . . . . . . . . . . . . . . . zero of a list ** ** 'ZeroList' is the extended dispatcher for the zero involving lists. That ** is, whenever zero for a list is called and 'ZeroSameMutFuncs' does not ** point to a special function, then 'ZeroList' is called. 'ZeroList' ** determines the extended type of the operand and then dispatches through ** 'ZeroSameMutFuncs' again. ** ** 'ZeroListDefault' is a generic function for the zero.
*/ static Obj ZeroListDefault(Obj list)
{
Obj res; // Obj elm; Int len; Int i;
// make the result list -- same mutability as argument
len = LEN_LIST( list ); if (len == 0) { return NEW_PLIST_WITH_MUTABILITY(IS_MUTABLE_OBJ(list), T_PLIST_EMPTY, 0);
}
res = NEW_PLIST_WITH_MUTABILITY( IS_MUTABLE_OBJ(list), T_PLIST, len );
SET_LEN_PLIST( res, len );
// enter zeroes everywhere // For now, lets just do the simplest and safest thing for (i = 1; i <= len; i++ )
{
Obj tmp = ELM0_LIST( list, i); if (tmp) {
tmp = ZERO_SAMEMUT(tmp);
SET_ELM_PLIST( res, i,tmp );
CHANGED_BAG( res);
}
} // Now adjust the result TNUM info
static Obj ZeroListMutDefault(Obj list)
{
Obj res; // Obj elm; Int len; Int i;
// make the result list -- always mutable
len = LEN_LIST( list ); if (len == 0) { return NewEmptyPlist();
}
res = NEW_PLIST( T_PLIST ,len );
SET_LEN_PLIST( res, len );
// enter zeroes everywhere // For now, lets just do the simplest and safest thing for (i = 1; i <= len; i++ )
{
Obj tmp = ELM0_LIST( list, i); if (tmp) {
tmp = ZERO_MUT(tmp);
SET_ELM_PLIST( res, i,tmp );
CHANGED_BAG( res);
}
} // Now adjust the result TNUM info
/* This is intended to be installed as a method for the Attribute Zero, for rectangular matrices rather than the Operation ZeroOp. This is useful because, knowing that we want an immutable result, we can (a) reuse a single row of zeros
(b) record that the result is a rectangular table */
static Obj FuncZERO_ATTR_MAT(Obj self, Obj mat)
{
Obj zrow;
UInt len;
UInt i;
Obj res;
len = LEN_LIST(mat); if (len == 0) return NewImmutableEmptyPlist();
zrow = ZERO_SAMEMUT(ELM_LIST(mat,1));
CheckedMakeImmutable(zrow);
res = NEW_PLIST_IMM(T_PLIST_TAB_RECT, len);
SET_LEN_PLIST(res,len); for (i = 1; i <= len; i++)
SET_ELM_PLIST(res,i,zrow); return res;
}
/**************************************************************************** ** *F AInvListDefault(<list>) . . . . . . . . . . . additive inverse of a list *F AInvMutListDefault ** ** 'AInvList' is the extended dispatcher for the additive inverse involving ** lists. That is, whenever the additive inverse for lists is called and ** 'AInvSameMutFuncs' does not point to a special function, then 'AInvList' ** is called.'AInvList' determines the extended type of the operand and then ** dispatches through 'AInvSameMutFuncs' again. ** ** 'AInvListDefault' is a generic function for the additive inverse.
*/
static Obj AInvMutListDefault(Obj list)
{
Obj res;
Obj elm; Int len; Int i;
/* make the result list -- always mutable, since this might be a method for
AdditiveInverseOp */
len = LEN_LIST( list ); if (len == 0) { return NewEmptyPlist();
}
res = NEW_PLIST( T_PLIST , len );
SET_LEN_PLIST( res, len );
// enter the additive inverses everywhere for ( i = 1; i <= len; i++ ) {
elm = ELM0_LIST( list, i ); if (elm) {
elm = AINV_MUT( elm );
SET_ELM_PLIST( res, i, elm );
CHANGED_BAG( res );
}
}
static Obj AInvListDefault(Obj list)
{
Obj res;
Obj elm; Int len; Int i;
// make the result list -- same mutability as input
len = LEN_LIST( list ); if (len == 0) { return NEW_PLIST_WITH_MUTABILITY(IS_MUTABLE_OBJ(list), T_PLIST_EMPTY, 0);
}
res = NEW_PLIST_WITH_MUTABILITY( IS_MUTABLE_OBJ(list), T_PLIST, len );
SET_LEN_PLIST( res, len );
// enter the additive inverses everywhere for ( i = 1; i <= len; i++ ) {
elm = ELM0_LIST( list, i ); if (elm) {
elm = AINV_SAMEMUT( elm );
SET_ELM_PLIST( res, i, elm );
CHANGED_BAG( res );
}
}
/**************************************************************************** ** *F DiffSclList(<listL>,<listR>) . . . . . difference of a scalar and a list *F DiffListScl(<listL>,<listR>) . . . . . difference of a list and a scalar *F DiffListList(<listL>,<listR>) . . . . . . . . . . difference of two lists ** ** 'DiffSclList' is a generic function for the first kind of difference, ** that of a scalar and a list. ** ** 'DiffListScl' is a generic function for the second kind of difference, ** that of a list and a scalar. ** ** 'DiffListList' is a generic function for the third kind of difference, ** that of two lists.
*/
Obj DiffSclList (
Obj listL,
Obj listR )
{
Obj listD; // difference, result
Obj elmD; // one element of difference list
Obj elmR; // one element of right operand Int len; // length Int i; // loop variable Int mut;
// make the result list
len = LEN_LIST( listR );
mut = IS_MUTABLE_OBJ(listL) || IS_MUTABLE_OBJ(listR); if (len == 0) { return NEW_PLIST_WITH_MUTABILITY(mut, T_PLIST_EMPTY, 0);
}
listD = NEW_PLIST_WITH_MUTABILITY(mut, T_PLIST, len);
SET_LEN_PLIST( listD, len );
// loop over the entries and subtract for ( i = 1; i <= len; i++ ) {
elmR = ELMV0_LIST( listR, i ); if (elmR)
{
elmD = DIFF( listL, elmR );
SET_ELM_PLIST( listD, i, elmD );
CHANGED_BAG( listD );
}
}
Obj DiffListScl (
Obj listL,
Obj listR )
{
Obj listD; // difference, result
Obj elmD; // one element of difference list
Obj elmL; // one element of left operand Int len; // length Int i; // loop variable Int mut;
// make the result list
len = LEN_LIST( listL );
mut = IS_MUTABLE_OBJ(listL) || IS_MUTABLE_OBJ(listR); if (len == 0) { return NEW_PLIST_WITH_MUTABILITY(mut, T_PLIST_EMPTY, 0);
}
listD = NEW_PLIST_WITH_MUTABILITY(mut, T_PLIST, len);
SET_LEN_PLIST( listD, len );
// loop over the entries and subtract for ( i = 1; i <= len; i++ ) {
elmL = ELMV0_LIST( listL, i ); if (elmL)
{
elmD = DIFF( elmL, listR );
SET_ELM_PLIST( listD, i, elmD );
CHANGED_BAG( listD );
}
Obj DiffListList (
Obj listL,
Obj listR )
{
Obj listD; // difference, result
Obj elmD; // one element of the difference
Obj elmL; // one element of the left list
Obj elmR; // one element of the right list Int i; // loop variable
UInt mutD; Int mut;
// Sort out mutability
mutD = 0; for (i = 1; i <= lenL; i++) if ((elmL = ELM0_LIST( listL, i)))
{
mutD = IS_MUTABLE_OBJ(elmL); break;
} for (i = 1; i <= lenR; i++) if ((elmR = ELM0_LIST( listR, i)))
{
mutD = mutD || IS_MUTABLE_OBJ(elmR); break;
}
// loop over the entries and subtract for ( i = 1; i <= lenD; i++ ) {
elmL = ELM0_LIST( listL, i );
elmR = ELM0_LIST( listR, i );
// Now compute the result 6 different cases! if (elmL)
{ if (elmR)
elmD= DIFF(elmL,elmR); elseif ( mutD)
elmD = SHALLOW_COPY_OBJ(elmL); else
elmD = elmL;
} elseif (elmR)
{ if (mutD)
elmD = AINV_MUT(elmR); else
elmD = AINV_SAMEMUT(elmR);
} else
elmD = 0;
if (elmD)
{
SET_ELM_PLIST( listD, i, elmD );
CHANGED_BAG( listD );
}
} // Now adjust the result TNUM info. There's not so much we can say // here with total reliability
/**************************************************************************** ** *F ProdSclList(<listL>,<listR>) . . . . . . product of a scalar and a list *F ProdListScl(<listL>,<listR>) . . . . . . product of a list and a scalar *F ProdListList(<listL>,<listR>) . . . . . . . . . . . product of two lists ** ** 'ProdSclList' is a generic function for the first kind of product, that ** of a scalar and a list. Note that this includes kind of product defines ** the product of a matrix with a list of matrices. ** ** 'ProdListScl' is a generic function for the second kind of product, that ** of a list and a scalar. Note that this kind of product defines the ** product of a matrix with a vector, the product of two matrices, and the ** product of a list of matrices and a matrix. ** ** 'ProdListList' is a generic function for the third kind of product, that ** of two lists. Note that this kind of product defines the product of two ** vectors, a vector and a matrix, and the product of a vector and a list of ** matrices.
*/
Obj ProdSclList (
Obj listL,
Obj listR )
{
Obj listP; // product, result
Obj elmP; // one element of product list
Obj elmR; // one element of right operand Int len; // length Int i; // loop variable Int mut;
// make the result list
len = LEN_LIST( listR );
mut = IS_MUTABLE_OBJ(listL) || IS_MUTABLE_OBJ(listR); if (len == 0) { return NEW_PLIST_WITH_MUTABILITY(mut, T_PLIST_EMPTY, 0);
}
listP = NEW_PLIST_WITH_MUTABILITY(mut, T_PLIST, len);
SET_LEN_PLIST( listP, len );
// loop over the entries and multiply for ( i = 1; i <= len; i++ ) {
elmR = ELMV0_LIST( listR, i ); if (elmR)
{
elmP = PROD( listL, elmR );
SET_ELM_PLIST( listP, i, elmP );
CHANGED_BAG( listP );
}
} if (IS_PLIST( listR ))
{ if (HAS_FILT_LIST(listR, FN_IS_DENSE))
SET_FILT_LIST( listP, FN_IS_DENSE ); elseif (HAS_FILT_LIST(listR, FN_IS_NDENSE))
SET_FILT_LIST( listP, FN_IS_NDENSE );
}
return listP;
}
Obj ProdListScl (
Obj listL,
Obj listR )
{
Obj listP; // product, result
Obj elmP; // one element of product list
Obj elmL; // one element of left operand Int len; // length Int i; // loop variable Int mut;
// make the result list
len = LEN_LIST( listL );
mut = IS_MUTABLE_OBJ(listL) || IS_MUTABLE_OBJ(listR); if (len == 0) { return NEW_PLIST_WITH_MUTABILITY(mut, T_PLIST_EMPTY, 0);
}
listP = NEW_PLIST_WITH_MUTABILITY(mut, T_PLIST, len);
SET_LEN_PLIST( listP, len );
// loop over the entries and multiply for ( i = 1; i <= len; i++ ) {
elmL = ELMV0_LIST( listL, i ); if (elmL) {
elmP = PROD( elmL, listR );
SET_ELM_PLIST( listP, i, elmP );
CHANGED_BAG( listP );
}
}
Obj ProdListList (
Obj listL,
Obj listR )
{
Obj listP; // product, result
Obj elmP; // one summand of the product
Obj elmL; // one element of the left list
Obj elmR; // one element of the right list Int lenL,lenR,len; // length Int i; // loop variable Int imm;
// get and check the length
lenL = LEN_LIST( listL );
lenR = LEN_LIST( listR );
len = (lenL < lenR) ? lenL : lenR; // loop over the entries and multiply and accumulate
listP = 0;
imm = 0; for (i = 1; i <= len; i++)
{
elmL = ELM0_LIST( listL, i );
elmR = ELM0_LIST( listR, i ); if (elmL && elmR)
{
elmP = PROD( elmL, elmR ); if (listP)
listP = SUM( listP, elmP ); else
{
listP = elmP;
imm = !IS_MUTABLE_OBJ(listP);
}
}
}
// TODO: This is possible expensive, we may be able to settle for // a cheaper check and call MakeImmutable() instead. if (imm && IS_MUTABLE_OBJ(listP))
CheckedMakeImmutable(listP);
if (!listP)
ErrorMayQuit("Inner product multiplication of lists: no summands", 0, 0);
// possibly adjust mutability if (!IS_MUTABLE_OBJ(prod)) switch (INT_INTOBJ(depthdiff)) { case -1: if (IS_MUTABLE_OBJ(listL))
prod = SHALLOW_COPY_OBJ(prod); break; case 1: if (IS_MUTABLE_OBJ(listR))
prod = SHALLOW_COPY_OBJ(prod); break; case 0: break; default:
ErrorMayQuit("PROD_LIST_LIST_DEFAULT: depth difference should be -1, " "0 or 1, not %d",
INT_INTOBJ(depthdiff), 0);
} return prod;
}
/**************************************************************************** ** *F OneMatrix(<list>, <mut>) . . . . . . . . . . . . . . . . . one of a list ** ** ** 'OneMatrix' is a generic function for the one. mut may be 0, for an ** immutable result, 1 for a result of the same mutability level as <list> ** and 2 for a fully mutable result.
*/
static Obj OneMatrix(Obj mat, UInt mut)
{
Obj res = 0; // one, result
Obj row; // one row of the result
Obj zero = 0; // zero element
Obj one = 0; // one element
UInt len; // length (and width) of matrix
UInt i, k; // loop variables BOOL rmut; // rows of result mutable? BOOL cmut; // result mutable?
GAP_ASSERT(mut <= 2);
// check that the operand is a *square* matrix
len = LEN_LIST( mat ); if ( len != LEN_LIST( ELM_LIST( mat, 1 ) ) ) {
ErrorMayQuit("Matrix ONE: must be square (not %d by %d)",
(Int)len, (Int)LEN_LIST(ELM_LIST(mat, 1)));
}
// get the zero and the one switch (mut) { case 0:
zero = ZERO_MUT( ELM_LIST( ELM_LIST( mat, 1 ), 1 ) );
one = ONE_SAMEMUT(zero);
CheckedMakeImmutable(zero);
CheckedMakeImmutable(one);
cmut = rmut = FALSE; break;
case 1:
zero = ZERO_MUT( ELM_LIST( ELM_LIST( mat, 1 ), 1 ) );
one = ONE_SAMEMUT(zero); if (IS_MUTABLE_OBJ(mat))
{
cmut = TRUE;
rmut = IS_MUTABLE_OBJ(ELM_LIST(mat, 1));
} else
cmut = rmut = FALSE; break;
case 2:
zero = ZERO_SAMEMUT( ELM_LIST( ELM_LIST( mat, 1 ), 1 ) );
one = ONE( zero );
cmut = rmut = TRUE; break;
}
// make the identity matrix
res = NEW_PLIST(T_PLIST, len);
SET_LEN_PLIST( res, len ); for ( i = 1; i <= len; i++ ) {
row = NEW_PLIST(T_PLIST, len);
SET_LEN_PLIST( row, len ); for ( k = 1; k <= len; k++ )
SET_ELM_PLIST( row, k, zero );
SET_ELM_PLIST( row, i, one ); if (!rmut)
MakeImmutableNoRecurse(row);
SET_ELM_PLIST( res, i, row );
CHANGED_BAG( res );
} if (!cmut)
MakeImmutableNoRecurse(res);
/**************************************************************************** ** *F InvList(<list>) . . . . . . . . . . . . . . . . . . . . inverse of a list *F InvMatrix(<list>) . . . . . . . . . . . . . . . . . . . inverse of a list ** ** 'InvList' is the extended dispatcher for inverses involving lists. That ** is, whenever inverse for a list is called and 'InvFuncs' does not point ** to a special function, then 'InvList' is called. 'InvList' determines ** the extended type of the operand and then dispatches through 'InvList' ** again. ** ** 'InvMatrix' is a generic function for the inverse. In nearly all ** circumstances, we should use a more efficient function based on ** calls to AddRowVector, etc.
*/
static Obj InvMatrix(Obj mat, UInt mut)
{
Obj res = 0; // power, result
Obj row; // one row of the matrix
Obj row2; // another row of the matrix
Obj elm; // one element of the matrix
Obj elm2; // another element of the matrix
Obj zero = 0; // zero element
Obj one = 0; // one element
UInt len; // length (and width) of matrix
UInt i, k, l; // loop variables BOOL rmut; // rows of result mutable? BOOL cmut; // result mutable?
GAP_ASSERT(mut <= 2);
// check that the operand is a *square* matrix
len = LEN_LIST( mat ); if ( len != LEN_LIST( ELM_LIST( mat, 1 ) ) ) {
ErrorMayQuit("Matrix INV: must be square (not %d by %d)",
(Int)len, (Int)LEN_LIST(ELM_LIST(mat, 1)));
}
// get the zero and the one switch (mut) { case 0:
zero = ZERO_MUT( ELM_LIST( ELM_LIST( mat, 1 ), 1 ) );
one = ONE_SAMEMUT(zero);
CheckedMakeImmutable(zero);
CheckedMakeImmutable(one);
cmut = rmut = FALSE; break;
case 1:
zero = ZERO_MUT( ELM_LIST( ELM_LIST( mat, 1 ), 1 ) );
one = ONE_SAMEMUT(zero); if (IS_MUTABLE_OBJ(mat)) {
cmut = TRUE;
rmut = IS_MUTABLE_OBJ(ELM_LIST(mat, 1));
} else
cmut = rmut = FALSE; break;
case 2:
zero = ZERO_SAMEMUT( ELM_LIST( ELM_LIST( mat, 1 ), 1 ) );
one = ONE( zero );
cmut = rmut = TRUE; break;
}
// make a matrix of the form $ ( Id_<len> | <mat> ) $
res = NEW_PLIST(T_PLIST, len);
SET_LEN_PLIST(res, len); for (i = 1; i <= len; i++) {
row = NEW_PLIST(T_PLIST, 2 * len);
SET_LEN_PLIST(row, 2 * len);
SET_ELM_PLIST(res, i, row);
CHANGED_BAG(res);
} for ( i = 1; i <= len; i++ ) {
row = ELM_PLIST( res, i ); for ( k = 1; k <= len; k++ )
SET_ELM_PLIST( row, k, zero );
SET_ELM_PLIST( row, i, one );
} for ( i = 1; i <= len; i++ ) {
row = ELM_PLIST( res, i );
row2 = ELM_LIST( mat, i ); for ( k = 1; k <= len; k++ ) {
SET_ELM_PLIST( row, k + len, ELM_LIST( row2, k ) );
CHANGED_BAG( row );
}
}
// make row operations to reach form $ ( <inv> | Id_<len> ) $ // loop over the columns of <mat> for ( k = len+1; k <= 2*len; k++ ) {
// find a nonzero entry in this column for ( i = k-len; i <= len; i++ ) { if ( ! EQ( ELM_PLIST( ELM_PLIST(res,i), k ), zero ) ) break;
} if ( len < i ) { return Fail;
}
// make the row the <k>-th row and normalize it
row = ELM_PLIST( res, i );
SET_ELM_PLIST( res, i, ELM_PLIST( res, k-len ) );
SET_ELM_PLIST( res, k-len, row ); if (mut < 2)
elm2 = INV_SAMEMUT( ELM_PLIST( row, k ) ); else
elm2 = INV( ELM_PLIST( row, k ) ); for ( l = 1; l <= 2*len; l++ ) {
elm = PROD( elm2, ELM_PLIST( row, l ) );
SET_ELM_PLIST( row, l, elm );
CHANGED_BAG( row );
}
// clear all entries in this column for ( i = 1; i <= len; i++ ) {
row2 = ELM_PLIST( res, i ); if (mut < 2)
elm = AINV_MUT( ELM_PLIST( row2, k ) ); else
elm = AINV_SAMEMUT( ELM_PLIST( row2, k ) ); if ( i != k-len && ! EQ(elm,zero) ) { for ( l = 1; l <= 2*len; l++ ) {
elm2 = PROD( elm, ELM_PLIST( row, l ) );
elm2 = SUM( ELM_PLIST( row2, l ), elm2 );
SET_ELM_PLIST( row2, l, elm2 );
CHANGED_BAG( row2 );
}
}
}
}
// throw away the right halves of each row for ( i = 1; i <= len; i++ ) {
row = ELM_PLIST(res, i);
SET_LEN_PLIST(row, len);
SHRINK_PLIST(row, len); if (!rmut)
MakeImmutableNoRecurse(row);
} if (!cmut)
MakeImmutableNoRecurse(res);
/**************************************************************************** ** *F FuncADD_ROW_VECTOR_5( <self>, <list1>, <list2>, <mult>, <from>, <to> ) ** ** This function adds <mult>*<list2>[i] destructively to <list1>[i] for ** each i in the range <from>..<to>. It does very little checking **
*/
// We need these to redispatch when the user has supplied a replacement value.
static Obj FuncADD_ROW_VECTOR_5(
Obj self, Obj list1, Obj list2, Obj mult, Obj from, Obj to)
{ Int ifrom = GetSmallInt("AddRowVector", from); Int ito = GetSmallInt("AddRowVector", to); if (ito > LEN_LIST(list1) || ito > LEN_LIST(list2))
ErrorMayQuit("AddRowVector: Upper limit too large", 0, 0); for (Int i = ifrom; i <= ito; i++) {
Obj el1 = ELM_LIST(list1, i);
Obj el2 = ELM_LIST(list2, i);
el2 = PROD(mult, el2);
el1 = SUM(el1, el2);
ASS_LIST(list1, i, el1);
CHANGED_BAG(list1);
} return 0;
}
/**************************************************************************** ** *F FuncADD_ROW_VECTOR_5_FAST(<self>,<list1>,<list2>,<mult>,<from>,<to>) ** ** This function adds <mult>*<list2>[i] destructively to <list1>[i] for ** each i in the range <from>..<to>. It does very little checking ** ** This version is specialised to the "fast" case where list1 and list2 are ** plain lists of cyclotomics and mult is a small integer.
*/ static Obj FuncADD_ROW_VECTOR_5_FAST(
Obj self, Obj list1, Obj list2, Obj mult, Obj from, Obj to)
{ Int ifrom = GetSmallInt("AddRowVector", from); Int ito = GetSmallInt("AddRowVector", to); if (ito > LEN_LIST(list1) || ito > LEN_LIST(list2))
ErrorMayQuit("AddRowVector: Upper limit too large", 0, 0);
Obj prd, sum; for (Int i = ifrom; i <= ito; i++) {
Obj e1 = ELM_PLIST(list1, i);
Obj e2 = ELM_PLIST(list2, i); if (!ARE_INTOBJS(e2, mult) || !PROD_INTOBJS(prd, e2, mult)) {
prd = PROD(e2, mult);
} if (!ARE_INTOBJS(e1, prd) || !SUM_INTOBJS(sum, e1, prd)) {
sum = SUM(e1, prd);
SET_ELM_PLIST(list1, i, sum);
CHANGED_BAG(list1);
} else
SET_ELM_PLIST(list1, i, sum);
} return 0;
}
/**************************************************************************** ** *F FuncADD_ROW_VECTOR_3( <self>, <list1>, <list2>, <mult> ) ** ** This function adds <mult>*<list2>[i] destructively to <list1>[i] for ** each i in the range 1..Length(<list1>). It does very little checking ** *T This could be speeded up still further by using special code for various ** types of list -- this version just uses generic list ops
*/ static Obj FuncADD_ROW_VECTOR_3(Obj self, Obj list1, Obj list2, Obj mult)
{
UInt i;
UInt len = LEN_LIST(list1);
Obj el1, el2;
RequireSameLength(SELF_NAME, list1, list2); for (i = 1; i <= len; i++)
{
el1 = ELMW_LIST(list1,i);
el2 = ELMW_LIST(list2,i);
el2 = PROD(mult, el2);
el1 = SUM(el1 , el2);
ASS_LIST(list1,i,el1);
CHANGED_BAG(list1);
} return 0;
}
/**************************************************************************** ** *F FuncADD_ROW_VECTOR_3_FAST( <self>, <list1>, <list2>, <mult> ) ** ** This function adds <mult>*<list2>[i] destructively to <list1>[i] for ** each i in the range 1..Length(<list1>). It does very little checking ** ** This version is specialised to the "fast" case where list1 and list2 are ** plain lists of cyclotomics and mult is a small integers
*/ static Obj FuncADD_ROW_VECTOR_3_FAST(Obj self, Obj list1, Obj list2, Obj mult)
{
UInt i;
Obj e1,e2, prd, sum;
UInt len = LEN_PLIST(list1);
RequireSameLength(SELF_NAME, list1, list2); for (i = 1; i <= len; i++)
{
e1 = ELM_PLIST(list1,i);
e2 = ELM_PLIST(list2,i); if ( !ARE_INTOBJS( e2, mult ) || !PROD_INTOBJS( prd, e2, mult ))
{
prd = PROD(e2,mult);
} if ( !ARE_INTOBJS(e1, prd) || !SUM_INTOBJS( sum, e1, prd) )
{
sum = SUM(e1,prd);
SET_ELM_PLIST(list1,i,sum);
CHANGED_BAG(list1);
} else
SET_ELM_PLIST(list1,i,sum);
} return 0;
}
/**************************************************************************** ** *F FuncADD_ROW_VECTOR_2( <self>, <list1>, <list2>) ** ** This function adds <list2>[i] destructively to <list1>[i] for ** each i in the range 1..Length(<list1>). It does very little checking ** *T This could be speeded up still further by using special code for various ** types of list -- this version just uses generic list ops
*/ static Obj FuncADD_ROW_VECTOR_2(Obj self, Obj list1, Obj list2)
{
UInt i;
Obj el1,el2;
UInt len = LEN_LIST(list1);
RequireSameLength(SELF_NAME, list1, list2); for (i = 1; i <= len; i++)
{
el1 = ELMW_LIST(list1,i);
el2 = ELMW_LIST(list2,i);
el1 = SUM(el1, el2 );
ASS_LIST(list1,i,el1);
CHANGED_BAG(list1);
} return 0;
}
/**************************************************************************** ** *F FuncADD_ROW_VECTOR_2_FAST( <self>, <list1>, <list2> ) ** ** This function adds <list2>[i] destructively to <list1>[i] for ** each i in the range 1..Length(<list1>). It does very little checking ** ** This version is specialised to the "fast" case where list1 and list2 are ** plain lists of cyclotomics
*/ static Obj FuncADD_ROW_VECTOR_2_FAST(Obj self, Obj list1, Obj list2)
{
UInt i;
Obj e1,e2, sum;
UInt len = LEN_PLIST(list1);
RequireSameLength(SELF_NAME, list1, list2); for (i = 1; i <= len; i++)
{
e1 = ELM_PLIST(list1,i);
e2 = ELM_PLIST(list2,i); if ( !ARE_INTOBJS(e1, e2) || !SUM_INTOBJS( sum, e1, e2) )
{
sum = SUM(e1,e2);
SET_ELM_PLIST(list1,i,sum);
CHANGED_BAG(list1);
} else
SET_ELM_PLIST(list1,i,sum);
} return 0;
}
/**************************************************************************** ** *F MULT_VECTOR_LEFT_RIGHT_2( <list>, <mult>, <left> ) ** ** This function destructively multiplies the entries of <list> by <mult>. ** It multiplies with <mult> from the left if <left> is not 0 and from the ** right otherwise. ** It does very little checking. **
*/ staticinline Obj MULT_VECTOR_LEFT_RIGHT_2(Obj list, Obj mult, UInt left)
{
UInt i;
Obj prd;
UInt len = LEN_LIST(list); if (left != 0) for (i = 1; i <= len; i++) {
prd = ELMW_LIST(list, i);
prd = PROD(mult, prd);
ASS_LIST(list, i, prd);
CHANGED_BAG(list);
} else for (i = 1; i <= len; i++) {
prd = ELMW_LIST(list, i);
prd = PROD(prd, mult);
ASS_LIST(list, i, prd);
CHANGED_BAG(list);
} return 0;
}
/**************************************************************************** ** *F FuncMULT_VECTOR_2_FAST( <self>, <list>, <mult> ) ** ** This function destructively multiplies the entries of <list> by <mult> ** It does very little checking ** ** This is the fast method for plain lists of cyclotomics and an integer ** multiplier
*/
static Obj FuncMULT_VECTOR_2_FAST(Obj self, Obj list, Obj mult)
{
UInt i;
Obj el,prd;
UInt len = LEN_PLIST(list); for (i = 1; i <= len; i++)
{
el = ELM_PLIST(list,i); if (!ARE_INTOBJS(el, mult) || !PROD_INTOBJS(prd,el,mult))
{
prd = PROD(el,mult);
SET_ELM_PLIST(list,i,prd);
CHANGED_BAG(list);
} else
SET_ELM_PLIST(list,i,prd);
} return 0;
}
/**************************************************************************** ** *F FuncPROD_VEC_MAT_DEFAULT( <self>, <vec>, <mat> ) ** ** This is a specialized version of PROD_LIST_LIST_DEFAULT, that uses ** AddRowVector rather than SUM and PROD.
*/
static Obj FuncPROD_VEC_MAT_DEFAULT(Obj self, Obj vec, Obj mat)
{
Obj res;
Obj elt;
Obj vecr;
UInt i,len;
Obj z;
res = (Obj) 0;
len = LEN_LIST(vec);
RequireSameLength(" * ", vec, mat);
elt = ELMW_LIST(vec,1);
z = ZERO_SAMEMUT(elt); for (i = 1; i <= len; i++)
{
elt = ELMW_LIST(vec,i); if (!EQ(elt,z))
{
vecr = ELMW_LIST(mat,i); if (res == (Obj)0)
{
res = SHALLOW_COPY_OBJ(vecr);
CALL_2ARGS(MultVectorLeftOp, res, elt);
} else
CALL_3ARGS(AddRowVectorOp, res, vecr, elt);
}
} if (res == (Obj)0)
res = ZERO_SAMEMUT(ELMW_LIST(mat,1)); if (!IS_MUTABLE_OBJ(vec) && !IS_MUTABLE_OBJ(mat))
CheckedMakeImmutable(res); return res;
}
/**************************************************************************** ** *F FuncINV_MAT_DEFAULT ** ** A faster version of InvMat for those matrices for whose rows AddRowVector ** and MultVectorLeft make sense (and might have fast kernel methods) **
*/
static Obj ConvertToMatrixRep;
static Obj InvMatWithRowVecs(Obj mat, UInt mut)
{
Obj res; // result
Obj matcopy; // copy of mat
Obj row = 0; // one row of matcopy
Obj row2; // corresponding row of res
Obj row3; // another row of matcopy
Obj x = 0; // one element of the matrix
Obj xi; // 1/x
Obj y; // another element of the matrix
Obj yi; // -y
Obj zero; // zero element
Obj zerov; // zero vector
Obj one; // one element
UInt len; // length (and width) of matrix
UInt i, k, j; // loop variables
// check that the operand is a *square* matrix
len = LEN_LIST( mat ); if ( len != LEN_LIST( ELM_LIST( mat, 1 ) ) ) {
ErrorMayQuit("Matrix INV: must be square (not %d by %d)",
(Int)len, (Int)LEN_LIST(ELM_LIST(mat, 1)));
}
// get the zero and the one
zerov = ZERO_SAMEMUT(ELMW_LIST(mat, 1));
zero = ZERO_SAMEMUT(ELMW_LIST(ELMW_LIST(mat, 1), 1));
one = ONE( zero );
// set up res (initially the identity) and matcopy
res = NEW_PLIST(T_PLIST,len);
matcopy = NEW_PLIST(T_PLIST,len);
SET_LEN_PLIST(res,len);
SET_LEN_PLIST(matcopy,len); for (i = 1; i <= len; i++)
{
row = SHALLOW_COPY_OBJ(zerov);
ASS_LIST(row,i,one);
SET_ELM_PLIST(res,i,row);
CHANGED_BAG(res);
row = SHALLOW_COPY_OBJ(ELM_LIST(mat,i));
SET_ELM_PLIST(matcopy,i,row);
CHANGED_BAG(matcopy);
}
/* Now to work, make matcopy an identity by row operations and
do the same row operations to res */
// outer loop over columns of matcopy for (i = 1; i <= len; i++)
{ // Find a non-zero leading entry that is in that column for (j = i; j <= len; j++)
{
row = ELM_PLIST(matcopy,j);
x = ELMW_LIST(row,i); if (!EQ(x,zero)) break;
}
// if there isn't one then the matrix is not invertible if (j > len) return Fail;
// Maybe swap two rows // But I will want this value anyway
row2 = ELM_PLIST(res,j); if (j != i)
{
SET_ELM_PLIST(matcopy,j,ELM_PLIST(matcopy,i));
SET_ELM_PLIST(res,j,ELM_PLIST(res,i));
SET_ELM_PLIST(matcopy,i,row);
SET_ELM_PLIST(res,i,row2);
}
/*Maybe rescale the row */ if (!EQ(x, one))
{
xi = INV(x);
CALL_2ARGS(MultVectorLeftOp, row, xi);
CALL_2ARGS(MultVectorLeftOp, row2, xi);
}
// Clear the entries. We know that we can ignore the entries in rows i..j for (k = 1; k < i; k++)
{
row3 = ELM_PLIST(matcopy,k);
y = ELMW_LIST(row3,i); if (!EQ(y,zero))
{
yi = AINV_SAMEMUT(y);
CALL_3ARGS(AddRowVectorOp, row3, row, yi);
CALL_3ARGS(AddRowVectorOp, ELM_PLIST(res,k), row2, yi);
}
} for (k = j+1; k <= len; k++)
{
row3 = ELM_PLIST(matcopy,k);
y = ELMW_LIST(row3,i); if (!EQ(y,zero))
{
yi = AINV_SAMEMUT(y);
CALL_3ARGS(AddRowVectorOp, row3, row, yi);
CALL_3ARGS(AddRowVectorOp, ELM_PLIST(res,k), row2, yi);
}
}
}
// Now we adjust mutability. Couldn't do it earlier, because AddRowVector, // etc. needs mutable target vectors switch (mut)
{ case 0:
CheckedMakeImmutable(res); break;
case 1: if (IS_MUTABLE_OBJ(mat))
{ if (!IS_MUTABLE_OBJ(ELM_LIST(mat,1))) for (i = 1; i <= len; i++)
CheckedMakeImmutable(ELM_LIST(res,i));
} else
CheckedMakeImmutable(res); break; case 2: break;
}
/**************************************************************************** ** *F FuncADD_TO_LIST_ENTRIES_PLIST_RANGE( <list>, <range>, <x> ) ** ** This is a method for the operation AddToListEntries, used mainly in ** the MPQS of Stefan Kohl. ** ** This method requires a plain list for the first argument, and a ** range for the second and is optimised for the case where x and the list ** entries are small integers, as are their sums
*/
static Obj
FuncADD_TO_LIST_ENTRIES_PLIST_RANGE(Obj self, Obj list, Obj range, Obj x)
{
UInt low, high, incr;
Obj y,z;
UInt i; if (!IS_INTOBJ(x)) return TRY_NEXT_METHOD;
low = GET_LOW_RANGE(range);
incr = GET_INC_RANGE(range);
high = low + incr*(GET_LEN_RANGE(range)-1); for (i = low; i <= high; i+= incr)
{
y = ELM_PLIST(list,i); if (!IS_INTOBJ(y) ||
!SUM_INTOBJS(z,x,y))
{
z = SUM(x,y);
SET_ELM_PLIST(list,i,z);
CHANGED_BAG(list);
} else
SET_ELM_PLIST(list,i,z);
} return (Obj) 0;
}
/**************************************************************************** ** *F MONOM_TOT_DEG_LEX( u, v ) . . . . . total degree lexicographical ordering ** ** This function implements the total degree plus lexicographical ordering ** for monomials of commuting indeterminates. It is in this file because ** monomials are presently implemented as lists of indeterminate-exponent ** pairs. Should there be more functions supporting polynomial arithmetic, ** then this function should go into a separate file. ** ** Examples: x^2y^3 < y^7, x^4 y^5 < x^3 y^6
*/ static Obj FuncMONOM_TOT_DEG_LEX ( Obj self, Obj u, Obj v ) {
Int4 i, lu, lv;
Obj total;
Obj lexico;
if (!IS_PLIST(u) || !IS_DENSE_LIST(u)) {
RequireArgument(SELF_NAME, u, "must be a dense plain list");
} if (!IS_PLIST(v) || !IS_DENSE_LIST(v)) {
RequireArgument(SELF_NAME, v, "must be a dense plain list");
}
lu = LEN_PLIST( u );
lv = LEN_PLIST( v );
// strip off common prefixes
i = 1; while ( i <= lu && i <= lv && EQ(ELM_PLIST( u,i ), ELM_PLIST( v,i )) )
i++;
// Is u a prefix of v ? Return true if u is a proper prefix. if ( i > lu ) return (lu < lv) ? True : False;
// Is v a prefix of u ? if ( i > lv ) returnFalse;
// Now determine the lexicographic order. The monomial is interpreted // as a string of indeterminates. if ( i % 2 == 1 ) { // The first difference between u and v is an indeterminate.
lexico = LT(ELM_PLIST( u, i ), ELM_PLIST( v, i )) ? True : False;
i++;
} else { // The first difference between u and v is an exponent.
lexico = LT(ELM_PLIST( v, i ), ELM_PLIST( u, i )) ? True : False;
}
// Now add up the remaining exponents in order to compare the total // degrees.
total = INTOBJ_INT(0); while ( i <= lu && i <= lv ) {
C_SUM_FIA( total, total, ELM_PLIST( u, i ) );
C_DIFF_FIA( total, total, ELM_PLIST( v, i ) );
i += 2;
}
// Only one of the following while loops is executed while ( i <= lu ) {
C_SUM_FIA( total, total, ELM_PLIST( u, i ) );
i += 2;
} while ( i <= lv ) {
C_DIFF_FIA( total, total, ELM_PLIST( v, i ) );
i += 2;
}
/**************************************************************************** ** *F MONOM_GRLEX( u, v ) . . . . . ``grlex'' ordering for internal monomials ** ** This function implements the ``grlex'' (degree, then lexicographic) ordering ** for monomials of commuting indeterminates with x_1>x_2>x_3 etc. (this ** is standard textbook usage). It is in this file because ** monomials are presently implemented as lists of indeterminate-exponent ** pairs. Should there be more functions supporting polynomial arithmetic, ** then this function should go into a separate file. ** ** Examples: x^2y^3 < y^7, x^4 y^5 < x^3 y^6
*/ static Obj FuncMONOM_GRLEX( Obj self, Obj u, Obj v ) {
Int4 i, lu, lv;
Obj total,ai,bi;
if (!IS_PLIST(u) || !IS_DENSE_LIST(u)) {
RequireArgument(SELF_NAME, u, "must be a dense plain list");
} if (!IS_PLIST(v) || !IS_DENSE_LIST(v)) {
RequireArgument(SELF_NAME, v, "must be a dense plain list");
}
lu = LEN_PLIST( u );
lv = LEN_PLIST( v );
// compare the total degrees
total = INTOBJ_INT(0); for (i=2;i<=lu;i+=2) {
C_SUM_FIA( total, total, ELM_PLIST( u, i ) );
}
for (i=2;i<=lv;i+=2) {
C_DIFF_FIA( total, total, ELM_PLIST( v, i ) );
}
if ( ! (EQ( total, INTOBJ_INT(0))) ) { // degrees differ, use these return LT( total, INTOBJ_INT(0)) ? True : False;
}
// now use lexicographic ordering
i=1; while (i<=lu && i<=lv) {
ai=ELM_PLIST(u,i);
bi=ELM_PLIST(v,i); if (LT(bi,ai)) { returnTrue;
} if (LT(ai,bi)) { returnFalse;
}
ai=ELM_PLIST(u,i+1);
bi=ELM_PLIST(v,i+1); if (LT(ai,bi)) { returnTrue;
} if (LT(bi,ai)) { returnFalse;
}
i+=2;
} if (i<lv) { returnTrue;
} returnFalse;
}
/**************************************************************************** ** *F ZIPPED_SUM_LISTS(z1,z2,zero,f) ** ** implements the `ZippedSum' function to add polynomials in external ** representation. This is time critical and thus in the kernel. ** the function assumes that all lists are plists.
*/ static Obj FuncZIPPED_SUM_LISTS( Obj self, Obj z1, Obj z2, Obj zero, Obj f ) {
Int l1,l2,i; Int i1,i2;
Obj sum,x,y;
Obj cmpfun,sumfun,a,b,c;
l1=LEN_LIST(z1);
l2=LEN_LIST(z2);
cmpfun=ELM_LIST(f,1);
sumfun=ELM_LIST(f,2);
sum=NEW_PLIST(T_PLIST,0);
i1=1;
i2=1; while ((i1<=l1) && (i2<=l2)) { // Pr("A= %d %d\n",i1,i2); // if z1[i1] = z2[i2] then
a=ELM_PLIST(z1,i1);
b=ELM_PLIST(z2,i2); if (EQ(a,b)) { // the entries are equal
x=ELM_PLIST(z1,i1+1);
y=ELM_PLIST(z2,i2+1); // Pr("EQ, %d %d\n",INT_INTOBJ(x),INT_INTOBJ(y));
c=CALL_2ARGS(sumfun,x,y); if (!(EQ(c,zero))) { // Pr("Added %d\n",INT_INTOBJ(c), 0);
AddList(sum,a);
AddList(sum,c);
}
i1=i1+2;
i2=i2+2;
} else { // compare
a=ELM_PLIST(z1,i1); // in case the EQ triggered a GC
b=ELM_PLIST(z2,i2); // Pr("B= %d %d\n",ELM_LIST(a,1),ELM_LIST(b,1));
c=CALL_2ARGS(cmpfun,a,b); // Pr("C= %d %d\n",c, 0);
if ( // this construct is taken from the compiler
(Obj)(UInt)(c != False) ) {
a=ELM_PLIST(z1,i1);
AddList(sum,a);
c=ELM_PLIST(z1,i1+1);
AddList(sum,c);
i1=i1+2;
} else {
b=ELM_PLIST(z2,i2);
AddList(sum,b);
c=ELM_PLIST(z2,i2+1);
AddList(sum,c);
i2=i2+2;
} // else
} /*else (elif)*/
} // while
for (i=i1;i<l1;i+=2) {
AddList(sum,ELM_PLIST(z1,i));
AddList(sum,ELM_PLIST(z1,i+1));
}
for (i=i2;i<l2;i+=2) {
AddList(sum,ELM_PLIST(z2,i));
AddList(sum,ELM_PLIST(z2,i+1));
} return sum;
}
/**************************************************************************** ** *F FuncMONOM_PROD(m1,m2) ** ** implements the multiplication of monomials. Both must be plain lists ** of integers.
*/ static Obj FuncMONOM_PROD( Obj self, Obj m1, Obj m2 ) {
// No kernel installations for One or Inverse any more
/* Sum. Here we can do list + non-list and non-list + list, we have to careful about list+list, though, as it might really be list+matrix or matrix+list
for the T_PLIST_CYC and T_PLIST_FFE cases, we know the nesting depth (1) and so we can do the cases of adding them to each other and to T_PLIST_TAB objects, which have at least nesting depth 2. Some of this will be overwritten in vector.c and vecffe.c
everything else needs to wait until the library */
