1 #include "misc/auxiliary.h"
2
3 #include "misc/intvec.h"
4 #include "misc/options.h"
5 #include "polys/monomials/p_polys.h"
6 #include "polys/matpol.h"
7 #include "polys/simpleideals.h"
8 #include "coeffs/longrat.h"
9 #include "Singular/feOpt.h"
10 #include "kernel/polys.h"
11 #include "kernel/mod2.h"
12
13 #ifdef HAVE_VSPACE
14
15 #define mpz_isNeg(A) ((A)->_mp_size<0)
16 number nlRInit (long i);
17
18 #include "kernel/oswrapper/vspace.h"
19 #include "kernel/ideals.h"
20 #include <sys/types.h>
21 #include <sys/wait.h>
22
23 // send a number via a string s
send_number(char * s,number n)24 static char *send_number(char *s, number n)
25 {
26 long *d=(long*)s;
27 if (SR_HDL(n) & SR_INT)
28 {
29 *d=(long)n;
30 s+=SIZEOF_LONG;
31 }
32 else
33 {
34 *d=n->s*2;/* n->s in 0..3: 0..6, use +8 for neg. numbers */
35 s+=SIZEOF_LONG;
36 if (mpz_isNeg(n->z)) { *d+=8; mpz_neg(n->z,n->z); }
37 size_t l;
38 d=(long*)s;
39 s+=SIZEOF_LONG;
40 mpz_export(s,&l,-1,sizeof(mp_limb_t),0,0,n->z);
41 *d=l;
42 s+=l*sizeof(mp_limb_t);
43 if (n->s!=3)
44 {
45 d=(long*)s;
46 s+=SIZEOF_LONG;
47 mpz_export(s,&l,-1,sizeof(mp_limb_t),0,0,n->n);
48 *d=l;
49 s+=l*sizeof(mp_limb_t);
50 }
51 }
52 return s;
53 }
54
get_number(char * s,number * n)55 static char * get_number(char *s, number *n)
56 {
57 // format: last bit 1: imm. number (long)
58 // otherwise: 0,2: size(long) mpz, size(long) mpz
59 // 6: size(long) mpz
60 // 8,10: size(long) -mpz, size(long) mpz
61 // 14: size(long) -mpz
62 long *d=(long*)s;
63 s+=SIZEOF_LONG;
64 if (((*d)&1)==1) // immidiate number
65 {
66 *n=(number)(*d);
67 }
68 else
69 {
70 *n=nlRInit(0);
71 BOOLEAN neg=(*d>=8);
72 if (neg) *d-=8;
73 (*n)->s=(*d)/2;
74 d=(long*)s;
75 s+=SIZEOF_LONG;
76 size_t l=*d;
77 mpz_realloc2((*n)->z,l*sizeof(mp_limb_t)*8);
78 mpz_import((*n)->z,l,-1,sizeof(mp_limb_t),0,0,s);
79 if (neg) mpz_neg((*n)->z,(*n)->z);
80 s+=l*sizeof(mp_limb_t);
81 if ((*n)->s!=3)
82 {
83 d=(long*)s;
84 s+=SIZEOF_LONG;
85 l=*d;
86 mpz_init2((*n)->n,l*sizeof(mp_limb_t)*8);
87 mpz_import((*n)->n,l,-1,sizeof(mp_limb_t),0,0,s);
88 s+=l*sizeof(mp_limb_t);
89 }
90 }
91 return s;
92 }
93
size_number(number n)94 static long size_number(number n)
95 {
96 long ll=SIZEOF_LONG;
97 if (SR_HDL(n) & SR_INT)
98 {
99 return SIZEOF_LONG;
100 }
101 else
102 {
103 if (n->s==3)
104 {
105 ll+=SIZEOF_LONG*2; /* n->s, mpz size */
106 long l=mpz_size1(n->z);
107 ll+=l*sizeof(mp_limb_t);
108 }
109 else
110 {
111 ll+=SIZEOF_LONG*3; /* n->s, mpz size(n->z) mpz size(n->n)*/
112 size_t l=mpz_size1(n->z);
113 ll+=l*sizeof(mp_limb_t);
114 l=mpz_size1(n->n);
115 ll+=l*sizeof(mp_limb_t);
116 }
117 }
118 return ll;
119 }
120
send_mon(char * s,poly m,const ring r)121 static char* send_mon(char *s, poly m, const ring r)
122 {
123 // format: number exp[0..r->ExpL_Size]
124 s=send_number(s,p_GetCoeff(m,r));
125 memcpy(s,m->exp,r->ExpL_Size*sizeof(long));
126 s+=r->ExpL_Size*sizeof(long);
127 return s;
128 }
129
get_mon(char * s,poly * m,const ring r)130 static char* get_mon(char *s, poly *m, const ring r)
131 {
132 (*m)=p_Init(r);
133 s=get_number(s,&p_GetCoeff(*m,r));
134 memcpy((*m)->exp,s,r->ExpL_Size*sizeof(long));
135 s+=r->ExpL_Size*sizeof(long);
136 return s;
137 }
138
size_mon(poly m,const ring r)139 static long size_mon(poly m, const ring r)
140 {
141 long ll=size_number(p_GetCoeff(m,r));
142 ll+=r->ExpL_Size*sizeof(long);
143 return ll;
144 }
145
send_poly(char * s,int ind,poly p,const ring r)146 static char* send_poly(char *s, int ind, poly p, const ring r)
147 {
148 // format: index(long) length(long) mon...
149 //p_Write(p,r);PrintLn();
150 long *d=(long*)s;
151 *d=ind;
152 s+=SIZEOF_LONG;
153 long l=pLength(p);
154 d=(long*)s;
155 *d=l;
156 s+=SIZEOF_LONG;
157 while(p!=NULL)
158 {
159 s=send_mon(s,p,r);
160 pIter(p);
161 }
162 return s;
163 }
164
get_poly(char * s,int & ind,poly * p,const ring r)165 static char* get_poly(char *s,int &ind, poly *p,const ring r)
166 {
167 long *d=(long*)s;
168 ind=*d;
169 s+=SIZEOF_LONG;
170 d=(long*)s;
171 long l=*d;
172 s+=SIZEOF_LONG;
173 for(long i=0;i<l;i++)
174 {
175 poly m;
176 s=get_mon(s,&m,r);
177 pNext(m)=*p;
178 *p=m;
179 }
180 *p=pReverse(*p);
181 return s;
182 }
183
size_poly(poly p,const ring r)184 static long size_poly(poly p, const ring r)
185 {
186 long l=SIZEOF_LONG*2;
187 while(p!=NULL)
188 {
189 l+=size_mon(p,r);
190 pIter(p);
191 }
192 return l;
193 }
194
id_ChineseRemainder_0(ideal * xx,number * q,int rl,const ring r)195 ideal id_ChineseRemainder_0(ideal *xx, number *q, int rl, const ring r)
196 {
197 int cnt=0;int rw=0; int cl=0;
198 // find max. size of xx[.]:
199 for(int j=rl-1;j>=0;j--)
200 {
201 int i=IDELEMS(xx[j])*xx[j]->nrows;
202 if (i>cnt) cnt=i;
203 if (xx[j]->nrows >rw) rw=xx[j]->nrows; // for lifting matrices
204 if (xx[j]->ncols >cl) cl=xx[j]->ncols; // for lifting matrices
205 }
206 if (rw*cl !=cnt)
207 {
208 WerrorS("format mismatch in CRT");
209 return NULL;
210 }
211 int cpus=(int)(long)feOptValue(FE_OPT_CPUS);
212 if ((cpus==1) || (2*cpus>=cnt))
213 /* at least 2 polys for each process, or switch to seriell version */
214 return id_ChineseRemainder(xx,q,rl,r);
215 ideal result=idInit(cnt,xx[0]->rank);
216 result->nrows=rw; // for lifting matrices
217 result->ncols=cl; // for lifting matrices
218 int parent_pid=getpid();
219 using namespace vspace;
220 vmem_init();
221 // Create a queue of int
222 VRef<Queue<int> > queue = vnew<Queue<int> >();
223 for(int i=cnt-1;i>=0; i--)
224 {
225 queue->enqueue(i); // the tasks: construct poly p[i]
226 }
227 for(int i=cpus;i>=0;i--)
228 {
229 queue->enqueue(-1); // stop sign, one for each child
230 }
231 // Create a queue of polys
232 VRef<Queue<VRef<VString> > > rqueue = vnew<Queue<VRef<VString> > >();
233 for (int i=0;i<cpus;i++)
234 {
235 int pid = fork_process();
236 if (pid==0) break; //child
237 }
238 if (parent_pid!=getpid()) // child ------------------------------------------
239 {
240 number *x=(number *)omAlloc(rl*sizeof(number));
241 poly *p=(poly *)omAlloc(rl*sizeof(poly));
242 CFArray inv_cache(rl);
243 EXTERN_VAR int n_SwitchChinRem;
244 n_SwitchChinRem=1;
245 loop
246 {
247 int ind=queue->dequeue();
248 if (ind== -1)
249 {
250 exit(0);
251 }
252
253 for(int j=rl-1;j>=0;j--)
254 {
255 if(ind>=IDELEMS(xx[j])*xx[j]->nrows) // out of range of this ideal
256 p[j]=NULL;
257 else
258 p[j]=xx[j]->m[ind];
259 }
260 poly res=p_ChineseRemainder(p,x,q,rl,inv_cache,r);
261 long l=size_poly(res,r);
262 //printf("size: %ld kB\n",(l+1023)/1024);
263 VRef<VString> msg = vstring(l+1);
264 char *s=(char*)msg->str();
265 send_poly(s,ind,res,r);
266 rqueue->enqueue(msg);
267 if (TEST_OPT_PROT) printf(".");
268 }
269 }
270 else // parent ---------------------------------------------------
271 {
272 if (TEST_OPT_PROT) printf("%d children created\n",cpus);
273 VRef<VString> msg;
274 while(cnt>0)
275 {
276 msg=rqueue->dequeue();
277 char *s=(char*)msg->str();
278 int ind;
279 poly p=NULL;
280 get_poly(s,ind,&p,r);
281 //printf("got res[%d]\n",ind);
282 result->m[ind]=p;
283 msg.free();
284 cnt--;
285 }
286 // removes queues
287 queue.free();
288 rqueue.free();
289 vmem_deinit();
290 }
291 return result;
292 }
293
id_Farey_0(ideal x,number N,const ring r)294 ideal id_Farey_0(ideal x, number N, const ring r)
295 {
296 int cnt=IDELEMS(x)*x->nrows;
297 int cpus=(int)(long)feOptValue(FE_OPT_CPUS);
298 if (2*cpus>=cnt) /* at least 2 polys for each process,
299 or switch to seriell version */
300 return id_Farey(x,N,r);
301 ideal result=idInit(cnt,x->rank);
302 result->nrows=x->nrows; // for lifting matrices
303 result->ncols=x->ncols; // for lifting matrices
304
305 int parent_pid=getpid();
306 using namespace vspace;
307 vmem_init();
308 // Create a queue of int
309 VRef<Queue<int> > queue = vnew<Queue<int> >();
310 for(int i=cnt-1;i>=0; i--)
311 {
312 queue->enqueue(i); // the tasks: construct poly p[i]
313 }
314 for(int i=cpus;i>=0;i--)
315 {
316 queue->enqueue(-1); // stop sign, one for each child
317 }
318 // Create a queue of polys
319 VRef<Queue<VRef<VString> > > rqueue = vnew<Queue<VRef<VString> > >();
320 for (int i=0;i<cpus;i++)
321 {
322 int pid = fork_process();
323 if (pid==0) break; //child
324 }
325 if (parent_pid!=getpid()) // child ------------------------------------------
326 {
327 loop
328 {
329 int ind=queue->dequeue();
330 if (ind== -1)
331 {
332 exit(0);
333 }
334
335 poly res=p_Farey(x->m[ind],N,r);
336 long l=size_poly(res,r);
337 VRef<VString> msg = vstring(l+1);
338 char *s=(char*)msg->str();
339 send_poly(s,ind,res,r);
340 rqueue->enqueue(msg);
341 if (TEST_OPT_PROT) printf(".");
342 }
343 }
344 else // parent ---------------------------------------------------
345 {
346 if (TEST_OPT_PROT) printf("%d children created\n",cpus);
347 VRef<VString> msg;
348 while(cnt>0)
349 {
350 msg=rqueue->dequeue();
351 char *s=(char*)msg->str();
352 int ind;
353 poly p=NULL;
354 get_poly(s,ind,&p,r);
355 //printf("got res[%d]\n",ind);
356 result->m[ind]=p;
357 msg.free();
358 cnt--;
359 }
360 // removes queues
361 queue.free();
362 rqueue.free();
363 vmem_deinit();
364 }
365 return result;
366 }
367
test_n(poly n)368 void test_n(poly n)
369 {
370 p_Write(n,currRing);
371 char *buf=(char*)omAlloc0(2048*1000);
372 int ll=size_poly(n,currRing);
373 printf("size: %d\n",ll);
374 char *s=send_poly(buf,12345,n,currRing);
375 printf("send len: %d\n",(int)(s-buf));
376 long *d=(long*)buf;
377 for(int i=0;i<=ll/SIZEOF_LONG;i++) printf("%ld ",d[i]);
378 printf("\n");
379 n=NULL;
380 s=get_poly(buf,ll,&n,currRing);
381 printf("read len: %d\n",(int)(s-buf));
382 Print(":index: %d\n",ll);
383 p_Write(n,currRing);
384 PrintLn();
385 omFree(buf);
386 }
387 #endif
388