Changeset fe9a47 in git


Ignore:
Timestamp:
Jul 6, 2017, 3:26:55 PM (7 years ago)
Author:
bendooru <bendooru@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
b7e650aa28103f9685a9dd9313ce231b375f1ac5
Parents:
f0fba1c6553d60f8414803dce6e04c16b8750b0e
git-author:
bendooru <bendooru@users.noreply.github.com>2017-07-06 15:26:55+02:00
git-committer:
bendooru <bendooru@users.noreply.github.com>2017-07-07 23:07:00+02:00
Message:
Interval/Box rig consistency

This improves the consistency of interval->R and box->R when
constructing boxes via the procedure `setRing(r)`. Fixes a problem where
calling `rootIsolation` over a real `basering` would fail (coefficient
fields between R[x1..xn] and R[x] differ)

Replaces `n*` procedures with `n_*` counterparts where using `currRing`
may make no sense.
Location:
Singular/dyn_modules/interval
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • Singular/dyn_modules/interval/interval.cc

    rf0fba1c rfe9a47  
    1212/* interval */
    1313
    14 interval::interval()
    15 {
    16     lower = nInit(0);
    17     upper = nInit(0);
    18     R = currRing;
     14interval::interval(ring r = currRing)
     15{
     16    lower = n_Init(0, r->cf);
     17    upper = n_Init(0, r->cf);
     18    R = r;
    1919    R->ref++;
    2020}
    2121
    22 interval::interval(number a)
    23 {
     22interval::interval(number a, ring r = currRing)
     23{
     24    // dangerous: check if a is in coefs r->cf
    2425    lower = a;
    25     upper = nCopy(a);
    26     R = currRing;
     26    upper = n_Copy(a, r->cf);
     27    R = r;
    2728    R->ref++;
    2829}
    2930
    30 interval::interval(number a, number b)
     31interval::interval(number a, number b, ring r = currRing)
    3132{
    3233    lower = a;
    3334    upper = b;
    34     R = currRing;
     35    R = r;
    3536    R->ref++;
    3637}
    … …  
    3839interval::interval(interval *I)
    3940{
    40     lower = nCopy(I->lower);
    41     upper = nCopy(I->upper);
     41    lower = n_Copy(I->lower, I->R->cf);
     42    upper = n_Copy(I->upper, I->R->cf);
    4243    R = I->R;
    4344    R->ref++;
    … …  
    4647interval::~interval()
    4748{
    48     nDelete(&lower);
    49     nDelete(&upper);
     49    n_Delete(&lower, R->cf);
     50    n_Delete(&upper, R->cf);
    5051    R->ref--;
     52}
     53
     54interval& interval::setRing(ring r)
     55{
     56    if (R != r)
     57    {
     58        if (R->cf != r->cf)
     59        {
     60            nMapFunc fun = n_SetMap(R->cf, r->cf);
     61            number lo = fun(lower, R->cf, r->cf),
     62                up = fun(upper, R->cf, r->cf);
     63            n_Delete(&lower, R->cf);
     64            n_Delete(&upper, R->cf);
     65            lower = lo;
     66            upper = up;
     67        }
     68        R->ref--;
     69        r->ref++;
     70        R = r;
     71    }
     72
     73    return *this;
    5174}
    5275
    … …  
    133156        interval* i = (interval*) d;
    134157
    135         // use nWrite since nothing better (?) exists
     158        // use n_Write since nothing better (?) exists
    136159        StringSetS("[");
    137         nWrite(i->lower);
     160        n_Write(i->lower, i->R->cf);
    138161        StringAppendS(", ");
    139         nWrite(i->upper);
     162        n_Write(i->upper, i->R->cf);
    140163        StringAppendS("]");
    141164
    … …  
    170193     */
    171194
    172     number n1, n2;
    173 
    174     if (args->Typ() == INT_CMD)
    175     {
    176         n1 = nInit((int)(long) args->Data());
    177     }
    178     else if (args->Typ() == NUMBER_CMD)
    179     {
    180         n1 = (number) args->CopyD();
    181     }
    182     else if (args->Typ() == intervalID)
    183     {
    184         interval *I = (interval*) args->Data();
    185         n1 = nCopy(I->lower);
    186         n2 = nCopy(I->upper);
     195    if (args->Typ() == intervalID)
     196    {
     197        RES = new interval((interval*) args->CopyD());
    187198    }
    188199    else
    189200    {
    190         WerrorS("Input not supported: first argument not int or number");
    191         return TRUE;
    192     }
    193 
    194     // check if second argument exists
    195     if (args->Typ() == intervalID)
    196     {
    197         RES = new interval(n1, n2);
    198     }
    199     else if (args->next == NULL)
    200     {
    201         RES = new interval(n1);
    202     }
    203     else
    204     {
    205         if (args->next->Typ() == INT_CMD)
    206         {
    207             n2 = nInit((int)(long) args->next->Data());
    208         }
    209         else if (args->next->Typ() == NUMBER_CMD)
    210         {
    211             n2 = (number) args->next->CopyD();
     201        number n1, n2;
     202
     203        if (args->Typ() == INT_CMD)
     204        {
     205            n1 = nInit((int)(long) args->Data());
     206        }
     207        else if (args->Typ() == NUMBER_CMD)
     208        {
     209            n1 = (number) args->CopyD();
    212210        }
    213211        else
    214212        {
    215             WerrorS("Input not supported: second argument not int or number");
     213            WerrorS("Input not supported: first argument not int or number");
    216214            return TRUE;
    217215        }
    218216
    219         RES = new interval(n1, n2);
     217        // check if second argument exists
     218        if (args->Typ() == intervalID)
     219        {
     220            RES = new interval(n1, n2);
     221        }
     222        else if (args->next == NULL)
     223        {
     224            RES = new interval(n1);
     225        }
     226        else
     227        {
     228            if (args->next->Typ() == INT_CMD)
     229            {
     230                n2 = nInit((int)(long) args->next->Data());
     231            }
     232            else if (args->next->Typ() == NUMBER_CMD)
     233            {
     234                n2 = (number) args->next->CopyD();
     235            }
     236            else
     237            {
     238                WerrorS("Input not supported: second argument not int or number");
     239                return TRUE;
     240            }
     241
     242            RES = new interval(n1, n2);
     243        }
    220244    }
    221245
    … …  
    246270        interval *I = (interval*) arg->Data();
    247271        result->rtyp = NUMBER_CMD;
    248         result->data = (void*) nSub(I->upper, I->lower);
     272        result->data = (void*) n_Sub(I->upper, I->lower, I->R->cf);
    249273        arg->CleanUp();
    250274        return FALSE;
    … …  
    259283static interval* intervalScalarMultiply(number a, interval *I)
    260284{
     285    // must assume a is in ring I->R
    261286    number lo, up;
    262287    if (nGreaterZero(a))
    263288    {
    264         lo = nMult(a, I->lower);
    265         up = nMult(a, I->upper);
     289        lo = n_Mult(a, I->lower, I->R->cf);
     290        up = n_Mult(a, I->upper, I->R->cf);
    266291    }
    267292    else
    268293    {
    269         lo = nMult(a, I->upper);
    270         up = nMult(a, I->lower);
    271     }
    272 
    273     nNormalize(lo);
    274     nNormalize(up);
    275 
    276     return new interval(lo, up);
     294        lo = n_Mult(a, I->upper, I->R->cf);
     295        up = n_Mult(a, I->lower, I->R->cf);
     296    }
     297
     298    n_Normalize(lo, I->R->cf);
     299    n_Normalize(up, I->R->cf);
     300
     301    return new interval(lo, up, I->R);
    277302}
    278303
    … …  
    281306    number lo, up;
    282307    number nums[4];
    283     nums[0] = nMult(I->lower, J->lower);
    284     nums[1] = nMult(I->lower, J->upper);
    285     nums[2] = nMult(I->upper, J->lower);
    286     nums[3] = nMult(I->upper, J->upper);
     308    nums[0] = n_Mult(I->lower, J->lower, I->R->cf);
     309    nums[1] = n_Mult(I->lower, J->upper, I->R->cf);
     310    nums[2] = n_Mult(I->upper, J->lower, I->R->cf);
     311    nums[3] = n_Mult(I->upper, J->upper, I->R->cf);
    287312
    288313    int i, imax = 0, imin = 0;
    289314    for (i = 1; i < 4; i++)
    290315    {
    291         if (nGreater(nums[i], nums[imax]))
     316        if (n_Greater(nums[i], nums[imax], I->R->cf))
    292317        {
    293318            imax = i;
    294319        }
    295         if (nGreater(nums[imin], nums[i]))
     320        if (n_Greater(nums[imin], nums[i], I->R->cf))
    296321        {
    297322            imin = i;
    … …  
    299324    }
    300325
    301     lo = nCopy(nums[imin]);
    302     up = nCopy(nums[imax]);
     326    lo = n_Copy(nums[imin], I->R->cf);
     327    up = n_Copy(nums[imax], I->R->cf);
    303328
    304329    // delete products
    305330    for (i = 0; i < 4; i++)
    306331    {
    307         nDelete(&nums[i]);
    308     }
    309 
    310     nNormalize(lo);
    311     nNormalize(up);
     332        n_Delete(&nums[i], I->R->cf);
     333    }
     334
     335    n_Normalize(lo, I->R->cf);
     336    n_Normalize(up, I->R->cf);
     337
     338    return new interval(lo, up, I->R);
     339}
     340
     341static interval* intervalAdd(interval *I, interval *J)
     342{
     343    number lo = n_Add(I->lower, J->lower, I->R->cf),
     344           up = n_Add(I->upper, J->upper, I->R->cf);
     345
     346    n_Normalize(lo, I->R->cf);
     347    n_Normalize(up, I->R->cf);
    312348
    313349    return new interval(lo, up);
    314350}
    315351
    316 static interval* intervalAdd(interval *I, interval *J)
    317 {
    318     number lo = nAdd(I->lower, J->lower),
    319            up = nAdd(I->upper, J->upper);
    320 
    321     nNormalize(lo);
    322     nNormalize(up);
    323 
    324     return new interval(lo, up);
    325 }
    326 
    327352static interval* intervalSubtract(interval *I, interval *J)
    328353{
    329     number lo = nSub(I->lower, J->upper),
    330            up = nSub(I->upper, J->lower);
    331 
    332     nNormalize(lo);
    333     nNormalize(up);
    334 
    335     return new interval(lo, up);
     354    number lo = n_Sub(I->lower, J->upper, I->R->cf),
     355           up = n_Sub(I->upper, J->lower, I->R->cf);
     356
     357    n_Normalize(lo, I->R->cf);
     358    n_Normalize(up, I->R->cf);
     359
     360    return new interval(lo, up, I->R);
    336361}
    337362
    338363static bool intervalEqual(interval *I, interval *J)
    339364{
    340     return nEqual(I->lower, J->lower) && nEqual(I->upper, J->upper);
     365    assume(I->R == J->R);
     366    return n_Equal(I->lower, J->lower, I->R->cf)
     367        && n_Equal(I->upper, J->upper, I->R->cf);
    341368}
    342369
    … …  
    344371static bool intervalContainsZero(interval *I)
    345372{
    346     number n = nMult(I->lower, I->upper);
    347     bool result = !nGreaterZero(n);
     373    number n = n_Mult(I->lower, I->upper, I->R->cf);
     374    bool result = !n_GreaterZero(n, I->R->cf);
    348375    // delete helper number
    349     nDelete(&n);
     376    n_Delete(&n, I->R->cf);
    350377
    351378    return result;
    … …  
    356383    if (p == 0)
    357384    {
    358         return new interval(nInit(1));
     385        return new interval(n_Init(1,I->R->cf), I->R);
    359386    }
    360387
    … …  
    362389    number lo, up;
    363390
    364     nPower(I->lower, p, &lo);
    365     nPower(I->upper, p, &up);
     391    n_Power(I->lower, p, &lo, I->R->cf);
     392    n_Power(I->upper, p, &up, I->R->cf);
    366393
    367394    // should work now
    368395    if (p % 2 == 1)
    369396    {
    370         return new interval(lo, up);
     397        return new interval(lo, up, I->R);
    371398    }
    372399    else
    … …  
    374401        // perform pointer swap if necessary
    375402        number tmp;
    376         if (nGreater(lo, up))
     403        if (n_Greater(lo, up, I->R->cf))
    377404        {
    378405            tmp = up;
    … …  
    383410        if (intervalContainsZero(I))
    384411        {
    385             nDelete(&lo);
    386             lo = nInit(0);
    387         }
    388         return new interval(lo, up);
     412            n_Delete(&lo, I->R->cf);
     413            lo = n_Init(0, I->R->cf);
     414        }
     415        return new interval(lo, up, I->R);
    389416    }
    390417}
    … …  
    428455            I1 = (interval*) i1->Data();
    429456            I2 = (interval*) i2->Data();
     457            if (I1->R != I2->R)
     458            {
     459                WerrorS("adding intervals defined in different rings not supported");
     460                return TRUE;
     461            }
    430462
    431463            RES = intervalAdd(I1, I2);
    … …  
    442474            I1 = (interval*) i1->Data();
    443475            I2 = (interval*) i2->Data();
     476            if (I1->R != I2->R)
     477            {
     478                WerrorS("subtracting intervals defined in different rings not supported");
     479                return TRUE;
     480            }
    444481
    445482            RES = intervalSubtract(I1, I2);
    … …  
    454491                I1 = (interval*) i1->Data();
    455492                I2 = (interval*) i2->Data();
     493                if (I1->R != I2->R)
     494                {
     495                    WerrorS("multiplying intervals defined in different rings not supported");
     496                    return TRUE;
     497                }
    456498
    457499                RES = intervalMultiply(I1, I2);
    … …  
    503545
    504546                // make sure I2 is invertible
    505                 if(intervalContainsZero(I2))
     547                if (intervalContainsZero(I2))
    506548                {
    507549                    WerrorS("second interval contains zero");
    … …  
    510552
    511553                number invlo, invup;
    512                 invlo = nInvers(I2->lower);
    513                 invup = nInvers(I2->upper);
     554                invlo = n_Invers(I2->lower, I2->R->cf);
     555                invup = n_Invers(I2->upper, I2->R->cf);
    514556
    515557                // inverse interval
    516                 interval *I2inv = new interval(invup, invlo);
     558                interval *I2inv = new interval(invup, invlo, I2->R);
    517559
    518560                if (i1->Typ() == intervalID)
    519561                {
    520562                    interval *I1 = (interval*) i1->Data();
     563                    if (I1->R != I2->R)
     564                    {
     565                        WerrorS("dividing intervals from different rings not supported");
     566                        delete I2inv;
     567                        return TRUE;
     568                    }
    521569                    RES = intervalMultiply(I1, I2inv);
    522570                }
    … …  
    815863            // make sure rings of boxes and their intervals are consistent
    816864            // this is important for serialization
    817             if (RES->R != RES->intervals[i]->R)
    818             {
    819                 if (RES->R->cf != RES->intervals[i]->R->cf)
    820                 {
    821                     WerrorS("Passing interval to ring with different coefficient field");
    822                     delete RES;
    823                     args->CleanUp();
    824                     return TRUE;
    825                 }
    826 
    827                 RES->intervals[i]->R->ref--;
    828                 RES->R->ref++;
    829                 RES->intervals[i]->R = RES->R;
    830             }
     865            RES->intervals[i]->setRing(RES->R);
    831866        }
    832867    }
    … …  
    909944            {
    910945                WerrorS("second argument not box");
    911             }
    912             if (result->Data() != NULL)
    913             {
    914                 delete (box*) result->Data();
     946                return TRUE;
    915947            }
    916948
    … …  
    918950            // maybe try to skip this initialisation
    919951            // copying def of box() results in segfault?
     952            if (B1->R != B2->R)
     953            {
     954                WerrorS("subtracting boxes from different rings not supported");
     955                return TRUE;
     956            }
    920957            RES = new box();
    921958            int i;
    … …  
    11171154
    11181155    RES->setInterval(i-1, new interval(I));
    1119 
    1120     // same as above, ensure ring consistency
    1121     if (RES->R != RES->intervals[i-1]->R)
    1122     {
    1123         if (RES->R->cf != RES->intervals[i-1]->R->cf)
    1124         {
    1125             WerrorS("Passing interval to ring with different coefficient field");
    1126             delete RES;
    1127             args->CleanUp();
    1128             return TRUE;
    1129         }
    1130 
    1131         RES->intervals[i]->R->ref--;
    1132         RES->R->ref++;
    1133         RES->intervals[i]->R = RES->R;
    1134     }
     1156    // ensure consistency
     1157    RES->intervals[i-1]->setRing(RES->R);
    11351158
    11361159    result->rtyp = boxID;

  • Singular/dyn_modules/interval/interval.h

    rf0fba1c rfe9a47  
    1010    ring R;
    1111
    12     interval();
    13     interval(number);
    14     interval(number, number);
     12    interval(ring);
     13    interval(number, ring);
     14    interval(number, number, ring);
    1515    interval(interval*);
    1616    ~interval();
     17
     18    interval& setRing(ring);
    1719};
    1820
Note: See TracChangeset for help on using the changeset viewer.