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