1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4
5 /*
6 * ABSTRACT:
7 */
8
9 #include "misc/auxiliary.h"
10 #include "omalloc/omalloc.h"
11
12 #include <math.h>
13 #include <string.h>
14
15 double wFunctionalMora(int *degw, int *lpol, int npol,
16 double *rel, double wx, double wwNsqr);
17 double wFunctionalBuch(int *degw, int *lpol, int npol,
18 double *rel, double wx, double wwNsqr);
19 void wAdd(int *A, int mons, int kn, int xx, int rvar);
20 void wNorm(int *degw, int *lpol, int npol, double *rel);
21 void wFirstSearch(int *A, int *x, int mons,
22 int *lpol, int npol, double *rel, double *fopt, double wNsqr, int rvar);
23 void wSecondSearch(int *A, int *x, int *lpol,
24 int npol, int mons, double *rel, double *fk, double wNsqr, int rvar);
25 void wGcd(int *x, int n);
26 /*0 implementation*/
27
28 VAR short * ecartWeights=NULL;
29
30 VAR double (*wFunctional)(int *degw, int *lpol, int npol,
31 double *rel, double wx, double wNsqr);
32
33
wFunctionalMora(int * degw,int * lpol,int npol,double * rel,double wx,double wNsqr)34 double wFunctionalMora(int *degw, int *lpol, int npol,
35 double *rel, double wx, double wNsqr)
36 {
37 int i, j, e1, ecu, ecl, ec;
38 int *ex;
39 double gfmax, gecart, ghom, pfmax;
40 double *r;
41
42 ex = degw;
43 r = rel;
44 gfmax = (double)0.0;
45 gecart = (double)0.4 + (double)npol;
46 ghom = (double)1.0;
47 for (i = 0; i < npol; i++)
48 {
49 ecl = ecu = e1 = *ex++;
50 for (j = lpol[i] - 1; j!=0; j--)
51 {
52 ec = *ex++;
53 if (ec > ecu)
54 ecu = ec;
55 else if (ec < ecl)
56 ecl = ec;
57 }
58 pfmax = (double)ecl / (double)ecu;
59 if (pfmax < ghom)
60 ghom = pfmax;
61 pfmax = (double)e1 / (double)ecu;
62 if (pfmax > (double)0.5)
63 gecart -= (pfmax * pfmax);
64 else
65 gecart -= (double)0.25;
66 ecu = 2 * ecu - ecl;
67 gfmax += (double)(ecu * ecu) * (*r++);
68 }
69 if (ghom > (double)0.8)
70 {
71 ghom *= (double)5.0;
72 gecart *= ((double)5.0 - ghom);
73 }
74 return (gfmax * gecart) / pow(wx, wNsqr);
75 }
76
77
wFunctionalBuch(int * degw,int * lpol,int npol,double * rel,double wx,double wNsqr)78 double wFunctionalBuch(int *degw, int *lpol, int npol,
79 double *rel, double wx, double wNsqr)
80 {
81 int i, j, ecl, ecu, ec;
82 int *ex;
83 double gfmax, ghom, pfmax;
84 double *r;
85
86 ex = degw;
87 r = rel;
88 gfmax = (double)0.0;
89 ghom = (double)1.0;
90 for (i = 0; i < npol; i++)
91 {
92 ecu = ecl = *ex++;
93 for (j = lpol[i] - 1; j!=0 ; j--)
94 {
95 ec = *ex++;
96 if (ec < ecl)
97 ecl = ec;
98 else if (ec > ecu)
99 ecu = ec;
100 }
101 pfmax = (double)ecl / (double)ecu;
102 if (pfmax < ghom)
103 ghom = pfmax;
104 gfmax += (double)(ecu * ecu) * (*r++);
105 }
106 if (ghom > (double)0.5)
107 gfmax *= ((double)1.0 - (ghom * ghom)) / (double)0.75;
108 return gfmax / pow(wx, wNsqr);
109 }
110
111
wSub(int * A,int mons,int kn,int xx,int rvar)112 static void wSub(int *A, int mons, int kn, int xx,int rvar)
113 {
114 int i, *B, *ex;
115
116 B = A + ((kn - 1) * mons);
117 ex = A + (rvar * mons);
118 i = mons;
119 if (xx == 1)
120 {
121 for (/* i=mons */; i!=0 ; i--)
122 *ex++ -= *B++;
123 }
124 else
125 {
126 for (/* i=mons */; i!=0 ; i--)
127 *ex++ -= (*B++) * xx;
128 }
129 }
130
131
wAdd(int * A,int mons,int kn,int xx,int rvar)132 void wAdd(int *A, int mons, int kn, int xx, int rvar)
133 {
134 int i, *B, *ex;
135
136 B = A + ((kn - 1) * mons);
137 ex = A + (rvar * mons);
138 i = mons;
139 if (xx == 1)
140 {
141 for (/* i=mons */; i!=0 ; i--)
142 *ex++ += *B++;
143 }
144 else
145 {
146 for (/* i=mons */; i!=0 ; i--)
147 *ex++ += (*B++) * xx;
148 }
149 }
150
151
wFirstSearch(int * A,int * x,int mons,int * lpol,int npol,double * rel,double * fopt,double wNsqr,int rvar)152 void wFirstSearch(int *A, int *x, int mons,
153 int *lpol, int npol, double *rel, double *fopt, double wNsqr, int rvar)
154 {
155 int a0, a, n, xn, t, xx, y1;
156 int *y, *degw, *xopt;
157 double fy, fmax, wx;
158 double *pr;
159 void *adr;
160
161 fy = *fopt;
162 n = rvar;
163 xn = n + 6 + (21 / n);
164 a0 = n * sizeof(double);
165 a = n * sizeof(int);
166 y = (int * )omAlloc((long)a);
167 adr = (void * )omAllocAligned((long)a0);
168 pr = adr;
169 *pr = (double)1.0;
170 *y = 0;
171 degw = A + (n * mons);
172 xopt = x + (n + 2);
173 t = 1;
174 loop
175 {
176 while (t < n)
177 {
178 xx = x[t] + 1;
179 wx = pr[t-1] * (double)xx;
180 y1 = y[t-1] + xx;
181 if ((y1 + n - t) <= xn)
182 {
183 pr[t] = wx;
184 y[t] = y1;
185 x[t] = xx;
186 if (xx > 1)
187 wAdd(A, mons, t, 1, rvar);
188 t++;
189 }
190 else
191 {
192 xx = x[t] - 1;
193 x[t] = 0;
194 if (xx!=0)
195 wSub(A, mons, t, xx, rvar);
196 t--;
197 if (t==0)
198 {
199 *fopt = fy;
200 omFreeSize((ADDRESS)y, (long)a);
201 omFreeSize((ADDRESS)adr, (long)a0);
202 return;
203 }
204 }
205 }
206 xx = xn - y[t-1];
207 wx = pr[t-1] * (double)xx;
208 x[t] = xx;
209 xx--;
210 if (xx!=0)
211 wAdd(A, mons, t, xx, rvar);
212 fmax = (*wFunctional)(degw, lpol, npol, rel, wx,wNsqr);
213 if (xx!=0)
214 wSub(A, mons, t, xx, rvar);
215 if (fmax < fy)
216 {
217 fy = fmax;
218 memcpy(xopt, x + 1, a);
219 }
220 t--;
221 } /* end loop */
222 }
223
224
wPrWeight(int * x,int n)225 static double wPrWeight(int *x, int n)
226 {
227 int i;
228 double y;
229
230 y = (double)x[n];
231 for (i = n - 1; i!=0 ; i--)
232 y *= (double)x[i];
233 return y;
234 }
235
236
wEstimate(int * A,int * x,int * lpol,int npol,int mons,double wx,double * rel,double * fopt,int * s0,int * s1,int * s2,double wNsqr,int rvar)237 static void wEstimate(int *A, int *x, int *lpol, int npol, int mons,
238 double wx, double *rel, double *fopt, int *s0, int *s1, int *s2, double wNsqr, int rvar)
239 {
240 int n, i1, i2, k0 = 0, k1 = 0, k2 = 0;
241 int *degw;
242 double fo1, fo2, fmax, wx1, wx2;
243
244 n = rvar;
245 degw = A + (n * mons);
246 fo2 = fo1 = (double)1.0e10;
247 for (i1 = n; i1!=0 ; i1--)
248 {
249 if (x[i1] > 1)
250 {
251 wSub(A, mons, i1, 1, rvar);
252 wx1 = wx - wx / (double)x[i1];
253 x[i1]--;
254 fmax = (*wFunctional)(degw, lpol, npol, rel, wx1,wNsqr);
255 if (fmax < fo1)
256 {
257 fo1 = fmax;
258 k0 = i1;
259 }
260 for (i2 = i1; i2!=0 ; i2--)
261 {
262 if (x[i2] > 1)
263 {
264 wSub(A, mons, i2, 1, rvar);
265 wx2 = wx1 - wx1 / (double)x[i2];
266 fmax = (*wFunctional)(degw, lpol, npol, rel, wx2, wNsqr);
267 if (fmax < fo2)
268 {
269 fo2 = fmax;
270 k1 = i1;
271 k2 = i2;
272 }
273 wAdd(A, mons, i2, 1, rvar);
274 }
275 }
276 wAdd(A, mons, i1, 1, rvar);
277 x[i1]++;
278 }
279 }
280 if (fo1 < fo2)
281 {
282 *fopt = fo1;
283 *s0 = k0;
284 }
285 else
286 {
287 *fopt = fo2;
288 *s0 = 0;
289 }
290 *s1 = k1;
291 *s2 = k2;
292 }
293
294
wSecondSearch(int * A,int * x,int * lpol,int npol,int mons,double * rel,double * fk,double wNsqr,int rvar)295 void wSecondSearch(int *A, int *x, int *lpol,
296 int npol, int mons, double *rel, double *fk, double wNsqr, int rvar)
297 {
298 int n, s0, s1, s2, *xopt;
299 double fx, fopt, wx;
300
301 n = rvar;
302 xopt = x + (n + 2);
303 fopt = *fk * (double)0.999999999999;
304 wx = wPrWeight(x, n);
305 loop
306 {
307 wEstimate(A, x, lpol, npol, mons, wx, rel, &fx, &s0, &s1, &s2, wNsqr, rvar);
308 if (fx > fopt)
309 {
310 if (s0!=0)
311 x[s0]--;
312 else if (s1!=0)
313 {
314 x[s1]--;
315 x[s2]--;
316 }
317 else
318 break;
319 }
320 else
321 {
322 fopt = fx;
323 if (s0!=0)
324 {
325 x[s0]--;
326 memcpy(xopt, x + 1, n * sizeof(int));
327 if (s1==0)
328 break;
329 }
330 else if (s1!=0)
331 {
332 x[s1]--;
333 x[s2]--;
334 memcpy(xopt, x + 1, n * sizeof(int));
335 }
336 else
337 break;
338 }
339 if (s0!=0)
340 wSub(A, mons, s0, 1, rvar);
341 else
342 {
343 wSub(A, mons, s1, 1, rvar);
344 wSub(A, mons, s2, 1, rvar);
345 }
346 wx = wPrWeight(x, n);
347 }
348 *fk = fopt;
349 }
350
351
wGcd(int * x,int n)352 void wGcd(int *x, int n)
353 {
354 int i, b, a, h;
355
356 i = n;
357 b = x[i];
358 loop
359 {
360 i--;
361 if (i==0)
362 break;
363 a = x[i];
364 if (a < b)
365 {
366 h = a;
367 a = b;
368 b = h;
369 }
370 do
371 {
372 h = a % b;
373 a = b;
374 b = h;
375 }
376 while (b!=0);
377 b = a;
378 if (b == 1)
379 return;
380 }
381 for (i = n; i!=0 ; i--)
382 x[i] /= b;
383 }
384
385
386 #if 0 /*currently unused*/
387 static void wSimple(int *x, int n)
388 {
389 int g, min, c, d, f, kopt, k, i;
390 int *xopt;
391 double sopt, s1, s2;
392
393 xopt = x + (n + 1);
394 kopt = k = g = 0;
395 min = 1000000;
396 for (i = n; i!=0 ; i--)
397 {
398 c = xopt[i];
399 if (c > 1)
400 {
401 if (c < min)
402 min = c;
403 if (c > k)
404 k = c;
405 }
406 else
407 g = 1;
408 }
409 k -= min;
410 if ((g==0) && (k < 4))
411 return;
412 if (k < min)
413 min = k+1;
414 sopt = (double)1.0e10;
415 for (k = min; k > 1; k--)
416 {
417 s2 = s1 = (double)0.0;
418 for(i = n; i!=0 ; i--)
419 {
420 c = xopt[i];
421 if (c > 1)
422 {
423 d = c / k;
424 d *= k;
425 f = d = c - d;
426 if (f!=0)
427 {
428 f = k - f;
429 if (f < d)
430 s2 += (double)f / (double)c;
431 else
432 s1 += (double)d / (double)c;
433 }
434 }
435 }
436 s1 += s2 + sqrt(s1 * s2);
437 s1 -= (double)0.01 * sqrt((double)k);
438 if (s1 < sopt)
439 {
440 sopt = s1;
441 kopt = k;
442 }
443 }
444 for(i = n; i!=0 ; i--)
445 {
446 x[i] = 1;
447 c = xopt[i];
448 if (c > 1)
449 {
450 d = c / kopt;
451 d *= kopt;
452 x[i] = d;
453 d = c - d;
454 if ((d!=0) && (kopt < 2 * d))
455 x[i] += kopt;
456 }
457 }
458 if (g==0)
459 wGcd(x, n);
460 }
461 #endif
462
wNorm(int * degw,int * lpol,int npol,double * rel)463 void wNorm(int *degw, int *lpol, int npol, double *rel)
464 {
465 int i, j, ecu, ec;
466 int *ex;
467 double *r;
468
469 ex = degw;
470 r = rel;
471 for (i = 0; i < npol; i++)
472 {
473 ecu = *ex++;
474 for (j = lpol[i] - 1; j!=0 ; j--)
475 {
476 ec = *ex++;
477 if (ec > ecu)
478 ecu = ec;
479 }
480 *r = (double)1.0 / (double)(ecu * ecu);
481 r++;
482 }
483 }
484