source: git/libpolys/polys/weight0.c @ 6ce030f

spielwiese
Last change on this file since 6ce030f was 6ce030f, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
removal of the $Id$ svn tag from everywhere NOTE: the git SHA1 may be used instead (only on special places) NOTE: the libraries Singular/LIB/*.lib still contain the marker due to our current use of svn
  • 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 <math.h>
10#include <string.h>
11#include <misc/auxiliary.h>
12#include <omalloc/omalloc.h>
13
14double wFunctionalMora(int *degw, int *lpol, int npol,
15       double *rel, double wx, double wwNsqr);
16double wFunctionalBuch(int *degw, int *lpol, int npol,
17       double *rel, double wx, double wwNsqr);
18void wAdd(int *A, int mons, int kn, int xx, int rvar);
19void wNorm(int *degw, int *lpol, int npol, double *rel);
20void wFirstSearch(int *A, int *x, int mons,
21        int *lpol, int npol, double *rel, double *fopt, double wNsqr, int rvar);
22void wSecondSearch(int *A, int *x, int *lpol,
23        int npol, int mons, double *rel, double *fk, double wNsqr, int rvar);
24void wGcd(int *x, int n);
25/*0 implementation*/
26
27short * ecartWeights=NULL;
28
29double (*wFunctional)(int *degw, int *lpol, int npol,
30       double *rel, double wx, double wNsqr);
31
32
33double wFunctionalMora(int *degw, int *lpol, int npol,
34       double *rel, double wx, double wNsqr)
35{
36  int  i, j, e1, ecu, ecl, ec;
37  int  *ex;
38  double gfmax, gecart, ghom, pfmax;
39  double *r;
40
41  ex = degw;
42  r = rel;
43  gfmax = (double)0.0;
44  gecart = (double)0.4 + (double)npol;
45  ghom = (double)1.0;
46  for (i = 0; i < npol; i++)
47  {
48    ecl = ecu = e1 = *ex++;
49    for (j = lpol[i] - 1; j!=0; j--)
50    {
51      ec = *ex++;
52      if (ec > ecu)
53        ecu = ec;
54      else if (ec < ecl)
55        ecl = ec;
56    }
57    pfmax = (double)ecl / (double)ecu;
58    if (pfmax < ghom)
59      ghom = pfmax;
60    pfmax = (double)e1 / (double)ecu;
61    if (pfmax > (double)0.5)
62      gecart -= (pfmax * pfmax);
63    else
64      gecart -= (double)0.25;
65    ecu = 2 * ecu - ecl;
66    gfmax += (double)(ecu * ecu) * (*r++);
67  }
68  if (ghom > (double)0.8)
69  {
70    ghom *= (double)5.0;
71    gecart *= ((double)5.0 - ghom);
72  }
73  return (gfmax * gecart) / pow(wx, wNsqr);
74}
75
76
77double wFunctionalBuch(int *degw, int *lpol, int npol,
78       double *rel, double wx, double wNsqr)
79{
80  int  i, j, ecl, ecu, ec;
81  int  *ex;
82  double gfmax, ghom, pfmax;
83  double *r;
84
85  ex = degw;
86  r = rel;
87  gfmax = (double)0.0;
88  ghom = (double)1.0;
89  for (i = 0; i < npol; i++)
90  {
91    ecu = ecl = *ex++;
92    for (j = lpol[i] - 1; j!=0 ; j--)
93    {
94      ec = *ex++;
95      if (ec < ecl)
96        ecl = ec;
97      else if (ec > ecu)
98        ecu = ec;
99    }
100    pfmax = (double)ecl / (double)ecu;
101    if (pfmax < ghom)
102      ghom = pfmax;
103    gfmax += (double)(ecu * ecu) * (*r++);
104  }
105  if (ghom > (double)0.5)
106    gfmax *= ((double)1.0 - (ghom * ghom)) / (double)0.75;
107  return gfmax / pow(wx, wNsqr);
108}
109
110
111static void wSub(int *A, int mons, int kn, int xx,int rvar)
112{
113  int  i, *B, *ex;
114
115  B = A + ((kn - 1) * mons);
116  ex = A + (rvar * mons);
117  i = mons;
118  if (xx == 1)
119  {
120    for (/* i=mons */; i!=0 ; i--)
121      *ex++ -= *B++;
122  }
123  else
124  {
125    for (/* i=mons */; i!=0 ; i--)
126      *ex++ -= (*B++) * xx;
127  }
128}
129
130
131void wAdd(int *A, int mons, int kn, int xx, int rvar)
132{
133  int  i, *B, *ex;
134
135  B = A + ((kn - 1) * mons);
136  ex = A + (rvar * mons);
137  i = mons;
138  if (xx == 1)
139  {
140    for (/* i=mons */; i!=0 ; i--)
141      *ex++ += *B++;
142  }
143  else
144  {
145    for (/* i=mons */; i!=0 ; i--)
146      *ex++ += (*B++) * xx;
147  }
148}
149
150
151void wFirstSearch(int *A, int *x, int mons,
152  int *lpol, int npol, double *rel, double *fopt, double wNsqr, int rvar)
153{
154  int  a0, a, n, xn, t, xx, y1;
155  int  *y, *degw, *xopt;
156  double  fy, fmax, wx;
157  double *pr;
158  void *adr;
159
160  fy = *fopt;
161  n = rvar;
162  xn = n + 6 + (21 / n);
163  a0 = n * sizeof(double);
164  a = n * sizeof(int);
165  y = (int * )omAlloc((long)a);
166  adr = (void * )omAllocAligned((long)a0);
167  pr = adr;
168  *pr = (double)1.0;
169  *y = 0;
170  degw = A + (n * mons);
171  xopt = x + (n + 2);
172  t = 1;
173  loop
174  {
175    while (t < n)
176    {
177      xx = x[t] + 1;
178      wx = pr[t-1] * (double)xx;
179      y1 = y[t-1] + xx;
180      if ((y1 + n - t) <= xn)
181      {
182        pr[t] = wx;
183        y[t] = y1;
184        x[t] = xx;
185        if (xx > 1)
186          wAdd(A, mons, t, 1, rvar);
187        t++;
188      }
189      else
190      {
191        xx = x[t] - 1;
192        x[t] = 0;
193        if (xx!=0)
194          wSub(A, mons, t, xx, rvar);
195        t--;
196        if (t==0)
197        {
198          *fopt = fy;
199          omFreeSize((ADDRESS)y, (long)a);
200          omFreeSize((ADDRESS)adr, (long)a0);
201          return;
202        }
203      }
204    }
205    xx = xn - y[t-1];
206    wx = pr[t-1] * (double)xx;
207    x[t] = xx;
208    xx--;
209    if (xx!=0)
210      wAdd(A, mons, t, xx, rvar);
211    fmax = (*wFunctional)(degw, lpol, npol, rel, wx,wNsqr);
212    if (xx!=0)
213      wSub(A, mons, t, xx, rvar);
214    if (fmax < fy)
215    {
216      fy = fmax;
217      memcpy(xopt, x + 1, a);
218    }
219    t--;
220  } /* end loop */
221}
222
223
224static double wPrWeight(int *x, int n)
225{
226  int i;
227  double y;
228
229  y = (double)x[n];
230  for (i = n - 1; i!=0 ; i--)
231    y *= (double)x[i];
232  return y;
233}
234
235
236static void wEstimate(int *A, int *x, int *lpol, int npol, int mons,
237double wx, double *rel, double *fopt, int *s0, int *s1, int *s2, double wNsqr, int rvar)
238{
239  int  n, i1, i2, k0 = 0, k1 = 0, k2 = 0;
240  int  *degw;
241  double fo1, fo2, fmax, wx1, wx2;
242
243  n = rvar;
244  degw = A + (n * mons);
245  fo2 = fo1 = (double)1.0e10;
246  for (i1 = n; i1!=0 ; i1--)
247  {
248    if (x[i1] > 1)
249    {
250      wSub(A, mons, i1, 1, rvar);
251      wx1 = wx - wx / (double)x[i1];
252      x[i1]--;
253      fmax = (*wFunctional)(degw, lpol, npol, rel, wx1,wNsqr);
254      if (fmax < fo1)
255      {
256        fo1 = fmax;
257        k0 = i1;
258      }
259      for (i2 = i1; i2!=0 ; i2--)
260      {
261        if (x[i2] > 1)
262        {
263          wSub(A, mons, i2, 1, rvar);
264          wx2 = wx1 - wx1 / (double)x[i2];
265          fmax = (*wFunctional)(degw, lpol, npol, rel, wx2, wNsqr);
266          if (fmax < fo2)
267          {
268            fo2 = fmax;
269            k1 = i1;
270            k2 = i2;
271          }
272          wAdd(A, mons, i2, 1, rvar);
273        }
274      }
275      wAdd(A, mons, i1, 1, rvar);
276      x[i1]++;
277    }
278  }
279  if (fo1 < fo2)
280  {
281    *fopt = fo1;
282    *s0 = k0;
283  }
284  else
285  {
286    *fopt = fo2;
287    *s0 = 0;
288  }
289  *s1 = k1;
290  *s2 = k2;
291}
292
293
294void wSecondSearch(int *A, int *x, int *lpol,
295int npol, int mons, double *rel, double *fk, double wNsqr, int rvar)
296{
297  int  n, s0, s1, s2, *xopt;
298  double  one, fx, fopt, wx;
299
300  n = rvar;
301  xopt = x + (n + 2);
302  fopt = *fk * (double)0.999999999999;
303  wx = wPrWeight(x, n);
304  one = (double)1.0;
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
386static void wSimple(int *x, int n)
387{
388  int g, min, c, d, f, kopt, k, i;
389  int *xopt;
390  double sopt, s1, s2;
391
392  xopt = x + (n + 1);
393  kopt = k = g = 0;
394  min = 1000000;
395  for (i = n; i!=0 ; i--)
396  {
397    c = xopt[i];
398    if (c > 1)
399    {
400      if (c < min)
401        min = c;
402      if (c > k)
403        k = c;
404    }
405    else
406      g = 1;
407  }
408  k -= min;
409  if ((g==0) && (k < 4))
410    return;
411  if (k < min)
412    min = k+1;
413  sopt = (double)1.0e10;
414  for (k = min; k > 1; k--)
415  {
416    s2 = s1 = (double)0.0;
417    for(i = n; i!=0 ; i--)
418    {
419      c = xopt[i];
420      if (c > 1)
421      {
422        d = c / k;
423        d *= k;
424        f = d = c - d;
425        if (f!=0)
426        {
427          f = k - f;
428          if (f < d)
429            s2 += (double)f / (double)c;
430          else
431            s1 += (double)d / (double)c;
432        }
433      }
434    }
435    s1 += s2 + sqrt(s1 * s2);
436    s1 -= (double)0.01 * sqrt((double)k);
437    if (s1 < sopt)
438    {
439      sopt = s1;
440      kopt = k;
441    }
442  }
443  for(i = n; i!=0 ; i--)
444  {
445    x[i] = 1;
446    c = xopt[i];
447    if (c > 1)
448    {
449      d = c / kopt;
450      d *= kopt;
451      x[i] = d;
452      d = c - d;
453      if ((d!=0) && (kopt < 2 * d))
454        x[i] += kopt;
455    }
456  }
457  if (g==0)
458    wGcd(x, n);
459}
460
461
462void wNorm(int *degw, int *lpol, int npol, double *rel)
463{
464  int  i, j, ecu, ec;
465  int  *ex;
466  double *r;
467
468  ex = degw;
469  r = rel;
470  for (i = 0; i < npol; i++)
471  {
472    ecu = *ex++;
473    for (j = lpol[i] - 1; j!=0 ; j--)
474    {
475      ec = *ex++;
476      if (ec > ecu)
477        ecu = ec;
478    }
479    *r = (double)1.0 / (double)(ecu * ecu);
480    r++;
481  }
482}
Note: See TracBrowser for help on using the repository browser.