Changeset dc79bd in git


Ignore:
Timestamp:
Oct 4, 2012, 8:16:12 PM (12 years ago)
Author:
Oleksandr Motsak <motsak@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '38dfc5131670d387a89455159ed1e071997eec94')
Children:
20c54056198f825ef62888512342886cf96c82a6
Parents:
9f6cc091f81a4d6c75b4696ae3e84bc53e419051
git-author:
Oleksandr Motsak <motsak@mathematik.uni-kl.de>2012-10-04 20:16:12+02:00
git-committer:
Oleksandr Motsak <motsak@mathematik.uni-kl.de>2012-10-05 18:15:05+02:00
Message:
Major Update for Enumerators + Fixes for Algrbraic & Transcendental extensions

General:
chg: cleanup + documentation + additional assumes

Enumerators:
chg: some polish on Enumerators
add: CRecursivePolyCoeffsEnumerator<ConverterPolicy> for recursive treatment of converted coeffs as Enumerators

Coeffs:
chg: use mpz_lcm for readability in nlClearDenominators + cleanup
add: nlClear*NoPositiveLead variants should not make LC positive
chg: all nlClear* are not static in order to be usable from alg / trans exts.
fix: fixed a bug in ndClearContent

AlgExt:
add: nCoeff_is_Q_algext for checking an alg. ext. of Q
add: naClear* for alg. ext. over Q
NOTE: coeffs are polynomials in Q[a] - one should simply consider each of them recursively as a collection of numbers...
NOTE: compute GCDs over Alg. Ext... + gcds of (int.) numbers!?
NOTE: trying to be conform with older Singular: no negative leading coeff. normalization
chg: Alg. Ext: use singclap_*gcd (instead of Frank's gcd-stuff)

p_poly:
add: p_Cleardenom_n/p_Cleardenom also clear content afterwards...
chg: major and minor changes to p_Content/p_Cleardenom_n/p_Cleardenom + cleanup
add: additionally trying to assure positive leading coeff. after p_Content/p_Cleardenom_n(/p_Cleardenom?)
NOTE: which should not be needed as n_ClearDenominators/n_ClearContent are supposed to assure that themselves!
add: more assumes to p_polys.cc
NOTE: massive usage of enumerators form p_* causes problems - only doing that for Q_a()!
NOTE: do -normalization over Q(x...)

TransExt:
add: ntClear* for trans. ext's
fix: correct ntGetDenom/ntGetNumerator (thanks to pSubstPar),
NOTE: no negative denominator out of ntGetNumerator/ntGetDenom!
add: first inefficient ntClearContent/Q  and ntClearDenominators/Q & F_p impl.
NOTE:  careful with the use of nlClear* ! (only over Q!)
add: added ntTest to transext.cc on most in/outs + use ntInit(poly)!

NOTE: trying to fix the monic-poly-gcd problem in ntClearDenominators!
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • Singular/maps_ip.cc

    r9f6cc0 rdc79bd  
    259259
    260260    number d = n_GetDenom(p_GetCoeff(p, currRing), currRing);
    261     if (!n_IsOne (d, currRing->cf)) WarnS("ignoring denominators of coefficients...");
     261    if (!n_IsOne (d, currRing->cf))
     262      WarnS("ignoring denominators of coefficients...");
    262263    n_Delete(&d, currRing);
    263264

  • libpolys/coeffs/Enumerator.h

    r9f6cc0 rdc79bd  
    5454    virtual void Reset() = 0;
    5555
    56     virtual ~IBaseEnumerator() {} // TODO: needed?
    57 
     56    /// Current position is inside the collection (not -1 or past the end)
     57    virtual bool IsValid() const = 0;
     58   
    5859  private:
     60    /// disable copy constructor and assigment operator
    5961    IBaseEnumerator(const IBaseEnumerator&);
    6062    void operator=(const IBaseEnumerator&);
    … …  
    6264  protected:
    6365    IBaseEnumerator(){}
    64 
    65     /// Current position is inside the collection (not -1 or past the end)
    66     virtual bool IsValid() const = 0;
     66    ~IBaseEnumerator() {} // TODO: needed?
     67
     68
    6769};
    6870
    … …  
    9092    virtual const_reference Current() const = 0;
    9193
    92     virtual ~IAccessor() {} // TODO: needed?
     94 protected:
     95    IAccessor(){}   
     96    ~IAccessor() {} // TODO: needed?
    9397 
    9498};

  • libpolys/coeffs/coeffs.h

    r9f6cc0 rdc79bd  
    774774{
    775775  assume(r != NULL);
    776   return ((!nCoeff_is_Ring(r)) && (n_GetChar(r) == 0) && nCoeff_is_Extension(r));
    777 }
     776  return ((n_GetChar(r) == 0) && nCoeff_is_Extension(r));
     777}
     778
     779
     780
    778781
    779782static inline BOOLEAN nCoeff_is_long_R(const coeffs r)
    … …  
    798801static inline BOOLEAN nCoeff_is_algExt(const coeffs r)
    799802{ assume(r != NULL); return (getCoeffType(r)==n_algExt); }
     803
     804/// is it an alg. ext. of Q?
     805static inline BOOLEAN nCoeff_is_Q_algext(const coeffs r)
     806{ assume(r != NULL); return ((n_GetChar(r) == 0) && nCoeff_is_algExt(r)); }
    800807
    801808/// TRUE iff r represents a transcendental extension field

  • libpolys/coeffs/longrat.cc

    r9f6cc0 rdc79bd  
    26412641}
    26422642
    2643 static void nlClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
     2643void nlClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
    26442644{
    26452645  assume(cf != NULL);
    … …  
    26612661  s=2147483647; // max. int
    26622662
    2663  
    2664   int lc_is_pos=nlGreaterZero(numberCollectionEnumerator.Current(),cf);
     2663  const BOOLEAN lc_is_pos=nlGreaterZero(numberCollectionEnumerator.Current(),cf);
    26652664
    26662665  int normalcount = 0;
    … …  
    27352734}
    27362735
    2737 static void nlClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
     2736void nlClearContentNoPositiveLead(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
     2737{
     2738  assume(cf != NULL);
     2739  assume(getCoeffType(cf) == ID);
     2740
     2741  numberCollectionEnumerator.Reset();
     2742
     2743  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
     2744  {
     2745    c = n_Init(1, cf);
     2746    return;
     2747  }
     2748
     2749  // all coeffs are given by integers!!!
     2750
     2751  // part 1, find a small candidate for gcd
     2752  number cand1,cand;
     2753  int s1,s;
     2754  s=2147483647; // max. int
     2755
     2756  const BOOLEAN lc_is_pos = TRUE; // nlGreaterZero(numberCollectionEnumerator.Current(),cf);
     2757
     2758  int normalcount = 0;
     2759  do
     2760  {
     2761    number& n = numberCollectionEnumerator.Current();
     2762    nlNormalize(n, cf); ++normalcount;
     2763    cand1 = n;
     2764
     2765    if (SR_HDL(cand1)&SR_INT) { cand=cand1; break; }
     2766    assume(cand1->s==3); // all coeffs should be integers // ==0?!! after printing
     2767    s1=mpz_size1(cand1->z);
     2768    if (s>s1)
     2769    {
     2770      cand=cand1;
     2771      s=s1;
     2772    }
     2773  } while (numberCollectionEnumerator.MoveNext() );
     2774
     2775//  assume( nlGreaterZero(cand,cf) ); // cand may be a negative integer!
     2776
     2777  cand=nlCopy(cand,cf);
     2778  // part 2: compute gcd(cand,all coeffs)
     2779
     2780  numberCollectionEnumerator.Reset();
     2781
     2782  while (numberCollectionEnumerator.MoveNext() )
     2783  {
     2784    number& n = numberCollectionEnumerator.Current();
     2785
     2786    if( (--normalcount) <= 0)
     2787      nlNormalize(n, cf);
     2788
     2789    nlInpGcd(cand, n, cf);
     2790
     2791    assume( nlGreaterZero(cand,cf) );
     2792
     2793    if(nlIsOne(cand,cf))
     2794    {
     2795      c = cand;
     2796
     2797      if(!lc_is_pos)
     2798      {
     2799        // make the leading coeff positive
     2800        c = nlNeg(c, cf);
     2801        numberCollectionEnumerator.Reset();
     2802
     2803        while (numberCollectionEnumerator.MoveNext() )
     2804        {
     2805          number& nn = numberCollectionEnumerator.Current();
     2806          nn = nlNeg(nn, cf);
     2807        }
     2808      }
     2809      return;
     2810    }
     2811  }
     2812
     2813  // part3: all coeffs = all coeffs / cand
     2814  if (!lc_is_pos)
     2815    cand = nlNeg(cand,cf);
     2816
     2817  c = cand;
     2818  numberCollectionEnumerator.Reset();
     2819
     2820  while (numberCollectionEnumerator.MoveNext() )
     2821  {
     2822    number& n = numberCollectionEnumerator.Current();
     2823    number t=nlIntDiv(n, cand, cf); // simple integer exact division, no ratios to remain
     2824    nlDelete(&n, cf);
     2825    n = t;
     2826  }
     2827}
     2828
     2829void nlClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
    27382830{
    27392831  assume(cf != NULL);
    … …  
    27452837  {
    27462838    c = n_Init(1, cf);
     2839//    assume( n_GreaterZero(c, cf) );
    27472840    return;
    27482841  }
    … …  
    27592852
    27602853  int s=0;
    2761 //  mpz_t tmp; mpz_init(tmp); // tmp = GMP int
    27622854 
    2763   int lc_is_pos=nlGreaterZero(numberCollectionEnumerator.Current(),cf);
     2855  const BOOLEAN lc_is_pos=nlGreaterZero(numberCollectionEnumerator.Current(),cf);
    27642856
    27652857  do
    … …  
    27812873        {
    27822874          mpz_lcm(cand->z, cand->z, cand1->n);
    2783 /*                 
    2784           mpz_gcd(tmp,cand->z,cand1->n); // tmp = GCD( cand->z, cand1->n )
    2785          
    2786           if (mpz_cmp_si(tmp,1)!=0)
    2787             mpz_divexact(cand->z,cand->z,tmp); // cand->z /= tmp
    2788          
    2789           mpz_mul(cand->z,cand->z,cand1->n); // cand->z *= cand1->n
    2790 */
    27912875        }
    27922876      }
    … …  
    28152899      }
    28162900    }
     2901//    assume( n_GreaterZero(c, cf) );
    28172902    return;
    28182903  }
    … …  
    28292914  c = cand;
    28302915 
     2916  while (numberCollectionEnumerator.MoveNext() )
     2917  {
     2918    number &n = numberCollectionEnumerator.Current();
     2919    n_InpMult(n, cand, cf);
     2920  }
     2921
     2922}
     2923
     2924
     2925void nlClearDenominatorsNoPositiveLead(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
     2926{
     2927  assume(cf != NULL);
     2928  assume(getCoeffType(cf) == ID);
     2929
     2930  numberCollectionEnumerator.Reset();
     2931
     2932  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
     2933  {
     2934    c = n_Init(1, cf);
     2935    return;
     2936  }
     2937
     2938  // all coeffs are given by integers after returning from this routine
     2939
     2940  // part 1, collect product of all denominators /gcds
     2941  number cand;
     2942  cand=ALLOC_RNUMBER();
     2943#if defined(LDEBUG)
     2944  cand->debug=123456;
     2945#endif
     2946  cand->s=3;
     2947
     2948  int s=0;
     2949
     2950  const BOOLEAN lc_is_pos = TRUE; // nlGreaterZero(numberCollectionEnumerator.Current(),cf);
     2951
     2952  do
     2953  {
     2954    number& cand1 = numberCollectionEnumerator.Current();
     2955
     2956    if (!(SR_HDL(cand1)&SR_INT))
     2957    {
     2958      nlNormalize(cand1, cf);
     2959      if ((!(SR_HDL(cand1)&SR_INT)) // not a short int
     2960          && (cand1->s==1))             // and is a normalised rational
     2961      {
     2962        if (s==0) // first denom, we meet
     2963        {
     2964          mpz_init_set(cand->z, cand1->n); // cand->z = cand1->n
     2965          s=1;
     2966        }
     2967        else // we have already something
     2968        {
     2969          mpz_lcm(cand->z, cand->z, cand1->n);
     2970        }
     2971      }
     2972    }
     2973  }
     2974  while (numberCollectionEnumerator.MoveNext() );
     2975
     2976
     2977  if (s==0) // nothing to do, all coeffs are already integers
     2978  {
     2979//    mpz_clear(tmp);
     2980    FREE_RNUMBER(cand);
     2981    if (lc_is_pos)
     2982      c=nlInit(1,cf);
     2983    else
     2984    {
     2985      // make the leading coeff positive
     2986      c=nlInit(-1,cf);
     2987
     2988      // TODO: incorporate the following into the loop below?
     2989      numberCollectionEnumerator.Reset();
     2990      while (numberCollectionEnumerator.MoveNext() )
     2991      {
     2992        number& n = numberCollectionEnumerator.Current();
     2993        n = nlNeg(n, cf);
     2994      }
     2995    }
     2996    return;
     2997  }
     2998
     2999  cand = nlShort3(cand);
     3000
     3001  // part2: all coeffs = all coeffs * cand
     3002  // make the lead coeff positive
     3003  numberCollectionEnumerator.Reset();
     3004
     3005  if (!lc_is_pos)
     3006    cand = nlNeg(cand, cf);
     3007
     3008  c = cand;
     3009
    28313010  while (numberCollectionEnumerator.MoveNext() )
    28323011  {

  • libpolys/coeffs/numbers.cc

    r9f6cc0 rdc79bd  
    129129
    130130  // no fractions
    131   assume(!(  nCoeff_is_Q(r) || nCoeff_is_Q_a(r) ));
     131  assume(!(  nCoeff_is_Q(r) ));
    132132  // all coeffs are given by integers!!!
    133133
    … …  
    172172
    173173  assume(!nCoeff_is_Ring(r));
    174   assume(nCoeff_is_Zp(r) || nCoeff_is_numeric(r) || nCoeff_is_GF(r) || nCoeff_is_Zp_a(r));
    175 
    176   c = curr;
    177  
    178   n_Normalize(c, r);
    179 
    180   if (!n_IsOne(c, r))
    181   {   
     174  assume(nCoeff_is_Zp(r) || nCoeff_is_numeric(r) || nCoeff_is_GF(r) || nCoeff_is_Zp_a(r) || nCoeff_is_Q_algext(r));
     175
     176  n_Normalize(curr, r); // Q: good/bad/ugly??
     177
     178  if (!n_IsOne(curr, r))
     179  {
     180    number t = curr; // takes over the curr! note: not a reference!!!
     181
    182182    curr = n_Init(1, r); // ???
    183183   
    184     number inv = n_Invers(c, r);
     184    number inv = n_Invers(t, r);
    185185   
    186186    while( numberCollectionEnumerator.MoveNext() )
    187187    {
    188188      number &n = numberCollectionEnumerator.Current();
    189       n_Normalize(n, r); // ?
    190189      n_InpMult(n, inv, r); // TODO: either this or directly divide!!!?
     190//      n_Normalize(n, r); // ?
    191191    }
    192192   
    193193    n_Delete(&inv, r);
    194   }
     194
     195    c = t;
     196  } else
     197    c = n_Copy(curr, r); // c == 1 and nothing else to do...
    195198}
    196199

  • libpolys/polys/PolyEnumerator.h

    r9f6cc0 rdc79bd  
    3333class CBasePolyEnumerator: public virtual IBaseEnumerator
    3434{
     35  template <class T>
     36  friend class CRecursivePolyCoeffsEnumerator;
    3537  private:
    36     const poly m_poly; ///< essentially immutable original iterable object
     38    poly m_poly; ///< essentially immutable original iterable object
    3739   
    3840    static const spolyrec m_prevposition_struct; ///< tag for "-1" position
    … …  
    4143    poly m_position; ///< current position in the iterable object
    4244
     45  public:
    4346    virtual bool IsValid() const
    4447    {
    … …  
    4649      return (m_position != NULL) && (m_position != &m_prevposition_struct);
    4750    }
    48    
    49   public:
    50     CBasePolyEnumerator(poly p):
     51
     52
     53    /// Reset this polynomial Enumerator with a different input polynomial
     54    void Reset(poly p)
     55    {
     56      m_poly = p;
     57      m_position = const_cast<poly>(&m_prevposition_struct);
     58      assume( !IsValid() );
     59    }
     60
     61    /// This enumerator is an empty polynomial by default
     62    CBasePolyEnumerator(poly p = NULL):
    5163        IBaseEnumerator(), m_poly(p), m_position(const_cast<poly>(&m_prevposition_struct))
    5264    {
    … …  
    6072      assume( !IsValid() );
    6173    }
     74
    6275
    6376    /// Advances the position to the next term of the polynomial.
    … …  
    117130{
    118131  public:
    119     CPolyCoeffsEnumerator(poly p): CBasePolyEnumerator(p) { assume(p != NULL); }
     132    CPolyCoeffsEnumerator(poly p): CBasePolyEnumerator(p) {}
    120133   
    121134    /// Gets the current element in the collection (read and write).
    … …  
    134147};
    135148
     149
     150struct NAConverter
     151{
     152  static inline poly convert(const number& n)
     153  {
     154    // suitable for alg. ext. numbers that are just polys actually
     155    return (poly)n;
     156  }
     157};
     158
     159/// go into polynomials over an alg. extension recursively
     160template <class ConverterPolicy>
     161class CRecursivePolyCoeffsEnumerator: public IPolyCoeffsEnumerator
     162{
     163  private:
     164    IPolyCoeffsEnumerator& m_global_enumerator; ///< iterates the input polynomial
     165    CBasePolyEnumerator m_local_enumerator; ///< iterates the current coeff. of m_global_enumerator
     166
     167  protected:
     168    virtual bool IsValid() const
     169    {
     170      return m_global_enumerator.IsValid() &&  m_local_enumerator.IsValid();
     171    }   
     172   
     173  public:
     174   
     175    /// NOTE: carefull: don't destruct the input enumerator before doing it with this one...
     176    /// this also changes the original IPolyCoeffsEnumerator& itr!
     177    CRecursivePolyCoeffsEnumerator(IPolyCoeffsEnumerator& itr): m_global_enumerator(itr), m_local_enumerator(NULL) {}
     178
     179    virtual bool MoveNext()
     180    {
     181      if( m_local_enumerator.MoveNext() )
     182        return true;
     183
     184      if( !m_global_enumerator.MoveNext() ) // at the end of the main input polynomial?
     185        return false;
     186
     187      // TODO: make the following changeable (metaprogramming: use policy?),
     188      // leave the following as default option...
     189      poly p = ConverterPolicy::convert(m_global_enumerator.Current()); // Assumes that these numbers are just polynomials!
     190      assume( p != NULL );
     191
     192      // the followig actually needs CPolyCoeffsEnumerator
     193      m_local_enumerator.Reset( p ); // -1 position in p :: to be skipped now!
     194
     195      if( m_local_enumerator.MoveNext() ) // should be true
     196        return true;
     197
     198      assume( FALSE ); return MoveNext(); // this should not happen as p should be non-zero, but just in case...
     199    }
     200   
     201    virtual void Reset()
     202    {
     203      m_global_enumerator.Reset();
     204      m_local_enumerator.Reset(NULL);
     205    }
     206
     207    /// Gets the current element in the collection (read and write).
     208    virtual IPolyCoeffsEnumerator::reference Current()
     209    {
     210      assume( IsValid() );
     211      return pGetCoeff(m_local_enumerator.m_position);
     212    }
     213
     214    /// Gets the current element in the collection (read only).
     215    virtual IPolyCoeffsEnumerator::const_reference Current() const
     216    {
     217      assume( IsValid() );
     218      return pGetCoeff(m_local_enumerator.m_position);
     219    }
     220};
     221
     222
    136223#endif
    137224/* #ifndef POLYENUMERATOR_H */

  • libpolys/polys/clapsing.h

    r9f6cc0 rdc79bd  
    2323//#include <kernel/longtrans.h>
    2424
    25 /* destroys f and g */
     25/// destroys f and g
    2626poly singclap_gcd ( poly f, poly g, const ring r );
    2727
    28 /*
     28
    2929// commented out!
    30 napoly singclap_alglcm ( napoly f, napoly g );
    31 void singclap_algdividecontent ( napoly f, napoly g, napoly &ff, napoly &gg );
    32 */
     30// poly singclap_alglcm ( poly f, poly g, const ring r );
     31// void singclap_algdividecontent ( napoly f, napoly g, napoly &ff, napoly &gg );
    3332
    3433poly singclap_resultant ( poly f, poly g , poly x, const ring r);

  • libpolys/polys/ext_fields/algext.cc

    r9f6cc0 rdc79bd  
    4444#include <polys/simpleideals.h>
    4545
     46#include <polys/PolyEnumerator.h>
     47
    4648#ifdef HAVE_FACTORY
     49#include <factory/factory.h>
    4750#include <polys/clapconv.h>
    48 #include <factory/factory.h>
     51#include <polys/clapsing.h>
    4952#endif
    5053
    51 #include "ext_fields/algext.h"
     54   
     55#include <polys/ext_fields/algext.h>
    5256#define TRANSEXT_PRIVATES 1
    53 #include "ext_fields/transext.h"
     57#include <polys/ext_fields/transext.h>
    5458
    5559#ifdef LDEBUG
    … …  
    103107void     naDelete(number *a, const coeffs cf);
    104108void     naCoeffWrite(const coeffs cf, BOOLEAN details);
    105 number   naIntDiv(number a, number b, const coeffs cf);
     109//number   naIntDiv(number a, number b, const coeffs cf);
    106110const char * naRead(const char *s, number *a, const coeffs cf);
    107111
    … …  
    652656  naTest(a); naTest(b);
    653657  if ((a == NULL) && (b == NULL)) WerrorS(nDivBy0);
    654   return (number)p_Gcd((poly)a, (poly)b, naRing);
     658  const ring R = naRing;
     659  return (number) singclap_gcd(p_Copy((poly)a, R), p_Copy((poly)b, R), R);
     660//  return (number)p_Gcd((poly)a, (poly)b, naRing);
    655661}
    656662
    … …  
    660666  if (a == NULL) WerrorS(nDivBy0);
    661667 
    662   poly aFactor = NULL; poly mFactor = NULL;
    663   poly theGcd = p_ExtGcd((poly)a, aFactor, naMinpoly, mFactor, naRing);
     668  poly aFactor = NULL; poly mFactor = NULL; poly theGcd = NULL;
     669// singclap_extgcd!
     670  const BOOLEAN ret = singclap_extgcd ((poly)a, naMinpoly, theGcd, aFactor, mFactor, naRing);
     671
     672  assume( !ret );
     673
     674//  if( ret ) theGcd = p_ExtGcd((poly)a, aFactor, naMinpoly, mFactor, naRing);
    664675 
    665676  naTest((number)theGcd); naTest((number)aFactor); naTest((number)mFactor);
    … …  
    845856}
    846857
    847 static void naClearContent(ICoeffsEnumerator& /*numberCollectionEnumerator*/, number& c, const coeffs cf)
     858
     859static void naClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
    848860{
    849861  assume(cf != NULL);
    850862  assume(getCoeffType(cf) == ID);
    851   assume(nCoeff_is_Q_a(cf)); // only over Q[a]/m(a), while the default impl. is used over Zp[a]/m(a) !
    852   // all coeffs are given by integers!!!
    853 
    854   c = n_Init(1, cf);
    855   assume(FALSE); // TODO: NOT YET IMPLEMENTED!!!
    856 
    857 //   numberCollectionEnumerator.Reset();
    858 //
    859 //   c = numberCollectionEnumerator.Current();
    860 //
    861 //   n_Normalize(c, r);
    862 //
    863 //   if (!n_IsOne(c, r))
    864 //   {   
    865 //     numberCollectionEnumerator.Current() = n_Init(1, r); // ???
    866 //
    867 //     number inv = n_Invers(c, r);
    868 //
    869 //     while( numberCollectionEnumerator.MoveNext() )
    870 //     {
    871 //       number &n = numberCollectionEnumerator.Current();
    872 //       n_Normalize(n, r); // ?
    873 //       n_InpMult(n, inv, r);
    874 //     }
    875 //
    876 //     n_Delete(&inv, r);
    877 //   }
     863  assume(nCoeff_is_Q_algext(cf)); // only over (Q[a]/m(a)), while the default impl. is used over Zp[a]/m(a) !
     864
     865  const ring   R = cf->extRing;
     866  assume(R != NULL);
     867  const coeffs Q = R->cf;
     868  assume(Q != NULL);
     869  assume(nCoeff_is_Q(Q)); 
     870
     871  numberCollectionEnumerator.Reset();
     872
     873  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
     874  {
     875    c = n_Init(1, cf);
     876    return;
     877  }
     878
     879  naTest(numberCollectionEnumerator.Current());
     880
     881  // part 1, find a small candidate for gcd
     882  int s1; int s=2147483647; // max. int
     883
     884  const BOOLEAN lc_is_pos=naGreaterZero(numberCollectionEnumerator.Current(),cf);
     885
     886  int normalcount = 0;
     887 
     888  poly cand1, cand;
     889 
     890  do
     891  {
     892    number& n = numberCollectionEnumerator.Current();
     893    naNormalize(n, cf); ++normalcount;
     894
     895    naTest(n);
     896   
     897    cand1 = (poly)n;
     898
     899    s1 = p_Deg(cand1, R); // naSize?
     900    if (s>s1)
     901    {
     902      cand = cand1;
     903      s = s1;
     904    }
     905  } while (numberCollectionEnumerator.MoveNext() );
     906
     907//  assume( nlGreaterZero(cand,cf) ); // cand may be a negative integer!
     908
     909  cand = p_Copy(cand, R);
     910  // part 2: compute gcd(cand,all coeffs)
     911
     912  numberCollectionEnumerator.Reset();
     913
     914  int length = 0;
     915  while (numberCollectionEnumerator.MoveNext() )
     916  {
     917    number& n = numberCollectionEnumerator.Current();
     918    ++length;
     919   
     920    if( (--normalcount) <= 0)
     921      naNormalize(n, cf);
     922   
     923    naTest(n);
     924
     925//    p_InpGcd(cand, (poly)n, R);
     926   
     927    cand = singclap_gcd(cand, p_Copy((poly)n, R), R);
     928
     929//    cand1 = p_Gcd(cand,(poly)n, R); p_Delete(&cand, R); cand = cand1;
     930   
     931    assume( naGreaterZero((number)cand, cf) ); // ???
     932/*
     933    if(p_IsConstant(cand,R))
     934    {
     935      c = cand;
     936
     937      if(!lc_is_pos)
     938      {
     939        // make the leading coeff positive
     940        c = nlNeg(c, cf);
     941        numberCollectionEnumerator.Reset();
     942
     943        while (numberCollectionEnumerator.MoveNext() )
     944        {
     945          number& nn = numberCollectionEnumerator.Current();
     946          nn = nlNeg(nn, cf);
     947        }
     948      }
     949      return;
     950    }
     951*/
     952   
     953  }
     954
     955  // part3: all coeffs = all coeffs / cand
     956  if (!lc_is_pos)
     957    cand = p_Neg(cand, R);
     958
     959  c = (number)cand; naTest(c); 
     960
     961  poly cInverse = (poly)naInvers(c, cf);
     962  assume(cInverse != NULL); // c is non-zero divisor!?
     963
     964 
     965  numberCollectionEnumerator.Reset();
     966
     967
     968  while (numberCollectionEnumerator.MoveNext() )
     969  {
     970    number& n = numberCollectionEnumerator.Current();
     971
     972    assume( length > 0 );
     973
     974    if( --length > 0 )
     975    {
     976      assume( cInverse != NULL );
     977      n = (number) p_Mult_q(p_Copy(cInverse, R), (poly)n, R);
     978    }
     979    else
     980    {
     981      n = (number) p_Mult_q(cInverse, (poly)n, R);
     982      cInverse = NULL;
     983      assume(length == 0);
     984    }
     985   
     986    definiteReduce((poly &)n, naMinpoly, cf);   
     987  }
     988 
     989  assume(length == 0);
     990  assume(cInverse == NULL); //   p_Delete(&cInverse, R);
     991
     992  // Quick and dirty fix for constant content clearing... !?
     993  CRecursivePolyCoeffsEnumerator<NAConverter> itr(numberCollectionEnumerator); // recursively treat the numbers as polys!
     994
     995  number cc;
     996
     997  extern void nlClearContentNoPositiveLead(ICoeffsEnumerator&, number&, const coeffs);
     998
     999  nlClearContentNoPositiveLead(itr, cc, Q); // TODO: get rid of (-LC) normalization!?
     1000
     1001  // over alg. ext. of Q // takes over the input number
     1002  c = (number) p_Mult_nn( (poly)c, cc, R);
     1003//      p_Mult_q(p_NSet(cc, R), , R);
     1004
     1005  n_Delete(&cc, Q);
     1006
     1007  // TODO: the above is not enough! need GCD's of polynomial coeffs...!
     1008/*
     1009  // old and wrong part of p_Content
     1010    if (rField_is_Q_a(r) && !CLEARENUMERATORS) // should not be used anymore if CLEARENUMERATORS is 1
     1011    {
     1012      // we only need special handling for alg. ext.
     1013      if (getCoeffType(r->cf)==n_algExt)
     1014      {
     1015        number hzz = n_Init(1, r->cf->extRing->cf);
     1016        p=ph;
     1017        while (p!=NULL)
     1018        { // each monom: coeff in Q_a
     1019          poly c_n_n=(poly)pGetCoeff(p);
     1020          poly c_n=c_n_n;
     1021          while (c_n!=NULL)
     1022          { // each monom: coeff in Q
     1023            d=n_Lcm(hzz,pGetCoeff(c_n),r->cf->extRing->cf);
     1024            n_Delete(&hzz,r->cf->extRing->cf);
     1025            hzz=d;
     1026            pIter(c_n);
     1027          }
     1028          pIter(p);
     1029        }
     1030        // hzz contains the 1/lcm of all denominators in c_n_n
     1031        h=n_Invers(hzz,r->cf->extRing->cf);
     1032        n_Delete(&hzz,r->cf->extRing->cf);
     1033        n_Normalize(h,r->cf->extRing->cf);
     1034        if(!n_IsOne(h,r->cf->extRing->cf))
     1035        {
     1036          p=ph;
     1037          while (p!=NULL)
     1038          { // each monom: coeff in Q_a
     1039            poly c_n=(poly)pGetCoeff(p);
     1040            while (c_n!=NULL)
     1041            { // each monom: coeff in Q
     1042              d=n_Mult(h,pGetCoeff(c_n),r->cf->extRing->cf);
     1043              n_Normalize(d,r->cf->extRing->cf);
     1044              n_Delete(&pGetCoeff(c_n),r->cf->extRing->cf);
     1045              pGetCoeff(c_n)=d;
     1046              pIter(c_n);
     1047            }
     1048            pIter(p);
     1049          }
     1050        }
     1051        n_Delete(&h,r->cf->extRing->cf);
     1052      }
     1053    }
     1054*/ 
     1055
     1056 
     1057//  c = n_Init(1, cf); assume(FALSE); // TODO: NOT YET IMPLEMENTED!!!
     1058}
     1059
     1060
     1061static void naClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
     1062{
     1063  assume(cf != NULL);
     1064  assume(getCoeffType(cf) == ID);
     1065  assume(nCoeff_is_Q_algext(cf)); // only over (Q[a]/m(a)), while the default impl. is used over Zp[a]/m(a) !
     1066
     1067  assume(cf->extRing != NULL);
     1068  const coeffs Q = cf->extRing->cf;
     1069  assume(Q != NULL);
     1070  assume(nCoeff_is_Q(Q)); 
     1071  number n;
     1072  CRecursivePolyCoeffsEnumerator<NAConverter> itr(numberCollectionEnumerator); // recursively treat the numbers as polys!
     1073
     1074  extern void nlClearDenominatorsNoPositiveLead(ICoeffsEnumerator&, number&, const coeffs);
     1075
     1076  nlClearDenominatorsNoPositiveLead(itr, n, Q); // this should probably be fine...
     1077  c = (number)p_NSet(n, cf->extRing); // over alg. ext. of Q // takes over the input number
    8781078}
    8791079
    … …  
    9511151  cf->nCoeffIsEqual  = naCoeffIsEqual;
    9521152  cf->cfInvers       = naInvers;
    953   cf->cfIntDiv       = naDiv;
     1153  cf->cfIntDiv       = naDiv; // ???
    9541154#ifdef HAVE_FACTORY
    9551155  cf->convFactoryNSingN=naConvFactoryNSingN;
    … …  
    9631163
    9641164  if( nCoeff_is_Q(R->cf) )
     1165  {
    9651166    cf->cfClearContent = naClearContent;
     1167    cf->cfClearDenominators = naClearDenominators;
     1168  }
    9661169 
    9671170  return FALSE;

  • libpolys/polys/ext_fields/transext.cc

    r9f6cc0 rdc79bd  
    3030*           is being replaced by a quotient of two polynomials over Z, or
    3131*           - if the denominator is a constant - by a polynomial over Q.
     32*
     33*           TODO: the description above needs a major update!!!
    3234*/
    3335#define TRANSEXT_PRIVATES
    … …  
    5456#endif
    5557
    56 #include "ext_fields/transext.h"
    57 #include "prCopy.h"
     58#include <polys/ext_fields/transext.h>
     59#include <polys/prCopy.h>
     60
     61#include <polys/PolyEnumerator.h>
     62
    5863
    5964/* constants for controlling the complexity of numbers */
    … …  
    6267#define BOUND_COMPLEXITY 10   /**< maximum complexity of a number */
    6368
    64 #define NUMIS1(f) (p_IsConstant(NUM(f), cf->extRing) && \
    65                    n_IsOne(p_GetCoeff(NUM(f), cf->extRing), \
    66                            cf->extRing->cf))
    67                    /**< TRUE iff num. represents 1 */
     69
     70static inline BOOLEAN p_IsOne(const poly p, const ring R)
     71{
     72  assume( p_Test(p, R) );
     73  return (p_IsConstant(p, R) && n_IsOne(p_GetCoeff(p, R), R->cf));
     74}
     75
     76/// TRUE iff num. represents 1
     77#define NUMIS1(f) (p_IsOne(NUM(f), cf->extRing))
     78
    6879#define COM(f) f->complexity
    6980
    7081
    7182#ifdef LDEBUG
    72 #define ntTest(a) ntDBTest(a,__FILE__,__LINE__,cf)
     83#define ntTest(a) assume(ntDBTest(a,__FILE__,__LINE__,cf))
    7384BOOLEAN  ntDBTest(number a, const char *f, const int l, const coeffs r);
    7485#else
    … …  
    89100#define ntCoeffs cf->extRing->cf
    90101
     102
     103extern void nlClearContent(ICoeffsEnumerator&, number&, const coeffs);
     104extern void nlClearContentNoPositiveLead(ICoeffsEnumerator&, number&, const coeffs);
     105
     106//extern void nlClearDenominators(ICoeffsEnumerator&, number&, const coeffs);
     107//extern void nlClearDenominatorsNoPositiveLead(ICoeffsEnumerator&, number&, const coeffs);
    91108
    92109
    … …  
    134151{
    135152  assume(getCoeffType(cf) == ID);
    136   fraction t = (fraction)a;
    137   if (IS0(t)) return TRUE;
    138   assume(NUM(t) != NULL);   /**< t != 0 ==> numerator(t) != 0 */
    139   p_Test(NUM(t), ntRing);
    140   if (!DENIS1(t))
    141   {
    142     p_Test(DEN(t), ntRing);
    143     if(p_IsConstant(DEN(t),ntRing) && (n_IsOne(pGetCoeff(DEN(t)),ntRing->cf)))
     153
     154  if (IS0(a)) return TRUE;
     155
     156  const fraction t = (fraction)a;
     157
     158  const poly num = NUM(t);
     159  assume(num != NULL);   /**< t != 0 ==> numerator(t) != 0 */
     160  assume( p_Test(num, ntRing) );
     161
     162  const poly den = DEN(t);
     163 
     164  if (den != NULL) // !DENIS1(f)
     165  {
     166    assume( p_Test(den, ntRing) );
     167     
     168    if(p_IsConstant(den, ntRing) && (n_IsOne(pGetCoeff(den), ntRing->cf)))
    144169    {
    145170      Print("?/1 in %s:%d\n",f,l);
    146        return FALSE;
    147     }
     171      return FALSE;
     172    }
     173       
     174    if( !n_GreaterZero(pGetCoeff(den), ntRing->cf) )
     175    {
     176      Print("negative sign of DEN. of a fraction in %s:%d\n",f,l);
     177      return FALSE;
     178    }
     179   
     180    // test that den is over integers!?
     181     
     182  } else
     183  {  // num != NULL // den == NULL
     184   
     185//    if( COM(t) != 0 )
     186//    {
     187//      Print("?//NULL with non-zero complexity: %d in %s:%d\n", COM(t), f, l);   
     188//      return FALSE;
     189//    }
     190    // test that nume is over integers!?
    148191  }
    149192  return TRUE;
    … …  
    171214BOOLEAN ntIsZero(number a, const coeffs cf)
    172215{
    173   ntTest(a);
     216  ntTest(a); // !!!
    174217  return (IS0(a));
    175218}
    … …  
    177220void ntDelete(number * a, const coeffs cf)
    178221{
     222  ntTest(*a); // !!!
    179223  fraction f = (fraction)(*a);
    180224  if (IS0(f)) return;
    … …  
    187231BOOLEAN ntEqual(number a, number b, const coeffs cf)
    188232{
    189   ntTest(a); ntTest(b);
     233  ntTest(a);
     234  ntTest(b);
    190235
    191236  /// simple tests
    … …  
    230275number ntCopy(number a, const coeffs cf)
    231276{
    232   ntTest(a);
     277  ntTest(a); // !!!
    233278  if (IS0(a)) return NULL;
    234279  fraction f = (fraction)a;
    … …  
    239284  DEN(result) = h;
    240285  COM(result) = COM(f);
     286  ntTest((number)result);
    241287  return (number)result;
    242288}
    243289
     290/// TODO: normalization of a!?
    244291number ntGetNumerator(number &a, const coeffs cf)
    245292{
    246293  ntTest(a);
    247294  definiteGcdCancellation(a, cf, FALSE);
     295 
    248296  if (IS0(a)) return NULL;
     297
    249298  fraction f = (fraction)a;
    250299  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
    251   BOOLEAN denis1= DENIS1 (f);
     300
     301  const BOOLEAN denis1= DENIS1 (f);
     302
    252303  if (getCoeffType (ntCoeffs) == n_Q && !denis1)
    253304    handleNestedFractionsOverQ (f, cf);
    254   NUM (result)= p_Copy (NUM (f), ntRing);
     305
     306  if (getCoeffType (ntCoeffs) == n_Q && denis1)
     307  {
     308    assume( DEN (f) == NULL );
     309   
     310    number g;
     311    // TODO/NOTE: the following should not be necessary (due to
     312    // Hannes!) as NUM (f) should be over Z!!!   
     313    CPolyCoeffsEnumerator itr(NUM(f));
     314
     315
     316    n_ClearDenominators(itr, g, ntRing->cf);
     317//    nlClearDenominators(itr, g, ntRing->cf);
     318
     319    if( !n_GreaterZero(g, ntRing->cf) )
     320    {
     321      NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
     322      g = n_Neg(g, ntRing->cf);
     323    }
     324
     325    // g should be a positive integer now!
     326    assume( n_GreaterZero(g, ntRing->cf) ); 
     327   
     328    if( !n_IsOne(g, ntRing->cf) )
     329    {
     330      DEN (f) = p_NSet(g, ntRing); // update COM(f)???
     331      COM (f) ++;
     332      assume( DEN (f) != NULL );           
     333    }
     334    else
     335      n_Delete(&g, ntRing->cf);
     336
     337    ntTest(a);
     338  }
     339
     340  // Call ntNormalize instead of above?!?
     341 
     342  NUM (result) = p_Copy (NUM (f), ntRing); // ???
    255343  DEN (result) = NULL;
    256344  COM (result) = 0;
    257   if (getCoeffType (ntCoeffs) == n_Q && denis1)
    258   {
    259     if (!p_IsConstant (NUM (result), ntRing) && pNext (NUM(result)) != NULL)
    260       p_Cleardenom (NUM(result), ntRing);
    261     else
    262     {
    263       number g= p_GetAllDenom (NUM (result), ntRing);
    264       NUM (result)= p_Mult_nn (NUM (result), g, ntRing);
    265     }
    266   }
     345 
     346  ntTest((number)result);
    267347  return (number)result;
    268348}
    269349
     350/// TODO: normalization of a!?
    270351number ntGetDenom(number &a, const coeffs cf)
    271352{
    … …  
    273354  definiteGcdCancellation(a, cf, FALSE);
    274355  fraction f = (fraction)a;
     356
    275357  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
    276   number g;
    277   if (IS0(f) || (DENIS1 (f) && getCoeffType (ntCoeffs) != n_Q))
     358  DEN (result)= NULL;
     359  COM (result)= 0;
     360
     361
     362  const BOOLEAN denis1 = DENIS1 (f);
     363 
     364  if( IS0(f) || (denis1 && getCoeffType (ntCoeffs) != n_Q) ) // */1 or 0
    278365  {
    279366    NUM (result)= p_One(ntRing);
    280     DEN (result)= NULL;
    281     COM (result)= 0;
    282   }
    283   else if (DENIS1 (f))
    284   {
    285     poly num= p_Copy (NUM (f), ntRing);
    286     if (!p_IsConstant (num, ntRing) && pNext(num) != NULL)
    287       p_Cleardenom_n (num, ntRing, g);
    288     else
    289       g= p_GetAllDenom (num, ntRing);
    290     result= (fraction) ntSetMap (ntRing->cf, cf) (g, ntRing->cf, cf);
    291   }
    292   else
    293   {
     367    ntTest((number)result);
     368    return (number)result;
     369  }
     370
     371  if (!denis1) // */* / Q
     372  {
     373    assume( DEN (f) != NULL );
     374
    294375    if (getCoeffType (ntCoeffs) == n_Q)
    295376      handleNestedFractionsOverQ (f, cf);
     377
     378    ntTest(a);
     379
     380    if( DEN (f) != NULL ) // is it ?? // 1 now???
     381    {
     382      assume( !p_IsOne(DEN (f), ntRing) );
     383
     384      NUM (result) = p_Copy (DEN (f), ntRing);
     385      ntTest((number)result);
     386      return (number)result;
     387    }
     388//    NUM (result) = p_One(ntRing); // NOTE: just in order to be sure...
     389  }
     390 
     391  // */1 / Q
     392  assume( getCoeffType (ntCoeffs) == n_Q );
     393  assume( DEN (f) == NULL );     
     394   
     395  number g;   
     396//    poly num= p_Copy (NUM (f), ntRing); // ???
     397
     398
     399  // TODO/NOTE: the following should not be necessary (due to
     400  // Hannes!) as NUM (f) should be over Z!!!
     401  CPolyCoeffsEnumerator itr(NUM(f));
     402   
     403  n_ClearDenominators(itr, g, ntRing->cf); // may return -1 :(((
     404//    nlClearDenominators(itr, g, ntRing->cf);
     405
     406   
     407  if( !n_GreaterZero(g, ntRing->cf) )
     408  {
     409//     NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
     410//     g = n_Neg(g, ntRing->cf);
     411    NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
     412    g = n_Neg(g, ntRing->cf);
     413  }
     414
     415  // g should be a positive integer now!
     416  assume( n_GreaterZero(g, ntRing->cf) ); 
     417
     418  if( !n_IsOne(g, ntRing->cf) )
     419  {
     420    assume( n_GreaterZero(g, ntRing->cf) );
     421    assume( !n_IsOne(g, ntRing->cf) );
     422   
     423    DEN (f) = p_NSet(g, ntRing); // update COM(f)???
     424    assume( DEN (f) != NULL );
     425    COM (f) ++;
     426     
    296427    NUM (result)= p_Copy (DEN (f), ntRing);
    297     DEN (result) = NULL;
    298     COM (result) = 0;
    299   }
     428  }
     429  else
     430  { // common denom == 1?
     431    NUM (result)= p_NSet(g, ntRing); // p_Copy (DEN (f), ntRing);
     432//  n_Delete(&g, ntRing->cf);
     433  }     
     434 
     435//    if (!p_IsConstant (num, ntRing) && pNext(num) != NULL)
     436//    else
     437//      g= p_GetAllDenom (num, ntRing);
     438//    result= (fraction) ntSetMap (ntRing->cf, cf) (g, ntRing->cf, cf);
     439
     440  ntTest((number)result);
    300441  return (number)result;
    301442}
    … …  
    303444BOOLEAN ntIsOne(number a, const coeffs cf)
    304445{
    305   ntTest(a);
     446  ntTest(a); // !!! 
    306447  definiteGcdCancellation(a, cf, FALSE);
    307448  fraction f = (fraction)a;
    … …  
    329470    NUM(f) = p_Neg(NUM(f), ntRing);
    330471  }
     472  ntTest(a);
    331473  return a;
    332474}
    … …  
    359501
    360502  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
    361 
    362   NUM(result) = p_NSet(n, A);
    363   DEN(result) = NULL;
     503   
     504  number den = n_GetDenom(n, C);
     505   
     506  assume( n_GreaterZero(den, C) );
     507   
     508  if( n_IsOne(den, C) )
     509  {
     510     NUM(result) = p_NSet(n, A);
     511     DEN(result) = NULL;       
     512     n_Delete(&den, C);
     513  } else
     514  {
     515     DEN(result) = p_NSet(den, A);     
     516     NUM(result) = p_NSet(n_GetNumerator(n, C), A);
     517     n_Delete(&n, C);     
     518  }
     519
    364520  COM(result) = 0;
     521
     522  ntTest((number)result);
     523
    365524  return (number)result;
    366525}
    … …  
    376535    //DEN(result) = NULL; // done by omAlloc0Bin
    377536    //COM(result) = 0; // done by omAlloc0Bin
     537    ntTest((number)result);
    378538    return (number)result;
    379539  }
    380540}
    381541
     542
     543/// takes over p!
    382544number ntInit(poly p, const coeffs cf)
    383545{
    384546  if (p == 0) return NULL;
    385   else
    386   {
    387     fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
    388     NUM(result) = p;
    389     DEN(result) = NULL;
    390     COM(result) = 0;
    391     return (number)result;
    392   }
     547   
     548    number g;
     549    // TODO/NOTE: the following should not be necessary (due to
     550    // Hannes!) as NUM (f) should be over Z!!!   
     551    CPolyCoeffsEnumerator itr(p);
     552
     553    n_ClearDenominators(itr, g, ntRing->cf);
     554//    nlClearDenominators(itr, g, ntRing->cf);
     555
     556    if( !n_GreaterZero(g, ntRing->cf) )
     557    {
     558      p = p_Neg(p, ntRing); // Ugly :(((
     559      g = n_Neg(g, ntRing->cf);
     560    }
     561
     562    // g should be a positive integer now!
     563    assume( n_GreaterZero(g, ntRing->cf) ); 
     564
     565    fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
     566     
     567    if( !n_IsOne(g, ntRing->cf) )
     568    {
     569      DEN (f) = p_NSet(g, ntRing);
     570//      COM (f) ++; // update COM(f)???
     571      assume( DEN (f) != NULL );           
     572    }
     573    else
     574    {
     575      DEN(f) = NULL;     
     576      n_Delete(&g, ntRing->cf);
     577    }   
     578   
     579    NUM(f) = p;
     580    COM(f) = 0;
     581
     582    ntTest((number)f);
     583    return (number)f;
    393584}
    394585
    … …  
    422613BOOLEAN ntGreater(number a, number b, const coeffs cf)
    423614{
    424   ntTest(a); ntTest(b);
     615  ntTest(a);
     616  ntTest(b);
    425617  number aNumCoeff = NULL; int aNumDeg = 0;
    426618  number bNumCoeff = NULL; int bNumDeg = 0;
    … …  
    497689number ntAdd(number a, number b, const coeffs cf)
    498690{
    499   ntTest(a); ntTest(b);
     691  ntTest(a);
     692  ntTest(b);
    500693  if (IS0(a)) return ntCopy(b, cf);
    501694  if (IS0(b)) return ntCopy(a, cf);
    … …  
    525718  COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
    526719  heuristicGcdCancellation((number)result, cf);
     720
     721//  ntTest((number)result);
     722   
    527723  return (number)result;
    528724}
    … …  
    530726number ntSub(number a, number b, const coeffs cf)
    531727{
    532   ntTest(a); ntTest(b);
     728  ntTest(a);
     729  ntTest(b);
    533730  if (IS0(a)) return ntNeg(ntCopy(b, cf), cf);
    534731  if (IS0(b)) return ntCopy(a, cf);
    … …  
    558755  COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
    559756  heuristicGcdCancellation((number)result, cf);
     757//  ntTest((number)result);
    560758  return (number)result;
    561759}
    … …  
    563761number ntMult(number a, number b, const coeffs cf)
    564762{
    565   ntTest(a); ntTest(b);
     763  ntTest(a); // !!!?
     764  ntTest(b); // !!!?
     765 
    566766  if (IS0(a) || IS0(b)) return NULL;
    567767
    … …  
    569769  fraction fb = (fraction)b;
    570770
    571   poly g = p_Copy(NUM(fa), ntRing);
    572   poly h = p_Copy(NUM(fb), ntRing);
    573   g = p_Mult_q(g, h, ntRing);
    574 
    575   if (g == NULL) return NULL;   /* may happen due to zero divisors */
    576 
    577   poly f;
    578   if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
    579   else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
    580   else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
    581   else /* both den's are != 1 */      f = p_Mult_q(p_Copy(DEN(fa), ntRing),
    582                                                    p_Copy(DEN(fb), ntRing),
    583                                                    ntRing);
     771  const poly g = pp_Mult_qq(NUM(fa), NUM(fb), ntRing);
     772
     773  if (g == NULL) return NULL; // may happen due to zero divisors???
    584774
    585775  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
     776
    586777  NUM(result) = g;
    587   DEN(result) = f;
    588   COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
    589   heuristicGcdCancellation((number)result, cf);
     778
     779  const poly da = DEN(fa);
     780  const poly db = DEN(fb);
     781
     782 
     783  if (db == NULL)
     784  {
     785    // b = ? // NULL
     786   
     787    if(da == NULL)
     788    { // both fa && fb are ?? // NULL!
     789      assume (da == NULL && db == NULL);
     790      DEN(result) = NULL;
     791      COM(result) = 0;
     792    }
     793    else
     794    {
     795      assume (da != NULL && db == NULL);
     796      DEN(result) = p_Copy(da, ntRing);
     797      COM(result) = COM(fa) + MULT_COMPLEXITY;
     798      heuristicGcdCancellation((number)result, cf);
     799    }
     800  } else
     801  { // b = ?? / ??
     802    if (da == NULL)
     803    { // a == ? // NULL
     804      assume( db != NULL && da == NULL);
     805      DEN(result) = p_Copy(db, ntRing);
     806      COM(result) = COM(fb) + MULT_COMPLEXITY;
     807      heuristicGcdCancellation((number)result, cf);
     808    }
     809    else /* both den's are != 1 */
     810    {
     811      assume (da != NULL && db != NULL);
     812      DEN(result) = pp_Mult_qq(da, db, ntRing);
     813      COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
     814      heuristicGcdCancellation((number)result, cf);
     815    }
     816  }
     817
     818//  ntTest((number)result);
     819 
    590820  return (number)result;
    591821}
    … …  
    593823number ntDiv(number a, number b, const coeffs cf)
    594824{
    595   ntTest(a); ntTest(b);
     825  ntTest(a);
     826  ntTest(b);
    596827  if (IS0(a)) return NULL;
    597828  if (IS0(b)) WerrorS(nDivBy0);
    … …  
    614845  COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
    615846  heuristicGcdCancellation((number)result, cf);
     847//  ntTest((number)result);
    616848  return (number)result;
    617849}
    … …  
    687919  }
    688920  *b = pow;
     921  ntTest(*b);
    689922}
    690923
    … …  
    8131046    p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
    8141047  }
     1048   
     1049  if( DEN(f) != NULL )
     1050    if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
     1051    {
     1052      NUM(f) = p_Neg(NUM(f), ntRing);
     1053      DEN(f) = p_Neg(DEN(f), ntRing);
     1054    } 
     1055   
     1056  ntTest((number)f); // TODO!
    8151057}
    8161058
    … …  
    8181060void heuristicGcdCancellation(number a, const coeffs cf)
    8191061{
    820   ntTest(a);
     1062//  ntTest(a); // !!!!????
    8211063  if (IS0(a)) return;
    8221064
    8231065  fraction f = (fraction)a;
    824   if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; return; }
     1066  if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; ntTest(a); return; }
     1067
     1068  assume( DEN(f) != NULL );
    8251069
    8261070  /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
    … …  
    8301074    p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
    8311075    COM(f) = 0;
    832     return;
    833   }
    834 
    835   if (COM(f) <= BOUND_COMPLEXITY) return;
    836   else definiteGcdCancellation(a, cf, TRUE);
    837 }
    838 
    839 /* modifies a */
     1076  } else
     1077  {
     1078    if (COM(f) > BOUND_COMPLEXITY)
     1079      definiteGcdCancellation(a, cf, TRUE);
     1080
     1081  // TODO: check if it is enough to put the following into definiteGcdCancellation?!
     1082  if( DEN(f) != NULL )
     1083    if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
     1084    {
     1085      NUM(f) = p_Neg(NUM(f), ntRing);
     1086      DEN(f) = p_Neg(DEN(f), ntRing);
     1087    } 
     1088  }
     1089 
     1090 
     1091  ntTest(a); // !!!!????
     1092}
     1093
     1094/// modifies a
    8401095void definiteGcdCancellation(number a, const coeffs cf,
    8411096                             BOOLEAN simpleTestsHaveAlreadyBeenPerformed)
    8421097{
    843   ntTest(a);
     1098  ntTest(a); // !!!!
    8441099
    8451100  fraction f = (fraction)a;
    … …  
    8561111      p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
    8571112      COM(f) = 0;
     1113      ntTest(a); // !!!!
    8581114      return;
    8591115    }
    … …  
    8611117
    8621118#ifdef HAVE_FACTORY
    863   /* singclap_gcd destroys its arguments; we hence need copies: */
    864   poly pNum = p_Copy(NUM(f), ntRing);
    865   poly pDen = p_Copy(DEN(f), ntRing);
    866 
    8671119  /* Note that, over Q, singclap_gcd will remove the denominators in all
    8681120     rational coefficients of pNum and pDen, before starting to compute
    8691121     the gcd. Thus, we do not need to ensure that the coefficients of
    8701122     pNum and pDen live in Z; they may well be elements of Q\Z. */
    871   poly pGcd = singclap_gcd(pNum, pDen, cf->extRing);
     1123  /* singclap_gcd destroys its arguments; we hence need copies: */
     1124  poly pGcd = singclap_gcd(p_Copy(NUM(f), ntRing), p_Copy(DEN(f), ntRing), cf->extRing);
    8721125  if (p_IsConstant(pGcd, ntRing) &&
    8731126      n_IsOne(p_GetCoeff(pGcd, ntRing), ntCoeffs))
    … …  
    9051158  COM(f) = 0;
    9061159  p_Delete(&pGcd, ntRing);
     1160
     1161  if( DEN(f) != NULL )
     1162    if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
     1163    {
     1164      NUM(f) = p_Neg(NUM(f), ntRing);
     1165      DEN(f) = p_Neg(DEN(f), ntRing);
     1166    }   
    9071167#endif /* HAVE_FACTORY */
     1168   
     1169  ntTest(a); // !!!!
    9081170}
    9091171
    … …  
    9321194    }   
    9331195  }
     1196  ntTest(a); // !!!!
    9341197}
    9351198
    … …  
    9581221    }
    9591222  }
     1223  ntTest(a);
    9601224}
    9611225
    … …  
    9641228  poly p;
    9651229  const char * result = p_Read(s, p, ntRing);
    966   if (p == NULL) { *a = NULL; return result; }
    967   else
    968   {
    969     fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
    970     NUM(f) = p;
    971     DEN(f) = NULL;
    972     COM(f) = 0;
    973     *a = (number)f;
    974     return result;
    975   }
     1230  if (p == NULL) *a = NULL;
     1231  else *a = ntInit(p, cf);
     1232  return result;
    9761233}
    9771234
    … …  
    9791236{
    9801237  definiteGcdCancellation(a, cf, FALSE);
     1238  ntTest(a); // !!!!
    9811239}
    9821240
    … …  
    10061264number ntLcm(number a, number b, const coeffs cf)
    10071265{
    1008   ntTest(a); ntTest(b);
     1266  ntTest(a);
     1267  ntTest(b);
    10091268  fraction fb = (fraction)b;
    10101269  if ((b==NULL)||(DEN(fb)==NULL)) return ntCopy(a,cf);
    … …  
    10261285    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
    10271286    NUM(result) = pp_Mult_qq(NUM(fa),DEN(fb),ntRing);
     1287
     1288    ntTest((number)result); // !!!!
     1289     
    10281290    return (number)result;
    10291291  }
    1030   else
    1031   { /* return pa*pb/gcd */
     1292   
     1293
     1294  /* return pa*pb/gcd */
    10321295    poly newNum = singclap_pdivide(NUM(fa), pGcd, ntRing);
    10331296    p_Delete(&pGcd,ntRing);
    10341297    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
    10351298    NUM(result) = p_Mult_q(p_Copy(DEN(fb),ntRing),newNum,ntRing);
     1299    ntTest((number)result); // !!!!
    10361300    return (number)result;
    1037   }
     1301 
    10381302#else
    10391303  Print("// factory needed: transext.cc:ntLcm\n");
    … …  
    10451309number ntGcd(number a, number b, const coeffs cf)
    10461310{
    1047   ntTest(a); ntTest(b);
     1311  ntTest(a);
     1312  ntTest(b);
    10481313  if (a==NULL) return ntCopy(b,cf);
    10491314  if (b==NULL) return ntCopy(a,cf);
    … …  
    10621327  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
    10631328  NUM(result) = pGcd;
     1329  ntTest((number)result); // !!!!
    10641330  return (number)result;
    10651331#else
    … …  
    11041370    }
    11051371  }
     1372  ntTest(a); // !!!!
    11061373  return numDegree + denDegree + noOfTerms;
    11071374}
    … …  
    11101377{
    11111378  ntTest(a);
    1112   if (IS0(a)) WerrorS(nDivBy0);
     1379  if (IS0(a))
     1380  {   
     1381    WerrorS(nDivBy0);
     1382    return NULL;
     1383  } 
    11131384  fraction f = (fraction)a;
     1385  assume( f != NULL );
     1386
    11141387  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
    1115   poly g;
    1116   if (DENIS1(f)) g = p_One(ntRing);
    1117   else           g = p_Copy(DEN(f), ntRing);
    1118   NUM(result) = g;
    1119   DEN(result) = p_Copy(NUM(f), ntRing);
    1120   COM(result) = COM(f);
     1388
     1389  assume( NUM(f) != NULL );
     1390  const poly den = DEN(f);
     1391 
     1392  if (den == NULL)
     1393    NUM(result) = p_One(ntRing);
     1394  else
     1395    NUM(result) = p_Copy(den, ntRing);
     1396
     1397  if( !NUMIS1(f) )
     1398  {
     1399    DEN(result) = p_Copy(NUM(f), ntRing);
     1400    COM(result) = COM(f);
     1401  }
     1402  else
     1403  {
     1404    DEN(result) = NULL;
     1405    COM(result) = 0;
     1406  }
     1407  ntTest((number)result); // !!!!
    11211408  return (number)result;
    11221409}
    … …  
    11261413{
    11271414  if (n_IsZero(a, src)) return NULL;
     1415  assume(n_Test(a, src));
    11281416  assume(src == dst->extRing->cf);
    1129   poly p = p_Init(dst->extRing);
    1130   number na=n_Copy(a, src);
    1131   n_Normalize(na, src);
    1132   p_SetCoeff0(p, na, dst->extRing);
     1417  return ntInit(p_NSet(n_Copy(a, src), dst->extRing), dst);
     1418}
     1419
     1420/* assumes that src = Z/p, dst = Q(t_1, ..., t_s) */
     1421number ntMapP0(number a, const coeffs src, const coeffs dst)
     1422{
     1423  if (n_IsZero(a, src)) return NULL;
     1424  assume(n_Test(a, src));
     1425  /* mapping via intermediate int: */
     1426  int n = n_Int(a, src);
     1427  number q = n_Init(n, dst->extRing->cf);
     1428  if (n_IsZero(q, dst->extRing->cf))
     1429  {
     1430    n_Delete(&q, dst->extRing->cf);
     1431    return NULL;
     1432  }
     1433  return ntInit(p_NSet(q, dst->extRing), dst);
     1434}
     1435
     1436/* assumes that either src = Q(t_1, ..., t_s), dst = Q(t_1, ..., t_s), or
     1437                       src = Z/p(t_1, ..., t_s), dst = Z/p(t_1, ..., t_s) */
     1438number ntCopyMap(number a, const coeffs cf, const coeffs dst)
     1439{
     1440//  if (n_IsZero(a, cf)) return NULL;
     1441   
     1442  ntTest(a);
     1443
     1444  if (IS0(a)) return NULL;
     1445   
     1446  const ring rSrc = cf->extRing;
     1447  const ring rDst = dst->extRing;
     1448 
     1449  if( rSrc == rDst )
     1450    return ntCopy(a, dst); // USUALLY WRONG!
     1451 
     1452  fraction f = (fraction)a;
     1453  poly g = prCopyR(NUM(f), rSrc, rDst);
     1454   
     1455  poly h = NULL;
     1456 
     1457  if (!DENIS1(f))
     1458     h = prCopyR(DEN(f), rSrc, rDst);
     1459   
     1460  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
     1461   
     1462  NUM(result) = g;
     1463  DEN(result) = h;
     1464  COM(result) = COM(f);
     1465  assume(n_Test((number)result, dst));
     1466  return (number)result; 
     1467}
     1468
     1469number ntCopyAlg(number a, const coeffs cf, const coeffs dst)
     1470{
     1471  assume( n_Test(a, cf) );
     1472  if (n_IsZero(a, cf)) return NULL;
     1473   
     1474  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
     1475  // DEN(f) = NULL; COM(f) = 0;
     1476  NUM(f) = prCopyR((poly)a, cf->extRing, dst->extRing);
     1477  assume(n_Test((number)f, dst));
     1478  return (number)f;
     1479}
     1480
     1481/* assumes that src = Q, dst = Z/p(t_1, ..., t_s) */
     1482number ntMap0P(number a, const coeffs src, const coeffs dst)
     1483{
     1484  assume( n_Test(a, src) );
     1485  if (n_IsZero(a, src)) return NULL;
     1486  int p = rChar(dst->extRing);
     1487  number q = nlModP(a, src, dst->extRing->cf);
     1488
     1489  if (n_IsZero(q, dst->extRing->cf))
     1490  {
     1491    n_Delete(&q, dst->extRing->cf);
     1492    return NULL;
     1493  }
     1494 
     1495  poly g = p_NSet(q, dst->extRing);
     1496
     1497  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
     1498  NUM(f) = g; // DEN(f) = NULL; COM(f) = 0;
     1499  assume(n_Test((number)f, dst));
     1500  return (number)f;
     1501}
     1502
     1503/* assumes that src = Z/p, dst = Z/p(t_1, ..., t_s) */
     1504number ntMapPP(number a, const coeffs src, const coeffs dst)
     1505{
     1506  assume( n_Test(a, src) );
     1507  if (n_IsZero(a, src)) return NULL;
     1508  assume(src == dst->extRing->cf);
     1509  poly p = p_One(dst->extRing);
     1510  p_SetCoeff(p, n_Copy(a, src), dst->extRing);
    11331511  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
    11341512  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
     1513  assume(n_Test((number)f, dst));
    11351514  return (number)f;
    11361515}
    11371516
    1138 /* assumes that src = Z/p, dst = Q(t_1, ..., t_s) */
    1139 number ntMapP0(number a, const coeffs src, const coeffs dst)
    1140 {
     1517/* assumes that src = Z/u, dst = Z/p(t_1, ..., t_s), where u != p */
     1518number ntMapUP(number a, const coeffs src, const coeffs dst)
     1519{
     1520  assume( n_Test(a, src) );
    11411521  if (n_IsZero(a, src)) return NULL;
    11421522  /* mapping via intermediate int: */
    … …  
    11531533  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
    11541534  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
    1155   return (number)f;
    1156 }
    1157 
    1158 /* assumes that either src = Q(t_1, ..., t_s), dst = Q(t_1, ..., t_s), or
    1159                        src = Z/p(t_1, ..., t_s), dst = Z/p(t_1, ..., t_s) */
    1160 number ntCopyMap(number a, const coeffs cf, const coeffs dst)
    1161 {
    1162 //  if (n_IsZero(a, cf)) return NULL;
    1163    
    1164   ntTest(a);
    1165 
    1166   if (IS0(a)) return NULL;
    1167    
    1168   const ring rSrc = cf->extRing;
    1169   const ring rDst = dst->extRing;
    1170  
    1171   if( rSrc == rDst )
    1172     return ntCopy(a, dst); // USUALLY WRONG!
    1173  
    1174   fraction f = (fraction)a;
    1175   poly g = prCopyR(NUM(f), rSrc, rDst);
    1176    
    1177   poly h = NULL;
    1178  
    1179   if (!DENIS1(f))
    1180      h = prCopyR(DEN(f), rSrc, rDst);
    1181    
    1182   fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
    1183    
    1184   NUM(result) = g;
    1185   DEN(result) = h;
    1186   COM(result) = COM(f);
    1187   return (number)result; 
    1188 }
    1189 
    1190 number ntCopyAlg(number a, const coeffs cf, const coeffs dst)
    1191 {
    1192   if (n_IsZero(a, cf)) return NULL;
    1193    
    1194   fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
    1195   // DEN(f) = NULL; COM(f) = 0;
    1196   NUM(f) = prCopyR((poly)a, cf->extRing, dst->extRing);
    1197   return (number)f;
    1198 }
    1199 
    1200 /* assumes that src = Q, dst = Z/p(t_1, ..., t_s) */
    1201 number ntMap0P(number a, const coeffs src, const coeffs dst)
    1202 {
    1203   if (n_IsZero(a, src)) return NULL;
    1204   int p = rChar(dst->extRing);
    1205   number q = nlModP(a, src, dst->extRing->cf);
    1206 
    1207   if (n_IsZero(q, dst->extRing->cf))
    1208   {
    1209     n_Delete(&q, dst->extRing->cf);
    1210     return NULL;
    1211   }
    1212  
    1213   poly g = p_NSet(q, dst->extRing);
    1214 
    1215   fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
    1216   NUM(f) = g; // DEN(f) = NULL; COM(f) = 0;
    1217   return (number)f;
    1218 }
    1219 
    1220 /* assumes that src = Z/p, dst = Z/p(t_1, ..., t_s) */
    1221 number ntMapPP(number a, const coeffs src, const coeffs dst)
    1222 {
    1223   if (n_IsZero(a, src)) return NULL;
    1224   assume(src == dst->extRing->cf);
    1225   poly p = p_One(dst->extRing);
    1226   p_SetCoeff(p, n_Copy(a, src), dst->extRing);
    1227   fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
    1228   NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
    1229   return (number)f;
    1230 }
    1231 
    1232 /* assumes that src = Z/u, dst = Z/p(t_1, ..., t_s), where u != p */
    1233 number ntMapUP(number a, const coeffs src, const coeffs dst)
    1234 {
    1235   if (n_IsZero(a, src)) return NULL;
    1236   /* mapping via intermediate int: */
    1237   int n = n_Int(a, src);
    1238   number q = n_Init(n, dst->extRing->cf);
    1239   poly p;
    1240   if (n_IsZero(q, dst->extRing->cf))
    1241   {
    1242     n_Delete(&q, dst->extRing->cf);
    1243     return NULL;
    1244   }
    1245   p = p_One(dst->extRing);
    1246   p_SetCoeff(p, q, dst->extRing);
    1247   fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
    1248   NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
     1535  assume(n_Test((number)f, dst));
    12491536  return (number)f;
    12501537}
    … …  
    13491636  //DEN(result) = NULL; // done by omAlloc0Bin
    13501637  //COM(result) = 0; // done by omAlloc0Bin
     1638  ntTest((number)result);
    13511639  return (number)result;
    13521640}
    … …  
    13631651static int ntParDeg(number a, const coeffs cf)
    13641652{
     1653  ntTest(a);
    13651654  if (IS0(a)) return -1;
    13661655  fraction fa = (fraction)a;
    … …  
    13861675  COM(f) = 0;
    13871676
     1677  ntTest((number)f);
     1678
    13881679  return (number)f;
    13891680}
    … …  
    13921683int ntIsParam(number m, const coeffs cf)
    13931684{
     1685  ntTest(m);
    13941686  assume(getCoeffType(cf) == ID);
    13951687
    … …  
    14051697}
    14061698
    1407 static void ntClearContent(ICoeffsEnumerator& /*numberCollectionEnumerator*/, number& c, const coeffs cf)
     1699struct NTNumConverter
     1700{
     1701  static inline poly convert(const number& n)
     1702  {
     1703    // suitable for trans. ext. numbers that are fractions of polys
     1704    return NUM((fraction)n); // return the numerator
     1705  }
     1706};
     1707
     1708
     1709static void ntClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
    14081710{
    14091711  assume(cf != NULL);
    14101712  assume(getCoeffType(cf) == ID);
    1411   assume(nCoeff_is_Q_a(cf)); // only over Q(a), while the default impl. is used over Zp(a) !
    1412   // all coeffs are given by integers!!!
    1413 
    1414   c = n_Init(1, cf);
    1415   assume(FALSE); // TODO: NOT YET IMPLEMENTED!!!
    1416 
    1417 //   numberCollectionEnumerator.Reset();
    1418 //
    1419 //   c = numberCollectionEnumerator.Current();
    1420 //
    1421 //   n_Normalize(c, r);
    1422 //
    1423 //   if (!n_IsOne(c, r))
    1424 //   {   
    1425 //     numberCollectionEnumerator.Current() = n_Init(1, r); // ???
    1426 //
    1427 //     number inv = n_Invers(c, r);
    1428 //
    1429 //     while( numberCollectionEnumerator.MoveNext() )
    1430 //     {
    1431 //       number &n = numberCollectionEnumerator.Current();
    1432 //       n_Normalize(n, r); // ?
    1433 //       n_InpMult(n, inv, r);
    1434 //     }
    1435 //
    1436 //     n_Delete(&inv, r);
    1437 //   }
    1438 }
    1439 
    1440 static void ntClearDenominators(ICoeffsEnumerator& /*numberCollectionEnumerator*/, number& c, const coeffs cf)
     1713  // all coeffs are given by fractions of polynomails over integers!!!
     1714  // without denominators!!!
     1715
     1716  const ring   R = cf->extRing;
     1717  assume(R != NULL);
     1718  const coeffs Q = R->cf;
     1719  assume(Q != NULL);
     1720  assume(nCoeff_is_Q(Q)); 
     1721
     1722 
     1723  numberCollectionEnumerator.Reset();
     1724
     1725  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
     1726  {
     1727    c = ntInit(1, cf);
     1728    return;
     1729  }
     1730
     1731  // all coeffs are given by integers after returning from this routine
     1732
     1733  // part 1, collect product of all denominators /gcds
     1734  poly cand = NULL;
     1735
     1736  do
     1737  {
     1738    number &n = numberCollectionEnumerator.Current();
     1739
     1740    ntNormalize(n, cf);
     1741
     1742    fraction f = (fraction)n;
     1743
     1744    assume( f != NULL );
     1745
     1746    const poly den = DEN(f);
     1747
     1748    assume( den == NULL ); // ?? / 1 ?
     1749
     1750    const poly num = NUM(f);
     1751
     1752    if( cand == NULL )
     1753      cand = p_Copy(num, R);
     1754    else
     1755      cand = singclap_gcd(cand, p_Copy(num, R), R); // gcd(cand, num)
     1756
     1757    if( p_IsConstant(cand, R) )
     1758      break;
     1759  }
     1760  while( numberCollectionEnumerator.MoveNext() ) ;
     1761
     1762   
     1763  // part2: all coeffs = all coeffs * cand
     1764  if( cand != NULL )
     1765  {
     1766  if( !p_IsConstant(cand, R) )
     1767  {
     1768    c = ntInit(cand, cf);   
     1769    numberCollectionEnumerator.Reset();
     1770    while (numberCollectionEnumerator.MoveNext() )
     1771    {
     1772      number &n = numberCollectionEnumerator.Current();
     1773      const number t = ntDiv(n, c, cf); // TODO: rewrite!?
     1774      ntDelete(&n, cf);
     1775      n = t;
     1776    }
     1777  } // else NUM (result) = p_One(R);
     1778  else { p_Delete(&cand, R); cand = NULL; }
     1779  }
     1780   
     1781  // Quick and dirty fix for constant content clearing: consider numerators???
     1782  CRecursivePolyCoeffsEnumerator<NTNumConverter> itr(numberCollectionEnumerator); // recursively treat the NUM(numbers) as polys!
     1783  number cc;
     1784   
     1785//  nlClearContentNoPositiveLead(itr, cc, Q); // TODO: get rid of (-LC) normalization!?
     1786  nlClearContent(itr, cc, Q);
     1787  number g = ntInit(p_NSet(cc, R), cf);
     1788   
     1789  if( cand != NULL )
     1790  {
     1791    number gg = ntMult(g, c, cf);
     1792    ntDelete(&g, cf);
     1793    ntDelete(&c, cf); c = gg;
     1794  } else
     1795    c = g;
     1796  ntTest(c);
     1797}
     1798
     1799static void ntClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
    14411800{
    14421801  assume(cf != NULL);
    14431802  assume(getCoeffType(cf) == ID); // both over Q(a) and Zp(a)!
    1444   // all coeffs are given by integers!!!
    1445 
    1446   c = n_Init(1, cf);
    1447   assume(FALSE); // TODO: NOT YET IMPLEMENTED!!!
     1803  // all coeffs are given by fractions of polynomails over integers!!!
     1804
     1805  numberCollectionEnumerator.Reset();
     1806
     1807  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
     1808  {
     1809    c = ntInit(1, cf);
     1810    return;
     1811  }
     1812
     1813  // all coeffs are given by integers after returning from this routine
     1814
     1815  // part 1, collect product of all denominators /gcds
     1816  poly cand = NULL;
     1817
     1818  const ring R = cf->extRing;
     1819  assume(R != NULL);
     1820
     1821  const coeffs Q = R->cf;
     1822  assume(Q != NULL);
     1823//  assume(nCoeff_is_Q(Q)); 
     1824
     1825  do
     1826  {
     1827    number &n = numberCollectionEnumerator.Current();
     1828   
     1829    ntNormalize(n, cf);
     1830
     1831    fraction f = (fraction)n;
     1832
     1833    assume( f != NULL );
     1834
     1835    const poly den = DEN(f);
     1836
     1837    if( den == NULL ) // ?? / 1 ?
     1838      continue;
     1839
     1840    if( cand == NULL )
     1841      cand = p_Copy(den, R);
     1842    else
     1843    {
     1844      // cand === LCM( cand, den )!!!!
     1845      // NOTE: maybe it's better to make the product and clearcontent afterwards!?
     1846      // TODO: move the following to factory?
     1847      poly gcd = singclap_gcd(p_Copy(cand, R), p_Copy(den, R), R); // gcd(cand, den) is monic no mater leading coeffs! :((((
     1848//      assume( n_IsOne(pGetCoeff(gcd), Q) ); // TODO: this may be wrong...
     1849      cand = p_Mult_q(cand, p_Copy(den, R), R); // cand *= den     
     1850      const poly t = singclap_pdivide( cand, gcd, R ); // cand' * den / gcd(cand', den)
     1851      p_Delete(&cand, R);
     1852      p_Delete(&gcd, R);
     1853      cand = t;
     1854    }
     1855  }
     1856  while( numberCollectionEnumerator.MoveNext() );
     1857
     1858  if( cand == NULL )
     1859  {
     1860    c = ntInit(1, cf);
     1861    return;
     1862  } 
     1863
     1864  c = ntInit(cand, cf);
     1865
     1866  numberCollectionEnumerator.Reset();
     1867   
     1868  number d = NULL;
     1869
     1870  while (numberCollectionEnumerator.MoveNext() )
     1871  {
     1872    number &n = numberCollectionEnumerator.Current();
     1873    number t = ntMult(n, c, cf); // TODO: rewrite!?
     1874    ntDelete(&n, cf);
     1875
     1876    ntNormalize(t, cf); // TODO: needed?
     1877    n = t;
     1878   
     1879    fraction f = (fraction)t;
     1880    assume( f != NULL );
     1881
     1882    const poly den = DEN(f);
     1883
     1884    if( den != NULL ) // ?? / ?? ?
     1885    {
     1886      assume( p_IsConstant(den, R) );
     1887      assume( pNext(den) == NULL );
     1888         
     1889      if( d == NULL )
     1890        d = n_Copy(pGetCoeff(den), Q);
     1891      else
     1892      {
     1893        number g = n_Lcm(d, pGetCoeff(den), Q);
     1894        n_Delete(&d, Q); d = g;
     1895      }
     1896    }
     1897  }
     1898   
     1899  if( d != NULL )
     1900  {
     1901    numberCollectionEnumerator.Reset();
     1902    while (numberCollectionEnumerator.MoveNext() )
     1903    {
     1904      number &n = numberCollectionEnumerator.Current();
     1905      fraction f = (fraction)n;
     1906
     1907      assume( f != NULL );
     1908
     1909      const poly den = DEN(f);
     1910
     1911      if( den == NULL ) // ?? / 1 ?
     1912        NUM(f) = p_Mult_nn(NUM(f), d, R);
     1913      else
     1914      {
     1915        assume( p_IsConstant(den, R) );
     1916        assume( pNext(den) == NULL );
     1917         
     1918        number ddd = n_Div(d, pGetCoeff(den), Q); // but be an integer now!!!
     1919        NUM(f) = p_Mult_nn(NUM(f), ddd, R);
     1920        n_Delete(&ddd, Q);
     1921       
     1922        p_Delete(&DEN(f), R);
     1923        DEN(f) = NULL; // TODO: check if this is needed!?
     1924      }     
     1925       
     1926      assume( DEN(f) == NULL );
     1927    }
     1928     
     1929    NUM(c) = p_Mult_nn(NUM(c), d, R);
     1930    n_Delete(&d, Q);
     1931  }
     1932   
     1933   
     1934  ntTest(c);
    14481935}
    14491936

  • libpolys/polys/monomials/p_polys.cc

    r9f6cc0 rdc79bd  
    20232023void p_Content(poly ph, const ring r)
    20242024{
     2025  assume( ph != NULL );
     2026
     2027  assume( r != NULL ); assume( r->cf != NULL ); const coeffs C = r->cf;
     2028
     2029
    20252030#if CLEARENUMERATORS
    2026   if( (ph != NULL) && nCoeff_is_Q(r->cf) )
    2027   {
     2031  if( 0 )
     2032  {
     2033      // experimentall (recursive enumerator treatment) of alg. Ext!
    20282034    CPolyCoeffsEnumerator itr(ph);
    20292035    n_ClearContent(itr, r->cf);
    20302036
    2031     assume( n_GreaterZero(pGetCoeff(ph),r->cf) );
    2032    
     2037    p_Test(ph, r); n_Test(pGetCoeff(ph), C);
     2038    assume(n_GreaterZero(pGetCoeff(ph), C)); // ??
     2039
     2040      // if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
    20332041    return;
    20342042  }
    20352043#endif
    2036 
     2044 
     2045 
    20372046#ifdef HAVE_RINGS
    20382047  if (rField_is_Ring(r))
    20392048  {
    2040     if ((ph!=NULL) && rField_has_Units(r))
     2049    if (rField_has_Units(r))
    20412050    {
    20422051      number k = n_GetUnit(pGetCoeff(ph),r->cf);
    … …  
    20532062          pIter(h);
    20542063        }
    2055         assume( n_GreaterZero(pGetCoeff(ph),r->cf) );
     2064//        assume( n_GreaterZero(pGetCoeff(ph),r->cf) );
     2065//        if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
    20562066      }
    20572067      n_Delete(&k,r->cf);
    … …  
    20702080  else
    20712081  {
     2082    assume( pNext(ph) != NULL );
    20722083#if CLEARENUMERATORS
    2073     if( (ph != NULL) && nCoeff_is_Q(r->cf) )
    2074     {
     2084    if( nCoeff_is_Q(r->cf) || nCoeff_is_Q_a(r->cf) )
     2085    {
     2086      // experimentall (recursive enumerator treatment) of alg. Ext!
    20752087      CPolyCoeffsEnumerator itr(ph);
    20762088      n_ClearContent(itr, r->cf);
    2077       assume( n_GreaterZero(pGetCoeff(ph),r->cf) );
     2089
     2090      p_Test(ph, r); n_Test(pGetCoeff(ph), C);
     2091      assume(n_GreaterZero(pGetCoeff(ph), C)); // ??
     2092     
     2093      // if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
    20782094      return;
    20792095    }
    … …  
    20822098    n_Normalize(pGetCoeff(ph),r->cf);
    20832099    if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
    2084     if (rField_is_Q(r))
     2100    if (rField_is_Q(r)) // should not be used anymore if CLEARENUMERATORS is 1
    20852101    {
    20862102      h=p_InitContent(ph,r);
    … …  
    21322148//    }
    21332149#endif
    2134     if (rField_is_Q_a(r))
     2150    if (rField_is_Q_a(r)) 
    21352151    {
    21362152      // we only need special handling for alg. ext.
    … …  
    24522468poly p_Cleardenom(poly ph, const ring r)
    24532469{
    2454   const coeffs C = r->cf;
     2470  if( ph == NULL )
     2471    return NULL;
     2472
     2473  assume( r != NULL ); assume( r->cf != NULL ); const coeffs C = r->cf;
     2474 
     2475#if CLEARENUMERATORS
     2476  if( 0 )
     2477  {
     2478    CPolyCoeffsEnumerator itr(ph);
     2479
     2480    n_ClearDenominators(itr, C);
     2481
     2482    n_ClearContent(itr, C); // divide out the content
     2483
     2484    p_Test(ph, r); n_Test(pGetCoeff(ph), C);
     2485    assume(n_GreaterZero(pGetCoeff(ph), C)); // ??
     2486//    if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r);
     2487
     2488    return ph;
     2489  }
     2490#endif
    24552491
    24562492  poly start=ph;
    24572493
    2458 #if CLEARENUMERATORS
    2459   if( (ph != NULL) && nCoeff_is_Q(C) )
    2460   {
    2461     CPolyCoeffsEnumerator itr(ph);
    2462     n_ClearDenominators(itr, C);
    2463     n_ClearContent(itr, C); // divide out the content
    2464 
    2465     assume( n_GreaterZero(pGetCoeff(ph),C) );
    2466    
    2467     return start;
    2468   }
    2469 #endif
    2470  
    24712494  number d, h;
    24722495  poly p;
    … …  
    24772500    p_Content(ph,r);
    24782501    assume( n_GreaterZero(pGetCoeff(ph),C) );
     2502    if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r);
    24792503    return start;
    24802504  }
    … …  
    24842508  {
    24852509    assume( n_GreaterZero(pGetCoeff(ph),C) );
     2510    if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r);
    24862511    return start;
    24872512  }
    … …  
    25072532
    25082533    assume( n_GreaterZero(pGetCoeff(ph),C) );
     2534    if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r);
    25092535   
    25102536    return start;
    25112537  }
    25122538
     2539  assume(pNext(p)!=NULL);
     2540
    25132541#if CLEARENUMERATORS
    2514   if( (ph != NULL) && nCoeff_is_Q(C) )
     2542  if( nCoeff_is_Q(C) || nCoeff_is_Q_a(C) )
    25152543  {
    25162544    CPolyCoeffsEnumerator itr(ph);
     2545
    25172546    n_ClearDenominators(itr, C);
     2547
    25182548    n_ClearContent(itr, C); // divide out the content
    25192549
    2520     assume( n_GreaterZero(pGetCoeff(ph),C) );
     2550    p_Test(ph, r); n_Test(pGetCoeff(ph), C);
     2551    assume(n_GreaterZero(pGetCoeff(ph), C)); // ??
     2552//    if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r);
    25212553   
    25222554    return start;
    … …  
    26172649
    26182650  assume( n_GreaterZero(pGetCoeff(ph),C) );
     2651  if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r);
    26192652 
    26202653  return start;
    … …  
    26252658  const coeffs C = r->cf;
    26262659  number d, h;
    2627   poly p;
    2628 
    2629   p = ph;
    2630 
    2631   assume(ph != NULL);
    2632 
    2633 
     2660
     2661  assume( ph != NULL );
     2662
     2663  poly p = ph;
    26342664
    26352665#if CLEARENUMERATORS
    2636   if( (ph != NULL) && nCoeff_is_Q(C) )
     2666  if( 0 )
    26372667  {
    26382668    CPolyCoeffsEnumerator itr(ph);
     2669
    26392670    n_ClearDenominators(itr, d, C); // multiply with common denom. d
    26402671    n_ClearContent(itr, h, C); // divide by the content h
    … …  
    26452676    n_Delete(&h, C);
    26462677
     2678    n_Test(c, C);
     2679
     2680    p_Test(ph, r); n_Test(pGetCoeff(ph), C);
     2681    assume(n_GreaterZero(pGetCoeff(ph), C)); // ??
     2682/*
     2683    if(!n_GreaterZero(pGetCoeff(ph),C))
     2684    {
     2685      ph = p_Neg(ph,r);
     2686      c = n_Neg(c, C);
     2687    }
     2688*/
     2689    return;
     2690  }
     2691#endif
     2692 
     2693 
     2694  if( pNext(p) == NULL )
     2695  {
     2696    c=n_Invers(pGetCoeff(p), C);
     2697    p_SetCoeff(p, n_Init(1, C), r);
     2698
    26472699    assume( n_GreaterZero(pGetCoeff(ph),C) );
     2700    if(!n_GreaterZero(pGetCoeff(ph),C))
     2701    {
     2702      ph = p_Neg(ph,r);
     2703      c = n_Neg(c, C);
     2704    }
    26482705   
    26492706    return;
    26502707  }
    2651 #endif
     2708
     2709  assume( pNext(p) != NULL );
    26522710 
    2653   if(pNext(p)==NULL)
    2654   {
    2655     c=n_Invers(pGetCoeff(p), C);
    2656     p_SetCoeff(p, n_Init(1, C), r);
    2657 
    2658     assume( n_GreaterZero(pGetCoeff(ph),C) );
     2711#if CLEARENUMERATORS
     2712  if( nCoeff_is_Q(C) || nCoeff_is_Q_a(C) )
     2713  {
     2714    CPolyCoeffsEnumerator itr(ph);
    26592715   
    2660     return;
    2661   }
    2662 
    2663 #if CLEARENUMERATORS
    2664   if( (ph != NULL) && nCoeff_is_Q(C) )
    2665   {
    2666     CPolyCoeffsEnumerator itr(ph);
    26672716    n_ClearDenominators(itr, d, C); // multiply with common denom. d
    26682717    n_ClearContent(itr, h, C); // divide by the content h
    … …  
    26732722    n_Delete(&h, C);
    26742723
    2675     assume( n_GreaterZero(pGetCoeff(ph),C) );
    2676    
     2724    n_Test(c, C);
     2725
     2726    p_Test(ph, r); n_Test(pGetCoeff(ph), C);
     2727    assume(n_GreaterZero(pGetCoeff(ph), C)); // ??
     2728/*
     2729    if(!n_GreaterZero(pGetCoeff(ph),C))
     2730    {
     2731      ph = p_Neg(ph,r);
     2732      c = n_Neg(c, C);
     2733    }
     2734*/
    26772735    return;
    26782736  }
    … …  
    27662824
    27672825  assume( n_GreaterZero(pGetCoeff(ph),C) );
     2826  if(!n_GreaterZero(pGetCoeff(ph),C))
     2827  {
     2828    ph = p_Neg(ph,r);
     2829    c = n_Neg(c, C);
     2830  }
     2831 
    27682832}
    27692833
Note: See TracChangeset for help on using the changeset viewer.