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
20extern "C"
21{
22#include "spasm.h"
23}
24
25spasm* 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
46spasm* 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
66matrix 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
84ideal 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);
98 M->m[Aj[px]]=p_Add_q(M->m[Aj[px]],p,R);
99 }
100 }
101 return M;
102}
103
104spasm* 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
140spasm* 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
149spasm* 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/*----------------------------------------------------------------*/
157VAR int SPASM_CMD;
158
159static void* sp_Init(blackbox* /*b*/)
160{
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}
174static void sp_destroy(blackbox* /*b*/, void *d)
175{
176 if (d!=NULL) spasm_csr_free((spasm*)d);
177 d=NULL;
178}
179static 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}
186static 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}
196static 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
236static 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}
247static 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}
258static 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}
269static 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}
280static 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
292extern "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
309extern "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