source: git/libpolys/polys/weight0.c @ 75e126

spielwiese
Last change on this file since 75e126 was 52e6ef, checked in by Hans Schoenemann <hannes@…>, 8 years ago
chg: use misc/auxiliary.h in libpolys and kernel/mod2.h in Singular for config
  • Property mode set to 100644
File size: 8.8 KB
Line 
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
15double wFunctionalMora(int *degw, int *lpol, int npol,
16       double *rel, double wx, double wwNsqr);
17double wFunctionalBuch(int *degw, int *lpol, int npol,
18       double *rel, double wx, double wwNsqr);
19void wAdd(int *A, int mons, int kn, int xx, int rvar);
20void wNorm(int *degw, int *lpol, int npol, double *rel);
21void wFirstSearch(int *A, int *x, int mons,
22        int *lpol, int npol, double *rel, double *fopt, double wNsqr, int rvar);
23void wSecondSearch(int *A, int *x, int *lpol,
24        int npol, int mons, double *rel, double *fk, double wNsqr, int rvar);
25void wGcd(int *x, int n);
26/*0 implementation*/
27
28short * ecartWeights=NULL;
29
30double (*wFunctional)(int *degw, int *lpol, int npol,
31       double *rel, double wx, double wNsqr);
32
33
34double 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
78double 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
112static 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
132void 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
152void 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
225static 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
237static void wEstimate(int *A, int *x, int *lpol, int npol, int mons,
238double 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
295void wSecondSearch(int *A, int *x, int *lpol,
296int 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
352void 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*/
387static 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
463void 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}
Note: See TracBrowser for help on using the repository browser.