source: git/kernel/weight0.c @ f3a8c2e

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