1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4
5 /*
6 * ABSTRACT:
7 */
8
9 #include "misc/auxiliary.h"
10
11 #include "misc/options.h"
12 #include "misc/intvec.h"
13
14 #include "polys/monomials/ring.h"
15 #include "polys/monomials/p_polys.h"
16
17 #include "polys/weight.h"
18
19 /*0 implementation*/
20 extern "C" THREAD_VAR double (*wFunctional)(int *degw, int *lpol, int npol,
21 double *rel, double wx, double wNsqr);
22 extern "C" double wFunctionalMora(int *degw, int *lpol, int npol,
23 double *rel, double wx, double wNsqr);
24 extern "C" double wFunctionalBuch(int *degw, int *lpol, int npol,
25 double *rel, double wx, double wNsqr);
26 extern "C" void wAdd(int *A, int mons, int kn, int xx, int rvar);
27 extern "C" void wNorm(int *degw, int *lpol, int npol, double *rel);
28 extern "C" void wFirstSearch(int *A, int *x, int mons,
29 int *lpol, int npol, double *rel, double *fopt, double wNsqr, int rvar);
30 extern "C" void wSecondSearch(int *A, int *x, int *lpol,
31 int npol, int mons, double *rel, double *fk, double wNsqr, int rvar);
32 extern "C" void wGcd(int *x, int n);
33
wDimensions(poly * s,int sl,int * lpol,int * npol,int * mons)34 static void wDimensions(poly* s, int sl, int *lpol, int *npol, int *mons)
35 {
36 int i, i1, j, k;
37 poly p, q;
38
39 i1 = j = 0;
40 for (i = 0; i <= sl; i++)
41 {
42 p = s[i];
43 if (p!=NULL)
44 {
45 k = 1;
46 q = pNext(p);
47 while (q!=NULL)
48 {
49 k++;
50 q = pNext(q);
51 }
52 if (k > 1)
53 {
54 lpol[i1++] = k;
55 j += k;
56 }
57 }
58 }
59 *npol = i1;
60 *mons = j;
61 }
62
63
wInit(poly * s,int sl,int mons,int * A,const ring R)64 static void wInit(poly* s, int sl, int mons, int *A, const ring R)
65 {
66 int n, a, i, j, *B, *C;
67 poly p, q;
68 int *pl;
69
70 B = A;
71 n = rVar(R);
72 a = (n + 1) * sizeof(int);
73 pl = (int *)omAlloc(a);
74 for (i = 0; i <= sl; i++)
75 {
76 p = s[i];
77 if (p!=NULL)
78 {
79 q = pNext(p);
80 if (q!=NULL)
81 {
82 C = B;
83 B++;
84 p_GetExpV(p, pl,R);
85 for (j = 0; j < n; j++)
86 {
87 *C = pl[j+1];
88 C += mons;
89 }
90 }
91 while (q!=NULL)
92 {
93 C = B;
94 B++;
95 p_GetExpV(q, pl,R);
96 for (j = 0; j < n; j++)
97 {
98 *C = pl[j+1];
99 C += mons;
100 }
101 pIter(q);
102 }
103 }
104 }
105 omFreeSize((ADDRESS)pl, a);
106 }
107
wCall(poly * s,int sl,int * x,double wNsqr,const ring R)108 void wCall(poly* s, int sl, int *x, double wNsqr, const ring R)
109 {
110 int n, q, npol, mons, i;
111 int *A, *xopt, *lpol, *degw;
112 double f1, fx, eps, *rel;
113 void *adr;
114
115 n = rVar(R);
116 lpol = (int * )omAlloc((sl + 1) * sizeof(int));
117 wDimensions(s, sl, lpol, &npol, &mons);
118 xopt = x + (n + 1);
119 for (i = n; i!=0; i--)
120 xopt[i] = 1;
121 if (mons==0)
122 {
123 omFreeSize((ADDRESS)lpol, (sl + 1) * sizeof(int));
124 return;
125 }
126 adr = (void * )omAllocAligned(npol * sizeof(double));
127 rel = (double*)adr;
128 q = (n + 1) * mons * sizeof(int);
129 A = (int * )omAlloc(q);
130 wInit(s, sl, mons, A, R);
131 degw = A + (n * mons);
132 memset(degw, 0, mons * sizeof(int));
133 for (i = n; i!=0; i--)
134 wAdd(A, mons, i, 1, rVar(R));
135 wNorm(degw, lpol, npol, rel);
136 f1 = (*wFunctional)(degw, lpol, npol, rel, (double)1.0, wNsqr);
137 if (TEST_OPT_PROT) Print("// %e\n",f1);
138 eps = f1;
139 fx = (double)2.0 * eps;
140 memset(x, 0, (n + 1) * sizeof(int));
141 wFirstSearch(A, x, mons, lpol, npol, rel, &fx, wNsqr, rVar(R));
142 if (TEST_OPT_PROT) Print("// %e\n",fx);
143 memcpy(x + 1, xopt + 1, n * sizeof(int));
144 memset(degw, 0, mons * sizeof(int));
145 for (i = n; i!=0; i--)
146 {
147 x[i] *= 16;
148 wAdd(A, mons, i, x[i], rVar(R));
149 }
150 wSecondSearch(A, x, lpol, npol, mons, rel, &fx, wNsqr, rVar(R));
151 if (TEST_OPT_PROT) Print("// %e\n",fx);
152 if (fx >= eps)
153 {
154 for (i = n; i!=0; i--)
155 xopt[i] = 1;
156 }
157 else
158 {
159 wGcd(xopt, n);
160 // if (BTEST1(22))
161 // {
162 // f1 = fx + (double)0.1 * (f1 - fx);
163 // wSimple(x, n);
164 // memset(degw, 0, mons * sizeof(int));
165 // for (i = n; i!=0; i--)
166 // wAdd(A, mons, i, x[i], rVar(R));
167 // eps = wPrWeight(x, n);
168 // fx = (*wFunctional)(degw, lpol, npol, rel, eps);
169 // if (fx < f1)
170 // {
171 // if (TEST_OPT_PROT) Print("// %e\n",fx);
172 // memcpy(xopt + 1, x + 1, n * sizeof(int));
173 // }
174 // }
175 }
176 omFreeSize((ADDRESS)A, q);
177 omFreeSize((ADDRESS)lpol, (sl + 1) * sizeof(int));
178 omFreeSize((ADDRESS)adr, npol * sizeof(double));
179 }
180
181
kEcartWeights(poly * s,int sl,short * eweight,const ring R)182 void kEcartWeights(poly* s, int sl, short *eweight, const ring R)
183 {
184 int n, i;
185 int *x;
186
187 *eweight = 0;
188 n = rVar(R);
189 if (rHasLocalOrMixedOrdering(R))
190 wFunctional = wFunctionalMora;
191 else
192 wFunctional = wFunctionalBuch;
193 x = (int * )omAlloc(2 * (n + 1) * sizeof(int));
194 wCall(s, sl, x, (double)2.0 / (double)n, R);
195 for (i = n; i!=0; i--)
196 eweight[i] = x[i + n + 1];
197 omFreeSize((ADDRESS)x, 2 * (n + 1) * sizeof(int));
198 }
199
iv2array(intvec * iv,const ring R)200 int * iv2array(intvec * iv, const ring R)
201 {
202 int *s=(int *)omAlloc0((rVar(R)+1)*sizeof(int));
203 int len=0;
204 if(iv!=NULL)
205 len=si_min(iv->length(),rVar(R)); // usually: rVar(R)==length()
206 int i;
207 //for(i=pVariables;i>len;i--) s[i]=1;
208 for(i=len;i>0;i--) s[i]=(*iv)[i-1];
209 return s;
210 }
211
212 /*2
213 *computes the degree of the leading term of the polynomial
214 *with respect to given ecartWeights
215 *used for Graebes method if BTEST1(31) is set
216 */
totaldegreeWecart(poly p,ring r)217 long totaldegreeWecart(poly p, ring r)
218 {
219 int i;
220 long j =0;
221
222 for (i=rVar(r); i>0; i--)
223 j += (int)(p_GetExp(p,i,r) * ecartWeights[i]);
224 return j;
225 }
226
227 /*2
228 *computes the degree of the leading term of the polynomial
229 *with respect to given weights
230 */
totaldegreeWecart_IV(poly p,ring r,const int * w)231 long totaldegreeWecart_IV(poly p, ring r, const int *w)
232 {
233 int i;
234 long j =0;
235
236 for (i=rVar(r); i>0; i--)
237 j += (long)((int)(p_GetExp(p,i,r) * w[i]));
238 return j;
239 }
240
241 /*2
242 *computes the maximal degree of all terms of the polynomial
243 *with respect to given ecartWeights and
244 *computes the length of the polynomial
245 *used for Graebes method if BTEST1(31) is set
246 */
maxdegreeWecart(poly p,int * l,ring r)247 long maxdegreeWecart(poly p,int *l, ring r)
248 {
249 short k=p_GetComp(p, r);
250 int ll=1;
251 long t,max;
252
253 max=totaldegreeWecart(p, r);
254 pIter(p);
255 while ((p!=NULL) && (p_GetComp(p, r)==k))
256 {
257 t=totaldegreeWecart(p, r);
258 if (t>max) max=t;
259 ll++;
260 pIter(p);
261 }
262 *l=ll;
263 return max;
264 }
265