Changeset 270e65 in git


Ignore:
Timestamp:
Sep 20, 2000, 3:25:42 PM (24 years ago)
Author:
Olaf Bachmann <obachman@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
822fc78417de82c543b235478d77e1cd5e3a4bf3
Parents:
44e09ae64325b228e21fa4d05f93f5a78c922bf6
Message:
* first try for FiledQ optimizations


git-svn-id: file:///usr/local/Singular/svn/trunk@4625 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
Singular
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • Singular/iparith.cc

    r44e09a r270e65  
    22*  Computer Algebra System SINGULAR     *
    33****************************************/
    4 /* $Id: iparith.cc,v 1.226 2000-09-19 12:43:28 Singular Exp $ */
     4/* $Id: iparith.cc,v 1.227 2000-09-20 13:25:39 obachman Exp $ */
    55
    66/*
    … …  
    33603360        case (int)jjidMinBase:    dArith1[i].p=(proc1)idMinBase; break;
    33613361        case (int)jjsyMinBase:    dArith1[i].p=(proc1)syMinBase; break;
    3362         case (int)jjpMaxComp:     dArith1[i].p=(proc1)pMaxComp; break;
     3362        case (int)jjpMaxComp:     dArith1[i].p=(proc1)pMaxCompProc; break;
    33633363        case (int)jjmpTrace:      dArith1[i].p=(proc1)mpTrace; break;
    33643364        case (int)jjmpTransp:     dArith1[i].p=(proc1)mpTransp; break;

  • Singular/longrat.cc

    r44e09a r270e65  
    22*  Computer Algebra System SINGULAR     *
    33****************************************/
    4 /* $Id: longrat.cc,v 1.31 2000-09-18 09:19:13 obachman Exp $ */
     4/* $Id: longrat.cc,v 1.32 2000-09-20 13:25:40 obachman Exp $ */
    55/*
    66* ABSTRACT: computation with long rational numbers (Hubert Grassmann)
    77*/
     8
     9#ifndef LONGRAT_CC
     10#define LONGRAT_CC
    811
    912#include <string.h>
    … …  
    2023#include "longrat.h"
    2124
    22 omBin rnumber_bin = omGetSpecBin(sizeof(rnumber));
    2325
    2426#ifndef BYTES_PER_MP_LIMB
    … …  
    2931#endif
    3032#endif
    31 
    32 static int nlPrimeM;
    33 static number nlMapP(number from)
    34 {
    35   number to;
    36   int save=npPrimeM;
    37   npPrimeM=nlPrimeM;
    38   to = nlInit(npInt(from));
    39   npPrimeM=save;
    40   return to;
    41 }
    42 
    43 //static number nlMapLongR(number from);
    44 static number nlMapR(number from);
    45 
    46 BOOLEAN nlSetMap(ring r)
    47 {
    48   if (rField_is_Q(r))
    49   {
    50     nMap = nlCopy;   /*Q -> Q*/
    51     return TRUE;
    52   }
    53   if (rField_is_Zp(r))
    54   {
    55     nlPrimeM=rChar(r);
    56     nMap = nlMapP; /* Z/p -> Q */
    57     return TRUE;
    58   }
    59   if (rField_is_R(r))
    60   {
    61     nMap = nlMapR; /* short R -> Q */
    62     return TRUE;
    63   }
    64 //  if (rField_is_long_R(r))
    65 //  {
    66 //    nMap = nlMapLongR; /* long R -> Q */
    67 //    return TRUE;
    68 //  }
    69   return FALSE;
    70 }
    7133
    7234/*-----------------------------------------------------------------*/
    … …  
    13193#endif
    13294
    133 
    13495#ifdef LDEBUG
    13596#define nlTest(a) nlDBTest(a,__FILE__,__LINE__)
     97BOOLEAN nlDBTest(number a, char *f,int l);
     98#else
     99#define nlTest(a) ((void)0)
     100#endif
     101
     102
     103/***************************************************************
     104 * 
     105 * Routines which are never inlined by p_Numbers.h
     106 *
     107 *******************************************************************/
     108#ifndef P_NUMBERS_H
     109
     110omBin rnumber_bin = omGetSpecBin(sizeof(rnumber));
     111
     112static int nlPrimeM;
     113static number nlMapP(number from)
     114{
     115  number to;
     116  int save=npPrimeM;
     117  npPrimeM=nlPrimeM;
     118  to = nlInit(npInt(from));
     119  npPrimeM=save;
     120  return to;
     121}
     122
     123//static number nlMapLongR(number from);
     124static number nlMapR(number from);
     125
     126BOOLEAN nlSetMap(ring r)
     127{
     128  if (rField_is_Q(r))
     129  {
     130    nMap = nlCopy;   /*Q -> Q*/
     131    return TRUE;
     132  }
     133  if (rField_is_Zp(r))
     134  {
     135    nlPrimeM=rChar(r);
     136    nMap = nlMapP; /* Z/p -> Q */
     137    return TRUE;
     138  }
     139  if (rField_is_R(r))
     140  {
     141    nMap = nlMapR; /* short R -> Q */
     142    return TRUE;
     143  }
     144//  if (rField_is_long_R(r))
     145//  {
     146//    nMap = nlMapLongR; /* long R -> Q */
     147//    return TRUE;
     148//  }
     149  return FALSE;
     150}
     151
     152#ifdef LDEBUG
    136153BOOLEAN nlDBTest(number a, char *f,int l)
    137154{
    … …  
    219236  return TRUE;
    220237}
    221 #else
    222 #define nlTest(x) (TRUE)
    223 #endif
    224 
    225 void nlNew (number * r)
    226 {
    227   *r=NULL;
    228 }
    229 
    230 /*2
    231 * z := i
    232 */
    233 static inline number nlRInit (int i)
    234 {
    235   number z=(number)omAllocBin(rnumber_bin);
    236 #if defined(LDEBUG)
    237   z->debug=123456;
    238 #endif
    239   mpz_init_set_si(&z->z,(long)i);
    240   z->s = 3;
    241   return z;
    242 }
    243 
    244 number nlInit (int i)
    245 {
    246   number n;
    247   if ( ((i << 3) >> 3) == i ) n=INT_TO_SR(i);
    248   else                        n=nlRInit(i);
    249 #ifdef LDEBUG
    250   nlTest(n);
    251 #endif
    252   return n;
    253 }
    254 
     238#endif
     239
     240
     241static inline number nlRInit (int i);
    255242
    256243number nlInit (number u)
    … …  
    427414}
    428415
    429 /*2
    430 * delete a
    431 */
    432 #ifdef LDEBUG
    433 void nlDBDelete (number * a,char *f, int l)
    434 #else
    435 void nlDelete (number * a)
    436 #endif
    437 {
    438   if (*a!=NULL)
    439   {
    440 #ifdef LDEBUG
    441     nlTest(*a);
    442 #endif
    443     if ((SR_HDL(*a) & SR_INT)==0)
    444     {
    445       switch ((*a)->s)
    446       {
    447         case 0:
    448         case 1:
    449           mpz_clear(&(*a)->n);
    450         case 3:
    451           mpz_clear(&(*a)->z);
    452 #ifdef LDEBUG
    453           (*a)->s=2;
    454 #endif
    455       }
    456       omFreeBin((ADDRESS) *a, rnumber_bin);
    457     }
    458     *a=NULL;
    459   }
    460 }
    461 
    462 /*2
    463 * copy a to b
    464 */
    465 number nlCopy(number a)
    466 {
    467   number b;
    468   if ((SR_HDL(a) & SR_INT)||(a==NULL))
    469   {
    470     return a;
    471   }
    472 #ifdef LDEBUG
    473   nlTest(a);
    474 #endif
    475   b=(number)omAllocBin(rnumber_bin);
    476 #if defined(LDEBUG)
    477   b->debug=123456;
    478 #endif
    479   switch (a->s)
    480   {
    481     case 0:
    482     case 1:
    483             nlGmpSimple(&a->n);
    484             mpz_init_set(&b->n,&a->n);
    485     case 3:
    486             nlGmpSimple(&a->z);
    487             mpz_init_set(&b->z,&a->z);
    488             break;
    489   }
    490   b->s = a->s;
    491 #ifdef LDEBUG
    492   nlTest(b);
    493 #endif
    494   return b;
    495 }
    496 
    497 /*2
    498 * za:= - za
    499 */
    500 number nlNeg (number a)
    501 {
    502 #ifdef LDEBUG
    503   nlTest(a);
    504 #endif
    505   if(SR_HDL(a) &SR_INT)
    506   {
    507     int r=SR_TO_INT(a);
    508     if (r==(-(1<<28))) a=nlRInit(1<<28);
    509     else               a=INT_TO_SR(-r);
    510   }
    511   else
    512   {
    513     mpz_neg(&a->z,&a->z);
    514     if ((a->s==3) && (mpz_size1(&a->z)<=MP_SMALL))
    515     {
    516       int ui=(int)mpz_get_si(&a->z);
    517       if ((((ui<<3)>>3)==ui)
    518       && (mpz_cmp_si(&a->z,(long)ui)==0))
    519       {
    520         mpz_clear(&a->z);
    521         omFreeBin((ADDRESS)a, rnumber_bin);
    522         a=INT_TO_SR(ui);
    523       }
    524     }
    525   }
    526 #ifdef LDEBUG
    527   nlTest(a);
    528 #endif
    529   return a;
    530 }
     416
    531417
    532418/*
    … …  
    627513}
    628514
     515
     516/*2
     517* u := a / b in Z, if b | a (else undefined)
     518*/
     519number   nlExactDiv(number a, number b)
     520{
     521  if (b==INT_TO_SR(0))
     522  {
     523    WerrorS("div. by 0");
     524    return INT_TO_SR(0);
     525  }
     526  if (a==INT_TO_SR(0))
     527    return INT_TO_SR(0);
     528  number u;
     529  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
     530  {
     531    /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
     532    if ((a==INT_TO_SR(-(1<<28)))&&(b==INT_TO_SR(-1)))
     533    {
     534      return nlRInit(1<<28);
     535    }
     536    int aa=SR_TO_INT(a);
     537    int bb=SR_TO_INT(b);
     538    return INT_TO_SR(aa/bb);
     539  }
     540  number bb=NULL;
     541  if (SR_HDL(b) & SR_INT)
     542  {
     543    bb=nlRInit(SR_TO_INT(b));
     544    b=bb;
     545  }
     546  u=(number)omAllocBin(rnumber_bin);
     547#if defined(LDEBUG)
     548  u->debug=123456;
     549#endif
     550  mpz_init(&u->z);
     551  /* u=a/b */
     552  u->s = 3;
     553  MPZ_EXACTDIV(&u->z,&a->z,&b->z);
     554  if (bb!=NULL)
     555  {
     556    mpz_clear(&bb->z);
     557    omFreeBin((ADDRESS)bb, rnumber_bin);
     558  }
     559  if (mpz_size1(&u->z)<=MP_SMALL)
     560  {
     561    int ui=(int)mpz_get_si(&u->z);
     562    if ((((ui<<3)>>3)==ui)
     563    && (mpz_cmp_si(&u->z,(long)ui)==0))
     564    {
     565      mpz_clear(&u->z);
     566      omFreeBin((ADDRESS)u, rnumber_bin);
     567      u=INT_TO_SR(ui);
     568    }
     569  }
     570#ifdef LDEBUG
     571  nlTest(u);
     572#endif
     573  return u;
     574}
     575
     576/*2
     577* u := a / b in Z
     578*/
     579number nlIntDiv (number a, number b)
     580{
     581  if (b==INT_TO_SR(0))
     582  {
     583    WerrorS("div. by 0");
     584    return INT_TO_SR(0);
     585  }
     586  if (a==INT_TO_SR(0))
     587    return INT_TO_SR(0);
     588  number u;
     589  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
     590  {
     591    /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
     592    if ((a==INT_TO_SR(-(1<<28)))&&(b==INT_TO_SR(-1)))
     593    {
     594      return nlRInit(1<<28);
     595    }
     596    /* consider the signs of a and b:
     597    *  + + -> u=a+b-1
     598    *  + - -> u=a-b-1
     599    *  - + -> u=a-b+1
     600    *  - - -> u=a+b+1
     601    */
     602    int aa=SR_TO_INT(a);
     603    int bb=SR_TO_INT(b);
     604    if (aa>0)
     605    {
     606      if (bb>0)
     607        return INT_TO_SR((aa+bb-1)/bb);
     608      else
     609        return INT_TO_SR((aa-bb-1)/bb);
     610    }
     611    else
     612    {
     613      if (bb>0)
     614      {
     615        return INT_TO_SR((aa-bb+1)/bb);
     616      }
     617      else
     618      {
     619        return INT_TO_SR((aa+bb+1)/bb);
     620      }
     621    }
     622  }
     623  if (SR_HDL(a) & SR_INT)
     624  {
     625    /* the small int -(1<<28) divided by 2^28 is 1   */
     626    if (a==INT_TO_SR(-(1<<28)))
     627    {
     628      if(mpz_cmp_si(&b->z,(long)(1<<28))==0)
     629      {
     630        return INT_TO_SR(-1);
     631      }
     632    }
     633    /* a is a small and b is a large int: -> 0 */
     634    return INT_TO_SR(0);
     635  }
     636  number bb=NULL;
     637  if (SR_HDL(b) & SR_INT)
     638  {
     639    bb=nlRInit(SR_TO_INT(b));
     640    b=bb;
     641  }
     642  u=(number)omAllocBin(rnumber_bin);
     643#if defined(LDEBUG)
     644  u->debug=123456;
     645#endif
     646  mpz_init_set(&u->z,&a->z);
     647  /* consider the signs of a and b:
     648  *  + + -> u=a+b-1
     649  *  + - -> u=a-b-1
     650  *  - + -> u=a-b+1
     651  *  - - -> u=a+b+1
     652  */
     653  if (mpz_isNeg(&a->z))
     654  {
     655    if (mpz_isNeg(&b->z))
     656    {
     657      mpz_add(&u->z,&u->z,&b->z);
     658    }
     659    else
     660    {
     661      mpz_sub(&u->z,&u->z,&b->z);
     662    }
     663    mpz_add_ui(&u->z,&u->z,1);
     664  }
     665  else
     666  {
     667    if (mpz_isNeg(&b->z))
     668    {
     669      mpz_sub(&u->z,&u->z,&b->z);
     670    }
     671    else
     672    {
     673      mpz_add(&u->z,&u->z,&b->z);
     674    }
     675    mpz_sub_ui(&u->z,&u->z,1);
     676  }
     677  /* u=u/b */
     678  u->s = 3;
     679  MPZ_DIV(&u->z,&u->z,&b->z);
     680  nlGmpSimple(&u->z);
     681  if (bb!=NULL)
     682  {
     683    mpz_clear(&bb->z);
     684    omFreeBin((ADDRESS)bb, rnumber_bin);
     685  }
     686  if (mpz_size1(&u->z)<=MP_SMALL)
     687  {
     688    int ui=(int)mpz_get_si(&u->z);
     689    if ((((ui<<3)>>3)==ui)
     690    && (mpz_cmp_si(&u->z,(long)ui)==0))
     691    {
     692      mpz_clear(&u->z);
     693      omFreeBin((ADDRESS)u, rnumber_bin);
     694      u=INT_TO_SR(ui);
     695    }
     696  }
     697#ifdef LDEBUG
     698  nlTest(u);
     699#endif
     700  return u;
     701}
     702
     703/*2
     704* u := a mod b in Z, u>=0
     705*/
     706number nlIntMod (number a, number b)
     707{
     708  if (b==INT_TO_SR(0))
     709  {
     710    WerrorS("div. by 0");
     711    return INT_TO_SR(0);
     712  }
     713  if (a==INT_TO_SR(0))
     714    return INT_TO_SR(0);
     715  number u;
     716  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
     717  {
     718    if ((int)a>0)
     719    {
     720      if ((int)b>0)
     721        return INT_TO_SR(SR_TO_INT(a)%SR_TO_INT(b));
     722      else
     723        return INT_TO_SR(SR_TO_INT(a)%(-SR_TO_INT(b)));
     724    }
     725    else
     726    {
     727      if ((int)b>0)
     728      {
     729        int i=(-SR_TO_INT(a))%SR_TO_INT(b);
     730        if ( i != 0 ) i = (SR_TO_INT(b))-i;
     731        return INT_TO_SR(i);
     732      }
     733      else
     734      {
     735        int i=(-SR_TO_INT(a))%(-SR_TO_INT(b));
     736        if ( i != 0 ) i = (-SR_TO_INT(b))-i;
     737        return INT_TO_SR(i);
     738      }
     739    }
     740  }
     741  if (SR_HDL(a) & SR_INT)
     742  {
     743    /* a is a small and b is a large int: -> a or (a+b) or (a-b) */
     744    if ((int)a<0)
     745    {
     746      if (mpz_isNeg(&b->z))
     747        return nlSub(a,b);
     748      /*else*/
     749        return nlAdd(a,b);
     750    }
     751    /*else*/
     752      return a;
     753  }
     754  number bb=NULL;
     755  if (SR_HDL(b) & SR_INT)
     756  {
     757    bb=nlRInit(SR_TO_INT(b));
     758    b=bb;
     759  }
     760  u=(number)omAllocBin(rnumber_bin);
     761#if defined(LDEBUG)
     762  u->debug=123456;
     763#endif
     764  mpz_init(&u->z);
     765  u->s = 3;
     766  mpz_mod(&u->z,&a->z,&b->z);
     767  if (bb!=NULL)
     768  {
     769    mpz_clear(&bb->z);
     770    omFreeBin((ADDRESS)bb, rnumber_bin);
     771  }
     772  if (mpz_isNeg(&u->z))
     773  {
     774    if (mpz_isNeg(&b->z))
     775      mpz_sub(&u->z,&u->z,&b->z);
     776    else
     777      mpz_add(&u->z,&u->z,&b->z);
     778  }
     779  nlGmpSimple(&u->z);
     780  if (mpz_size1(&u->z)<=MP_SMALL)
     781  {
     782    int ui=(int)mpz_get_si(&u->z);
     783    if ((((ui<<3)>>3)==ui)
     784    && (mpz_cmp_si(&u->z,(long)ui)==0))
     785    {
     786      mpz_clear(&u->z);
     787      omFreeBin((ADDRESS)u, rnumber_bin);
     788      u=INT_TO_SR(ui);
     789    }
     790  }
     791#ifdef LDEBUG
     792  nlTest(u);
     793#endif
     794  return u;
     795}
     796
     797/*2
     798* u := a / b
     799*/
     800number nlDiv (number a, number b)
     801{
     802  number u;
     803  if (nlIsZero(b))
     804  {
     805    WerrorS("div. by 0");
     806    return INT_TO_SR(0);
     807  }
     808  u=(number)omAllocBin(rnumber_bin);
     809  u->s=0;
     810#if defined(LDEBUG)
     811  u->debug=123456;
     812#endif
     813// ---------- short / short ------------------------------------
     814  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
     815  {
     816    int i=SR_TO_INT(a);
     817    int j=SR_TO_INT(b);
     818    int r=i%j;
     819    if (r==0)
     820    {
     821      omFreeBin((ADDRESS)u, rnumber_bin);
     822      return INT_TO_SR(i/j);
     823    }
     824    mpz_init_set_si(&u->z,(long)i);
     825    mpz_init_set_si(&u->n,(long)j);
     826  }
     827  else
     828  {
     829    mpz_init(&u->z);
     830// ---------- short / long ------------------------------------
     831    if (SR_HDL(a) & SR_INT)
     832    {
     833      // short a / (z/n) -> (a*n)/z
     834      if (b->s<2)
     835      {
     836        if ((int)a>=0)
     837          mpz_mul_ui(&u->z,&b->n,SR_TO_INT(a));
     838        else
     839        {
     840          mpz_mul_ui(&u->z,&b->n,-SR_TO_INT(a));
     841          mpz_neg(&u->z,&u->z);
     842        }
     843      }
     844      else
     845      // short a / long z -> a/z
     846      {
     847        mpz_set_si(&u->z,SR_TO_INT(a));
     848      }
     849      if (mpz_cmp(&u->z,&b->z)==0)
     850      {
     851        mpz_clear(&u->z);
     852        omFreeBin((ADDRESS)u, rnumber_bin);
     853        return INT_TO_SR(1);
     854      }
     855      mpz_init_set(&u->n,&b->z);
     856    }
     857// ---------- long / short ------------------------------------
     858    else if (SR_HDL(b) & SR_INT)
     859    {
     860      mpz_set(&u->z,&a->z);
     861      // (z/n) / b -> z/(n*b)
     862      if (a->s<2)
     863      {
     864        mpz_init_set(&u->n,&a->n);
     865        if ((int)b>0)
     866          mpz_mul_ui(&u->n,&u->n,SR_TO_INT(b));
     867        else
     868        {
     869          mpz_mul_ui(&u->n,&u->n,-SR_TO_INT(b));
     870          mpz_neg(&u->z,&u->z);
     871        }
     872      }
     873      else
     874      // long z / short b -> z/b
     875      {
     876        //mpz_set(&u->z,&a->z);
     877        mpz_init_set_si(&u->n,SR_TO_INT(b));
     878      }
     879    }
     880// ---------- long / long ------------------------------------
     881    else
     882    {
     883      //u->s=0;
     884      mpz_set(&u->z,&a->z);
     885      mpz_init_set(&u->n,&b->z);
     886      if (a->s<2) mpz_mul(&u->n,&u->n,&a->n);
     887      if (b->s<2) mpz_mul(&u->z,&u->z,&b->n);
     888    }
     889  }
     890  if (mpz_isNeg(&u->n))
     891  {
     892    mpz_neg(&u->z,&u->z);
     893    mpz_neg(&u->n,&u->n);
     894  }
     895  if (mpz_cmp_si(&u->n,(long)1)==0)
     896  {
     897    mpz_clear(&u->n);
     898    if (mpz_size1(&u->z)<=MP_SMALL)
     899    {
     900      int ui=(int)mpz_get_si(&u->z);
     901      if ((((ui<<3)>>3)==ui)
     902      && (mpz_cmp_si(&u->z,(long)ui)==0))
     903      {
     904        mpz_clear(&u->z);
     905        omFreeBin((ADDRESS)u, rnumber_bin);
     906        return INT_TO_SR(ui);
     907      }
     908    }
     909    u->s=3;
     910  }
     911#ifdef LDEBUG
     912  nlTest(u);
     913#endif
     914  return u;
     915}
     916
     917#if (defined(i386)) && (defined(HAVE_LIBGMP1))
     918/*
     919* compute x^exp for x in Z
     920* there seems to be an bug in mpz_pow_ui for i386
     921*/
     922static inline void nlPow (MP_INT * res,MP_INT * x,int exp)
     923{
     924  if (exp==0)
     925  {
     926    mpz_set_si(res,(long)1);
     927  }
     928  else
     929  {
     930    MP_INT xh;
     931    mpz_init(&xh);
     932    mpz_set(res,x);
     933    exp--;
     934    while (exp!=0)
     935    {
     936      mpz_mul(&xh,res,x);
     937      mpz_set(res,&xh);
     938      exp--;
     939    }
     940    mpz_clear(&xh);
     941  }
     942}
     943#endif
     944
     945/*2
     946* u:= x ^ exp
     947*/
     948void nlPower (number x,int exp,number * u)
     949{
     950  if (!nlIsZero(x))
     951  {
     952#ifdef LDEBUG
     953    nlTest(x);
     954#endif
     955    number aa=NULL;
     956    if (SR_HDL(x) & SR_INT)
     957    {
     958      aa=nlRInit(SR_TO_INT(x));
     959      x=aa;
     960    }
     961    *u=(number)omAllocBin(rnumber_bin);
     962#if defined(LDEBUG)
     963    (*u)->debug=123456;
     964#endif
     965    mpz_init(&(*u)->z);
     966    (*u)->s = x->s;
     967#if (!defined(i386)) || (defined(HAVE_LIBGMP2))
     968    mpz_pow_ui(&(*u)->z,&x->z,(unsigned long)exp);
     969#else
     970    nlPow(&(*u)->z,&x->z,exp);
     971#endif
     972    if (x->s<2)
     973    {
     974      mpz_init(&(*u)->n);
     975#if (!defined(i386)) || (defined(HAVE_LIBGMP2))
     976      mpz_pow_ui(&(*u)->n,&x->n,(unsigned long)exp);
     977#else
     978      nlPow(&(*u)->n,&x->n,exp);
     979#endif
     980    }
     981    if (aa!=NULL)
     982    {
     983      mpz_clear(&aa->z);
     984      omFreeBin((ADDRESS)aa, rnumber_bin);
     985    }
     986    if (((*u)->s<=2) && (mpz_cmp_si(&(*u)->n,(long)1)==0))
     987    {
     988      mpz_clear(&(*u)->n);
     989      (*u)->s=3;
     990    }
     991    if (((*u)->s==3) && (mpz_size1(&(*u)->z)<=MP_SMALL))
     992    {
     993      int ui=(int)mpz_get_si(&(*u)->z);
     994      if ((((ui<<3)>>3)==ui)
     995      && (mpz_cmp_si(&(*u)->z,(long)ui)==0))
     996      {
     997        mpz_clear(&(*u)->z);
     998        omFreeBin((ADDRESS)*u, rnumber_bin);
     999        *u=INT_TO_SR(ui);
     1000      }
     1001    }
     1002  }
     1003  else
     1004    *u = INT_TO_SR(0);
     1005#ifdef LDEBUG
     1006  nlTest(*u);
     1007#endif
     1008}
     1009
     1010
     1011/*2
     1012* za >= 0 ?
     1013*/
     1014BOOLEAN nlGreaterZero (number a)
     1015{
     1016#ifdef LDEBUG
     1017  nlTest(a);
     1018#endif
     1019  if (SR_HDL(a) & SR_INT) return SR_HDL(a)>=0;
     1020  return (!mpz_isNeg(&a->z));
     1021}
     1022
     1023/*2
     1024* a > b ?
     1025*/
     1026BOOLEAN nlGreater (number a, number b)
     1027{
     1028#ifdef LDEBUG
     1029  nlTest(a);
     1030  nlTest(b);
     1031#endif
     1032  number r;
     1033  BOOLEAN rr;
     1034  r=nlSub(a,b);
     1035  rr=(!nlIsZero(r)) && (nlGreaterZero(r));
     1036  nlDelete(&r);
     1037  return rr;
     1038}
     1039
     1040/*2
     1041* a == -1 ?
     1042*/
     1043BOOLEAN nlIsMOne (number a)
     1044{
     1045#ifdef LDEBUG
     1046  if (a==NULL) return FALSE;
     1047  nlTest(a);
     1048#endif
     1049  if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(-1));
     1050  return FALSE;
     1051}
     1052
     1053/*2
     1054* result =gcd(a,b)
     1055*/
     1056number nlGcd(number a, number b)
     1057{
     1058  number result;
     1059#ifdef LDEBUG
     1060  nlTest(a);
     1061  nlTest(b);
     1062#endif
     1063  //nlNormalize(a);
     1064  //nlNormalize(b);
     1065  if ((SR_HDL(a)==5)||(a==INT_TO_SR(-1))
     1066  ||  (SR_HDL(b)==5)||(b==INT_TO_SR(-1))) return INT_TO_SR(1);
     1067  if (SR_HDL(a) & SR_HDL(b) & SR_INT)
     1068  {
     1069    int i=SR_TO_INT(a);
     1070    int j=SR_TO_INT(b);
     1071    if((i==0)||(j==0))
     1072      return INT_TO_SR(1);
     1073    int l;
     1074    i=ABS(i);
     1075    j=ABS(j);
     1076    do
     1077    {
     1078      l=i%j;
     1079      i=j;
     1080      j=l;
     1081    } while (l!=0);
     1082    return INT_TO_SR(i);
     1083  }
     1084  if (((!(SR_HDL(a) & SR_INT))&&(a->s<2))
     1085  ||  ((!(SR_HDL(b) & SR_INT))&&(b->s<2))) return INT_TO_SR(1);
     1086  number aa=NULL;
     1087  if (SR_HDL(a) & SR_INT)
     1088  {
     1089    aa=nlRInit(SR_TO_INT(a));
     1090    a=aa;
     1091  }
     1092  else
     1093  if (SR_HDL(b) & SR_INT)
     1094  {
     1095    aa=nlRInit(SR_TO_INT(b));
     1096    b=aa;
     1097  }
     1098  result=(number)omAllocBin(rnumber_bin);
     1099#if defined(LDEBUG)
     1100  result->debug=123456;
     1101#endif
     1102  mpz_init(&result->z);
     1103  mpz_gcd(&result->z,&a->z,&b->z);
     1104  nlGmpSimple(&result->z);
     1105  result->s = 3;
     1106  if (mpz_size1(&result->z)<=MP_SMALL)
     1107  {
     1108    int ui=(int)mpz_get_si(&result->z);
     1109    if ((((ui<<3)>>3)==ui)
     1110    && (mpz_cmp_si(&result->z,(long)ui)==0))
     1111    {
     1112      mpz_clear(&result->z);
     1113      omFreeBin((ADDRESS)result, rnumber_bin);
     1114      result=INT_TO_SR(ui);
     1115    }
     1116  }
     1117  if (aa!=NULL)
     1118  {
     1119    mpz_clear(&aa->z);
     1120    omFreeBin((ADDRESS)aa, rnumber_bin);
     1121  }
     1122#ifdef LDEBUG
     1123  nlTest(result);
     1124#endif
     1125  return result;
     1126}
     1127
     1128/*2
     1129* simplify x
     1130*/
     1131void nlNormalize (number &x)
     1132{
     1133  if ((SR_HDL(x) & SR_INT) ||(x==NULL))
     1134    return;
     1135#ifdef LDEBUG
     1136  if (!nlTest(x)) return;
     1137#endif
     1138  if (x->s==3)
     1139  {
     1140    if (mpz_cmp_ui(&x->z,(long)0)==0)
     1141    {
     1142      nlDelete(&x);
     1143      x=nlInit(0);
     1144      return;
     1145    }
     1146    if (mpz_size1(&x->z)<=MP_SMALL)
     1147    {
     1148      int ui=(int)mpz_get_si(&x->z);
     1149      if ((((ui<<3)>>3)==ui)
     1150      && (mpz_cmp_si(&x->z,(long)ui)==0))
     1151      {
     1152        mpz_clear(&x->z);
     1153        omFreeBin((ADDRESS)x, rnumber_bin);
     1154        x=INT_TO_SR(ui);
     1155        return;
     1156      }
     1157    }
     1158  }
     1159  if (x->s!=1) nlGmpSimple(&x->z);
     1160  if (x->s==0)
     1161  {
     1162    BOOLEAN divided=FALSE;
     1163    MP_INT gcd;
     1164    mpz_init(&gcd);
     1165    mpz_gcd(&gcd,&x->z,&x->n);
     1166    x->s=1;
     1167    if (mpz_cmp_si(&gcd,(long)1)!=0)
     1168    {
     1169      MP_INT r;
     1170      mpz_init(&r);
     1171      MPZ_EXACTDIV(&r,&x->z,&gcd);
     1172      mpz_set(&x->z,&r);
     1173      MPZ_EXACTDIV(&r,&x->n,&gcd);
     1174      mpz_set(&x->n,&r);
     1175      mpz_clear(&r);
     1176      nlGmpSimple(&x->z);
     1177      divided=TRUE;
     1178    }
     1179    mpz_clear(&gcd);
     1180    nlGmpSimple(&x->n);
     1181    if (mpz_cmp_si(&x->n,(long)1)==0)
     1182    {
     1183      mpz_clear(&x->n);
     1184      if (mpz_size1(&x->z)<=MP_SMALL)
     1185      {
     1186        int ui=(int)mpz_get_si(&x->z);
     1187        if ((((ui<<3)>>3)==ui)
     1188        && (mpz_cmp_si(&x->z,(long)ui)==0))
     1189        {
     1190          mpz_clear(&x->z);
     1191#if defined(LDEBUG)
     1192          x->debug=654324;
     1193#endif
     1194          omFreeBin((ADDRESS)x, rnumber_bin);
     1195          x=INT_TO_SR(ui);
     1196          return;
     1197        }
     1198      }
     1199      x->s=3;
     1200    }
     1201    else if (divided)
     1202    {
     1203      _mpz_realloc(&x->n,mpz_size1(&x->n));
     1204    }
     1205    if (divided) _mpz_realloc(&x->z,mpz_size1(&x->z));
     1206  }
     1207#ifdef LDEBUG
     1208  nlTest(x);
     1209#endif
     1210}
     1211
     1212/*2
     1213* returns in result->z the lcm(a->z,b->n)
     1214*/
     1215number nlLcm(number a, number b)
     1216{
     1217  number result;
     1218#ifdef LDEBUG
     1219  nlTest(a);
     1220  nlTest(b);
     1221#endif
     1222  if ((SR_HDL(b) & SR_INT)
     1223  || (b->s==3))
     1224  {
     1225    // b is 1/(b->n) => b->n is 1 => result is a
     1226    return nlCopy(a);
     1227  }
     1228  number aa=NULL;
     1229  if (SR_HDL(a) & SR_INT)
     1230  {
     1231    aa=nlRInit(SR_TO_INT(a));
     1232    a=aa;
     1233  }
     1234  result=(number)omAllocBin(rnumber_bin);
     1235#if defined(LDEBUG)
     1236  result->debug=123456;
     1237#endif
     1238  result->s=3;
     1239  MP_INT gcd;
     1240  mpz_init(&gcd);
     1241  mpz_init(&result->z);
     1242  mpz_gcd(&gcd,&a->z,&b->n);
     1243  if (mpz_cmp_si(&gcd,(long)1)!=0)
     1244  {
     1245    MP_INT bt;
     1246    mpz_init_set(&bt,&b->n);
     1247    MPZ_EXACTDIV(&bt,&bt,&gcd);
     1248    mpz_mul(&result->z,&a->z,&bt);
     1249    mpz_clear(&bt);
     1250  }
     1251  else
     1252    mpz_mul(&result->z,&a->z,&b->n);
     1253  mpz_clear(&gcd);
     1254  if (aa!=NULL)
     1255  {
     1256    mpz_clear(&aa->z);
     1257    omFreeBin((ADDRESS)aa, rnumber_bin);
     1258  }
     1259  nlGmpSimple(&result->z);
     1260  if (mpz_size1(&result->z)<=MP_SMALL)
     1261  {
     1262    int ui=(int)mpz_get_si(&result->z);
     1263    if ((((ui<<3)>>3)==ui)
     1264    && (mpz_cmp_si(&result->z,(long)ui)==0))
     1265    {
     1266      mpz_clear(&result->z);
     1267      omFreeBin((ADDRESS)result, rnumber_bin);
     1268      return INT_TO_SR(ui);
     1269    }
     1270  }
     1271#ifdef LDEBUG
     1272  nlTest(result);
     1273#endif
     1274  return result;
     1275}
     1276
     1277int nlModP(number n, int p)
     1278{
     1279  if (SR_HDL(n) & SR_INT)
     1280  {
     1281    int i=SR_TO_INT(n);
     1282    if (i<0) return (p-((-i)%p));
     1283    return i%p;
     1284  }
     1285  int iz=mpz_mmod_ui(NULL,&n->z,(unsigned long)p);
     1286  if (n->s!=3)
     1287  {
     1288    int in=mpz_mmod_ui(NULL,&n->n,(unsigned long)p);
     1289    return (int)npDiv((number)iz,(number)in);
     1290  }
     1291  return iz;
     1292}
     1293
     1294/*2
     1295* acces to denominator, other 1 for integers
     1296*/
     1297number   nlGetDenom(number &n)
     1298{
     1299  if (!(SR_HDL(n) & SR_INT))
     1300  {
     1301    if (n->s==0)
     1302    {
     1303      nlNormalize(n);
     1304    }
     1305    if (!(SR_HDL(n) & SR_INT))
     1306    {
     1307      if (n->s!=3)
     1308      {
     1309        number u=(number)omAllocBin(rnumber_bin);
     1310        u->s=3;
     1311#if defined(LDEBUG)
     1312        u->debug=123456;
     1313#endif
     1314
     1315        mpz_init_set(&u->z,&n->n);
     1316        {
     1317          int ui=(int)mpz_get_si(&u->z);
     1318          if ((((ui<<3)>>3)==ui)
     1319          && (mpz_cmp_si(&u->z,(long)ui)==0))
     1320          {
     1321            mpz_clear(&u->z);
     1322            omFreeBin((ADDRESS)u, rnumber_bin);
     1323            return INT_TO_SR(ui);
     1324          }
     1325        }
     1326        return u;
     1327      }
     1328    }
     1329  }
     1330  return INT_TO_SR(1);
     1331}
     1332
     1333#endif
     1334
     1335
     1336/***************************************************************
     1337 * 
     1338 * Routines which might be inlined by p_Numbers.h
     1339 *
     1340 *******************************************************************/
     1341#if defined(DO_LINLINE) || !defined(P_NUMBERS_H)
     1342
     1343// routines which are always inlined/static
     1344/*2
     1345* z := i
     1346*/
     1347static inline number nlRInit (int i)
     1348{
     1349  number z=(number)omAllocBin(rnumber_bin);
     1350#if defined(LDEBUG)
     1351  z->debug=123456;
     1352#endif
     1353  mpz_init_set_si(&z->z,(long)i);
     1354  z->s = 3;
     1355  return z;
     1356}
     1357
     1358/*2
     1359* a = b ?
     1360*/
     1361LINLINE BOOLEAN nlEqual (number a, number b)
     1362{
     1363#ifdef LDEBUG
     1364  nlTest(a);
     1365  nlTest(b);
     1366#endif
     1367// short - short
     1368  if (SR_HDL(a) & SR_HDL(b) & SR_INT) return a==b;
     1369//  long - short
     1370  BOOLEAN bo;
     1371  if (SR_HDL(b) & SR_INT)
     1372  {
     1373    if (a->s!=0) return FALSE;
     1374    number n=b; b=a; a=n;
     1375  }
     1376//  short - long
     1377  if (SR_HDL(a) & SR_INT)
     1378  {
     1379    if (b->s!=0)
     1380      return FALSE;
     1381    if (((int)a > 0) && (mpz_isNeg(&b->z)))
     1382      return FALSE;
     1383    if (((int)a < 0) && (!mpz_isNeg(&b->z)))
     1384      return FALSE;
     1385    MP_INT  bb;
     1386    mpz_init_set(&bb,&b->n);
     1387    if ((int)a<0)
     1388    {
     1389      mpz_neg(&bb,&bb);
     1390      mpz_mul_ui(&bb,&bb,(long)-SR_TO_INT(a));
     1391    }
     1392    else
     1393    {
     1394      mpz_mul_ui(&bb,&bb,(long)SR_TO_INT(a));
     1395    }
     1396    bo=(mpz_cmp(&bb,&b->z)==0);
     1397    mpz_clear(&bb);
     1398    return bo;
     1399  }
     1400// long - long
     1401  if (((a->s==1) && (b->s==3))
     1402  ||  ((b->s==1) && (a->s==3)))
     1403    return FALSE;
     1404  if (mpz_isNeg(&a->z)&&(!mpz_isNeg(&b->z)))
     1405    return FALSE;
     1406  if (mpz_isNeg(&b->z)&&(!mpz_isNeg(&a->z)))
     1407    return FALSE;
     1408  nlGmpSimple(&a->z);
     1409  nlGmpSimple(&b->z);
     1410  if (a->s<2)
     1411    nlGmpSimple(&a->n);
     1412  if (b->s<2)
     1413    nlGmpSimple(&b->n);
     1414  MP_INT  aa;
     1415  MP_INT  bb;
     1416  mpz_init_set(&aa,&a->z);
     1417  mpz_init_set(&bb,&b->z);
     1418  if (a->s<2) mpz_mul(&bb,&bb,&a->n);
     1419  if (b->s<2) mpz_mul(&aa,&aa,&b->n);
     1420  bo=(mpz_cmp(&aa,&bb)==0);
     1421  mpz_clear(&aa);
     1422  mpz_clear(&bb);
     1423  return bo;
     1424}
     1425
     1426LINLINE number nlInit (int i)
     1427{
     1428  number n;
     1429  if ( ((i << 3) >> 3) == i ) n=INT_TO_SR(i);
     1430  else                        n=nlRInit(i);
     1431#ifdef LDEBUG
     1432  nlTest(n);
     1433#endif
     1434  return n;
     1435}
     1436
     1437
     1438/*2
     1439* a == 1 ?
     1440*/
     1441LINLINE BOOLEAN nlIsOne (number a)
     1442{
     1443#ifdef LDEBUG
     1444  if (a==NULL) return FALSE;
     1445  nlTest(a);
     1446#endif
     1447  if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(1));
     1448  return FALSE;
     1449}
     1450
     1451LINLINE BOOLEAN nlIsZero (number a)
     1452{
     1453// #ifdef LDEBUG
     1454//   nlTest(a);
     1455// #endif
     1456  //if (a==INT_TO_SR(0)) return TRUE;
     1457  //if (SR_HDL(a) & SR_INT) return FALSE;
     1458  //number aa=nlCopy(a);
     1459  //nlNormalize(aa);
     1460  return (a==INT_TO_SR(0));
     1461}
     1462
     1463
     1464/*2
     1465* copy a to b
     1466*/
     1467LINLINE number nlCopy(number a)
     1468{
     1469  number b;
     1470  if ((SR_HDL(a) & SR_INT)||(a==NULL))
     1471  {
     1472    return a;
     1473  }
     1474#ifdef LDEBUG
     1475  nlTest(a);
     1476#endif
     1477  b=(number)omAllocBin(rnumber_bin);
     1478#if defined(LDEBUG)
     1479  b->debug=123456;
     1480#endif
     1481  switch (a->s)
     1482  {
     1483    case 0:
     1484    case 1:
     1485            nlGmpSimple(&a->n);
     1486            mpz_init_set(&b->n,&a->n);
     1487    case 3:
     1488            nlGmpSimple(&a->z);
     1489            mpz_init_set(&b->z,&a->z);
     1490            break;
     1491  }
     1492  b->s = a->s;
     1493#ifdef LDEBUG
     1494  nlTest(b);
     1495#endif
     1496  return b;
     1497}
     1498
     1499LINLINE void nlNew (number * r)
     1500{
     1501  *r=NULL;
     1502}
     1503
     1504/*2
     1505* delete a
     1506*/
     1507#ifdef LDEBUG
     1508void nlDBDelete (number * a,char *f, int l)
     1509#else
     1510LINLINE void nlDelete (number * a)
     1511#endif
     1512{
     1513  if (*a!=NULL)
     1514  {
     1515#ifdef LDEBUG
     1516    nlTest(*a);
     1517#endif
     1518    if ((SR_HDL(*a) & SR_INT)==0)
     1519    {
     1520      switch ((*a)->s)
     1521      {
     1522        case 0:
     1523        case 1:
     1524          mpz_clear(&(*a)->n);
     1525        case 3:
     1526          mpz_clear(&(*a)->z);
     1527#ifdef LDEBUG
     1528          (*a)->s=2;
     1529#endif
     1530      }
     1531      omFreeBin((ADDRESS) *a, rnumber_bin);
     1532    }
     1533    *a=NULL;
     1534  }
     1535}
     1536
     1537/*2
     1538* za:= - za
     1539*/
     1540LINLINE number nlNeg (number a)
     1541{
     1542#ifdef LDEBUG
     1543  nlTest(a);
     1544#endif
     1545  if(SR_HDL(a) &SR_INT)
     1546  {
     1547    int r=SR_TO_INT(a);
     1548    if (r==(-(1<<28))) a=nlRInit(1<<28);
     1549    else               a=INT_TO_SR(-r);
     1550  }
     1551  else
     1552  {
     1553    mpz_neg(&a->z,&a->z);
     1554    if ((a->s==3) && (mpz_size1(&a->z)<=MP_SMALL))
     1555    {
     1556      int ui=(int)mpz_get_si(&a->z);
     1557      if ((((ui<<3)>>3)==ui)
     1558      && (mpz_cmp_si(&a->z,(long)ui)==0))
     1559      {
     1560        mpz_clear(&a->z);
     1561        omFreeBin((ADDRESS)a, rnumber_bin);
     1562        a=INT_TO_SR(ui);
     1563      }
     1564    }
     1565  }
     1566#ifdef LDEBUG
     1567  nlTest(a);
     1568#endif
     1569  return a;
     1570}
     1571
    6291572/*2
    6301573* u:= a + b
    6311574*/
    632 number nlAdd (number a, number b)
     1575LINLINE number nlAdd (number a, number b)
    6331576{
    6341577  number u;
    … …  
    8621805* u:= a - b
    8631806*/
    864 number nlSub (number a, number b)
     1807LINLINE number nlSub (number a, number b)
    8651808{
    8661809  number u;
    … …  
    11602103* u := a * b
    11612104*/
    1162 number nlMult (number a, number b)
     2105LINLINE number nlMult (number a, number b)
    11632106{
    11642107  number u;
    … …  
    13192262}
    13202263
    1321 /*2
    1322 * u := a / b in Z, if b | a (else undefined)
    1323 */
    1324 number   nlExactDiv(number a, number b)
    1325 {
    1326   if (b==INT_TO_SR(0))
    1327   {
    1328     WerrorS("div. by 0");
    1329     return INT_TO_SR(0);
    1330   }
    1331   if (a==INT_TO_SR(0))
    1332     return INT_TO_SR(0);
    1333   number u;
    1334   if (SR_HDL(a) & SR_HDL(b) & SR_INT)
    1335   {
    1336     /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
    1337     if ((a==INT_TO_SR(-(1<<28)))&&(b==INT_TO_SR(-1)))
    1338     {
    1339       return nlRInit(1<<28);
    1340     }
    1341     int aa=SR_TO_INT(a);
    1342     int bb=SR_TO_INT(b);
    1343     return INT_TO_SR(aa/bb);
    1344   }
    1345   number bb=NULL;
    1346   if (SR_HDL(b) & SR_INT)
    1347   {
    1348     bb=nlRInit(SR_TO_INT(b));
    1349     b=bb;
    1350   }
    1351   u=(number)omAllocBin(rnumber_bin);
    1352 #if defined(LDEBUG)
    1353   u->debug=123456;
    1354 #endif
    1355   mpz_init(&u->z);
    1356   /* u=a/b */
    1357   u->s = 3;
    1358   MPZ_EXACTDIV(&u->z,&a->z,&b->z);
    1359   if (bb!=NULL)
    1360   {
    1361     mpz_clear(&bb->z);
    1362     omFreeBin((ADDRESS)bb, rnumber_bin);
    1363   }
    1364   if (mpz_size1(&u->z)<=MP_SMALL)
    1365   {
    1366     int ui=(int)mpz_get_si(&u->z);
    1367     if ((((ui<<3)>>3)==ui)
    1368     && (mpz_cmp_si(&u->z,(long)ui)==0))
    1369     {
    1370       mpz_clear(&u->z);
    1371       omFreeBin((ADDRESS)u, rnumber_bin);
    1372       u=INT_TO_SR(ui);
    1373     }
    1374   }
    1375 #ifdef LDEBUG
    1376   nlTest(u);
    1377 #endif
    1378   return u;
    1379 }
    1380 
    1381 /*2
    1382 * u := a / b in Z
    1383 */
    1384 number nlIntDiv (number a, number b)
    1385 {
    1386   if (b==INT_TO_SR(0))
    1387   {
    1388     WerrorS("div. by 0");
    1389     return INT_TO_SR(0);
    1390   }
    1391   if (a==INT_TO_SR(0))
    1392     return INT_TO_SR(0);
    1393   number u;
    1394   if (SR_HDL(a) & SR_HDL(b) & SR_INT)
    1395   {
    1396     /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
    1397     if ((a==INT_TO_SR(-(1<<28)))&&(b==INT_TO_SR(-1)))
    1398     {
    1399       return nlRInit(1<<28);
    1400     }
    1401     /* consider the signs of a and b:
    1402     *  + + -> u=a+b-1
    1403     *  + - -> u=a-b-1
    1404     *  - + -> u=a-b+1
    1405     *  - - -> u=a+b+1
    1406     */
    1407     int aa=SR_TO_INT(a);
    1408     int bb=SR_TO_INT(b);
    1409     if (aa>0)
    1410     {
    1411       if (bb>0)
    1412         return INT_TO_SR((aa+bb-1)/bb);
    1413       else
    1414         return INT_TO_SR((aa-bb-1)/bb);
    1415     }
    1416     else
    1417     {
    1418       if (bb>0)
    1419       {
    1420         return INT_TO_SR((aa-bb+1)/bb);
    1421       }
    1422       else
    1423       {
    1424         return INT_TO_SR((aa+bb+1)/bb);
    1425       }
    1426     }
    1427   }
    1428   if (SR_HDL(a) & SR_INT)
    1429   {
    1430     /* the small int -(1<<28) divided by 2^28 is 1   */
    1431     if (a==INT_TO_SR(-(1<<28)))
    1432     {
    1433       if(mpz_cmp_si(&b->z,(long)(1<<28))==0)
    1434       {
    1435         return INT_TO_SR(-1);
    1436       }
    1437     }
    1438     /* a is a small and b is a large int: -> 0 */
    1439     return INT_TO_SR(0);
    1440   }
    1441   number bb=NULL;
    1442   if (SR_HDL(b) & SR_INT)
    1443   {
    1444     bb=nlRInit(SR_TO_INT(b));
    1445     b=bb;
    1446   }
    1447   u=(number)omAllocBin(rnumber_bin);
    1448 #if defined(LDEBUG)
    1449   u->debug=123456;
    1450 #endif
    1451   mpz_init_set(&u->z,&a->z);
    1452   /* consider the signs of a and b:
    1453   *  + + -> u=a+b-1
    1454   *  + - -> u=a-b-1
    1455   *  - + -> u=a-b+1
    1456   *  - - -> u=a+b+1
    1457   */
    1458   if (mpz_isNeg(&a->z))
    1459   {
    1460     if (mpz_isNeg(&b->z))
    1461     {
    1462       mpz_add(&u->z,&u->z,&b->z);
    1463     }
    1464     else
    1465     {
    1466       mpz_sub(&u->z,&u->z,&b->z);
    1467     }
    1468     mpz_add_ui(&u->z,&u->z,1);
    1469   }
    1470   else
    1471   {
    1472     if (mpz_isNeg(&b->z))
    1473     {
    1474       mpz_sub(&u->z,&u->z,&b->z);
    1475     }
    1476     else
    1477     {
    1478       mpz_add(&u->z,&u->z,&b->z);
    1479     }
    1480     mpz_sub_ui(&u->z,&u->z,1);
    1481   }
    1482   /* u=u/b */
    1483   u->s = 3;
    1484   MPZ_DIV(&u->z,&u->z,&b->z);
    1485   nlGmpSimple(&u->z);
    1486   if (bb!=NULL)
    1487   {
    1488     mpz_clear(&bb->z);
    1489     omFreeBin((ADDRESS)bb, rnumber_bin);
    1490   }
    1491   if (mpz_size1(&u->z)<=MP_SMALL)
    1492   {
    1493     int ui=(int)mpz_get_si(&u->z);
    1494     if ((((ui<<3)>>3)==ui)
    1495     && (mpz_cmp_si(&u->z,(long)ui)==0))
    1496     {
    1497       mpz_clear(&u->z);
    1498       omFreeBin((ADDRESS)u, rnumber_bin);
    1499       u=INT_TO_SR(ui);
    1500     }
    1501   }
    1502 #ifdef LDEBUG
    1503   nlTest(u);
    1504 #endif
    1505   return u;
    1506 }
    1507 
    1508 /*2
    1509 * u := a mod b in Z, u>=0
    1510 */
    1511 number nlIntMod (number a, number b)
    1512 {
    1513   if (b==INT_TO_SR(0))
    1514   {
    1515     WerrorS("div. by 0");
    1516     return INT_TO_SR(0);
    1517   }
    1518   if (a==INT_TO_SR(0))
    1519     return INT_TO_SR(0);
    1520   number u;
    1521   if (SR_HDL(a) & SR_HDL(b) & SR_INT)
    1522   {
    1523     if ((int)a>0)
    1524     {
    1525       if ((int)b>0)
    1526         return INT_TO_SR(SR_TO_INT(a)%SR_TO_INT(b));
    1527       else
    1528         return INT_TO_SR(SR_TO_INT(a)%(-SR_TO_INT(b)));
    1529     }
    1530     else
    1531     {
    1532       if ((int)b>0)
    1533       {
    1534         int i=(-SR_TO_INT(a))%SR_TO_INT(b);
    1535         if ( i != 0 ) i = (SR_TO_INT(b))-i;
    1536         return INT_TO_SR(i);
    1537       }
    1538       else
    1539       {
    1540         int i=(-SR_TO_INT(a))%(-SR_TO_INT(b));
    1541         if ( i != 0 ) i = (-SR_TO_INT(b))-i;
    1542         return INT_TO_SR(i);
    1543       }
    1544     }
    1545   }
    1546   if (SR_HDL(a) & SR_INT)
    1547   {
    1548     /* a is a small and b is a large int: -> a or (a+b) or (a-b) */
    1549     if ((int)a<0)
    1550     {
    1551       if (mpz_isNeg(&b->z))
    1552         return nlSub(a,b);
    1553       /*else*/
    1554         return nlAdd(a,b);
    1555     }
    1556     /*else*/
    1557       return a;
    1558   }
    1559   number bb=NULL;
    1560   if (SR_HDL(b) & SR_INT)
    1561   {
    1562     bb=nlRInit(SR_TO_INT(b));
    1563     b=bb;
    1564   }
    1565   u=(number)omAllocBin(rnumber_bin);
    1566 #if defined(LDEBUG)
    1567   u->debug=123456;
    1568 #endif
    1569   mpz_init(&u->z);
    1570   u->s = 3;
    1571   mpz_mod(&u->z,&a->z,&b->z);
    1572   if (bb!=NULL)
    1573   {
    1574     mpz_clear(&bb->z);
    1575     omFreeBin((ADDRESS)bb, rnumber_bin);
    1576   }
    1577   if (mpz_isNeg(&u->z))
    1578   {
    1579     if (mpz_isNeg(&b->z))
    1580       mpz_sub(&u->z,&u->z,&b->z);
    1581     else
    1582       mpz_add(&u->z,&u->z,&b->z);
    1583   }
    1584   nlGmpSimple(&u->z);
    1585   if (mpz_size1(&u->z)<=MP_SMALL)
    1586   {
    1587     int ui=(int)mpz_get_si(&u->z);
    1588     if ((((ui<<3)>>3)==ui)
    1589     && (mpz_cmp_si(&u->z,(long)ui)==0))
    1590     {
    1591       mpz_clear(&u->z);
    1592       omFreeBin((ADDRESS)u, rnumber_bin);
    1593       u=INT_TO_SR(ui);
    1594     }
    1595   }
    1596 #ifdef LDEBUG
    1597   nlTest(u);
    1598 #endif
    1599   return u;
    1600 }
    1601 
    1602 /*2
    1603 * u := a / b
    1604 */
    1605 number nlDiv (number a, number b)
    1606 {
    1607   number u;
    1608   if (nlIsZero(b))
    1609   {
    1610     WerrorS("div. by 0");
    1611     return INT_TO_SR(0);
    1612   }
    1613   u=(number)omAllocBin(rnumber_bin);
    1614   u->s=0;
    1615 #if defined(LDEBUG)
    1616   u->debug=123456;
    1617 #endif
    1618 // ---------- short / short ------------------------------------
    1619   if (SR_HDL(a) & SR_HDL(b) & SR_INT)
    1620   {
    1621     int i=SR_TO_INT(a);
    1622     int j=SR_TO_INT(b);
    1623     int r=i%j;
    1624     if (r==0)
    1625     {
    1626       omFreeBin((ADDRESS)u, rnumber_bin);
    1627       return INT_TO_SR(i/j);
    1628     }
    1629     mpz_init_set_si(&u->z,(long)i);
    1630     mpz_init_set_si(&u->n,(long)j);
    1631   }
    1632   else
    1633   {
    1634     mpz_init(&u->z);
    1635 // ---------- short / long ------------------------------------
    1636     if (SR_HDL(a) & SR_INT)
    1637     {
    1638       // short a / (z/n) -> (a*n)/z
    1639       if (b->s<2)
    1640       {
    1641         if ((int)a>=0)
    1642           mpz_mul_ui(&u->z,&b->n,SR_TO_INT(a));
    1643         else
    1644         {
    1645           mpz_mul_ui(&u->z,&b->n,-SR_TO_INT(a));
    1646           mpz_neg(&u->z,&u->z);
    1647         }
    1648       }
    1649       else
    1650       // short a / long z -> a/z
    1651       {
    1652         mpz_set_si(&u->z,SR_TO_INT(a));
    1653       }
    1654       if (mpz_cmp(&u->z,&b->z)==0)
    1655       {
    1656         mpz_clear(&u->z);
    1657         omFreeBin((ADDRESS)u, rnumber_bin);
    1658         return INT_TO_SR(1);
    1659       }
    1660       mpz_init_set(&u->n,&b->z);
    1661     }
    1662 // ---------- long / short ------------------------------------
    1663     else if (SR_HDL(b) & SR_INT)
    1664     {
    1665       mpz_set(&u->z,&a->z);
    1666       // (z/n) / b -> z/(n*b)
    1667       if (a->s<2)
    1668       {
    1669         mpz_init_set(&u->n,&a->n);
    1670         if ((int)b>0)
    1671           mpz_mul_ui(&u->n,&u->n,SR_TO_INT(b));
    1672         else
    1673         {
    1674           mpz_mul_ui(&u->n,&u->n,-SR_TO_INT(b));
    1675           mpz_neg(&u->z,&u->z);
    1676         }
    1677       }
    1678       else
    1679       // long z / short b -> z/b
    1680       {
    1681         //mpz_set(&u->z,&a->z);
    1682         mpz_init_set_si(&u->n,SR_TO_INT(b));
    1683       }
    1684     }
    1685 // ---------- long / long ------------------------------------
    1686     else
    1687     {
    1688       //u->s=0;
    1689       mpz_set(&u->z,&a->z);
    1690       mpz_init_set(&u->n,&b->z);
    1691       if (a->s<2) mpz_mul(&u->n,&u->n,&a->n);
    1692       if (b->s<2) mpz_mul(&u->z,&u->z,&b->n);
    1693     }
    1694   }
    1695   if (mpz_isNeg(&u->n))
    1696   {
    1697     mpz_neg(&u->z,&u->z);
    1698     mpz_neg(&u->n,&u->n);
    1699   }
    1700   if (mpz_cmp_si(&u->n,(long)1)==0)
    1701   {
    1702     mpz_clear(&u->n);
    1703     if (mpz_size1(&u->z)<=MP_SMALL)
    1704     {
    1705       int ui=(int)mpz_get_si(&u->z);
    1706       if ((((ui<<3)>>3)==ui)
    1707       && (mpz_cmp_si(&u->z,(long)ui)==0))
    1708       {
    1709         mpz_clear(&u->z);
    1710         omFreeBin((ADDRESS)u, rnumber_bin);
    1711         return INT_TO_SR(ui);
    1712       }
    1713     }
    1714     u->s=3;
    1715   }
    1716 #ifdef LDEBUG
    1717   nlTest(u);
    1718 #endif
    1719   return u;
    1720 }
    1721 
    1722 #if (defined(i386)) && (defined(HAVE_LIBGMP1))
    1723 /*
    1724 * compute x^exp for x in Z
    1725 * there seems to be an bug in mpz_pow_ui for i386
    1726 */
    1727 static inline void nlPow (MP_INT * res,MP_INT * x,int exp)
    1728 {
    1729   if (exp==0)
    1730   {
    1731     mpz_set_si(res,(long)1);
    1732   }
    1733   else
    1734   {
    1735     MP_INT xh;
    1736     mpz_init(&xh);
    1737     mpz_set(res,x);
    1738     exp--;
    1739     while (exp!=0)
    1740     {
    1741       mpz_mul(&xh,res,x);
    1742       mpz_set(res,&xh);
    1743       exp--;
    1744     }
    1745     mpz_clear(&xh);
    1746   }
    1747 }
    1748 #endif
    1749 
    1750 /*2
    1751 * u:= x ^ exp
    1752 */
    1753 void nlPower (number x,int exp,number * u)
    1754 {
    1755   if (!nlIsZero(x))
    1756   {
    1757 #ifdef LDEBUG
    1758     nlTest(x);
    1759 #endif
    1760     number aa=NULL;
    1761     if (SR_HDL(x) & SR_INT)
    1762     {
    1763       aa=nlRInit(SR_TO_INT(x));
    1764       x=aa;
    1765     }
    1766     *u=(number)omAllocBin(rnumber_bin);
    1767 #if defined(LDEBUG)
    1768     (*u)->debug=123456;
    1769 #endif
    1770     mpz_init(&(*u)->z);
    1771     (*u)->s = x->s;
    1772 #if (!defined(i386)) || (defined(HAVE_LIBGMP2))
    1773     mpz_pow_ui(&(*u)->z,&x->z,(unsigned long)exp);
    1774 #else
    1775     nlPow(&(*u)->z,&x->z,exp);
    1776 #endif
    1777     if (x->s<2)
    1778     {
    1779       mpz_init(&(*u)->n);
    1780 #if (!defined(i386)) || (defined(HAVE_LIBGMP2))
    1781       mpz_pow_ui(&(*u)->n,&x->n,(unsigned long)exp);
    1782 #else
    1783       nlPow(&(*u)->n,&x->n,exp);
    1784 #endif
    1785     }
    1786     if (aa!=NULL)
    1787     {
    1788       mpz_clear(&aa->z);
    1789       omFreeBin((ADDRESS)aa, rnumber_bin);
    1790     }
    1791     if (((*u)->s<=2) && (mpz_cmp_si(&(*u)->n,(long)1)==0))
    1792     {
    1793       mpz_clear(&(*u)->n);
    1794       (*u)->s=3;
    1795     }
    1796     if (((*u)->s==3) && (mpz_size1(&(*u)->z)<=MP_SMALL))
    1797     {
    1798       int ui=(int)mpz_get_si(&(*u)->z);
    1799       if ((((ui<<3)>>3)==ui)
    1800       && (mpz_cmp_si(&(*u)->z,(long)ui)==0))
    1801       {
    1802         mpz_clear(&(*u)->z);
    1803         omFreeBin((ADDRESS)*u, rnumber_bin);
    1804         *u=INT_TO_SR(ui);
    1805       }
    1806     }
    1807   }
    1808   else
    1809     *u = INT_TO_SR(0);
    1810 #ifdef LDEBUG
    1811   nlTest(*u);
    1812 #endif
    1813 }
    1814 
    1815 BOOLEAN nlIsZero (number a)
    1816 {
    1817 // #ifdef LDEBUG
    1818 //   nlTest(a);
    1819 // #endif
    1820   //if (a==INT_TO_SR(0)) return TRUE;
    1821   //if (SR_HDL(a) & SR_INT) return FALSE;
    1822   //number aa=nlCopy(a);
    1823   //nlNormalize(aa);
    1824   return (a==INT_TO_SR(0));
    1825 }
    1826 
    1827 /*2
    1828 * za >= 0 ?
    1829 */
    1830 BOOLEAN nlGreaterZero (number a)
    1831 {
    1832 #ifdef LDEBUG
    1833   nlTest(a);
    1834 #endif
    1835   if (SR_HDL(a) & SR_INT) return SR_HDL(a)>=0;
    1836   return (!mpz_isNeg(&a->z));
    1837 }
    1838 
    1839 /*2
    1840 * a > b ?
    1841 */
    1842 BOOLEAN nlGreater (number a, number b)
    1843 {
    1844 #ifdef LDEBUG
    1845   nlTest(a);
    1846   nlTest(b);
    1847 #endif
    1848   number r;
    1849   BOOLEAN rr;
    1850   r=nlSub(a,b);
    1851   rr=(!nlIsZero(r)) && (nlGreaterZero(r));
    1852   nlDelete(&r);
    1853   return rr;
    1854 }
    1855 
    1856 /*2
    1857 * a = b ?
    1858 */
    1859 BOOLEAN nlEqual (number a, number b)
    1860 {
    1861 #ifdef LDEBUG
    1862   nlTest(a);
    1863   nlTest(b);
    1864 #endif
    1865 // short - short
    1866   if (SR_HDL(a) & SR_HDL(b) & SR_INT) return a==b;
    1867 //  long - short
    1868   BOOLEAN bo;
    1869   if (SR_HDL(b) & SR_INT)
    1870   {
    1871     if (a->s!=0) return FALSE;
    1872     number n=b; b=a; a=n;
    1873   }
    1874 //  short - long
    1875   if (SR_HDL(a) & SR_INT)
    1876   {
    1877     if (b->s!=0)
    1878       return FALSE;
    1879     if (((int)a > 0) && (mpz_isNeg(&b->z)))
    1880       return FALSE;
    1881     if (((int)a < 0) && (!mpz_isNeg(&b->z)))
    1882       return FALSE;
    1883     MP_INT  bb;
    1884     mpz_init_set(&bb,&b->n);
    1885     if ((int)a<0)
    1886     {
    1887       mpz_neg(&bb,&bb);
    1888       mpz_mul_ui(&bb,&bb,(long)-SR_TO_INT(a));
    1889     }
    1890     else
    1891     {
    1892       mpz_mul_ui(&bb,&bb,(long)SR_TO_INT(a));
    1893     }
    1894     bo=(mpz_cmp(&bb,&b->z)==0);
    1895     mpz_clear(&bb);
    1896     return bo;
    1897   }
    1898 // long - long
    1899   if (((a->s==1) && (b->s==3))
    1900   ||  ((b->s==1) && (a->s==3)))
    1901     return FALSE;
    1902   if (mpz_isNeg(&a->z)&&(!mpz_isNeg(&b->z)))
    1903     return FALSE;
    1904   if (mpz_isNeg(&b->z)&&(!mpz_isNeg(&a->z)))
    1905     return FALSE;
    1906   nlGmpSimple(&a->z);
    1907   nlGmpSimple(&b->z);
    1908   if (a->s<2)
    1909     nlGmpSimple(&a->n);
    1910   if (b->s<2)
    1911     nlGmpSimple(&b->n);
    1912   MP_INT  aa;
    1913   MP_INT  bb;
    1914   mpz_init_set(&aa,&a->z);
    1915   mpz_init_set(&bb,&b->z);
    1916   if (a->s<2) mpz_mul(&bb,&bb,&a->n);
    1917   if (b->s<2) mpz_mul(&aa,&aa,&b->n);
    1918   bo=(mpz_cmp(&aa,&bb)==0);
    1919   mpz_clear(&aa);
    1920   mpz_clear(&bb);
    1921   return bo;
    1922 }
    1923 
    1924 /*2
    1925 * a == 1 ?
    1926 */
    1927 BOOLEAN nlIsOne (number a)
    1928 {
    1929 #ifdef LDEBUG
    1930   if (a==NULL) return FALSE;
    1931   nlTest(a);
    1932 #endif
    1933   if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(1));
    1934   return FALSE;
    1935 }
    1936 
    1937 /*2
    1938 * a == -1 ?
    1939 */
    1940 BOOLEAN nlIsMOne (number a)
    1941 {
    1942 #ifdef LDEBUG
    1943   if (a==NULL) return FALSE;
    1944   nlTest(a);
    1945 #endif
    1946   if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(-1));
    1947   return FALSE;
    1948 }
    1949 
    1950 /*2
    1951 * result =gcd(a,b)
    1952 */
    1953 number nlGcd(number a, number b)
    1954 {
    1955   number result;
    1956 #ifdef LDEBUG
    1957   nlTest(a);
    1958   nlTest(b);
    1959 #endif
    1960   //nlNormalize(a);
    1961   //nlNormalize(b);
    1962   if ((SR_HDL(a)==5)||(a==INT_TO_SR(-1))
    1963   ||  (SR_HDL(b)==5)||(b==INT_TO_SR(-1))) return INT_TO_SR(1);
    1964   if (SR_HDL(a) & SR_HDL(b) & SR_INT)
    1965   {
    1966     int i=SR_TO_INT(a);
    1967     int j=SR_TO_INT(b);
    1968     if((i==0)||(j==0))
    1969       return INT_TO_SR(1);
    1970     int l;
    1971     i=ABS(i);
    1972     j=ABS(j);
    1973     do
    1974     {
    1975       l=i%j;
    1976       i=j;
    1977       j=l;
    1978     } while (l!=0);
    1979     return INT_TO_SR(i);
    1980   }
    1981   if (((!(SR_HDL(a) & SR_INT))&&(a->s<2))
    1982   ||  ((!(SR_HDL(b) & SR_INT))&&(b->s<2))) return INT_TO_SR(1);
    1983   number aa=NULL;
    1984   if (SR_HDL(a) & SR_INT)
    1985   {
    1986     aa=nlRInit(SR_TO_INT(a));
    1987     a=aa;
    1988   }
    1989   else
    1990   if (SR_HDL(b) & SR_INT)
    1991   {
    1992     aa=nlRInit(SR_TO_INT(b));
    1993     b=aa;
    1994   }
    1995   result=(number)omAllocBin(rnumber_bin);
    1996 #if defined(LDEBUG)
    1997   result->debug=123456;
    1998 #endif
    1999   mpz_init(&result->z);
    2000   mpz_gcd(&result->z,&a->z,&b->z);
    2001   nlGmpSimple(&result->z);
    2002   result->s = 3;
    2003   if (mpz_size1(&result->z)<=MP_SMALL)
    2004   {
    2005     int ui=(int)mpz_get_si(&result->z);
    2006     if ((((ui<<3)>>3)==ui)
    2007     && (mpz_cmp_si(&result->z,(long)ui)==0))
    2008     {
    2009       mpz_clear(&result->z);
    2010       omFreeBin((ADDRESS)result, rnumber_bin);
    2011       result=INT_TO_SR(ui);
    2012     }
    2013   }
    2014   if (aa!=NULL)
    2015   {
    2016     mpz_clear(&aa->z);
    2017     omFreeBin((ADDRESS)aa, rnumber_bin);
    2018   }
    2019 #ifdef LDEBUG
    2020   nlTest(result);
    2021 #endif
    2022   return result;
    2023 }
    2024 
    2025 /*2
    2026 * simplify x
    2027 */
    2028 void nlNormalize (number &x)
    2029 {
    2030   if ((SR_HDL(x) & SR_INT) ||(x==NULL))
    2031     return;
    2032 #ifdef LDEBUG
    2033   if (!nlTest(x)) return;
    2034 #endif
    2035   if (x->s==3)
    2036   {
    2037     if (mpz_cmp_ui(&x->z,(long)0)==0)
    2038     {
    2039       nlDelete(&x);
    2040       x=nlInit(0);
    2041       return;
    2042     }
    2043     if (mpz_size1(&x->z)<=MP_SMALL)
    2044     {
    2045       int ui=(int)mpz_get_si(&x->z);
    2046       if ((((ui<<3)>>3)==ui)
    2047       && (mpz_cmp_si(&x->z,(long)ui)==0))
    2048       {
    2049         mpz_clear(&x->z);
    2050         omFreeBin((ADDRESS)x, rnumber_bin);
    2051         x=INT_TO_SR(ui);
    2052         return;
    2053       }
    2054     }
    2055   }
    2056   if (x->s!=1) nlGmpSimple(&x->z);
    2057   if (x->s==0)
    2058   {
    2059     BOOLEAN divided=FALSE;
    2060     MP_INT gcd;
    2061     mpz_init(&gcd);
    2062     mpz_gcd(&gcd,&x->z,&x->n);
    2063     x->s=1;
    2064     if (mpz_cmp_si(&gcd,(long)1)!=0)
    2065     {
    2066       MP_INT r;
    2067       mpz_init(&r);
    2068       MPZ_EXACTDIV(&r,&x->z,&gcd);
    2069       mpz_set(&x->z,&r);
    2070       MPZ_EXACTDIV(&r,&x->n,&gcd);
    2071       mpz_set(&x->n,&r);
    2072       mpz_clear(&r);
    2073       nlGmpSimple(&x->z);
    2074       divided=TRUE;
    2075     }
    2076     mpz_clear(&gcd);
    2077     nlGmpSimple(&x->n);
    2078     if (mpz_cmp_si(&x->n,(long)1)==0)
    2079     {
    2080       mpz_clear(&x->n);
    2081       if (mpz_size1(&x->z)<=MP_SMALL)
    2082       {
    2083         int ui=(int)mpz_get_si(&x->z);
    2084         if ((((ui<<3)>>3)==ui)
    2085         && (mpz_cmp_si(&x->z,(long)ui)==0))
    2086         {
    2087           mpz_clear(&x->z);
    2088 #if defined(LDEBUG)
    2089           x->debug=654324;
    2090 #endif
    2091           omFreeBin((ADDRESS)x, rnumber_bin);
    2092           x=INT_TO_SR(ui);
    2093           return;
    2094         }
    2095       }
    2096       x->s=3;
    2097     }
    2098     else if (divided)
    2099     {
    2100       _mpz_realloc(&x->n,mpz_size1(&x->n));
    2101     }
    2102     if (divided) _mpz_realloc(&x->z,mpz_size1(&x->z));
    2103   }
    2104 #ifdef LDEBUG
    2105   nlTest(x);
    2106 #endif
    2107 }
    2108 
    2109 /*2
    2110 * returns in result->z the lcm(a->z,b->n)
    2111 */
    2112 number nlLcm(number a, number b)
    2113 {
    2114   number result;
    2115 #ifdef LDEBUG
    2116   nlTest(a);
    2117   nlTest(b);
    2118 #endif
    2119   if ((SR_HDL(b) & SR_INT)
    2120   || (b->s==3))
    2121   {
    2122     // b is 1/(b->n) => b->n is 1 => result is a
    2123     return nlCopy(a);
    2124   }
    2125   number aa=NULL;
    2126   if (SR_HDL(a) & SR_INT)
    2127   {
    2128     aa=nlRInit(SR_TO_INT(a));
    2129     a=aa;
    2130   }
    2131   result=(number)omAllocBin(rnumber_bin);
    2132 #if defined(LDEBUG)
    2133   result->debug=123456;
    2134 #endif
    2135   result->s=3;
    2136   MP_INT gcd;
    2137   mpz_init(&gcd);
    2138   mpz_init(&result->z);
    2139   mpz_gcd(&gcd,&a->z,&b->n);
    2140   if (mpz_cmp_si(&gcd,(long)1)!=0)
    2141   {
    2142     MP_INT bt;
    2143     mpz_init_set(&bt,&b->n);
    2144     MPZ_EXACTDIV(&bt,&bt,&gcd);
    2145     mpz_mul(&result->z,&a->z,&bt);
    2146     mpz_clear(&bt);
    2147   }
    2148   else
    2149     mpz_mul(&result->z,&a->z,&b->n);
    2150   mpz_clear(&gcd);
    2151   if (aa!=NULL)
    2152   {
    2153     mpz_clear(&aa->z);
    2154     omFreeBin((ADDRESS)aa, rnumber_bin);
    2155   }
    2156   nlGmpSimple(&result->z);
    2157   if (mpz_size1(&result->z)<=MP_SMALL)
    2158   {
    2159     int ui=(int)mpz_get_si(&result->z);
    2160     if ((((ui<<3)>>3)==ui)
    2161     && (mpz_cmp_si(&result->z,(long)ui)==0))
    2162     {
    2163       mpz_clear(&result->z);
    2164       omFreeBin((ADDRESS)result, rnumber_bin);
    2165       return INT_TO_SR(ui);
    2166     }
    2167   }
    2168 #ifdef LDEBUG
    2169   nlTest(result);
    2170 #endif
    2171   return result;
    2172 }
    2173 
    2174 int nlModP(number n, int p)
    2175 {
    2176   if (SR_HDL(n) & SR_INT)
    2177   {
    2178     int i=SR_TO_INT(n);
    2179     if (i<0) return (p-((-i)%p));
    2180     return i%p;
    2181   }
    2182   int iz=mpz_mmod_ui(NULL,&n->z,(unsigned long)p);
    2183   if (n->s!=3)
    2184   {
    2185     int in=mpz_mmod_ui(NULL,&n->n,(unsigned long)p);
    2186     return (int)npDiv((number)iz,(number)in);
    2187   }
    2188   return iz;
    2189 }
    2190 
    2191 /*2
    2192 * acces to denominator, other 1 for integers
    2193 */
    2194 number   nlGetDenom(number &n)
    2195 {
    2196   if (!(SR_HDL(n) & SR_INT))
    2197   {
    2198     if (n->s==0)
    2199     {
    2200       nlNormalize(n);
    2201     }
    2202     if (!(SR_HDL(n) & SR_INT))
    2203     {
    2204       if (n->s!=3)
    2205       {
    2206         number u=(number)omAllocBin(rnumber_bin);
    2207         u->s=3;
    2208 #if defined(LDEBUG)
    2209         u->debug=123456;
    2210 #endif
    2211 
    2212         mpz_init_set(&u->z,&n->n);
    2213         {
    2214           int ui=(int)mpz_get_si(&u->z);
    2215           if ((((ui<<3)>>3)==ui)
    2216           && (mpz_cmp_si(&u->z,(long)ui)==0))
    2217           {
    2218             mpz_clear(&u->z);
    2219             omFreeBin((ADDRESS)u, rnumber_bin);
    2220             return INT_TO_SR(ui);
    2221           }
    2222         }
    2223         return u;
    2224       }
    2225     }
    2226   }
    2227   return INT_TO_SR(1);
    2228 }
     2264#endif // DO_LINLINE
     2265
     2266#endif // LONGRAT_CC

  • Singular/longrat.h

    r44e09a r270e65  
    44*  Computer Algebra System SINGULAR     *
    55****************************************/
    6 /* $Id: longrat.h,v 1.17 2000-08-14 12:56:36 obachman Exp $ */
     6/* $Id: longrat.h,v 1.18 2000-09-20 13:25:41 obachman Exp $ */
    77/*
    88* ABSTRACT: computation with long rational numbers
    … …  
    5252};
    5353
     54// allow inlining only from p_Numbers.h and if ! LDEBUG
     55
     56#if defined(DO_LINLINE) && defined(P_NUMBERS_H) && !defined(LDEBUG)
     57#define LINLINE static inline
     58#else
     59#define LINLINE
     60#undef DO_LINLINE
     61#endif // DO_LINLINE
     62
     63LINLINE BOOLEAN  nlEqual(number a, number b);
     64LINLINE number   nlInit(int i);
     65LINLINE BOOLEAN  nlIsOne(number a);
     66LINLINE BOOLEAN  nlIsZero(number za);
     67LINLINE number   nlCopy(number a);
     68LINLINE void     nlNew(number *r);
     69#ifndef LDEBUG
     70LINLINE void     nlDelete(number *a);
     71#endif
     72LINLINE number   nlNeg(number za);
     73LINLINE number   nlAdd(number la, number li);
     74LINLINE number   nlSub(number la, number li);
     75LINLINE number   nlMult(number a, number b);
     76
    5477number   nlGcd(number a, number b);
    5578number   nlLcm(number a, number b);   /*special routine !*/
    5679BOOLEAN  nlGreater(number a, number b);
    57 BOOLEAN  nlEqual(number a, number b);
    58 BOOLEAN  nlIsOne(number a);
    5980BOOLEAN  nlIsMOne(number a);
    60 void     nlNew(number *r);
    61 number   nlInit(int i);
    6281number   nlInit(number i);
    6382int      nlInt(number &n);
    64 BOOLEAN  nlIsZero(number za);
    6583BOOLEAN  nlGreaterZero(number za);
    66 number   nlNeg(number za);
    6784number   nlInvers(number a);
    6885void     nlNormalize(number &x);
    69 number   nlAdd(number la, number li);
    70 number   nlSub(number la, number li);
    71 number   nlMult(number a, number b);
    7286number   nlDiv(number a, number b);
    7387number   nlExactDiv(number a, number b);
    … …  
    7589number   nlIntMod(number a, number b);
    7690void     nlPower(number x, int exp, number *lu);
    77 number   nlCopy(number a);
    7891char *   nlRead (char *s, number *a);
    7992void     nlWrite(number &a);
    … …  
    8598void     nlDBDelete(number *a, char *f, int l);
    8699#define  nlDelete(A) nlDBDelete(A,__FILE__,__LINE__)
    87 #else
    88 void     nlDelete(number *a);
    89100#endif
    90101

  • Singular/p_Numbers.h

    r44e09a r270e65  
    77 *  Author:  obachman (Olaf Bachmann)
    88 *  Created: 8/00
    9  *  Version: $Id: p_Numbers.h,v 1.3 2000-09-18 09:19:27 obachman Exp $
     9 *  Version: $Id: p_Numbers.h,v 1.4 2000-09-20 13:25:41 obachman Exp $
    1010 *******************************************************************/
    1111#ifndef P_NUMBERS_H
    … …  
    3232#define n_Sub_FieldZp(n1, n2, r)    npSubM(n1, n2)
    3333
     34#define DO_LINLINE
     35#include "longrat.cc"
     36#define n_Copy_FieldQ(n, r)        nlCopy(n)
     37#define n_Delete_FieldQ(n, r)      nlDelete(n)
     38#define n_Mult_FieldQ(n1, n2, r)   nlMult(n1,n2)
     39#define n_Add_FieldQ(n1, n2, r)    nlAdd(n1, n2)
     40#define n_IsZero_FieldQ(n, r)      nlIsZero(n)
     41#define n_Equal_FieldQ(n1, n2, r)  nlEqual(n1, n2)
     42#define n_Neg_FieldQ(n, r)             nlNeg(n)
     43#define n_Sub_FieldQ(n1, n2, r)    nlSub(n1, n2)
    3444#endif
    3545

  • Singular/p_Procs.cc

    r44e09a r270e65  
    77 *  Author:  obachman (Olaf Bachmann)
    88 *  Created: 8/00
    9  *  Version: $Id: p_Procs.cc,v 1.9 2000-09-18 10:26:13 obachman Exp $
     9 *  Version: $Id: p_Procs.cc,v 1.10 2000-09-20 13:25:41 obachman Exp $
    1010 *******************************************************************/
    1111#include <string.h>
     12
     13#define PDEBUG 2
    1214#include "mod2.h"
    1315
    … …  
    1820 *
    1921 *******************************************************************/
    20 // define to enable/disable ptest in p_Procs
    2122
    2223/***************************************************************
    … …  
    4243//   1 -- plus FieldZp_Length*_OrdGeneral procs
    4344//   2 -- plus FieldZp_Length*_Ord* procs
    44 //   3 -- plus Field*_Length*_OrdGeneral procs
    45 //   4 -- all Field*_Length*_Ord* procs
     45//   3 -- plus FieldQ_Length*_Ord*
     46//   4 -- plus FieldGeneral_Length*_OrdGeneral procs
     47//   5 -- all Field*_Length*_Ord* procs
    4648#ifdef NDEBUG
    47 const int HAVE_FAST_P_PROCS = 4;
     49const int HAVE_FAST_P_PROCS = 5;
    4850#else
    4951const int HAVE_FAST_P_PROCS = 0;
    … …  
    5355//   0 -- only FieldGeneral
    5456//   1 -- special cases for FieldZp
     57//   2 -- plus special cases for FieldQ
    5558//   nothing else is implemented, yet
    5659const int HAVE_FAST_FIELD = 1;
    … …  
    101104  FieldR,
    102105  FieldGF,
     106  FieldQ,
    103107#if HAVE_MORE_FIELDS_IMPLEMENTED
    104   FieldQ,
    105108  FieldLong_R,
    106109  FieldLong_C,
    … …  
    192195      case FieldR: return "FieldR";
    193196      case FieldGF: return "FieldGF";
     197      case FieldQ: return "FieldQ";
    194198#if HAVE_MORE_FIELDS_IMPLEMENTED
    195       case FieldQ: return "FieldQ";
    196199      case FieldLong_R: return "FieldLong_R";
    197200      case FieldLong_C: return "FieldLong_C";
    … …  
    331334static inline void FastP_ProcsFilter(p_Field &field, p_Length &length, p_Ord &ord)
    332335{
    333   if (HAVE_FAST_P_PROCS >= 4) return;
     336  if (HAVE_FAST_P_PROCS >= 5) return;
     337 
     338  if (HAVE_FAST_P_PROCS < 3 && field == FieldQ)
     339    field = FieldGeneral;
    334340 
    335341  if ((HAVE_FAST_P_PROCS == 0) ||
    336       (HAVE_FAST_P_PROCS <= 3 && field != FieldZp))
     342      (HAVE_FAST_P_PROCS <= 4 && field != FieldZp && field != FieldQ))
    337343  {
    338344    field = FieldGeneral;
    … …  
    341347    return;
    342348  }
    343   if (HAVE_FAST_P_PROCS == 1 || (HAVE_FAST_P_PROCS == 3 && field != FieldZp))
     349  if (HAVE_FAST_P_PROCS == 1 || (HAVE_FAST_P_PROCS == 4 && field != FieldZp))
    344350    ord = OrdGeneral;
    345351}
    … …  
    347353static inline void FastFieldFilter(p_Field &field)
    348354{
    349   if (HAVE_FAST_FIELD == 0 || field != FieldZp)
     355  if (HAVE_FAST_FIELD <= 0 ||
     356      (HAVE_FAST_FIELD == 1 && field != FieldZp) ||
     357      (field != FieldZp && field != FieldQ))
    350358    field = FieldGeneral;
    351359}
    … …  
    520528  if (rField_is_R(r)) return FieldR;
    521529  if (rField_is_GF(r)) return FieldGF;
     530  if (rField_is_Q(r)) return FieldQ;
    522531#ifdef HAVE_MORE_FIELDS_IMPLEMENTED
    523   if (rField_is_Q(r)) return FieldQ;
    524532  if (rField_is_long_R(r)) return FieldLong_R;
    525533  if (rField_is_long_C(r)) return FieldLong_C;

  • Singular/pp_Mult_mm__Template.cc

    r44e09a r270e65  
    77 *  Author:  obachman (Olaf Bachmann)
    88 *  Created: 8/00
    9  *  Version: $Id: pp_Mult_mm__Template.cc,v 1.3 2000-09-18 09:19:31 obachman Exp $
     9 *  Version: $Id: pp_Mult_mm__Template.cc,v 1.4 2000-09-20 13:25:42 obachman Exp $
    1010 *******************************************************************/
    1111
    … …  
    2121{
    2222  p_Test(p, ri);
    23   p_LmTest(p, ri);
     23  p_LmTest(m, ri);
    2424  if (p == NULL) return NULL;
    2525  spolyrec rp;
    … …  
    2929  DECLARE_LENGTH(const unsigned long length = ri->ExpLSize);
    3030  const unsigned long* m_e = m->exp;
    31   pAssume(!n_IsZero(ln,r));
     31  pAssume(!n_IsZero(ln,ri));
     32  pAssume1(p_GetComp(m, ri) == 0 || p_MaxComp(p, ri) == 0);
     33
    3234
    3335  if (spNoether == NULL)
Note: See TracChangeset for help on using the changeset viewer.