My Project
sispasm.cc
Go to the documentation of this file.
1 /*
2  * provides (after defintion of a ring with coeffs in Z/p)
3  * - type spasm
4  * - assignment smatrix ->spasm, matrix ->spasm
5  * - printing/string(spasm)
6  * - transpose(spasm) -> spasm
7  * - to_matrix(spams) -> matrix
8  * - to_smatrix(spasm) -> smatrix
9  * - spasm_kernel(spasm)->spasm
10  * - spasm_rref(spasm) -> spasm
11 */
12 #include "singularconfig.h"
14 #include "kernel/ideals.h"
15 #include "Singular/ipid.h"
16 #include "Singular/ipshell.h"
17 #include "Singular/blackbox.h"
18 #include "Singular/mod_lib.h"
19 #ifdef HAVE_SPASM_H
20 extern "C"
21 {
22 #include "spasm.h"
23 }
24 
25 spasm* conv_matrix2spasm(matrix M, const ring R)
26 {
27  int i=MATROWS(M);
28  int j=MATCOLS(M);
29  spasm_triplet *T = spasm_triplet_alloc(i, j, 1, R->cf->ch, 1);
30  for (int ii=1;ii<=i;ii++)
31  {
32  for(int jj=1;jj<=j;jj++)
33  {
34  poly p;
35  if ((p=MATELEM(M,ii,jj))!=NULL)
36  {
37  spasm_add_entry(T,ii-1,jj-1,(spasm_GFp)(long)pGetCoeff(p));
38  }
39  }
40  }
41  spasm* A=spasm_compress(T);
42  spasm_triplet_free(T);
43  return A;
44 }
45 
46 spasm* conv_smatrix2spasm(ideal M, const ring R)
47 {
48  int i=MATROWS((matrix)M);
49  int j=MATCOLS((matrix)M);
50  spasm_triplet *T = spasm_triplet_alloc(i, j, 1, R->cf->ch, 1);
51  for(int jj=0;jj<j;jj++)
52  {
53  poly p=M->m[jj];
54  while (p!=NULL)
55  {
56  int ii=p_GetComp(p,R);
57  spasm_add_entry(T,ii-1,jj,(spasm_GFp)(long)pGetCoeff(p));
58  pIter(p);
59  }
60  }
61  spasm* A=spasm_compress(T);
62  spasm_triplet_free(T);
63  return A;
64 }
65 
66 matrix conv_spasm2matrix(spasm *A, const ring R)
67 {
68  matrix M=mpNew(A->n,A->m);
69  int n=A->n;
70  int *Aj = A->j;
71  int *Ap = A->p;
72  spasm_GFp *Ax = A->x;
73  for (int i = 0; i < n; i++)
74  {
75  for (int px = Ap[i]; px < Ap[i + 1]; px++)
76  {
77  spasm_GFp x = (Ax != NULL) ? Ax[px] : 1;
78  MATELEM(M,i+1,Aj[px] + 1)=p_ISet(x,R);
79  }
80  }
81  return M;
82 }
83 
84 ideal conv_spasm2smatrix(spasm *A, const ring R)
85 {
86  ideal M=idInit(A->m,A->n);
87  int n=A->n;
88  int *Aj = A->j;
89  int *Ap = A->p;
90  spasm_GFp *Ax = A->x;
91  for (int i = 0; i < n; i++)
92  {
93  for (int px = Ap[i]; px < Ap[i + 1]; px++)
94  {
95  spasm_GFp x = (Ax != NULL) ? Ax[px] : 1;
96  poly p=p_ISet(x,R);
97  p_SetComp(p,i+1,R);p_SetmComp(p,R);
98  M->m[Aj[px]]=p_Add_q(M->m[Aj[px]],p,R);
99  }
100  }
101  return M;
102 }
103 
104 spasm* sp_kernel(spasm* A, const ring R)
105 { /*from kernel.c*/
106  int n = A->n;
107  int m = A->m;
108  int* p = (int*)spasm_malloc(n * sizeof(int));
109  int * qinv = (int*)spasm_malloc(m * sizeof(int));
110  spasm_find_pivots(A, p, qinv); /* this does some useless stuff, but
111  * pushes zero rows to the bottom */
112  spasm* A_clean = spasm_permute(A, p, SPASM_IDENTITY_PERMUTATION, SPASM_WITH_NUMERICAL_VALUES);
113  spasm_csr_free(A);
114  A = A_clean;
115  for (int i = 0; i < n; i++)
116  {
117  if (spasm_row_weight(A, i) == 0)
118  {
119  //fprintf(stderr, "[kernel] ignoring %d empty rows\n", n - i);
120  A->n = i;
121  n = i;
122  break;
123  }
124  }
125 
126  spasm* A_t = spasm_transpose(A, SPASM_WITH_NUMERICAL_VALUES);
127  spasm_find_pivots(A_t, qinv, p);
128 
129  spasm* K = spasm_kernel(A_t, qinv);
130  free(p);
131  free(qinv);
132  spasm_csr_free(A_t);
133  return K;
134 }
135 
136 spasm* sp_rref(spasm* A)
137 { /* from rref_gplu.c: compute an echelonized form, WITHOUT COLUMN PERMUTATION */
138  spasm_lu *LU = spasm_LU(A, SPASM_IDENTITY_PERMUTATION, 1);
139  spasm *U = spasm_transpose(LU->L, 1);
140  spasm_make_pivots_unitary(U, SPASM_IDENTITY_PERMUTATION, U->n);
141  spasm_free_LU(LU);
142  return U;
143 }
144 
145 spasm* sp_Mult_v(spasm* A, int *v)
146 {
147  int *y=(int*)omAlloc0(A->n*sizeof(int));
148  spasm *AA=spasm_submatrix(A,0,A->n,0,A->m,1); /*copy A*/
149  spasm_gaxpy(AA,v,y);
150  return AA;
151 }
152 /*----------------------------------------------------------------*/
153 VAR int SPASM_CMD;
154 
155 static void* sp_Init(blackbox* /*b*/)
156 {
157  if ((currRing!=NULL)&&(rField_is_Zp(currRing)))
158  {
159  spasm_triplet *T = spasm_triplet_alloc(0, 0, 1, currRing->cf->ch, 1);
160  spasm* A=spasm_compress(T);
161  spasm_triplet_free(T);
162  return (void*)A;
163  }
164  else
165  {
166  WerrorS("ring with Z/p coeffs required");
167  return NULL;
168  }
169 }
170 static void sp_destroy(blackbox* /*b*/, void *d)
171 {
172  if (d!=NULL) spasm_csr_free((spasm*)d);
173  d=NULL;
174 }
175 static char* sp_String(blackbox* /*b*/, void *d)
176 {
177  char buf[30];
178  spasm* A=(spasm*)d;
179  sprintf(buf,"spasm matrix %dx%d",A->n,A->m);
180  return omStrDup(buf);
181 }
182 static void* sp_Copy(blackbox* /*b*/, void *d)
183 {
184  if (d!=NULL)
185  {
186  spasm* A=(spasm*)d;
187  spasm* B=spasm_submatrix(A,0,A->n,0,A->m,1);
188  return (void*)B;
189  }
190  return NULL;
191 }
192 static BOOLEAN sp_Assign(leftv l, leftv r)
193 {
194  spasm* A;
195  void*d=l->Data();
196  if (d!=NULL) spasm_csr_free((spasm*)d);
197  if (l->rtyp==IDHDL)
198  {
199  IDDATA((idhdl)l->data) = (char*)NULL;
200  }
201  else
202  {
203  l->data = (void*)NULL;
204  }
205 
206  if (r->Typ()==l->Typ())
207  {
208  A=(spasm*)sp_Copy(NULL,r->Data());
209  }
210  else if (r->Typ()==SMATRIX_CMD)
211  {
212  A=conv_smatrix2spasm((ideal)r->Data(),currRing);
213  }
214  else if (r->Typ()==MATRIX_CMD)
215  {
216  A=conv_matrix2spasm((matrix)r->Data(),currRing);
217  }
218  else
219  return TRUE;
220 
221  if (l->rtyp==IDHDL)
222  {
223  IDDATA((idhdl)l->data) = (char*)A;
224  }
225  else
226  {
227  l->data = (void*)A;
228  }
229  return FALSE;
230 }
231 
232 static BOOLEAN to_smatrix(leftv res, leftv args)
233 {
234  leftv u = args;
235  if ((u!=NULL) && (u->Typ()==SPASM_CMD))
236  {
237  res->rtyp=SMATRIX_CMD;
238  res->data=(void*)conv_spasm2smatrix((spasm*)u->Data(),currRing);
239  return FALSE;
240  }
241  return TRUE;
242 }
243 static BOOLEAN to_matrix(leftv res, leftv args)
244 {
245  leftv u = args;
246  if ((u!=NULL) && (u->Typ()==SPASM_CMD))
247  {
248  res->rtyp=MATRIX_CMD;
249  res->data=(void*)conv_spasm2matrix((spasm*)u->Data(),currRing);
250  return FALSE;
251  }
252  return TRUE;
253 }
254 static BOOLEAN kernel(leftv res, leftv args)
255 {
256  leftv u = args;
257  if ((u!=NULL) && (u->Typ()==SPASM_CMD))
258  {
259  res->rtyp=SPASM_CMD;
260  res->data=(void*)sp_kernel((spasm*)u->Data(),currRing);
261  return FALSE;
262  }
263  return TRUE;
264 }
265 static BOOLEAN rref(leftv res, leftv args)
266 {
267  leftv u = args;
268  if ((u!=NULL) && (u->Typ()==SPASM_CMD))
269  {
270  res->rtyp=SPASM_CMD;
271  res->data=(void*)sp_rref((spasm*)u->Data());
272  return FALSE;
273  }
274  return TRUE;
275 }
276 static BOOLEAN sp_Op1(int op,leftv l, leftv r)
277 {
278  if(op==TRANSPOSE_CMD)
279  {
280  l->rtyp=r->Typ();
281  l->data=(void*)spasm_transpose((spasm*)r->Data(),SPASM_WITH_NUMERICAL_VALUES);
282  return FALSE;
283  }
284  return blackboxDefaultOp1(op,l,r);
285 }
286 /*----------------------------------------------------------------*/
287 // initialisation of the module
288 extern "C" int SI_MOD_INIT(sispasm)(SModulFunctions* p)
289 {
290  blackbox *b=(blackbox*)omAlloc0(sizeof(blackbox));
291  b->blackbox_destroy=sp_destroy;
292  b->blackbox_String=sp_String;
293  b->blackbox_Init=sp_Init;
294  b->blackbox_Copy=sp_Copy;
295  b->blackbox_Assign=sp_Assign;
296  b->blackbox_Op1=sp_Op1;
297  SPASM_CMD=setBlackboxStuff(b,"spasm");
298  p->iiAddCproc("spasm.so","spasm_kernel",FALSE,kernel);
299  p->iiAddCproc("spasm.so","spasm_rref",FALSE,rref);
300  p->iiAddCproc("spasm.so","to_smatrix",FALSE,to_smatrix);
301  p->iiAddCproc("spasm.so","to_matrix",FALSE,to_matrix);
302  return (MAX_TOK);
303 }
304 #else
305 extern "C" int SI_MOD_INIT(sispasm)(SModulFunctions* psModulFunctions)
306 {
307  PrintS("no spasm support\n");
308  return MAX_TOK;
309 }
310 #endif
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
int setBlackboxStuff(blackbox *bb, const char *n)
define a new type
Definition: blackbox.cc:142
BOOLEAN blackboxDefaultOp1(int op, leftv l, leftv r)
default procedure blackboxDefaultOp1, to be called as "default:" branch
Definition: blackbox.cc:78
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
Variable x
Definition: cfModGcd.cc:4082
int p
Definition: cfModGcd.cc:4078
CanonicalForm b
Definition: cfModGcd.cc:4103
Definition: idrec.h:35
Class used for (list of) interpreter objects.
Definition: subexpr.h:83
int Typ()
Definition: subexpr.cc:1011
void * Data()
Definition: subexpr.cc:1154
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:53
CanonicalForm res
Definition: facAbsFact.cc:60
b *CanonicalForm B
Definition: facBivar.cc:52
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
int j
Definition: facHensel.cc:110
void WerrorS(const char *s)
Definition: feFopen.cc:24
#define VAR
Definition: globaldefs.h:5
@ MATRIX_CMD
Definition: grammar.cc:286
@ SMATRIX_CMD
Definition: grammar.cc:291
#define IDDATA(a)
Definition: ipid.h:126
STATIC_VAR jList * T
Definition: janet.cc:30
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:37
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
#define MATROWS(i)
Definition: matpol.h:26
#define MATCOLS(i)
Definition: matpol.h:27
#define p_GetComp(p, r)
Definition: monomials.h:64
#define pIter(p)
Definition: monomials.h:37
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
#define omStrDup(s)
Definition: omAllocDecl.h:263
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define free
Definition: omAllocFunc.c:14
#define NULL
Definition: omList.c:12
poly p_ISet(long i, const ring r)
returns the poly representing the integer i
Definition: p_polys.cc:1297
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:938
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:249
#define p_SetmComp
Definition: p_polys.h:246
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
void PrintS(const char *s)
Definition: reporter.cc:284
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:501
int status int void * buf
Definition: si_signals.h:59
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
#define R
Definition: sirandom.c:27
#define A
Definition: sirandom.c:24
#define M
Definition: sirandom.c:25
int SI_MOD_INIT() sispasm(SModulFunctions *psModulFunctions)
Definition: sispasm.cc:305
#define IDHDL
Definition: tok.h:31
@ TRANSPOSE_CMD
Definition: tok.h:191
@ MAX_TOK
Definition: tok.h:218