source: git/libpolys/polys/weight.cc @ a44bcf

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