source: git/libpolys/polys/weight.cc @ 975db18

jengelh-datetimespielwiese
Last change on this file since 975db18 was 0033111, checked in by Hans Schoenemann <hannes@…>, 11 years ago
syntax fix
  • Property mode set to 100644
File size: 6.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5
6/*
7* ABSTRACT:
8*/
9
10#include <math.h>
11//#include <kernel/mod2.h>
12#include <omalloc/omalloc.h>
13#include <misc/options.h>
14#include <polys/monomials/p_polys.h>
15#include <misc/intvec.h>
16//#include <kernel/febase.h>
17//#include <kernel/ideals.h>
18#include <polys/monomials/ring.h>
19#include <polys/weight.h>
20
21/*0 implementation*/
22extern "C" double (*wFunctional)(int *degw, int *lpol, int npol,
23       double *rel, double wx, double wNsqr);
24extern "C" double wFunctionalMora(int *degw, int *lpol, int npol,
25       double *rel, double wx, double wNsqr);
26extern "C" double wFunctionalBuch(int *degw, int *lpol, int npol,
27       double *rel, double wx, double wNsqr);
28extern "C" void wAdd(int *A, int mons, int kn, int xx, int rvar);
29extern "C" void wNorm(int *degw, int *lpol, int npol, double *rel);
30extern "C" void wFirstSearch(int *A, int *x, int mons,
31        int *lpol, int npol, double *rel, double *fopt, double wNsqr, int rvar);
32extern "C" void wSecondSearch(int *A, int *x, int *lpol,
33        int npol, int mons, double *rel, double *fk, double wNsqr, int rvar);
34extern "C" void wGcd(int *x, int n);
35
36static void wDimensions(poly* s, int sl, int *lpol, int *npol, int *mons)
37{
38  int  i, i1, j, k;
39  poly p, q;
40
41  i1 = j = 0;
42  for (i = 0; i <= sl; i++)
43  {
44    p = s[i];
45    if (p!=NULL)
46    {
47      k = 1;
48      q = pNext(p);
49      while (q!=NULL)
50      {
51        k++;
52        q = pNext(q);
53      }
54      if (k > 1)
55      {
56        lpol[i1++] = k;
57        j += k;
58      }
59    }
60  }
61  *npol = i1;
62  *mons = j;
63}
64
65
66static void wInit(poly* s, int sl, int mons, int *A, const ring R)
67{
68  int  n, a, i, j, *B, *C;
69  poly p, q;
70  int *pl;
71
72  B = A;
73  n = rVar(R);
74  a = (n + 1) * sizeof(int);
75  pl = (int *)omAlloc(a);
76  for (i = 0; i <= sl; i++)
77  {
78    p = s[i];
79    if (p!=NULL)
80    {
81      q = pNext(p);
82      if (q!=NULL)
83      {
84        C = B;
85        B++;
86        p_GetExpV(p, pl,R);
87        for (j = 0; j < n; j++)
88        {
89          *C = pl[j+1];
90          C += mons;
91        }
92      }
93      while (q!=NULL)
94      {
95        C = B;
96        B++;
97        p_GetExpV(q, pl,R);
98        for (j = 0; j < n; j++)
99        {
100          *C = pl[j+1];
101          C += mons;
102        }
103        pIter(q);
104      }
105    }
106  }
107  omFreeSize((ADDRESS)pl, a);
108}
109
110void wCall(poly* s, int sl, int *x, double wNsqr, const ring R)
111{
112  int  n, q, npol, mons, i;
113  int  *A, *xopt, *lpol, *degw;
114  double  f1, fx, eps, *rel;
115  void *adr;
116
117  n = rVar(R);
118  lpol = (int * )omAlloc((sl + 1) * sizeof(int));
119  wDimensions(s, sl, lpol, &npol, &mons);
120  xopt = x + (n + 1);
121  for (i = n; i!=0; i--)
122    xopt[i] = 1;
123  if (mons==0)
124  {
125    omFreeSize((ADDRESS)lpol, (sl + 1) * sizeof(int));
126    return;
127  }
128  adr = (void * )omAllocAligned(npol * sizeof(double));
129  rel = (double*)adr;
130  q = (n + 1) * mons * sizeof(int);
131  A = (int * )omAlloc(q);
132  wInit(s, sl, mons, A, R);
133  degw = A + (n * mons);
134  memset(degw, 0, mons * sizeof(int));
135  for (i = n; i!=0; i--)
136    wAdd(A, mons, i, 1, rVar(R));
137  wNorm(degw, lpol, npol, rel);
138  f1 = (*wFunctional)(degw, lpol, npol, rel, (double)1.0, wNsqr);
139  if (TEST_OPT_PROT) Print("// %e\n",f1);
140  eps = f1;
141  fx = (double)2.0 * eps;
142  memset(x, 0, (n + 1) * sizeof(int));
143  wFirstSearch(A, x, mons, lpol, npol, rel, &fx, wNsqr, rVar(R));
144  if (TEST_OPT_PROT) Print("// %e\n",fx);
145  memcpy(x + 1, xopt + 1, n * sizeof(int));
146  memset(degw, 0, mons * sizeof(int));
147  for (i = n; i!=0; i--)
148  {
149    x[i] *= 16;
150    wAdd(A, mons, i, x[i], rVar(R));
151  }
152  wSecondSearch(A, x, lpol, npol, mons, rel, &fx, wNsqr, rVar(R));
153  if (TEST_OPT_PROT) Print("// %e\n",fx);
154  if (fx >= eps)
155  {
156    for (i = n; i!=0; i--)
157      xopt[i] = 1;
158  }
159  else
160  {
161    wGcd(xopt, n);
162//    if (BTEST1(22))
163//    {
164//      f1 = fx + (double)0.1 * (f1 - fx);
165//      wSimple(x, n);
166//      memset(degw, 0, mons * sizeof(int));
167//      for (i = n; i!=0; i--)
168//        wAdd(A, mons, i, x[i], rVar(R));
169//      eps = wPrWeight(x, n);
170//      fx = (*wFunctional)(degw, lpol, npol, rel, eps);
171//      if (fx < f1)
172//      {
173//        if (TEST_OPT_PROT) Print("// %e\n",fx);
174//        memcpy(xopt + 1, x + 1, n * sizeof(int));
175//      }
176//    }
177  }
178  omFreeSize((ADDRESS)A, q);
179  omFreeSize((ADDRESS)lpol, (sl + 1) * sizeof(int));
180  omFreeSize((ADDRESS)adr, npol * sizeof(double));
181}
182
183
184void kEcartWeights(poly* s, int sl, short *eweight, const ring R)
185{
186  int  n, i;
187  int  *x;
188
189  *eweight = 0;
190  n = rVar(R);
191  if (rHasLocalOrMixedOrdering(R))
192    wFunctional = wFunctionalMora;
193  else
194    wFunctional = wFunctionalBuch;
195  x = (int * )omAlloc(2 * (n + 1) * sizeof(int));
196  wCall(s, sl, x, (double)2.0 / (double)n, R);
197  for (i = n; i!=0; i--)
198    eweight[i] = x[i + n + 1];
199  omFreeSize((ADDRESS)x, 2 * (n + 1) * sizeof(int));
200}
201
202short * iv2array(intvec * iv, const ring R)
203{
204  short *s=(short *)omAlloc0((rVar(R)+1)*sizeof(short));
205  int len=0;
206  if(iv!=NULL)
207    len=si_min(iv->length(),rVar(R)); // usually: rVar(R)==length()
208  int i;
209  //for(i=pVariables;i>len;i--) s[i]=1;
210  for(i=len;i>0;i--)               s[i]=(*iv)[i-1];
211  return s;
212}
213
214/*2
215*computes the degree of the leading term of the polynomial
216*with respect to given ecartWeights
217*used for Graebes method if BTEST1(31) is set
218*/
219long totaldegreeWecart(poly p, ring r)
220{
221  int i;
222  long j =0;
223
224  for (i=rVar(r); i>0; i--)
225    j += (int)(p_GetExp(p,i,r) * ecartWeights[i]);
226  return  j;
227}
228
229/*2
230*computes the degree of the leading term of the polynomial
231*with respect to given weights
232*/
233long totaldegreeWecart_IV(poly p, ring r, const short *w)
234{
235  int i;
236  long j =0;
237
238  for (i=rVar(r); i>0; i--)
239    j += (long)((int)(p_GetExp(p,i,r) * w[i]));
240  return  j;
241}
242
243/*2
244*computes the maximal degree of all terms of the polynomial
245*with respect to given ecartWeights and
246*computes the length of the polynomial
247*used for Graebes method if BTEST1(31) is set
248*/
249long maxdegreeWecart(poly p,int *l, ring r)
250{
251  short k=p_GetComp(p, r);
252  int ll=1;
253  long t,max;
254
255  max=totaldegreeWecart(p, r);
256  pIter(p);
257  while ((p!=NULL) && (p_GetComp(p, r)==k))
258  {
259    t=totaldegreeWecart(p, r);
260    if (t>max) max=t;
261    ll++;
262    pIter(p);
263  }
264  *l=ll;
265  return max;
266}
Note: See TracBrowser for help on using the repository browser.