Changeset 2c7f28 in git


Ignore:
Timestamp:
May 27, 2011, 4:25:28 PM (13 years ago)
Author:
Frank Seelisch <seelisch@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
600ea94eb3125257cff2947d86f88bde3f272148
Parents:
9bb545791780a99344a26d6177decb3870f6da87
git-author:
Frank Seelisch <seelisch@mathematik.uni-kl.de>2011-05-27 16:25:28+02:00
git-committer:
Mohamed Barakat <mohamed.barakat@rwth-aachen.de>2011-11-09 12:39:06+01:00
Message:
impl. of rat. funct. fields (no tests yet, no cancellation of gcd's yet)
Location:
libpolys/polys
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • libpolys/polys/Makefile.am

    r9bb5457 r2c7f28  
    4949        kbuckets.cc sbuckets.cc weight.cc weight0.c simpleideals.cc matpol.cc \
    5050        ${USE_P_PROCS_STATIC_CC} ${USE_P_PROCS_DYNAMIC_CC} templates/mod_raw.cc \
    51         ext_fields/algext.cc clapsing.cc clapconv.cc
     51        ext_fields/algext.cc ext_fields/transext.cc clapsing.cc clapconv.cc
    5252
    5353BUILT_SOURCES = templates/p_Procs.inc
  • libpolys/polys/ext_fields/algext.cc

    r9bb5457 r2c7f28  
    8888  if (a == NULL) return TRUE;
    8989  p_Test((poly)a, naRing);
    90   assume(p_Deg((poly)a, naRing) <= p_Deg(naMinpoly, naRing));
     90  assume(p_Totaldegree((poly)a, naRing) <= p_Totaldegree(naMinpoly, naRing));
    9191  return TRUE;
    9292}
     
    128128}
    129129
    130 BOOLEAN naEqual (number a, number b, const coeffs cf)
     130BOOLEAN naEqual(number a, number b, const coeffs cf)
    131131{
    132132  naTest(a); naTest(b);
     
    139139  /// deg test
    140140  int aDeg = 0;
    141   if (a != NULL) aDeg = p_Deg((poly)a, naRing);
     141  if (a != NULL) aDeg = p_Totaldegree((poly)a, naRing);
    142142  int bDeg = 0;
    143   if (b != NULL) bDeg = p_Deg((poly)b, naRing);
     143  if (b != NULL) bDeg = p_Totaldegree((poly)b, naRing);
    144144  if (aDeg != bDeg) return FALSE;
    145145 
     
    172172{
    173173  naTest(a);
    174   if (p_GetExp((poly)a, 1, naRing) != 0) return FALSE;
    175   return n_IsOne(p_GetCoeff((poly)a, naRing), naCoeffs);
     174  poly aAsPoly = (poly)a;
     175  if (!p_IsConstant(aAsPoly, naRing)) return FALSE;
     176  return n_IsOne(p_GetCoeff(aAsPoly, naRing), naCoeffs);
    176177}
    177178
     
    179180{
    180181  naTest(a);
    181   if (p_GetExp((poly)a, 1, naRing) != 0) return FALSE;
    182   return n_IsMOne(p_GetCoeff((poly)a, naRing), naCoeffs);
     182  poly aAsPoly = (poly)a;
     183  if (!p_IsConstant(aAsPoly, naRing)) return FALSE;
     184  return n_IsMOne(p_GetCoeff(aAsPoly, naRing), naCoeffs);
    183185}
    184186
     
    206208{
    207209  naTest(a);
    208   if (p_GetExp((poly)a, 1, naRing) != 0) return 0;
    209   return n_Int(p_GetCoeff((poly)a, naRing), naCoeffs);
     210  poly aAsPoly = (poly)a;
     211  if (!p_IsConstant(aAsPoly, naRing)) return 0;
     212  return n_Int(p_GetCoeff(aAsPoly, naRing), naCoeffs);
    210213}
    211214
     
    213216BOOLEAN naGreater(number a, number b, const coeffs cf)
    214217{
     218  naTest(a); naTest(b);
    215219  if (naIsZero(a, cf)) return FALSE;
    216220  if (naIsZero(b, cf)) return TRUE;
    217221  int aDeg = 0;
    218   if (a != NULL) aDeg = p_Deg((poly)a, naRing);
     222  if (a != NULL) aDeg = p_Totaldegree((poly)a, naRing);
    219223  int bDeg = 0;
    220   if (b != NULL) bDeg = p_Deg((poly)b, naRing);
     224  if (b != NULL) bDeg = p_Totaldegree((poly)b, naRing);
    221225  return (aDeg > bDeg);
    222226}
     
    228232  if (a == NULL)                                            return FALSE;
    229233  if (n_GreaterZero(p_GetCoeff((poly)a, naRing), naCoeffs)) return TRUE;
    230   if (p_Deg((poly)a, naRing) > 0)                           return TRUE;
     234  if (p_Totaldegree((poly)a, naRing) > 0)                   return TRUE;
    231235  return FALSE;
    232236}
     
    298302   for |exp| >= 8 compute power along binary presentation of |exp|, e.g.
    299303      p^13 = p^1 * p^4 * p^8, where we utilise that
    300       p^(2^(k+1)) = p^(2^k) * p^(2^k)
     304      p^(2^(k+1)) = p^(2^k) * p^(2^k);
    301305   intermediate reduction modulo the minimal polynomial is controlled by
    302306   the in-place method heuristicReduce(poly, poly, coeffs); see there.
     
    318322  int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
    319323 
    320   /* now compute 'a' to the 'expAbs'-th power */
     324  /* now compute a^expAbs */
    321325  poly pow; poly aAsPoly = (poly)a;
    322326  if (expAbs <= 7)
     
    363367}
    364368
    365 /* may reduce p module the reducer by calling definiteReduce;
     369/* may reduce p modulo the reducer by calling definiteReduce;
    366370   the decision is made based on the following heuristic
    367371   (which should also only be changed here in this method):
     
    374378  p_Test((poly)reducer, naRing);
    375379  #endif
    376   if (p_Deg(p, naRing) > 10 * p_Deg(reducer, naRing))
     380  if (p_Totaldegree(p, naRing) > 10 * p_Totaldegree(reducer, naRing))
    377381    definiteReduce(p, reducer, cf);
    378382}
     
    424428  AlgExtInfo *e = (AlgExtInfo *)param;
    425429  /* for extension coefficient fields we expect the underlying
    426      polynomials rings to be IDENTICAL, i.e. the SAME OBJECT;
     430     polynomial rings to be IDENTICAL, i.e. the SAME OBJECT;
    427431     this expectation is based on the assumption that we have properly
    428432     registered cf and perform reference counting rather than creating
     
    445449  {
    446450    noOfTerms++;
    447     int d = 0;
    448     for (int i = 1; i <= rVar(naRing); i++)
    449       d += p_GetExp(aAsPoly, i, naRing);
     451    int d = p_GetExp(aAsPoly, 1, naRing);
    450452    if (d > theDegree) theDegree = d;
    451453    pIter(aAsPoly);
     
    496498number naMap00(number a, const coeffs src, const coeffs dst)
    497499{
     500  if (n_IsZero(a, src)) return NULL;
    498501  assume(src == dst->extRing->cf);
    499502  poly result = p_One(dst->extRing);
     
    505508number naMapP0(number a, const coeffs src, const coeffs dst)
    506509{
     510  if (n_IsZero(a, src)) return NULL;
    507511  /* mapping via intermediate int: */
    508512  int n = n_Int(a, src);
     
    514518
    515519/* assumes that either src = Q(a), dst = Q(a), or
    516                        src = Zp(a), dst = Zp(a)     */
     520                       src = Z/p(a), dst = Z/p(a)     */
    517521number naCopyMap(number a, const coeffs src, const coeffs dst)
    518522{
     
    523527number naMap0P(number a, const coeffs src, const coeffs dst)
    524528{
     529  if (n_IsZero(a, src)) return NULL;
    525530  int p = rChar(dst->extRing);
    526531  int n = nlModP(a, p, src);
     
    534539number naMapPP(number a, const coeffs src, const coeffs dst)
    535540{
     541  if (n_IsZero(a, src)) return NULL;
    536542  assume(src == dst->extRing->cf);
    537543  poly result = p_One(dst->extRing);
     
    543549number naMapUP(number a, const coeffs src, const coeffs dst)
    544550{
     551  if (n_IsZero(a, src)) return NULL;
    545552  /* mapping via intermediate int: */
    546553  int n = n_Int(a, src);
     
    554561{
    555562  /* dst is expected to be an algebraic field extension */
    556   assume(getCoeffType(dst) == n_algExt);
     563  assume(getCoeffType(dst) == naID);
    557564 
    558565  int h = 0; /* the height of the extension tower given by dst */
     
    585592  if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst))
    586593  {
    587     if (strcmp(rParameter(src->extRing)[0],
    588                rParameter(dst->extRing)[0]) == 0)
     594    if (strcmp(rRingVar(0, src->extRing),
     595               rRingVar(0, dst->extRing)) == 0)
    589596      return naCopyMap;                                  /// Q(a)   --> Q(a)
    590597    else
     
    605612
    606613BOOLEAN naInitChar(coeffs cf, void * infoStruct)
    607 {
    608   assume( getCoeffType(cf) == naID );
    609  
     614
    610615  AlgExtInfo *e = (AlgExtInfo *)infoStruct;
    611616  /// first check whether cf->extRing != NULL and delete old ring???
  • libpolys/polys/ext_fields/transext.cc

    r9bb5457 r2c7f28  
    44/* $Id$ */
    55/*
    6 * ABSTRACT: numbers in a transcendental extension field K(t_1, .., t_s)
    7 *           Assuming that we have a coeffs object cf, then these numbers
    8 *           are quotients of polynomials in the polynomial ring
    9 *           K[t_1, .., t_s] represented by cf->algring.
     6* ABSTRACT: numbers in a rational function field K(t_1, .., t_s) with
     7*           transcendental variables t_1, ..., t_s, where s >= 1.
     8*           Denoting the implemented coeffs object by cf, then these numbers
     9*           are represented as quotients of polynomials in the polynomial
     10*           ring K[t_1, .., t_s] represented by cf->algring.
    1011*/
    1112
     
    5556number   ntGcd(number a, number b, const coeffs cf);
    5657number   ntLcm(number a, number b, const coeffs cf);
    57 number   ntSize(number a, const coeffs cf);
     58int      ntSize(number a, const coeffs cf);
    5859void     ntDelete(number * a, const coeffs cf);
    5960void     ntCoeffWrite(const coeffs cf);
     
    6263static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param);
    6364
     65void heuristicGcdCancellation(number a, const coeffs cf);
     66void definiteGcdCancellation(number a, const coeffs cf);
     67
    6468#ifdef LDEBUG
    6569BOOLEAN ntDBTest(number a, const char *f, const int l, const coeffs cf)
    6670{
    6771  assume(getCoeffType(cf) == ntID);
    68   fraction f = (fraction)a;
    69   if (is0(f)) return TRUE;
    70   assume(num(f) != NULL);   /**< f != 0 ==> numerator(f) != 0 */
    71   p_Test(num(f), ntRing);
    72   if (!denIs1(f)) p_Test(den(f), ntRing);
     72  fraction t = (fraction)a;
     73  if (is0(t)) return TRUE;
     74  assume(num(t) != NULL);   /**< t != 0 ==> numerator(t) != 0 */
     75  p_Test(num(t), ntRing);
     76  if (!denIs1(t)) p_Test(den(t), ntRing);
    7377  return TRUE;
    7478}
    7579#endif
    76 
    77 void heuristicReduce(poly &p, poly reducer, const coeffs cf);
    78 void definiteReduce(poly &p, poly reducer, const coeffs cf);
    7980
    8081/* returns the bottom field in this field extension tower; if the tower
     
    9697}
    9798
    98 BOOLEAN naIsZero(number a, const coeffs cf)
    99 {
    100   naTest(a);
    101   return (a == NULL);
    102 }
    103 
    104 void naDelete(number * a, const coeffs cf)
    105 {
    106   if (*a == NULL) return;
    107   poly aAsPoly = (poly)(*a);
    108   p_Delete(&aAsPoly, naRing);
     99BOOLEAN ntIsZero(number a, const coeffs cf)
     100{
     101  ntTest(a);
     102  return (is0(a));
     103}
     104
     105void ntDelete(number * a, const coeffs cf)
     106{
     107  fraction f = (fraction)(*a);
     108  if (is0(f)) return;
     109  p_Delete(&num(f), ntRing);
     110  if (!denIs1(f)) p_Delete(&den(f), ntRing);
     111  omFreeBin((ADDRESS)f, fractionObjectBin);
    109112  *a = NULL;
    110113}
    111114
    112 BOOLEAN naEqual (number a, number b, const coeffs cf)
    113 {
    114   naTest(a); naTest(b);
     115BOOLEAN ntEqual(number a, number b, const coeffs cf)
     116{
     117  ntTest(a); ntTest(b);
    115118 
    116119  /// simple tests
    117120  if (a == b) return TRUE;
    118   if ((a == NULL) && (b != NULL)) return FALSE;
    119   if ((b == NULL) && (a != NULL)) return FALSE;
    120 
    121   /// deg test
    122   int aDeg = 0;
    123   if (a != NULL) aDeg = p_Deg((poly)a, naRing);
    124   int bDeg = 0;
    125   if (b != NULL) bDeg = p_Deg((poly)b, naRing);
    126   if (aDeg != bDeg) return FALSE;
    127  
    128   /// subtraction test
    129   number c = naSub(a, b, cf);
    130   BOOLEAN result = naIsZero(c, cf);
    131   naDelete(&c, naCoeffs);
    132   return result;
    133 }
    134 
    135 number naCopy(number a, const coeffs cf)
    136 {
    137   naTest(a);
    138   if (a == NULL) return NULL;
    139   return (number)p_Copy((poly)a, naRing);
    140 }
    141 
    142 number naGetNumerator(number &a, const coeffs cf)
    143 {
    144   return naCopy(a, cf);
    145 }
    146 
    147 number naGetDenom(number &a, const coeffs cf)
    148 {
    149   naTest(a);
    150   return naInit(1, cf);
    151 }
    152 
    153 BOOLEAN naIsOne(number a, const coeffs cf)
    154 {
    155   naTest(a);
    156   if (p_GetExp((poly)a, 1, naRing) != 0) return FALSE;
    157   return n_IsOne(p_GetCoeff((poly)a, naRing), naCoeffs);
    158 }
    159 
    160 BOOLEAN naIsMOne(number a, const coeffs cf)
    161 {
    162   naTest(a);
    163   if (p_GetExp((poly)a, 1, naRing) != 0) return FALSE;
    164   return n_IsMOne(p_GetCoeff((poly)a, naRing), naCoeffs);
     121  if ((is0(a)) && (!is0(b))) return FALSE;
     122  if ((is0(b)) && (!is0(a))) return FALSE;
     123 
     124  /// cheap test if gcd's have been cancelled in both numbers
     125  fraction fa = (fraction)a;
     126  fraction fb = (fraction)b;
     127  if ((c(fa) == 1) && (c(fb) == 1))
     128  {
     129    poly f = p_Add_q(p_Copy(num(fa), ntRing),
     130                     p_Neg(p_Copy(num(fb), ntRing), ntRing),
     131                     ntRing);
     132    if (f != NULL) { p_Delete(&f, ntRing); return FALSE; }
     133    if (denIs1(fa) && denIs1(fb))  return TRUE;
     134    if (denIs1(fa) && !denIs1(fb)) return FALSE;
     135    if (!denIs1(fa) && denIs1(fb)) return FALSE;
     136    f = p_Add_q(p_Copy(den(fa), ntRing),
     137                p_Neg(p_Copy(den(fb), ntRing), ntRing),
     138                ntRing);
     139    if (f != NULL) { p_Delete(&f, ntRing); return FALSE; }
     140    return TRUE;
     141  }
     142 
     143  /* default: the more expensive multiplication test
     144              a/b = c/d  <==>  a*d = b*c */
     145  poly f = p_Copy(num(fa), ntRing);
     146  if (!denIs1(fb)) f = p_Mult_q(f, p_Copy(den(fb), ntRing), ntRing);
     147  poly g = p_Copy(num(fb), ntRing);
     148  if (!denIs1(fa)) g = p_Mult_q(g, p_Copy(den(fa), ntRing), ntRing);
     149  poly h = p_Add_q(f, p_Neg(g, ntRing), ntRing);
     150  if (h == NULL) return TRUE;
     151  else
     152  {
     153    p_Delete(&h, ntRing);
     154    return FALSE;
     155  }
     156}
     157
     158number ntCopy(number a, const coeffs cf)
     159{
     160  ntTest(a);
     161  if (is0(a)) return NULL;
     162  fraction f = (fraction)a;
     163  poly g = p_Copy(num(f), ntRing);
     164  poly h = NULL; if (!denIs1(f)) h = p_Copy(den(f), ntRing);
     165  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
     166  num(result) = g;
     167  den(result) = h;
     168  c(result) = c(f);
     169  return (number)result;
     170}
     171
     172number ntGetNumerator(number &a, const coeffs cf)
     173{
     174  ntTest(a);
     175  definiteGcdCancellation(a, cf);
     176  if (is0(a)) return NULL;
     177  fraction f = (fraction)a;
     178  poly g = p_Copy(num(f), ntRing);
     179  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
     180  num(result) = g;
     181  den(result) = NULL;
     182  c(result) = 0;
     183  return (number)result;
     184}
     185
     186number ntGetDenom(number &a, const coeffs cf)
     187{
     188  ntTest(a);
     189  definiteGcdCancellation(a, cf);
     190  fraction f = (fraction)a;
     191  poly g;
     192  if (is0(f) || denIs1(f)) g = p_One(ntRing);
     193  else g = p_Copy(den(f), ntRing);
     194  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
     195  num(result) = g;
     196  den(result) = NULL;
     197  c(result) = 0;
     198  return (number)result;
     199}
     200
     201BOOLEAN ntIsOne(number a, const coeffs cf)
     202{
     203  ntTest(a);
     204  definiteGcdCancellation(a, cf);
     205  fraction f = (fraction)a;
     206  if (!denIs1(f)) return FALSE;
     207  poly g = num(f);
     208  if (!p_IsConstant(g, ntRing)) return FALSE;
     209  return n_IsOne(p_GetCoeff(g, ntRing), ntCoeffs);
     210}
     211
     212BOOLEAN ntIsMOne(number a, const coeffs cf)
     213{
     214  ntTest(a);
     215  definiteGcdCancellation(a, cf);
     216  fraction f = (fraction)a;
     217  if (!denIs1(f)) return FALSE;
     218  poly g = num(f);
     219  if (!p_IsConstant(g, ntRing)) return FALSE;
     220  return n_IsMOne(p_GetCoeff(g, ntRing), ntCoeffs);
    165221}
    166222
    167223/// this is in-place, modifies a
    168 number naNeg(number a, const coeffs cf)
    169 {
    170   naTest(a);
    171   if (a != NULL) a = (number)p_Neg((poly)a, naRing);
     224number ntNeg(number a, const coeffs cf)
     225{
     226  ntTest(a);
     227  if (!is0(a))
     228  {
     229    fraction f = (fraction)a;
     230    num(f) = p_Neg(num(f), ntRing);
     231  }
    172232  return a;
    173233}
    174234
    175 number naImPart(number a, const coeffs cf)
    176 {
    177   naTest(a);
     235number ntImPart(number a, const coeffs cf)
     236{
     237  ntTest(a);
    178238  return NULL;
    179239}
    180240
    181 number naInit(int i, const coeffs cf)
     241number ntInit(int i, const coeffs cf)
    182242{
    183243  if (i == 0) return NULL;
    184   else        return (number)p_ISet(i, naRing);
    185 }
    186 
    187 int naInt(number &a, const coeffs cf)
    188 {
    189   naTest(a);
    190   if (p_GetExp((poly)a, 1, naRing) != 0) return 0;
    191   return n_Int(p_GetCoeff((poly)a, naRing), naCoeffs);
    192 }
    193 
    194 /* TRUE iff (a != 0 and (b == 0 or deg(a) > deg(b))) */
    195 BOOLEAN naGreater(number a, number b, const coeffs cf)
    196 {
    197   if (naIsZero(a, cf)) return FALSE;
    198   if (naIsZero(b, cf)) return TRUE;
    199   int aDeg = 0;
    200   if (a != NULL) aDeg = p_Deg((poly)a, naRing);
    201   int bDeg = 0;
    202   if (b != NULL) bDeg = p_Deg((poly)b, naRing);
    203   return (aDeg > bDeg);
    204 }
    205 
    206 /* TRUE iff a != 0 and (LC(a) > 0 or deg(a) > 0) */
    207 BOOLEAN naGreaterZero(number a, const coeffs cf)
    208 {
    209   naTest(a);
    210   if (a == NULL)                                            return FALSE;
    211   if (n_GreaterZero(p_GetCoeff((poly)a, naRing), naCoeffs)) return TRUE;
    212   if (p_Deg((poly)a, naRing) > 0)                           return TRUE;
    213   return FALSE;
    214 }
    215 
    216 void naCoeffWrite(const coeffs cf)
    217 {
    218   char *x = rRingVar(0, naRing);
    219   Print("//   Coefficients live in the extension field K[%s]/<f(%s)>\n", x, x);
    220   Print("//   with the minimal polynomial f(%s) = %s\n", x,
    221         p_String(naMinpoly, naRing));
    222   PrintS("//   and K: "); n_CoeffWrite(cf->extRing->cf);
    223 }
    224 
    225 number naPar(int i, const coeffs cf)
    226 {
    227   assume(i == 1);   // there is only one parameter in this extension field
    228   poly p = p_ISet(1, naRing);   // p = 1
    229   p_SetExp(p, 1, 1, naRing);    // p = the sole extension variable
    230   p_Setm(p, naRing);
    231   return (number)p;
    232 }
    233 
    234 number naAdd(number a, number b, const coeffs cf)
    235 {
    236   naTest(a); naTest(b);
    237   if (a == NULL) return naCopy(b, cf);
    238   if (b == NULL) return naCopy(a, cf);
    239   poly aPlusB = p_Add_q(p_Copy((poly)a, naRing),
    240                         p_Copy((poly)b, naRing), naRing);
    241   definiteReduce(aPlusB, naMinpoly, cf);
    242   return (number)aPlusB;
    243 }
    244 
    245 number naSub(number a, number b, const coeffs cf)
    246 {
    247   naTest(a); naTest(b);
    248   if (b == NULL) return naCopy(a, cf);
    249   poly minusB = p_Neg(p_Copy((poly)b, naRing), naRing);
    250   if (a == NULL) return (number)minusB;
    251   poly aMinusB = p_Add_q(p_Copy((poly)a, naRing), minusB, naRing);
    252   definiteReduce(aMinusB, naMinpoly, cf);
    253   return (number)aMinusB;
    254 }
    255 
    256 number naMult(number a, number b, const coeffs cf)
    257 {
    258   naTest(a); naTest(b);
    259   if (a == NULL) return NULL;
    260   if (b == NULL) return NULL;
    261   poly aTimesB = p_Mult_q(p_Copy((poly)a, naRing),
    262                           p_Copy((poly)b, naRing), naRing);
    263   definiteReduce(aTimesB, naMinpoly, cf);
    264   return (number)aTimesB;
    265 }
    266 
    267 number naDiv(number a, number b, const coeffs cf)
    268 {
    269   naTest(a); naTest(b);
    270   if (b == NULL) WerrorS(nDivBy0);
    271   if (a == NULL) return NULL;
    272   poly bInverse = (poly)naInvers(b, cf);
    273   poly aDivB = p_Mult_q(p_Copy((poly)a, naRing), bInverse, naRing);
    274   definiteReduce(aDivB, naMinpoly, cf);
    275   return (number)aDivB;
     244  else
     245  {
     246    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
     247    num(result) = p_ISet(i, ntRing);
     248    den(result) = NULL;
     249    c(result) = 0;
     250    return (number)result;
     251  }
     252}
     253
     254int ntInt(number &a, const coeffs cf)
     255{
     256  ntTest(a);
     257  if (is0(a)) return 0;
     258  definiteGcdCancellation(a, cf);
     259  fraction f = (fraction)a;
     260  if (!denIs1(f)) return 0;
     261  if (!p_IsConstant(num(f), ntRing)) return 0;
     262  return n_Int(p_GetCoeff(num(f), ntRing), ntCoeffs);
     263}
     264
     265/* This method will only consider the numerators of a and b, without
     266   cancelling gcd's before.
     267   Moreover it may return TRUE only if one or both numerators
     268   are zero or if their degrees are equal. Then TRUE is returned iff
     269   coeff(numerator(a)) > coeff(numerator(b));
     270   In all other cases, FALSE will be returned. */
     271BOOLEAN ntGreater(number a, number b, const coeffs cf)
     272{
     273  ntTest(a); ntTest(b);
     274  number aNumCoeff = NULL; int aNumDeg = 0;
     275  number bNumCoeff = NULL; int bNumDeg = 0;
     276  if (!is0(a))
     277  {
     278    fraction fa = (fraction)a;
     279    aNumDeg = p_Totaldegree(num(fa), ntRing);
     280    aNumCoeff = p_GetCoeff(num(fa), ntRing);
     281  }
     282  if (!is0(b))
     283  {
     284    fraction fb = (fraction)b;
     285    bNumDeg = p_Totaldegree(num(fb), ntRing);
     286    bNumCoeff = p_GetCoeff(num(fb), ntRing);
     287  }
     288  if (aNumDeg != bNumDeg) return FALSE;
     289  else return n_Greater(aNumCoeff, bNumCoeff, ntCoeffs);
     290}
     291
     292/* this method will only consider the numerator of a, without cancelling
     293   the gcd before;
     294   returns TRUE iff the leading coefficient of the numerator of a is > 0
     295                    or the leading term of the numerator of a is not a
     296                    constant */
     297BOOLEAN ntGreaterZero(number a, const coeffs cf)
     298{
     299  ntTest(a);
     300  if (is0(a)) return FALSE;
     301  fraction f = (fraction)a;
     302  poly g = num(f);
     303  return (n_GreaterZero(p_GetCoeff(g, ntRing), ntCoeffs) ||
     304          (!p_LmIsConstant(g, ntRing)));
     305}
     306
     307void ntCoeffWrite(const coeffs cf)
     308{
     309  PrintS("//   Coefficients live in the rational function field\n");
     310  Print("//   K(");
     311  for (int i = 0; i < rVar(ntRing); i++)
     312  {
     313    if (i > 0) PrintS(", ");
     314    Print("%s", rRingVar(i, ntRing));
     315  }
     316  PrintS(") with\n");
     317  PrintS("//   K: "); n_CoeffWrite(cf->extRing->cf);
     318}
     319
     320/* the i-th parameter */
     321number ntPar(int i, const coeffs cf)
     322{
     323  assume((1 <= i) && (i <= rVar(ntRing)));
     324  poly p = p_ISet(1, ntRing);
     325  p_SetExp(p, i, 1, ntRing);
     326  p_Setm(p, ntRing);
     327  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
     328  num(result) = p;
     329  den(result) = NULL;
     330  c(result) = 0;
     331  return (number)result;
     332}
     333
     334number ntAdd(number a, number b, const coeffs cf)
     335{
     336  ntTest(a); ntTest(b);
     337  if (is0(a)) return ntCopy(b, cf);
     338  if (is0(b)) return ntCopy(a, cf);
     339 
     340  fraction fa = (fraction)a;
     341  fraction fb = (fraction)b;
     342  poly f;
     343  if      (denIs1(fa) && denIs1(fb))  f = NULL;
     344  else if (!denIs1(fa) && denIs1(fb)) f = p_Copy(den(fa), ntRing);
     345  else if (denIs1(fa) && !denIs1(fb)) f = p_Copy(den(fb), ntRing);
     346  else /* both denom's are != 1 */    f = p_Mult_q(p_Copy(den(fa), ntRing),
     347                                                   p_Copy(den(fb), ntRing),
     348                                                   ntRing);
     349  poly g = p_Copy(num(fa), ntRing);
     350  if (!denIs1(fb)) g = p_Mult_q(g, p_Copy(den(fb), ntRing), ntRing);
     351  poly h = p_Copy(num(fb), ntRing);
     352  if (!denIs1(fa)) h = p_Mult_q(h, p_Copy(den(fa), ntRing), ntRing);
     353  g = p_Add_q(g, h, ntRing);
     354  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
     355  num(result) = g;
     356  den(result) = f;
     357  c(result) = c(fa) + c(fb) + ADD_COMPLEXITY;
     358  heuristicGcdCancellation((number)result, cf);
     359  return (number)result;
     360}
     361
     362number ntSub(number a, number b, const coeffs cf)
     363{
     364  ntTest(a); ntTest(b);
     365  if (is0(a)) return ntNeg(ntCopy(b, cf), cf);
     366  if (is0(b)) return ntCopy(a, cf);
     367 
     368  fraction fa = (fraction)a;
     369  fraction fb = (fraction)b;
     370  poly f;
     371  if      (denIs1(fa) && denIs1(fb))  f = NULL;
     372  else if (!denIs1(fa) && denIs1(fb)) f = p_Copy(den(fa), ntRing);
     373  else if (denIs1(fa) && !denIs1(fb)) f = p_Copy(den(fb), ntRing);
     374  else /* both den's are != 1 */      f = p_Mult_q(p_Copy(den(fa), ntRing),
     375                                                   p_Copy(den(fb), ntRing),
     376                                                   ntRing);
     377  poly g = p_Copy(num(fa), ntRing);
     378  if (!denIs1(fb)) g = p_Mult_q(g, p_Copy(den(fb), ntRing), ntRing);
     379  poly h = p_Copy(num(fb), ntRing);
     380  if (!denIs1(fa)) h = p_Mult_q(h, p_Copy(den(fa), ntRing), ntRing);
     381  g = p_Add_q(g, p_Neg(h, ntRing), ntRing);
     382  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
     383  num(result) = g;
     384  den(result) = f;
     385  c(result) = c(fa) + c(fb) + ADD_COMPLEXITY;
     386  heuristicGcdCancellation((number)result, cf);
     387  return (number)result;
     388}
     389
     390number ntMult(number a, number b, const coeffs cf)
     391{
     392  ntTest(a); ntTest(b);
     393  if (is0(a) || is0(b)) return NULL;
     394 
     395  fraction fa = (fraction)a;
     396  fraction fb = (fraction)b;
     397  poly f;
     398  if      (denIs1(fa) && denIs1(fb))  f = NULL;
     399  else if (!denIs1(fa) && denIs1(fb)) f = p_Copy(den(fa), ntRing);
     400  else if (denIs1(fa) && !denIs1(fb)) f = p_Copy(den(fb), ntRing);
     401  else /* both den's are != 1 */      f = p_Mult_q(p_Copy(den(fa), ntRing),
     402                                                   p_Copy(den(fb), ntRing),
     403                                                   ntRing);
     404  poly g = p_Copy(num(fa), ntRing);
     405  poly h = p_Copy(num(fb), ntRing);
     406  g = p_Mult_q(g, h, ntRing);
     407  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
     408  num(result) = g;
     409  den(result) = f;
     410  c(result) = c(fa) + c(fb) + MULT_COMPLEXITY;
     411  heuristicGcdCancellation((number)result, cf);
     412  return (number)result;
     413}
     414
     415number ntDiv(number a, number b, const coeffs cf)
     416{
     417  ntTest(a); ntTest(b);
     418  if (is0(a)) return NULL;
     419  if (is0(b)) WerrorS(nDivBy0);
     420 
     421  fraction fa = (fraction)a;
     422  fraction fb = (fraction)b;
     423  poly f = p_Copy(num(fb), ntRing);
     424  if (!denIs1(fa)) f = p_Mult_q(f, p_Copy(den(fa), ntRing), ntRing);
     425  poly g = p_Copy(num(fa), ntRing);
     426  if (!denIs1(fb)) g = p_Mult_q(g, p_Copy(den(fb), ntRing), ntRing);
     427  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
     428  num(result) = g;
     429  den(result) = f;
     430  c(result) = c(fa) + c(fb) + MULT_COMPLEXITY;
     431  heuristicGcdCancellation((number)result, cf);
     432  return (number)result;
    276433}
    277434
     
    280437   for |exp| >= 8 compute power along binary presentation of |exp|, e.g.
    281438      p^13 = p^1 * p^4 * p^8, where we utilise that
    282       p^(2^(k+1)) = p^(2^k) * p^(2^k)
    283    intermediate reduction modulo the minimal polynomial is controlled by
    284    the in-place method heuristicReduce(poly, poly, coeffs); see there.
     439      p^(2^(k+1)) = p^(2^k) * p^(2^k);
     440   intermediate cancellation is controlled by the in-place method
     441   heuristicGcdCancellation; see there.
    285442*/
    286 void naPower(number a, int exp, number *b, const coeffs cf)
    287 {
    288   naTest(a);
     443void ntPower(number a, int exp, number *b, const coeffs cf)
     444{
     445  ntTest(a);
    289446 
    290447  /* special cases first */
    291   if (a == NULL)
     448  if (is0(a))
    292449  {
    293450    if (exp >= 0) *b = NULL;
    294451    else          WerrorS(nDivBy0);
    295452  }
    296   else if (exp ==  0) *b = naInit(1, cf);
    297   else if (exp ==  1) *b = naCopy(a, cf);
    298   else if (exp == -1) *b = naInvers(a, cf);
     453  else if (exp ==  0) *b = ntInit(1, cf);
     454  else if (exp ==  1) *b = ntCopy(a, cf);
     455  else if (exp == -1) *b = ntInvers(a, cf);
    299456 
    300457  int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
    301458 
    302   /* now compute 'a' to the 'expAbs'-th power */
    303   poly pow; poly aAsPoly = (poly)a;
     459  /* now compute a^expAbs */
     460  number pow; number t;
    304461  if (expAbs <= 7)
    305462  {
    306     pow = p_Copy(aAsPoly, naRing);
     463    pow = ntCopy(a, cf);
    307464    for (int i = 2; i <= expAbs; i++)
    308465    {
    309       pow = p_Mult_q(pow, p_Copy(aAsPoly, naRing), naRing);
    310       heuristicReduce(pow, naMinpoly, cf);
     466      t = ntMult(pow, a, cf);
     467      ntDelete(&pow, cf);
     468      pow = t;
     469      heuristicGcdCancellation(pow, cf);
    311470    }
    312     definiteReduce(pow, naMinpoly, cf);
    313471  }
    314472  else
    315473  {
    316     pow = p_ISet(1, naRing);
    317     poly factor = p_Copy(aAsPoly, naRing);
     474    pow = ntInit(1, cf);
     475    number factor = ntCopy(a, cf);
    318476    while (expAbs != 0)
    319477    {
    320478      if (expAbs & 1)
    321479      {
    322         pow = p_Mult_q(pow, p_Copy(factor, naRing), naRing);
    323         heuristicReduce(pow, naMinpoly, cf);
     480        t = ntMult(pow, factor, cf);
     481        ntDelete(&pow, cf);
     482        pow = t;
     483        heuristicGcdCancellation(pow, cf);
    324484      }
    325485      expAbs = expAbs / 2;
    326486      if (expAbs != 0)
    327487      {
    328         factor = p_Mult_q(factor, factor, naRing);
    329         heuristicReduce(factor, naMinpoly, cf);
     488        t = ntMult(factor, factor, cf);
     489        ntDelete(&factor, cf);
     490        factor = t;
     491        heuristicGcdCancellation(factor, cf);
    330492      }
    331493    }
    332     p_Delete(&factor, naRing);
    333     definiteReduce(pow, naMinpoly, cf);
     494    ntDelete(&factor, cf);
    334495  }
    335496 
    336497  /* invert if original exponent was negative */
    337   number n = (number)pow;
    338498  if (exp < 0)
    339499  {
    340     number m = naInvers(n, cf);
    341     naDelete(&n, cf);
    342     n = m;
    343   }
    344   *b = n;
    345 }
    346 
    347 /* may reduce p module the reducer by calling definiteReduce;
    348    the decision is made based on the following heuristic
    349    (which should also only be changed here in this method):
    350       if (deg(p) > 10*deg(reducer) then perform reduction;
    351    modifies p */
    352 void heuristicReduce(poly &p, poly reducer, const coeffs cf)
    353 {
    354   #ifdef LDEBUG
    355   p_Test((poly)p, naRing);
    356   p_Test((poly)reducer, naRing);
    357   #endif
    358   if (p_Deg(p, naRing) > 10 * p_Deg(reducer, naRing))
    359     definiteReduce(p, reducer, cf);
    360 }
    361 
    362 void naWrite(number &a, const coeffs cf)
    363 {
    364   naTest(a);
    365   if (a == NULL)
     500    t = ntInvers(pow, cf);
     501    ntDelete(&pow, cf);
     502    pow = t;
     503  }
     504  *b = pow;
     505}
     506
     507/* modifies a */
     508void heuristicGcdCancellation(number a, const coeffs cf)
     509{
     510  ntTest(a);
     511  if (is0(a)) return;
     512  fraction f = (fraction)a;
     513  if (denIs1(f) || n_IsOne(p_GetCoeff(num(f), ntRing), ntCoeffs))
     514  { c(f) = 0; return; }
     515  if (c(f) <= BOUND_COMPLEXITY) return;
     516  else definiteGcdCancellation(a, cf);
     517}
     518
     519/* modifies a */
     520void definiteGcdCancellation(number a, const coeffs cf)
     521{
     522  ntTest(a);
     523  if (is0(a)) return;
     524  fraction f = (fraction)a;
     525  if (denIs1(f) || n_IsOne(p_GetCoeff(num(f), ntRing), ntCoeffs))
     526  { c(f) = 0; return; }
     527  if (c(f) > BOUND_COMPLEXITY)
     528  {
     529    /* TO BE IMPLEMENTED!
     530       for the time, cancellation of gcd's does not take place */
     531    Print("// TO BE IMPLEMENTED: transext.cc:definiteGcdCancellation\n");
     532    Print("// (complexity of number = %d exceeds the bound = %d\n",
     533          c(f), BOUND_COMPLEXITY);
     534  }
     535}
     536
     537void ntWrite(number &a, const coeffs cf)
     538{
     539  ntTest(a);
     540  definiteGcdCancellation(a, cf);
     541  if (is0(a))
    366542    StringAppendS("0");
    367543  else
    368544  {
    369     poly aAsPoly = (poly)a;
    370     /* basically, just write aAsPoly using p_Write,
    371        but use brackets around the output, if a is not
    372        a constant living in naCoeffs = cf->extRing->cf */
    373     BOOLEAN useBrackets = !(p_IsConstant(aAsPoly, naRing));
     545    fraction f = (fraction)a;
     546    BOOLEAN useBrackets = !(p_IsConstant(num(f), ntRing)) && (!denIs1(f));
    374547    if (useBrackets) StringAppendS("(");
    375     p_String0(aAsPoly, naRing, naRing);
     548    p_String0(num(f), ntRing, ntRing);
    376549    if (useBrackets) StringAppendS(")");
    377   }
    378 }
    379 
    380 const char * naRead(const char *s, number *a, const coeffs cf)
    381 {
    382   poly aAsPoly;
    383   const char * result = p_Read(s, aAsPoly, naRing);
    384   *a = (number)aAsPoly;
    385   return result;
    386 }
    387 
    388 /* implemented by the rule lcm(a, b) = a * b / gcd(a, b) */
    389 number naLcm(number a, number b, const coeffs cf)
    390 {
    391   naTest(a); naTest(b);
    392   if (a == NULL) return NULL;
    393   if (b == NULL) return NULL;
    394   number theProduct = (number)p_Mult_q(p_Copy((poly)a, naRing),
    395                                        p_Copy((poly)b, naRing), naRing);
    396   /* note that theProduct needs not be reduced w.r.t. naMinpoly;
    397      but the final division will take care of the necessary reduction */
    398   number theGcd = naGcd(a, b, cf);
    399   return naDiv(theProduct, theGcd, cf);
    400 }
    401 
    402 /* expects *param to be castable to AlgExtInfo */
    403 static BOOLEAN naCoeffIsEqual(const coeffs cf, n_coeffType n, void * param)
    404 {
    405   if (naID != n) return FALSE;
    406   AlgExtInfo *e = (AlgExtInfo *)param;
    407   /* for extension coefficient fields we expect the underlying
    408      polynomials rings to be IDENTICAL, i.e. the SAME OBJECT;
     550    if (!denIs1(f))
     551    {
     552      StringAppendS("/");
     553      useBrackets = !p_IsConstant(den(f), ntRing);
     554      if (useBrackets) StringAppendS("(");
     555      p_String0(den(f), ntRing, ntRing);
     556      if (useBrackets) StringAppendS(")");
     557    }
     558  }
     559}
     560
     561/* this just reads polynomials;
     562   the rest will be taken care of by the interpreter: When it encounters
     563   a "/", then it will impose a division that will finally lead to a
     564   rational function */
     565const char * ntRead(const char *s, number *a, const coeffs cf)
     566{
     567  poly p;
     568  const char * result = p_Read(s, p, ntRing);
     569  if (p == NULL) { *a = NULL; return result; }
     570  else
     571  {
     572    fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
     573    num(f) = p;
     574    den(f) = NULL;
     575    c(f) = 0;
     576    *a = (number)f;
     577    return result;
     578  }
     579}
     580
     581/* expects *param to be castable to TransExtInfo */
     582static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param)
     583{
     584  if (ntID != n) return FALSE;
     585  TransExtInfo *e = (TransExtInfo *)param;
     586  /* for rational function fields we expect the underlying
     587     polynomial rings to be IDENTICAL, i.e. the SAME OBJECT;
    409588     this expectation is based on the assumption that we have properly
    410589     registered cf and perform reference counting rather than creating
    411590     multiple copies of the same coefficient field/domain/ring */
    412   return (naRing == e->r);
    413   /* (Note that then also the minimal ideals will necessarily be
    414      the same, as they are attached to the ring.) */
    415 }
    416 
    417 int naSize(number a, const coeffs cf)
    418 {
    419   if (a == NULL) return -1;
     591  return (ntRing == e->r);
     592}
     593
     594number ntLcm(number a, number b, const coeffs cf)
     595{
     596  ntTest(a); ntTest(b);
     597  /* TO BE IMPLEMENTED!
     598     for the time, we simply return NULL, representing the number zero */
     599  Print("// TO BE IMPLEMENTED: transext.cc:ntLcm\n");
     600  return NULL;
     601}
     602
     603number ntGcd(number a, number b, const coeffs cf)
     604{
     605  ntTest(a); ntTest(b);
     606  /* TO BE IMPLEMENTED!
     607     for the time, we simply return NULL, representing the number zero */
     608  Print("// TO BE IMPLEMENTED: transext.cc:ntGcd\n");
     609  return NULL;
     610}
     611
     612int ntSize(number a, const coeffs cf)
     613{
     614  ntTest(a);
     615  if (is0(a)) return -1;
    420616  /* this has been taken from the old implementation of field extensions,
    421      where we computed the sum of the degree and the number of terms in
    422      (poly)a; so we leave it at that, for the time being;
    423      maybe, the number of terms alone is a better measure? */
    424   poly aAsPoly = (poly)a;
    425   int theDegree = 0; int noOfTerms = 0;
    426   while (aAsPoly != NULL)
     617     where we computed the sum of the degrees and the numbers of terms in
     618     the numerator and denominator of a; so we leave it at that, for the
     619     time being */
     620  fraction f = (fraction)a;
     621  poly p = num(f);
     622  int noOfTerms = 0;
     623  int numDegree = 0;
     624  while (p != NULL)
    427625  {
    428626    noOfTerms++;
    429627    int d = 0;
    430     for (int i = 1; i <= rVar(naRing); i++)
    431       d += p_GetExp(aAsPoly, i, naRing);
    432     if (d > theDegree) theDegree = d;
    433     pIter(aAsPoly);
    434   }
    435   return theDegree + noOfTerms;
    436 }
    437 
    438 /* performs polynomial division and overrides p by the remainder
    439    of division of p by the reducer;
    440    modifies p */
    441 void definiteReduce(poly &p, poly reducer, const coeffs cf)
    442 {
    443   #ifdef LDEBUG
    444   p_Test((poly)p, naRing);
    445   p_Test((poly)reducer, naRing);
    446   #endif
    447   p_PolyDiv(p, reducer, FALSE, naRing);
    448 }
    449 
    450 /* IMPORTANT NOTE: Since an algebraic field extension is again a field,
    451                    the gcd of two elements is not very interesting. (It
    452                    is actually any unit in the field, i.e., any non-
    453                    zero element.) Note that the below method does not operate
    454                    in this strong sense but rather computes the gcd of
    455                    two given elements in the underlying polynomial ring. */
    456 number naGcd(number a, number b, const coeffs cf)
    457 {
    458   naTest(a); naTest(b);
    459   if ((a == NULL) && (b == NULL)) WerrorS(nDivBy0);
    460   return (number)p_Gcd((poly)a, (poly)b, naRing);
    461 }
    462 
    463 number naInvers(number a, const coeffs cf)
    464 {
    465   naTest(a);
    466   if (a == NULL) WerrorS(nDivBy0);
    467   poly aFactor = NULL; poly mFactor = NULL;
    468   poly theGcd = p_ExtGcd((poly)a, aFactor, naMinpoly, mFactor, naRing);
    469   naTest((number)theGcd); naTest((number)aFactor); naTest((number)mFactor);
    470   /* the gcd must be 1 since naMinpoly is irreducible and a != NULL: */
    471   assume(naIsOne((number)theGcd, cf));     
    472   p_Delete(&theGcd, naRing);
    473   p_Delete(&mFactor, naRing);
    474   return (number)(aFactor);
    475 }
    476 
    477 /* assumes that src = Q, dst = Q(a) */
    478 number naMap00(number a, const coeffs src, const coeffs dst)
    479 {
     628    for (int i = 1; i <= rVar(ntRing); i++)
     629      d += p_GetExp(p, i, ntRing);
     630    if (d > numDegree) numDegree = d;
     631    pIter(p);
     632  }
     633  int denDegree = 0;
     634  if (!denIs1(f))
     635  {
     636    p = den(f);
     637    while (p != NULL)
     638    {
     639      noOfTerms++;
     640      int d = 0;
     641      for (int i = 1; i <= rVar(ntRing); i++)
     642        d += p_GetExp(p, i, ntRing);
     643      if (d > denDegree) denDegree = d;
     644      pIter(p);
     645    }
     646  }
     647  return numDegree + denDegree + noOfTerms;
     648}
     649
     650number ntInvers(number a, const coeffs cf)
     651{
     652  ntTest(a);
     653  if (is0(a)) WerrorS(nDivBy0);
     654  fraction f = (fraction)a;
     655  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
     656  poly g;
     657  if (denIs1(f)) g = p_One(ntRing);
     658  else           g = p_Copy(den(f), ntRing);
     659  num(result) = g;
     660  den(result) = p_Copy(num(f), ntRing);
     661  c(result) = c(f);
     662  return (number)result;
     663}
     664
     665/* assumes that src = Q, dst = Q(t_1, ..., t_s) */
     666number ntMap00(number a, const coeffs src, const coeffs dst)
     667{
     668  if (n_IsZero(a, src)) return NULL;
    480669  assume(src == dst->extRing->cf);
    481   poly result = p_One(dst->extRing);
    482   p_SetCoeff(result, naCopy(a, src), dst->extRing);
    483   return (number)result;
    484 }
    485 
    486 /* assumes that src = Z/p, dst = Q(a) */
    487 number naMapP0(number a, const coeffs src, const coeffs dst)
    488 {
     670  poly p = p_One(dst->extRing);
     671  p_SetCoeff(p, ntCopy(a, src), dst->extRing);
     672  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
     673  num(f) = p; den(f) = NULL; c(f) = 0;
     674  return (number)f;
     675}
     676
     677/* assumes that src = Z/p, dst = Q(t_1, ..., t_s) */
     678number ntMapP0(number a, const coeffs src, const coeffs dst)
     679{
     680  if (n_IsZero(a, src)) return NULL;
    489681  /* mapping via intermediate int: */
    490682  int n = n_Int(a, src);
    491683  number q = n_Init(n, dst->extRing->cf);
    492   poly result = p_One(dst->extRing);
    493   p_SetCoeff(result, q, dst->extRing);
    494   return (number)result;
    495 }
    496 
    497 /* assumes that either src = Q(a), dst = Q(a), or
    498                        src = Zp(a), dst = Zp(a)     */
    499 number naCopyMap(number a, const coeffs src, const coeffs dst)
    500 {
    501   return naCopy(a, dst);
    502 }
    503 
    504 /* assumes that src = Q, dst = Z/p(a) */
    505 number naMap0P(number a, const coeffs src, const coeffs dst)
    506 {
     684  poly p;
     685  if (n_IsZero(q, dst->extRing->cf))
     686  {
     687    n_Delete(&q, dst->extRing->cf);
     688    return NULL;
     689  }
     690  p = p_One(dst->extRing);
     691  p_SetCoeff(p, q, dst->extRing);
     692  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
     693  num(f) = p; den(f) = NULL; c(f) = 0;
     694  return (number)f;
     695}
     696
     697/* assumes that either src = Q(t_1, ..., t_s), dst = Q(t_1, ..., t_s), or
     698                       src = Z/p(t_1, ..., t_s), dst = Z/p(t_1, ..., t_s) */
     699number ntCopyMap(number a, const coeffs src, const coeffs dst)
     700{
     701  return ntCopy(a, dst);
     702}
     703
     704/* assumes that src = Q, dst = Z/p(t_1, ..., t_s) */
     705number ntMap0P(number a, const coeffs src, const coeffs dst)
     706{
     707  if (n_IsZero(a, src)) return NULL;
    507708  int p = rChar(dst->extRing);
    508709  int n = nlModP(a, p, src);
    509710  number q = n_Init(n, dst->extRing->cf);
    510   poly result = p_One(dst->extRing);
    511   p_SetCoeff(result, q, dst->extRing);
    512   return (number)result;
    513 }
    514 
    515 /* assumes that src = Z/p, dst = Z/p(a) */
    516 number naMapPP(number a, const coeffs src, const coeffs dst)
    517 {
     711  poly g;
     712  if (n_IsZero(q, dst->extRing->cf))
     713  {
     714    n_Delete(&q, dst->extRing->cf);
     715    return NULL;
     716  }
     717  g = p_One(dst->extRing);
     718  p_SetCoeff(g, q, dst->extRing);
     719  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
     720  num(f) = g; den(f) = NULL; c(f) = 0;
     721  return (number)f;
     722}
     723
     724/* assumes that src = Z/p, dst = Z/p(t_1, ..., t_s) */
     725number ntMapPP(number a, const coeffs src, const coeffs dst)
     726{
     727  if (n_IsZero(a, src)) return NULL;
    518728  assume(src == dst->extRing->cf);
    519   poly result = p_One(dst->extRing);
    520   p_SetCoeff(result, naCopy(a, src), dst->extRing);
    521   return (number)result;
    522 }
    523 
    524 /* assumes that src = Z/u, dst = Z/p(a), where u != p */
    525 number naMapUP(number a, const coeffs src, const coeffs dst)
    526 {
     729  poly p = p_One(dst->extRing);
     730  p_SetCoeff(p, ntCopy(a, src), dst->extRing);
     731  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
     732  num(f) = p; den(f) = NULL; c(f) = 0;
     733  return (number)f;
     734}
     735
     736/* assumes that src = Z/u, dst = Z/p(t_1, ..., t_s), where u != p */
     737number ntMapUP(number a, const coeffs src, const coeffs dst)
     738{
     739  if (n_IsZero(a, src)) return NULL;
    527740  /* mapping via intermediate int: */
    528741  int n = n_Int(a, src);
    529742  number q = n_Init(n, dst->extRing->cf);
    530   poly result = p_One(dst->extRing);
    531   p_SetCoeff(result, q, dst->extRing);
    532   return (number)result;
    533 }
    534 
    535 nMapFunc naSetMap(const coeffs src, const coeffs dst)
    536 {
    537   /* dst is expected to be an algebraic field extension */
    538   assume(getCoeffType(dst) == n_algExt);
     743  poly p;
     744  if (n_IsZero(q, dst->extRing->cf))
     745  {
     746    n_Delete(&q, dst->extRing->cf);
     747    return NULL;
     748  }
     749  p = p_One(dst->extRing);
     750  p_SetCoeff(p, q, dst->extRing);
     751  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
     752  num(f) = p; den(f) = NULL; c(f) = 0;
     753  return (number)f;
     754}
     755
     756nMapFunc ntSetMap(const coeffs src, const coeffs dst)
     757{
     758  /* dst is expected to be a rational function field */
     759  assume(getCoeffType(dst) == ntID);
    539760 
    540761  int h = 0; /* the height of the extension tower given by dst */
     
    546767  if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
    547768 
     769  /* Let T denote the sequence of transcendental extension variables, i.e.,
     770     K[t_1, ..., t_s] =: K[T];
     771     Let moreover, for any such sequence T, T' denote any subsequence of T
     772     of the form t_1, ..., t_w with w <= s. */
     773 
    548774  if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst))
    549     return naMap00;                                      /// Q     -->  Q(a)
     775    return ntMap00;                                      /// Q       -->  Q(T)
    550776 
    551777  if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
    552     return naMapP0;                                      /// Z/p   -->  Q(a)
     778    return ntMapP0;                                      /// Z/p     -->  Q(T)
    553779 
    554780  if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst))
    555     return naMap0P;                                      /// Q      --> Z/p(a)
     781    return ntMap0P;                                      /// Q       --> Z/p(T)
    556782 
    557783  if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
    558784  {
    559     if (src->ch == dst->ch) return naMapPP;              /// Z/p    --> Z/p(a)
    560     else return naMapUP;                                 /// Z/u    --> Z/p(a)
     785    if (src->ch == dst->ch) return ntMapPP;              /// Z/p     --> Z/p(T)
     786    else return ntMapUP;                                 /// Z/u     --> Z/p(T)
    561787  }
    562788 
     
    567793  if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst))
    568794  {
    569     if (strcmp(rParameter(src->extRing)[0],
    570                rParameter(dst->extRing)[0]) == 0)
    571       return naCopyMap;                                  /// Q(a)   --> Q(a)
    572     else
    573       return NULL;                                       /// Q(b)   --> Q(a)
     795    if (rVar(src->extRing) > rVar(dst->extRing)) return NULL;
     796    for (int i = 0; i < rVar(src->extRing); i++)
     797      if (strcmp(rRingVar(i, src->extRing),
     798                 rRingVar(i, dst->extRing)) != 0) return NULL;
     799      return ntCopyMap;                                  /// Q(T')   --> Q(T)
    574800  }
    575801 
    576802  if (nCoeff_is_Zp(bSrc) && nCoeff_is_Zp(bDst))
    577803  {
    578     if (strcmp(rParameter(src->extRing)[0],
    579                rParameter(dst->extRing)[0]) == 0)
    580       return naCopyMap;                                  /// Z/p(a) --> Z/p(a)
    581     else
    582       return NULL;                                       /// Z/p(b) --> Z/p(a)
     804    if (rVar(src->extRing) > rVar(dst->extRing)) return NULL;
     805    for (int i = 0; i < rVar(src->extRing); i++)
     806      if (strcmp(rRingVar(i, src->extRing),
     807                 rRingVar(i, dst->extRing)) != 0) return NULL;
     808      return ntCopyMap;                                  /// Z/p(T') --> Z/p(T)
    583809  }
    584810 
     
    586812}
    587813
    588 BOOLEAN naInitChar(coeffs cf, void * infoStruct)
    589 {
    590   assume( getCoeffType(cf) == naID );
    591  
    592   AlgExtInfo *e = (AlgExtInfo *)infoStruct;
     814BOOLEAN ntInitChar(coeffs cf, void * infoStruct)
     815
     816  TransExtInfo *e = (TransExtInfo *)infoStruct;
    593817  /// first check whether cf->extRing != NULL and delete old ring???
    594818  cf->extRing           = e->r;
    595   cf->extRing->minideal = e->i;
    596 
    597   assume(cf->extRing                     != NULL);      // extRing;
    598   assume((cf->extRing->minideal          != NULL) &&    // minideal has one
    599          (IDELEMS(cf->extRing->minideal) != 0)    &&    // non-zero generator
    600          (cf->extRing->minideal->m[0]    != NULL)    ); // at m[0];
    601   assume(cf->extRing->cf                 != NULL);      // extRing->cf;
    602   assume(getCoeffType(cf) == naID);                     // coeff type;
     819  cf->extRing->minideal = NULL;
     820
     821  assume(cf->extRing                != NULL);      // extRing;
     822  assume(cf->extRing->cf            != NULL);      // extRing->cf;
     823  assume(getCoeffType(cf) == ntID);                // coeff type;
    603824 
    604825  /* propagate characteristic up so that it becomes
     
    606827  cf->ch = cf->extRing->cf->ch;
    607828 
    608   #ifdef LDEBUG
    609   p_Test((poly)naMinpoly, naRing);
    610   #endif
    611  
    612   cf->cfGreaterZero  = naGreaterZero;
    613   cf->cfGreater      = naGreater;
    614   cf->cfEqual        = naEqual;
    615   cf->cfIsZero       = naIsZero;
    616   cf->cfIsOne        = naIsOne;
    617   cf->cfIsMOne       = naIsMOne;
    618   cf->cfInit         = naInit;
    619   cf->cfInt          = naInt;
    620   cf->cfNeg          = naNeg;
    621   cf->cfPar          = naPar;
    622   cf->cfAdd          = naAdd;
    623   cf->cfSub          = naSub;
    624   cf->cfMult         = naMult;
    625   cf->cfDiv          = naDiv;
    626   cf->cfExactDiv     = naDiv;
    627   cf->cfPower        = naPower;
    628   cf->cfCopy         = naCopy;
    629   cf->cfWrite        = naWrite;
    630   cf->cfRead         = naRead;
    631   cf->cfDelete       = naDelete;
    632   cf->cfSetMap       = naSetMap;
    633   cf->cfGetDenom     = naGetDenom;
    634   cf->cfGetNumerator = naGetNumerator;
    635   cf->cfRePart       = naCopy;
    636   cf->cfImPart       = naImPart;
    637   cf->cfCoeffWrite   = naCoeffWrite;
    638   cf->cfDBTest       = naDBTest;
    639   cf->cfGcd          = naGcd;
    640   cf->cfLcm          = naLcm;
    641   cf->cfSize         = naSize;
    642   cf->nCoeffIsEqual  = naCoeffIsEqual;
    643   cf->cfInvers       = naInvers;
    644   cf->cfIntDiv       = naDiv;
     829  cf->cfGreaterZero  = ntGreaterZero;
     830  cf->cfGreater      = ntGreater;
     831  cf->cfEqual        = ntEqual;
     832  cf->cfIsZero       = ntIsZero;
     833  cf->cfIsOne        = ntIsOne;
     834  cf->cfIsMOne       = ntIsMOne;
     835  cf->cfInit         = ntInit;
     836  cf->cfInt          = ntInt;
     837  cf->cfNeg          = ntNeg;
     838  cf->cfPar          = ntPar;
     839  cf->cfAdd          = ntAdd;
     840  cf->cfSub          = ntSub;
     841  cf->cfMult         = ntMult;
     842  cf->cfDiv          = ntDiv;
     843  cf->cfExactDiv     = ntDiv;
     844  cf->cfPower        = ntPower;
     845  cf->cfCopy         = ntCopy;
     846  cf->cfWrite        = ntWrite;
     847  cf->cfRead         = ntRead;
     848  cf->cfDelete       = ntDelete;
     849  cf->cfSetMap       = ntSetMap;
     850  cf->cfGetDenom     = ntGetDenom;
     851  cf->cfGetNumerator = ntGetNumerator;
     852  cf->cfRePart       = ntCopy;
     853  cf->cfImPart       = ntImPart;
     854  cf->cfCoeffWrite   = ntCoeffWrite;
     855  cf->cfDBTest       = ntDBTest;
     856  cf->cfGcd          = ntGcd;
     857  cf->cfLcm          = ntLcm;
     858  cf->cfSize         = ntSize;
     859  cf->nCoeffIsEqual  = ntCoeffIsEqual;
     860  cf->cfInvers       = ntInvers;
     861  cf->cfIntDiv       = ntDiv;
    645862 
    646863  return FALSE;
  • libpolys/polys/ext_fields/transext.h

    r9bb5457 r2c7f28  
    66/* $Id$ */
    77/*
    8 * ABSTRACT: numbers in a transcendental extension field K(t_1, .., t_s)
    9 *           Assuming that we have a coeffs object cf, then these numbers
    10 *           are quotients of polynomials in the polynomial ring
    11 *           K[t_1, .., t_s] represented by cf->algring.
     8* ABSTRACT: numbers in a rational function field K(t_1, .., t_s) with
     9*           transcendental variables t_1, ..., t_s, where s >= 1.
     10*           Denoting the implemented coeffs object by cf, then these numbers
     11*           are represented as quotients of polynomials in the polynomial
     12*           ring K[t_1, .., t_s] represented by cf->algring.
    1213*/
    1314
     
    2021typedef struct sip_sideal * ideal;
    2122
     23struct spolyrec;
     24typedef struct spolyrec polyrec;
     25typedef polyrec * poly;
     26
    2227/// struct for passing initialization parameters to naInitChar
    23 typedef struct { ring r; } transExtInfo;
     28typedef struct { ring r; } TransExtInfo;
    2429
    2530/* a number in K(t_1, .., t_s) is represented by either NULL
     
    2833   if the denominator is 1, the member 'denominator' is NULL;
    2934   as a consequence of the above we get: if some number n is not NULL, then
    30    n->numerator cannot be NULL */
    31 struct { poly numerator; poly denominator; } ip_fraction;
    32 typedef struct ip_fraction * fraction;
    33 /* some useful accessors: */
    34 #define is0(n) (n == NULL) /**< TRUE iff n represents 0 in K(t_1, .., t_s) */
    35 #define num(n) n->numerator
    36 #define den(n) n->denominator
    37 #define denIs1(n) (n->denominator == NULL) /**< TRUE iff the denom. is 1 */
     35   n->numerator cannot be NULL;
     36   The member 'complexity' attempts to capture the complexity of any given
     37   number n, i.e., starting with a bunch of numbers n_i that have their gcd's
     38   cancelled out, n may be constructed from the n_i's by using field
     39   arithmetics (+, -, *, /). If we never cancel out gcd's during this process,
     40   n will become rather complex. The larger the attribute 'complexity' of n
     41   is, the more likely it is that n contains some non-trivial gcd. Thus, this
     42   attribute will be used by a heuristic method to cancel out gcd's from time
     43   to time. (This heuristic may be set up such that cancellation can be
     44   enforced after each arithmetic operation, or such that it will never take
     45   place.) Moreover, the 'complexity' of n is zero iff the gcd in n (that is,
     46   the gcd of its numerator and denominator) is trivial. */
     47struct fractionObject
     48{
     49  poly numerator;
     50  poly denominator;
     51  int complexity;
     52};
     53typedef struct fractionObject * fraction;
     54
     55omBin fractionObjectBin = omGetSpecBin(sizeof(fractionObject));
     56
     57/* constants for controlling the complexity of numbers */
     58#define ADD_COMPLEXITY 1   /**< complexity increase due to + and - */
     59#define MULT_COMPLEXITY 2   /**< complexity increase due to * and / */
     60#define BOUND_COMPLEXITY 10   /**< maximum complexity of a number */
     61
     62/* some useful accessors for fractions: */
     63#define is0(f) (f == NULL) /**< TRUE iff n represents 0 in K(t_1, .., t_s) */
     64#define num(f) f->numerator
     65#define den(f) f->denominator
     66#define denIs1(f) (f->denominator == NULL) /**< TRUE iff den. represents 1 */
     67#define c(f) f->complexity
    3868
    3969/// Get a mapping function from src into the domain of this type (n_transExt)
     
    6898number   ntGcd(number a, number b, const coeffs cf);
    6999number   ntLcm(number a, number b, const coeffs cf);
    70 number   ntSize(number a, const coeffs cf);
     100int      ntSize(number a, const coeffs cf);
    71101void     ntDelete(number * a, const coeffs cf);
    72102void     ntCoeffWrite(const coeffs cf);
Note: See TracChangeset for help on using the changeset viewer.