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,module ->spasm, matrix ->spasm
5  * - printing/string(spasm)
6  * - transpose(spasm) -> spasm
7  * - nrows(spasm) -> int
8  * - ncols(spasm) -> int
9  * - to_matrix(spams) -> matrix
10  * - to_smatrix(spasm) -> smatrix
11  * - spasm_kernel(spasm)->spasm
12  * - spasm_rref(spasm) -> spasm
13  * - <spasm>[<int>,<int>] -> number: reading an entry (get_spasm_entry)
14 */
15 #include "singularconfig.h"
17 #include "kernel/ideals.h"
18 #include "Singular/ipid.h"
19 #include "Singular/ipshell.h"
20 #include "Singular/blackbox.h"
21 #include "Singular/mod_lib.h"
22 #ifdef HAVE_SPASM_H
23 extern "C"
24 {
25 #include "spasm.h"
26 }
27 
28 static spasm* conv_matrix2spasm(matrix M, const ring R)
29 {
30  int i=MATROWS(M);
31  int j=MATCOLS(M);
32  spasm_triplet *T = spasm_triplet_alloc(i, j, 1, R->cf->ch, 1);
33  for (int ii=1;ii<=i;ii++)
34  {
35  for(int jj=1;jj<=j;jj++)
36  {
37  poly p;
38  if ((p=MATELEM(M,ii,jj))!=NULL)
39  {
40  spasm_add_entry(T,ii-1,jj-1,(spasm_GFp)(long)pGetCoeff(p));
41  }
42  }
43  }
44  spasm* A=spasm_compress(T);
45  spasm_triplet_free(T);
46  return A;
47 }
48 
49 static spasm* conv_smatrix2spasm(ideal M, const ring R)
50 {
51  int i=MATROWS((matrix)M);
52  int j=MATCOLS((matrix)M);
53  spasm_triplet *T = spasm_triplet_alloc(i, j, 1, R->cf->ch, 1);
54  for(int jj=0;jj<j;jj++)
55  {
56  poly p=M->m[jj];
57  while (p!=NULL)
58  {
59  int ii=p_GetComp(p,R);
60  spasm_add_entry(T,ii-1,jj,(spasm_GFp)(long)pGetCoeff(p));
61  pIter(p);
62  }
63  }
64  spasm* A=spasm_compress(T);
65  spasm_triplet_free(T);
66  return A;
67 }
68 
69 static matrix conv_spasm2matrix(spasm *A, const ring R)
70 {
71  matrix M=mpNew(A->n,A->m);
72  int n=A->n;
73  int *Aj = A->j;
74  int *Ap = A->p;
75  spasm_GFp *Ax = A->x;
76  for (int i = 0; i < n; i++)
77  {
78  for (int px = Ap[i]; px < Ap[i + 1]; px++)
79  {
80  spasm_GFp x = (Ax != NULL) ? Ax[px] : 1;
81  MATELEM(M,i+1,Aj[px] + 1)=p_ISet(x,R);
82  }
83  }
84  return M;
85 }
86 
87 static ideal conv_spasm2smatrix(spasm *A, const ring R)
88 {
89  ideal M=idInit(A->m,A->n);
90  int n=A->n;
91  int *Aj = A->j;
92  int *Ap = A->p;
93  spasm_GFp *Ax = A->x;
94  for (int i = 0; i < n; i++)
95  {
96  for (int px = Ap[i]; px < Ap[i + 1]; px++)
97  {
98  spasm_GFp x = (Ax != NULL) ? Ax[px] : 1;
99  poly p=p_ISet(x,R);
100  p_SetComp(p,i+1,R);p_SetmComp(p,R);
101  M->m[Aj[px]]=p_Add_q(M->m[Aj[px]],p,R);
102  }
103  }
104  return M;
105 }
106 
107 static number get_spasm_entry(spasm *A, int i, int j, const ring R)
108 {
109  matrix M=mpNew(A->n,A->m);
110  int n=A->n;
111  int *Aj = A->j;
112  int *Ap = A->p;
113  i--;j--;
114  spasm_GFp *Ax = A->x;
115  if (i<n)
116  {
117  for (int px = Ap[i]; px < Ap[i + 1]; px++)
118  {
119  spasm_GFp x = (Ax != NULL) ? Ax[px] : 1;
120  if (j==Aj[px])
121  return n_Init(x,R->cf);
122  }
123  }
124  return n_Init(0,R->cf);;
125 }
126 
127 static spasm* sp_kernel(spasm* A, const ring R)
128 {
129  int n = A->n;
130  int m = A->m;
131  int* p = (int*)spasm_malloc(n * sizeof(int));
132  int * qinv = (int*)spasm_malloc(m * sizeof(int));
133  spasm_find_pivots(A, p, qinv); /* this does some useless stuff, but
134  * pushes zero rows to the bottom */
135 #if 0
136  /*from kernel.c*/
137  spasm* A_clean = spasm_permute(A, p, SPASM_IDENTITY_PERMUTATION, SPASM_WITH_NUMERICAL_VALUES);
138  A = A_clean;
139  for (int i = 0; i < n; i++)
140  {
141  if (spasm_row_weight(A, i) == 0)
142  {
143  //fprintf(stderr, "[kernel] ignoring %d empty rows\n", n - i);
144  A->n = i;
145  n = i;
146  break;
147  }
148  }
149 
150  spasm* A_t = spasm_transpose(A, SPASM_WITH_NUMERICAL_VALUES);
151  spasm_find_pivots(A_t, qinv, p);
152 
153  spasm* K = spasm_kernel(A_t, qinv);
154  spasm_csr_free(A_t);
155 #else
156  spasm* K = spasm_kernel(A, p);
157 #endif
158  free(p);
159  free(qinv);
160  return K;
161 }
162 
163 static spasm* sp_rref(spasm* A)
164 { /* from rref_gplu.c: compute an echelonized form, WITHOUT COLUMN PERMUTATION */
165  spasm_lu *LU = spasm_LU(A, SPASM_IDENTITY_PERMUTATION, 1);
166  spasm *U = spasm_transpose(LU->L, 1);
167  spasm_make_pivots_unitary(U, SPASM_IDENTITY_PERMUTATION, U->n);
168  spasm_free_LU(LU);
169  return U;
170 }
171 
172 static spasm* sp_Mult_v(spasm* A, int *v)
173 {
174  int *y=(int*)omAlloc0(A->n*sizeof(int));
175  spasm *AA=spasm_submatrix(A,0,A->n,0,A->m,1); /*copy A*/
176  spasm_gaxpy(AA,v,y);
177  return AA;
178 }
179 /*----------------------------------------------------------------*/
180 VAR int SPASM_CMD;
181 
182 static void* sp_Init(blackbox* /*b*/)
183 {
184  if ((currRing!=NULL)&&(rField_is_Zp(currRing)))
185  {
186  spasm_triplet *T = spasm_triplet_alloc(0, 0, 1, currRing->cf->ch, 1);
187  spasm* A=spasm_compress(T);
188  spasm_triplet_free(T);
189  return (void*)A;
190  }
191  else
192  {
193  WerrorS("ring with Z/p coeffs required");
194  return NULL;
195  }
196 }
197 static void sp_destroy(blackbox* /*b*/, void *d)
198 {
199  if (d!=NULL) spasm_csr_free((spasm*)d);
200  d=NULL;
201 }
202 static char* sp_String(blackbox* /*b*/, void *d)
203 {
204  char buf[30];
205  spasm* A=(spasm*)d;
206  sprintf(buf,"spasm matrix %dx%d",A->n,A->m);
207  return omStrDup(buf);
208 }
209 static void* sp_Copy(blackbox* /*b*/, void *d)
210 {
211  if (d!=NULL)
212  {
213  spasm* A=(spasm*)d;
214  spasm* B=spasm_submatrix(A,0,A->n,0,A->m,1);
215  return (void*)B;
216  }
217  return NULL;
218 }
219 static BOOLEAN sp_Assign(leftv l, leftv r)
220 {
221  spasm* A;
222  void*d=l->Data();
223  if (d!=NULL) spasm_csr_free((spasm*)d);
224  if (l->rtyp==IDHDL)
225  {
226  IDDATA((idhdl)l->data) = (char*)NULL;
227  }
228  else
229  {
230  l->data = (void*)NULL;
231  }
232  int rt=r->Typ();
233 
234  if (rt==l->Typ())
235  {
236  A=(spasm*)sp_Copy(NULL,r->Data());
237  }
238  else if ((rt==SMATRIX_CMD)||(rt==MODUL_CMD))
239  {
240  A=conv_smatrix2spasm((ideal)r->Data(),currRing);
241  }
242  else if (rt==MATRIX_CMD)
243  {
244  A=conv_matrix2spasm((matrix)r->Data(),currRing);
245  }
246  else
247  return TRUE;
248 
249  if (l->rtyp==IDHDL)
250  {
251  IDDATA((idhdl)l->data) = (char*)A;
252  }
253  else
254  {
255  l->data = (void*)A;
256  }
257  return FALSE;
258 }
259 
260 static BOOLEAN to_smatrix(leftv res, leftv args)
261 {
262  leftv u = args;
263  if ((u!=NULL) && (u->Typ()==SPASM_CMD))
264  {
265  res->rtyp=SMATRIX_CMD;
266  res->data=(void*)conv_spasm2smatrix((spasm*)u->Data(),currRing);
267  return FALSE;
268  }
269  return TRUE;
270 }
271 static BOOLEAN to_matrix(leftv res, leftv args)
272 {
273  leftv u = args;
274  if ((u!=NULL) && (u->Typ()==SPASM_CMD))
275  {
276  res->rtyp=MATRIX_CMD;
277  res->data=(void*)conv_spasm2matrix((spasm*)u->Data(),currRing);
278  return FALSE;
279  }
280  return TRUE;
281 }
282 static BOOLEAN kernel(leftv res, leftv args)
283 {
284  leftv u = args;
285  if ((u!=NULL) && (u->Typ()==SPASM_CMD))
286  {
287  res->rtyp=SPASM_CMD;
288  res->data=(void*)sp_kernel((spasm*)u->Data(),currRing);
289  return FALSE;
290  }
291  return TRUE;
292 }
293 static BOOLEAN rref(leftv res, leftv args)
294 {
295  leftv u = args;
296  if ((u!=NULL) && (u->Typ()==SPASM_CMD))
297  {
298  res->rtyp=SPASM_CMD;
299  res->data=(void*)sp_rref((spasm*)u->Data());
300  return FALSE;
301  }
302  return TRUE;
303 }
304 static BOOLEAN sp_Op1(int op,leftv res, leftv arg)
305 {
306  if(op==TRANSPOSE_CMD)
307  {
308  res->rtyp=arg->Typ();
309  res->data=(void*)spasm_transpose((spasm*)arg->Data(),SPASM_WITH_NUMERICAL_VALUES);
310  return FALSE;
311  }
312  else if (op==COLS_CMD)
313  {
314  spasm* A=(spasm*)arg->Data();
315  res->rtyp=INT_CMD;
316  res->data=(void*)(long)A->m;
317  return FALSE;
318  }
319  else if (op==ROWS_CMD)
320  {
321  spasm* A=(spasm*)arg->Data();
322  res->rtyp=INT_CMD;
323  res->data=(void*)(long)A->n;
324  return FALSE;
325  }
326  return blackboxDefaultOp1(op,res,arg);
327 }
328 static BOOLEAN sp_Op3(int op,leftv res, leftv a1, leftv a2, leftv a3)
329 {
330  if ((op=='[')
331  && (a2->Typ()==INT_CMD)
332  && (a3->Typ()==INT_CMD))
333  {
334  res->rtyp=NUMBER_CMD;
335  res->data=(char*)get_spasm_entry((spasm*)a1->Data(),(int)(long)a2->Data(),
336  (int)(long)a3->Data(),currRing);
337  return FALSE;
338  }
339  return blackboxDefaultOp3(op,res,a1,a2,a3);
340 }
341 /*----------------------------------------------------------------*/
342 // initialisation of the module
343 extern "C" int SI_MOD_INIT(sispasm)(SModulFunctions* p)
344 {
345  blackbox *b=(blackbox*)omAlloc0(sizeof(blackbox));
346  b->blackbox_destroy=sp_destroy;
347  b->blackbox_String=sp_String;
348  b->blackbox_Init=sp_Init;
349  b->blackbox_Copy=sp_Copy;
350  b->blackbox_Assign=sp_Assign;
351  b->blackbox_Op1=sp_Op1;
352  b->blackbox_Op3=sp_Op3;
353  SPASM_CMD=setBlackboxStuff(b,"spasm");
354  p->iiAddCproc("spasm.so","spasm_kernel",FALSE,kernel);
355  p->iiAddCproc("spasm.so","spasm_rref",FALSE,rref);
356  p->iiAddCproc("spasm.so","to_smatrix",FALSE,to_smatrix);
357  p->iiAddCproc("spasm.so","to_matrix",FALSE,to_matrix);
358  return (MAX_TOK);
359 }
360 #else
361 extern "C" int SI_MOD_INIT(sispasm)(SModulFunctions* psModulFunctions)
362 {
363  PrintS("no spasm support\n");
364  return MAX_TOK;
365 }
366 #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 blackboxDefaultOp3(int, leftv, leftv, leftv, leftv)
default procedure blackboxDefaultOp3, to be called as "default:" branch
Definition: blackbox.cc:102
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
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition: coeffs.h:535
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
@ MODUL_CMD
Definition: grammar.cc:287
@ SMATRIX_CMD
Definition: grammar.cc:291
@ NUMBER_CMD
Definition: grammar.cc:288
#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:934
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:245
#define p_SetmComp
Definition: p_polys.h:242
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:500
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:361
#define IDHDL
Definition: tok.h:31
@ TRANSPOSE_CMD
Definition: tok.h:191
@ INT_CMD
Definition: tok.h:96
@ ROWS_CMD
Definition: tok.h:175
@ MAX_TOK
Definition: tok.h:218
@ COLS_CMD
Definition: tok.h:52