Changeset 4bc0478 in git


Ignore:
Timestamp:
Jul 2, 1998, 3:28:27 PM (26 years ago)
Author:
Wilfred Pohl <pohl@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
911089175d64029950da36384bcbb7951882ca63
Parents:
2302f759d08efc14b09786b9bdd4caea94a48333
Message:
CodeWarrior: #include "weight0.c" in weight.cc


git-svn-id: file:///usr/local/Singular/svn/trunk@2276 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/weight.cc

    r2302f75 r4bc0478  
    22*  Computer Algebra System SINGULAR     *
    33****************************************/
    4 /* $Id: weight.cc,v 1.9 1998-06-13 14:20:07 Singular Exp $ */
     4/* $Id: weight.cc,v 1.10 1998-07-02 13:28:27 pohl Exp $ */
    55
    66/*
    … …  
    4747extern double wNsqr;
    4848#else
    49 short * ecartWeights=NULL;
    50 
    51 double wNsqr;
    52 double (*wFunctional)(int *degw, int *lpol, int npol,
    53        double *rel, double wx);
    54 
    55 
    56 static double wFunctionalMora(int *degw, int *lpol, int npol,
    57        double *rel, double wx)
    58 {
    59   int  i, j, e1, ecu, ecl, ec;
    60   int  *ex;
    61   double gfmax, gecart, ghom, pfmax;
    62   double *r;
    63 
    64   ex = degw;
    65   r = rel;
    66   gfmax = (double)0.0;
    67   gecart = (double)0.4 + (double)npol;
    68   ghom = (double)1.0;
    69   for (i = 0; i < npol; i++)
    70   {
    71     ecl = ecu = e1 = *ex++;
    72     for (j = lpol[i] - 1; j!=0; j--)
    73     {
    74       ec = *ex++;
    75       if (ec > ecu)
    76         ecu = ec;
    77       else if (ec < ecl)
    78         ecl = ec;
    79     }
    80     pfmax = (double)ecl / (double)ecu;
    81     if (pfmax < ghom)
    82       ghom = pfmax;
    83     pfmax = (double)e1 / (double)ecu;
    84     if (pfmax > (double)0.5)
    85       gecart -= (pfmax * pfmax);
    86     else
    87       gecart -= (double)0.25;
    88     ecu = 2 * ecu - ecl;
    89     gfmax += (double)(ecu * ecu) * (*r++);
    90   }
    91   if (ghom > (double)0.8)
    92   {
    93     ghom *= (double)5.0;
    94     gecart *= ((double)5.0 - ghom);
    95   }
    96   return (gfmax * gecart) / pow(wx, wNsqr);
    97 }
    98 
    99 
    100 static double wFunctionalBuch(int *degw, int *lpol, int npol,
    101        double *rel, double wx)
    102 {
    103   int  i, j, ecl, ecu, ec;
    104   int  *ex;
    105   double gfmax, ghom, pfmax;
    106   double *r;
    107 
    108   ex = degw;
    109   r = rel;
    110   gfmax = (double)0.0;
    111   ghom = (double)1.0;
    112   for (i = 0; i < npol; i++)
    113   {
    114     ecu = ecl = *ex++;
    115     for (j = lpol[i] - 1; j!=0; j--)
    116     {
    117       ec = *ex++;
    118       if (ec < ecl)
    119         ecl = ec;
    120       else if (ec > ecu)
    121         ecu = ec;
    122     }
    123     pfmax = (double)ecl / (double)ecu;
    124     if (pfmax < ghom)
    125       ghom = pfmax;
    126     gfmax += (double)(ecu * ecu) * (*r++);
    127   }
    128   if (ghom > (double)0.5)
    129     gfmax *= ((double)1.0 - (ghom * ghom)) / (double)0.75;
    130   return gfmax / pow(wx, wNsqr);
    131 }
    132 
    133 
    134 static void wSub(int *A, int mons, int kn, int xx)
    135 {
    136   int  i, *B, *ex;
    137 
    138   B = A + ((kn - 1) * mons);
    139   ex = A + (pVariables * mons);
    140   if (xx == 1)
    141   {
    142     for (i = mons; i!=0; i--)
    143       *ex++ -= *B++;
    144   }
    145   else
    146   {
    147     for (i = mons; i!=0; i--)
    148       *ex++ -= (*B++) * xx;
    149   }
    150 }
    151 
    152 
    153 static void wAdd(int *A, int mons, int kn, int xx)
    154 {
    155   int  i, *B, *ex;
    156 
    157   B = A + ((kn - 1) * mons);
    158   ex = A + (pVariables * mons);
    159   if (xx == 1)
    160   {
    161     for (i = mons; i!=0; i--)
    162       *ex++ += *B++;
    163   }
    164   else
    165   {
    166     for (i = mons; i!=0; i--)
    167       *ex++ += (*B++) * xx;
    168   }
    169 }
    170 
    171 
    172 static void wFirstSearch(int *A, int *x, int mons,
    173   int *lpol, int npol, double *rel, double *fopt)
    174 {
    175   int  a0, a, n, xn, t, xx, y1;
    176   int  *y, *degw, *xopt;
    177   double  fy, fmax, wx;
    178   double *pr;
    179   void *adr;
    180 
    181   fy = *fopt;
    182   n = pVariables;
    183   xn = n + 6 + (21 / n);
    184   a0 = n * sizeof(double) + 8;
    185   a = n * sizeof(int);
    186   y = (int * )Alloc(a);
    187   adr = (void * )Alloc(a0);
    188   pr = wDouble(adr);
    189   *pr = (double)1.0;
    190   *y = 0;
    191   degw = A + (n * mons);
    192   xopt = x + (n + 2);
    193   t = 1;
    194   loop
    195   {
    196     while (t < n)
    197     {
    198       xx = x[t] + 1;
    199       wx = pr[t-1] * (double)xx;
    200       y1 = y[t-1] + xx;
    201       if ((y1 + n - t) <= xn)
    202       {
    203         pr[t] = wx;
    204         y[t] = y1;
    205         x[t] = xx;
    206         if (xx > 1)
    207           wAdd(A, mons, t, 1);
    208         t++;
    209       }
    210       else
    211       {
    212         xx = x[t] - 1;
    213         x[t] = 0;
    214         if (xx!=0)
    215           wSub(A, mons, t, xx);
    216         t--;
    217         if (t==0)
    218         {
    219           *fopt = fy;
    220           Free((ADDRESS)y, a);
    221           Free((ADDRESS)adr, a0);
    222           return;
    223         }
    224       }
    225     }
    226     xx = xn - y[t-1];
    227     wx = pr[t-1] * (double)xx;
    228     x[t] = xx;
    229     xx--;
    230     if (xx!=0)
    231       wAdd(A, mons, t, xx);
    232     fmax = (*wFunctional)(degw, lpol, npol, rel, wx);
    233     if (xx!=0)
    234       wSub(A, mons, t, xx);
    235     if (fmax < fy)
    236     {
    237       fy = fmax;
    238       memcpy(xopt, x + 1, a);
    239     }
    240     t--;
    241   } /* end loop */
    242 }
    243 
    244 
    245 static double wPrWeight(int *x, int n)
    246 {
    247   int i;
    248   double y;
    249 
    250   y = (double)x[n];
    251   for (i = n - 1; i!=0; i--)
    252     y *= (double)x[i];
    253   return y;
    254 }
    255 
    256 
    257 static void wEstimate(int *A, int *x, int *lpol, int npol, int mons,
    258 double wx, double *rel, double *fopt, int *s0, int *s1, int *s2)
    259 {
    260   int  n, i1, i2, k0 = 0, k1 = 0, k2 = 0;
    261   int  *degw;
    262   double fo1, fo2, fmax, wx1, wx2;
    263 
    264   n = pVariables;
    265   degw = A + (n * mons);
    266   fo2 = fo1 = (double)1.0e10;
    267   for (i1 = n; i1!=0; i1--)
    268   {
    269     if (x[i1] > 1)
    270     {
    271       wSub(A, mons, i1, 1);
    272       wx1 = wx - wx / (double)x[i1];
    273       x[i1]--;
    274       fmax = (*wFunctional)(degw, lpol, npol, rel, wx1);
    275       if (fmax < fo1)
    276       {
    277         fo1 = fmax;
    278         k0 = i1;
    279       }
    280       for (i2 = i1; i2!=0; i2--)
    281       {
    282         if (x[i2] > 1)
    283         {
    284           wSub(A, mons, i2, 1);
    285           wx2 = wx1 - wx1 / (double)x[i2];
    286           fmax = (*wFunctional)(degw, lpol, npol, rel, wx2);
    287           if (fmax < fo2)
    288           {
    289             fo2 = fmax;
    290             k1 = i1;
    291             k2 = i2;
    292           }
    293           wAdd(A, mons, i2, 1);
    294         }
    295       }
    296       wAdd(A, mons, i1, 1);
    297       x[i1]++;
    298     }
    299   }
    300   if (fo1 < fo2)
    301   {
    302     *fopt = fo1;
    303     *s0 = k0;
    304   }
    305   else
    306   {
    307     *fopt = fo2;
    308     *s0 = 0;
    309   }
    310   *s1 = k1;
    311   *s2 = k2;
    312 }
    313 
    314 
    315 static void wSecondSearch(int *A, int *x, int *lpol,
    316 int npol, int mons, double *rel, double *fk)
    317 {
    318   int  n, s0, s1, s2, *xopt;
    319   double  one, fx, fopt, wx;
    320 
    321   n = pVariables;
    322   xopt = x + (n + 2);
    323   fopt = *fk * (double)0.999999999999;
    324   wx = wPrWeight(x, n);
    325   one = (double)1.0;
    326   loop
    327   {
    328     wEstimate(A, x, lpol, npol, mons, wx, rel, &fx, &s0, &s1, &s2);
    329     if (fx > fopt)
    330     {
    331       if (s0!=0)
    332         x[s0]--;
    333       else if (s1!=0)
    334       {
    335         x[s1]--;
    336         x[s2]--;
    337       }
    338       else
    339         break;
    340     }
    341     else
    342     {
    343       fopt = fx;
    344       if (s0!=0)
    345       {
    346         x[s0]--;
    347         memcpy(xopt, x + 1, n * sizeof(int));
    348         if (s1==0)
    349           break;
    350       }
    351       else if (s1!=0)
    352       {
    353         x[s1]--;
    354         x[s2]--;
    355         memcpy(xopt, x + 1, n * sizeof(int));
    356       }
    357       else
    358         break;
    359     }
    360     if (s0!=0)
    361       wSub(A, mons, s0, 1);
    362     else
    363     {
    364       wSub(A, mons, s1, 1);
    365       wSub(A, mons, s2, 1);
    366     }
    367     wx = wPrWeight(x, n);
    368   }
    369   *fk = fopt;
    370 }
    371 
    372 
    373 static void wGcd(int *x, int n)
    374 {
    375   int i, b, a, h;
    376 
    377   i = n;
    378   b = x[i];
    379   loop
    380   {
    381     i--;
    382     if (i==0)
    383       break;
    384     a = x[i];
    385     if (a < b)
    386     {
    387       h = a;
    388       a = b;
    389       b = h;
    390     }
    391     do
    392     {
    393       h = a % b;
    394       a = b;
    395       b = h;
    396     }
    397     while (b!=0);
    398     b = a;
    399     if (b == 1)
    400       return;
    401   }
    402   for (i = n; i; i--)
    403     x[i] /= b;
    404 }
    405 
    406 
    407 static void wSimple(int *x, int n)
    408 {
    409   int g, min, c, d, f, kopt, k, i;
    410   int *xopt;
    411   double sopt, s1, s2;
    412 
    413   xopt = x + (n + 1);
    414   kopt = k = g = 0;
    415   min = 1000000;
    416   for (i = n; i!=0; i--)
    417   {
    418     c = xopt[i];
    419     if (c > 1)
    420     {
    421       if (c < min)
    422         min = c;
    423       if (c > k)
    424         k = c;
    425     }
    426     else
    427       g = 1;
    428   }
    429   k -= min;
    430   if ((g==0) && (k < 4))
    431     return;
    432   if (k < min)
    433     min = k+1;
    434   sopt = (double)1.0e10;
    435   for (k = min; k > 1; k--)
    436   {
    437     s2 = s1 = (double)0.0;
    438     for(i = n; i!=0; i--)
    439     {
    440       c = xopt[i];
    441       if (c > 1)
    442       {
    443         d = c / k;
    444         d *= k;
    445         f = d = c - d;
    446         if (f!=0)
    447         {
    448           f = k - f;
    449           if (f < d)
    450             s2 += (double)f / (double)c;
    451           else
    452             s1 += (double)d / (double)c;
    453         }
    454       }
    455     }
    456     s1 += s2 + sqrt(s1 * s2);
    457     s1 -= (double)0.01 * sqrt((double)k);
    458     if (s1 < sopt)
    459     {
    460       sopt = s1;
    461       kopt = k;
    462     }
    463   }
    464   for(i = n; i!=0; i--)
    465   {
    466     x[i] = 1;
    467     c = xopt[i];
    468     if (c > 1)
    469     {
    470       d = c / kopt;
    471       d *= kopt;
    472       x[i] = d;
    473       d = c - d;
    474       if ((d!=0) && (kopt < 2 * d))
    475         x[i] += kopt;
    476     }
    477   }
    478   if (g==0)
    479     wGcd(x, n);
    480 }
    481 
    482 
    483 static void wNorm(int *degw, int *lpol, int npol, double *rel)
    484 {
    485   int  i, j, ecu, ec;
    486   int  *ex;
    487   double *r;
    488 
    489   ex = degw;
    490   r = rel;
    491   for (i = 0; i < npol; i++)
    492   {
    493     ecu = *ex++;
    494     for (j = lpol[i] - 1; j!=0; j--)
    495     {
    496       ec = *ex++;
    497       if (ec > ecu)
    498         ecu = ec;
    499     }
    500     *r = (double)1.0 / (double)(ecu * ecu);
    501     r++;
    502   }
    503 }
    504 #endif /* __MWERKS */
    505 
     49#include "weight0.c"
     50#endif /* __MWERKS__ */
    50651
    50752static void wDimensions(polyset s, int sl, int *lpol, int *npol, int *mons)
Note: See TracChangeset for help on using the changeset viewer.