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 {
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 #if 0
113  /*from kernel.c*/
114  spasm* A_clean = spasm_permute(A, p, SPASM_IDENTITY_PERMUTATION, SPASM_WITH_NUMERICAL_VALUES);
115  A = A_clean;
116  for (int i = 0; i < n; i++)
117  {
118  if (spasm_row_weight(A, i) == 0)
119  {
120  //fprintf(stderr, "[kernel] ignoring %d empty rows\n", n - i);
121  A->n = i;
122  n = i;
123  break;
124  }
125  }
126 
127  spasm* A_t = spasm_transpose(A, SPASM_WITH_NUMERICAL_VALUES);
128  spasm_find_pivots(A_t, qinv, p);
129 
130  spasm* K = spasm_kernel(A_t, qinv);
131  spasm_csr_free(A_t);
132 #else
133  spasm* K = spasm_kernel(A, p);
134 #endif
135  free(p);
136  free(qinv);
137  return K;
138 }
139 
140 spasm* sp_rref(spasm* A)
141 { /* from rref_gplu.c: compute an echelonized form, WITHOUT COLUMN PERMUTATION */
142  spasm_lu *LU = spasm_LU(A, SPASM_IDENTITY_PERMUTATION, 1);
143  spasm *U = spasm_transpose(LU->L, 1);
144  spasm_make_pivots_unitary(U, SPASM_IDENTITY_PERMUTATION, U->n);
145  spasm_free_LU(LU);
146  return U;
147 }
148 
149 spasm* sp_Mult_v(spasm* A, int *v)
150 {
151  int *y=(int*)omAlloc0(A->n*sizeof(int));
152  spasm *AA=spasm_submatrix(A,0,A->n,0,A->m,1); /*copy A*/
153  spasm_gaxpy(AA,v,y);
154  return AA;
155 }
156 /*----------------------------------------------------------------*/
157 VAR int SPASM_CMD;
158 
159 static void* sp_Init(blackbox* /*b*/)
160 {
161  if ((currRing!=NULL)&&(rField_is_Zp(currRing)))
162  {
163  spasm_triplet *T = spasm_triplet_alloc(0, 0, 1, currRing->cf->ch, 1);
164  spasm* A=spasm_compress(T);
165  spasm_triplet_free(T);
166  return (void*)A;
167  }
168  else
169  {
170  WerrorS("ring with Z/p coeffs required");
171  return NULL;
172  }
173 }
174 static void sp_destroy(blackbox* /*b*/, void *d)
175 {
176  if (d!=NULL) spasm_csr_free((spasm*)d);
177  d=NULL;
178 }
179 static char* sp_String(blackbox* /*b*/, void *d)
180 {
181  char buf[30];
182  spasm* A=(spasm*)d;
183  sprintf(buf,"spasm matrix %dx%d",A->n,A->m);
184  return omStrDup(buf);
185 }
186 static void* sp_Copy(blackbox* /*b*/, void *d)
187 {
188  if (d!=NULL)
189  {
190  spasm* A=(spasm*)d;
191  spasm* B=spasm_submatrix(A,0,A->n,0,A->m,1);
192  return (void*)B;
193  }
194  return NULL;
195 }
196 static BOOLEAN sp_Assign(leftv l, leftv r)
197 {
198  spasm* A;
199  void*d=l->Data();
200  if (d!=NULL) spasm_csr_free((spasm*)d);
201  if (l->rtyp==IDHDL)
202  {
203  IDDATA((idhdl)l->data) = (char*)NULL;
204  }
205  else
206  {
207  l->data = (void*)NULL;
208  }
209 
210  if (r->Typ()==l->Typ())
211  {
212  A=(spasm*)sp_Copy(NULL,r->Data());
213  }
214  else if (r->Typ()==SMATRIX_CMD)
215  {
216  A=conv_smatrix2spasm((ideal)r->Data(),currRing);
217  }
218  else if (r->Typ()==MATRIX_CMD)
219  {
220  A=conv_matrix2spasm((matrix)r->Data(),currRing);
221  }
222  else
223  return TRUE;
224 
225  if (l->rtyp==IDHDL)
226  {
227  IDDATA((idhdl)l->data) = (char*)A;
228  }
229  else
230  {
231  l->data = (void*)A;
232  }
233  return FALSE;
234 }
235 
236 static BOOLEAN to_smatrix(leftv res, leftv args)
237 {
238  leftv u = args;
239  if ((u!=NULL) && (u->Typ()==SPASM_CMD))
240  {
241  res->rtyp=SMATRIX_CMD;
242  res->data=(void*)conv_spasm2smatrix((spasm*)u->Data(),currRing);
243  return FALSE;
244  }
245  return TRUE;
246 }
247 static BOOLEAN to_matrix(leftv res, leftv args)
248 {
249  leftv u = args;
250  if ((u!=NULL) && (u->Typ()==SPASM_CMD))
251  {
252  res->rtyp=MATRIX_CMD;
253  res->data=(void*)conv_spasm2matrix((spasm*)u->Data(),currRing);
254  return FALSE;
255  }
256  return TRUE;
257 }
258 static BOOLEAN kernel(leftv res, leftv args)
259 {
260  leftv u = args;
261  if ((u!=NULL) && (u->Typ()==SPASM_CMD))
262  {
263  res->rtyp=SPASM_CMD;
264  res->data=(void*)sp_kernel((spasm*)u->Data(),currRing);
265  return FALSE;
266  }
267  return TRUE;
268 }
269 static BOOLEAN rref(leftv res, leftv args)
270 {
271  leftv u = args;
272  if ((u!=NULL) && (u->Typ()==SPASM_CMD))
273  {
274  res->rtyp=SPASM_CMD;
275  res->data=(void*)sp_rref((spasm*)u->Data());
276  return FALSE;
277  }
278  return TRUE;
279 }
280 static BOOLEAN sp_Op1(int op,leftv l, leftv r)
281 {
282  if(op==TRANSPOSE_CMD)
283  {
284  l->rtyp=r->Typ();
285  l->data=(void*)spasm_transpose((spasm*)r->Data(),SPASM_WITH_NUMERICAL_VALUES);
286  return FALSE;
287  }
288  return blackboxDefaultOp1(op,l,r);
289 }
290 /*----------------------------------------------------------------*/
291 // initialisation of the module
292 extern "C" int SI_MOD_INIT(sispasm)(SModulFunctions* p)
293 {
294  blackbox *b=(blackbox*)omAlloc0(sizeof(blackbox));
295  b->blackbox_destroy=sp_destroy;
296  b->blackbox_String=sp_String;
297  b->blackbox_Init=sp_Init;
298  b->blackbox_Copy=sp_Copy;
299  b->blackbox_Assign=sp_Assign;
300  b->blackbox_Op1=sp_Op1;
301  SPASM_CMD=setBlackboxStuff(b,"spasm");
302  p->iiAddCproc("spasm.so","spasm_kernel",FALSE,kernel);
303  p->iiAddCproc("spasm.so","spasm_rref",FALSE,rref);
304  p->iiAddCproc("spasm.so","to_smatrix",FALSE,to_smatrix);
305  p->iiAddCproc("spasm.so","to_matrix",FALSE,to_matrix);
306  return (MAX_TOK);
307 }
308 #else
309 extern "C" int SI_MOD_INIT(sispasm)(SModulFunctions* psModulFunctions)
310 {
311  PrintS("no spasm support\n");
312  return MAX_TOK;
313 }
314 #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:309
#define IDHDL
Definition: tok.h:31
@ TRANSPOSE_CMD
Definition: tok.h:191
@ MAX_TOK
Definition: tok.h:218