Changeset 06df101 in git for libpolys


Ignore:
Timestamp:
Jul 6, 2011, 4:41:07 PM (12 years ago)
Author:
Frank Seelisch <seelisch@…>
Branches:
(u'spielwiese', '91e5db82acc17434e4062bcfa44e6efa7d41fd30')
Children:
6c19d854882ab99fd19b9b972a32606106c2f2ee
Parents:
2d3091c7bb239cd4787574560ab475ac88621da7
git-author:
Frank Seelisch <seelisch@mathematik.uni-kl.de>2011-07-06 16:41:07+02:00
git-committer:
Mohamed Barakat <mohamed.barakat@rwth-aachen.de>2011-11-09 12:39:12+01:00
Message:
handled nested fractions (when over Q), found bug in longrat.cc (see coeffs_test.h:507)
Location:
libpolys
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • libpolys/polys/ext_fields/transext.cc

    r2d3091c r06df101  
    77*           transcendental variables t_1, ..., t_s, where s >= 1.
    88*           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->extring.
     9*           are represented as quotients of polynomials living in the
     10*           polynomial ring K[t_1, .., t_s] represented by cf->extring.
     11*
     12*           An element of K(t_1, .., t_s) may have numerous representations,
     13*           due to the possibility of common polynomial factors in the
     14*           numerator and denominator. This problem is handled by a
     15*           cancellation heuristic: Each number "knows" its complexity
     16*           which is 0 if and only if common factors have definitely been
     17*           cancelled, and some positive integer otherwise.
     18*           Each arithmetic operation of two numbers with complexities c1
     19*           and c2 will result in a number of complexity c1 + c2 + some
     20*           penalty (specific for each arithmetic operation; see constants
     21*           in the *.h file). Whenever the resulting complexity exceeds a
     22*           certain threshold (see constant in the *.h file), then the
     23*           cancellation heuristic will call 'factory' to compute the gcd
     24*           and cancel it out in the given number. (This definite cancel-
     25*           lation will also be performed at the beginning of ntWrite,
     26*           ensuring that any output is free of common factors.
     27*           For the special case of K = Q (i.e., when computing over the
     28*           rationals), this definite cancellation procedure will also take
     29*           care of nested fractions: If there are fractional coefficients
     30*           in the numerator or denominator of a number, then this number
     31*           is being replaced by a quotient of two polynomials over Z, or
     32*           - if the denominator is a constant - by a polynomial over Q.
    1133*/
    1234
    … …  
    7193void heuristicGcdCancellation(number a, const coeffs cf);
    7294void definiteGcdCancellation(number a, const coeffs cf,
    73                              BOOLEAN skipSimpleTests);
     95                             BOOLEAN simpleTestsHaveAlreadyBeenPerformed);
     96void handleNestedFractionsOverQ(fraction f, const coeffs cf);
    7497
    7598#ifdef LDEBUG
    … …  
    529552}
    530553
     554/* assumes that cf represents the rationals, i.e. Q, and will only
     555   be called in that case;
     556   assumes furthermore that f != NULL and that the denominator of f != 1;
     557   generally speaking, this method removes denominators in the rational
     558   coefficients of the numerator and denominator of 'a';
     559   more concretely, the following normalizations will be performed,
     560   where t^alpha denotes a monomial in the transcendental variables t_k
     561   (1) if 'a' is of the form
     562          (sum_alpha a_alpha/b_alpha * t^alpha)
     563          -------------------------------------
     564            (sum_beta c_beta/d_beta * t^beta)
     565       with integers a_alpha, b_alpha, c_beta, d_beta, then both the
     566       numerator and the denominator will be multiplied by the LCM of
     567       the b_alpha's and the d_beta's (if this LCM is != 1),
     568   (2) if 'a' is - e.g. after having performed step (1) - of the form
     569          (sum_alpha a_alpha * t^alpha)
     570          -----------------------------
     571            (sum_beta c_beta * t^beta)
     572       with integers a_alpha, c_beta, and with a non-constant denominator,
     573       then both the numerator and the denominator will be divided by the
     574       GCD of the a_alpha's and the c_beta's (if this GCD is != 1),
     575   (3) if 'a' is - e.g. after having performed steps (1) and (2) - of the
     576       form
     577          (sum_alpha a_alpha * t^alpha)
     578          -----------------------------
     579                        c
     580       with integers a_alpha, and c != 1, then 'a' will be replaced by
     581       (sum_alpha a_alpha/c * t^alpha);
     582   this procedure does not alter COM(f) (this has to be done by the
     583   calling procedure);
     584   modifies f */
     585void handleNestedFractionsOverQ(fraction f, const coeffs cf)
     586{
     587  assume(nCoeff_is_Q(ntCoeffs));
     588  assume(!IS0(f));
     589  assume(!DENIS1(f));
     590 
     591  if (!p_IsConstant(DEN(f), ntRing))
     592  { /* step (1); see documentation of this procedure above */
     593    number lcmOfDenominators = n_Init(1, ntCoeffs);
     594    number c; number tmp;
     595    poly p = NUM(f);
     596    /* careful when using n_Lcm!!! It computes the lcm of the numerator
     597       of the 1st argument and the denominator of the 2nd!!! */
     598    while (p != NULL)
     599    {
     600      c = p_GetCoeff(p, ntRing);
     601      tmp = n_Lcm(lcmOfDenominators, c, ntCoeffs);
     602      n_Delete(&lcmOfDenominators, ntCoeffs);
     603      lcmOfDenominators = tmp;
     604      pIter(p);
     605    }
     606    p = DEN(f);
     607    while (p != NULL)
     608    {
     609      c = p_GetCoeff(p, ntRing);
     610      tmp = n_Lcm(lcmOfDenominators, c, ntCoeffs);
     611      n_Delete(&lcmOfDenominators, ntCoeffs);
     612      lcmOfDenominators = tmp;
     613      pIter(p);
     614    }
     615    if (!n_IsOne(lcmOfDenominators, ntCoeffs))
     616    { /* multiply NUM(f) and DEN(f) with lcmOfDenominators */
     617      NUM(f) = p_Mult_nn(NUM(f), lcmOfDenominators, ntRing);
     618      DEN(f) = p_Mult_nn(DEN(f), lcmOfDenominators, ntRing);
     619    }
     620    n_Delete(&lcmOfDenominators, ntCoeffs);
     621    if (!p_IsConstant(DEN(f), ntRing))
     622    { /* step (2); see documentation of this procedure above */
     623      p = NUM(f);
     624      number gcdOfCoefficients = n_Copy(p_GetCoeff(p, ntRing), ntCoeffs);
     625      pIter(p);
     626      while ((p != NULL) && (!n_IsOne(gcdOfCoefficients, ntCoeffs)))
     627      {
     628        c = p_GetCoeff(p, ntRing);
     629        tmp = n_Gcd(c, gcdOfCoefficients, ntCoeffs);
     630        n_Delete(&gcdOfCoefficients, ntCoeffs);
     631        gcdOfCoefficients = tmp;
     632        pIter(p);
     633      }
     634      p = DEN(f);
     635      while ((p != NULL) && (!n_IsOne(gcdOfCoefficients, ntCoeffs)))
     636      {
     637        c = p_GetCoeff(p, ntRing);
     638        tmp = n_Gcd(c, gcdOfCoefficients, ntCoeffs);
     639        n_Delete(&gcdOfCoefficients, ntCoeffs);
     640        gcdOfCoefficients = tmp;
     641        pIter(p);
     642      }
     643      if (!n_IsOne(gcdOfCoefficients, ntCoeffs))
     644      { /* divide NUM(f) and DEN(f) by gcdOfCoefficients */
     645        number inverseOfGcdOfCoefficients = n_Invers(gcdOfCoefficients,
     646                                                     ntCoeffs);
     647        NUM(f) = p_Mult_nn(NUM(f), inverseOfGcdOfCoefficients, ntRing);
     648        DEN(f) = p_Mult_nn(DEN(f), inverseOfGcdOfCoefficients, ntRing);
     649        n_Delete(&inverseOfGcdOfCoefficients, ntCoeffs);
     650      }
     651      n_Delete(&gcdOfCoefficients, ntCoeffs);
     652    }
     653  }
     654  if (p_IsConstant(DEN(f), ntRing) &&
     655      (!n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs)))
     656  { /* step (3); see documentation of this procedure above */
     657    number inverseOfDen = n_Invers(p_GetCoeff(DEN(f), ntRing), ntCoeffs);
     658    NUM(f) = p_Mult_nn(NUM(f), inverseOfDen, ntRing);
     659    n_Delete(&inverseOfDen, ntCoeffs);
     660    p_Delete(&DEN(f), ntRing);
     661    DEN(f) = NULL;
     662  }
     663 
     664  /* Now, due to the above computations, DEN(f) may have become the
     665     1-polynomial which needs to be represented by NULL: */
     666  if ((DEN(f) != NULL) &&
     667      p_IsConstant(DEN(f), ntRing) &&
     668      n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs))
     669  {
     670    p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
     671  }
     672}
     673
    531674/* modifies a */
    532675void heuristicGcdCancellation(number a, const coeffs cf)
    … …  
    553696/* modifies a */
    554697void definiteGcdCancellation(number a, const coeffs cf,
    555                              BOOLEAN skipSimpleTests)
     698                             BOOLEAN simpleTestsHaveAlreadyBeenPerformed)
    556699{
    557700  ntTest(a);
    … …  
    559702  fraction f = (fraction)a;
    560703 
    561   if (!skipSimpleTests)
     704  if (!simpleTestsHaveAlreadyBeenPerformed)
    562705  {
    563706    if (IS0(a)) return;
    … …  
    574717  }
    575718 
    576 #ifdef HAVE_FACTORY
    577   /* singclap_gcd destroys its arguments. We hence need copies: */
     719#ifdef HAVE_FACTORY 
     720  /* singclap_gcd destroys its arguments; we hence need copies: */
    578721  poly pNum = p_Copy(NUM(f), ntRing);
    579722  poly pDen = p_Copy(DEN(f), ntRing);
     723 
     724  /* Note that, over Q, singclap_gcd will remove the denominators in all
     725     rational coefficients of pNum and pDen, before starting to compute
     726     the gcd. Thus, we do not need to ensure that the coefficients of
     727     pNum and pDen live in Z; they may well be elements of Q\Z. */
    580728  poly pGcd = singclap_gcd(pNum, pDen, cf->extRing);
    581729  if (p_IsConstant(pGcd, ntRing) &&
    582730      n_IsOne(p_GetCoeff(pGcd, ntRing), ntCoeffs))
    583   {
    584       /* gcd = 1; nothing to cancel */
     731  { /* gcd = 1; nothing to cancel;
     732       Suppose the given rational function field is over Q. Although the
     733       gcd is 1, we may have produced fractional coefficients in NUM(f),
     734       DEN(f), or both, due to previous arithmetics. The next call will
     735       remove those nested fractions, in case there are any. */
     736    if (nCoeff_is_Q(ntCoeffs)) handleNestedFractionsOverQ(f, cf);
    585737  }
    586738  else
    587   {
    588       poly newNum = singclap_pdivide(NUM(f), pGcd, ntRing);
    589       p_Delete(&NUM(f), ntRing);
    590       NUM(f) = newNum;
    591       poly newDen = singclap_pdivide(DEN(f), pGcd, ntRing);
     739  { /* We divide both NUM(f) and DEN(f) by the gcd which is known
     740       to be != 1. */
     741    poly newNum = singclap_pdivide(NUM(f), pGcd, ntRing);
     742    p_Delete(&NUM(f), ntRing);
     743    NUM(f) = newNum;
     744    poly newDen = singclap_pdivide(DEN(f), pGcd, ntRing);
     745    p_Delete(&DEN(f), ntRing);
     746    DEN(f) = newDen;
     747    if (p_IsConstant(DEN(f), ntRing) &&
     748        n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs))
     749    {
     750      /* DEN(f) = 1 needs to be represented by NULL! */
    592751      p_Delete(&DEN(f), ntRing);
    593       if (p_IsConstant(newDen, ntRing) &&
    594           n_IsOne(p_GetCoeff(newDen, ntRing), ntCoeffs))
    595       {
    596         /* newDen = 1 needs to be represented by NULL */
    597         p_Delete(&newDen, ntRing);
    598         newDen = NULL;
    599       }
    600       DEN(f) = newDen;
     752      newDen = NULL;
     753    }
     754    else
     755    { /* Note that over Q, by cancelling the gcd, we may have produced
     756         fractional coefficients in NUM(f), DEN(f), or both. The next
     757         call will remove those nested fractions, in case there are
     758         any. */
     759      if (nCoeff_is_Q(ntCoeffs)) handleNestedFractionsOverQ(f, cf);
     760    }
    601761  }
    602762  COM(f) = 0;

  • libpolys/polys/ext_fields/transext.h

    r2d3091c r06df101  
    99*           transcendental variables t_1, ..., t_s, where s >= 1.
    1010*           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->extring.
     11*           are represented as quotients of polynomials living in the
     12*           polynomial ring K[t_1, .., t_s] represented by cf->extring.
     13*
     14*           An element of K(t_1, .., t_s) may have numerous representations,
     15*           due to the possibility of common polynomial factors in the
     16*           numerator and denominator. This problem is handled by a
     17*           cancellation heuristic: Each number "knows" its complexity
     18*           which is 0 if and only if common factors have definitely been
     19*           cancelled, and some positive integer otherwise.
     20*           Each arithmetic operation of two numbers with complexities c1
     21*           and c2 will result in a number of complexity c1 + c2 + some
     22*           penalty (specific for each arithmetic operation; see constants
     23*           in the *.h file). Whenever the resulting complexity exceeds a
     24*           certain threshold (see constant in the *.h file), then the
     25*           cancellation heuristic will call 'factory' to compute the gcd
     26*           and cancel it out in the given number. (This definite cancel-
     27*           lation will also be performed at the beginning of ntWrite,
     28*           ensuring that any output is free of common factors.
     29*           For the special case of K = Q (i.e., when computing over the
     30*           rationals), this definite cancellation procedure will also take
     31*           care of nested fractions: If there are fractional coefficients
     32*           in the numerator or denominator of a number, then this number
     33*           is being replaced by a quotient of two polynomials over Z, or
     34*           - if the denominator is a constant - by a polynomial over Q.
    1335*/
    1436

  • libpolys/tests/coeffs_test.h

    r2d3091c r06df101  
    481481     setGMPFloatDigits( 10, 5 ); // Init global variables in mpr_complex.cc for gmp_float's... // Note that this seems also to be required for Z_2^m (and Zn?)!????
    482482     simple(n_long_C);
    483    } 
    484 
    485    
    486    // polynomial rings needed for extentions: n_algExt, n_transExt !
    487    
     483   }
     484   
     485   void test_Q_special()
     486   {
     487     const coeffs cf = nInitChar(n_Q, NULLp);
     488
     489     if (cf == NULLp)
     490       clog << ( "Test: could not get this coeff. domain" );
     491
     492     TS_ASSERT_DIFFERS(cf->cfCoeffWrite, NULLp);
     493
     494     if (cf->cfCoeffWrite != NULL )
     495     {
     496       clog << "Coeff-domain: "  << endl;
     497       n_CoeffWrite(cf); PrintLn();
     498     }
     499     
     500     number q1 = n_Init(21, cf);
     501     number q2 = n_Init(2, cf);
     502     number q3 = n_Div(q1, q2, cf);
     503     number q4 = n_Init(30, cf);
     504     number q5 = n_Mult(q3, q4, cf);
     505     TS_ASSERT(n_Test(q5, cf));
     506     Print("21/2 * 30 = %d\n", n_Int(q5, cf));
     507     TS_ASSERT(n_Test(q5, cf));
     508     n_Delete(&q1, cf);
     509     n_Delete(&q2, cf);
     510     n_Delete(&q3, cf);
     511     n_Delete(&q4, cf);
     512     n_Delete(&q5, cf);
     513   }
    488514};
    489515

  • libpolys/tests/polys_test.h

    r2d3091c r06df101  
    220220  {
    221221    poly t = p_ISet(c, r);
     222    if (exp > 0) { p_SetExp(t, i, exp, r); p_Setm(t, r); }
     223    p = p_Add_q(p, t, r);
     224  }
     225  /* assumes that r is over Q;
     226     replaces p by p + c1 / c2 * var(i)^exp (in-place);
     227     expects exp >= 0, and c2 != 0 */
     228  void plusTermOverQ(poly &p, int c1, int c2, int i, int exp, const ring r)
     229  {
     230    number c1AsN = n_Init(c1, r->cf);
     231    number c2AsN = n_Init(c2, r->cf);
     232    number c =  n_Div(c1AsN, c2AsN, r->cf);
     233    poly t = p_ISet(1, r); p_SetCoeff(t, c, r);
    222234    if (exp > 0) { p_SetExp(t, i, exp, r); p_Setm(t, r); }
    223235    p = p_Add_q(p, t, r);
    … …  
    24932505    poly theProduct = p_Mult_q(entry, factor, s);
    24942506    p_Write(theProduct, s);
    2495     clog << "ending multiplication..." << endl;
     2507    clog << "...ending multiplication" << endl;
    24962508    n_Delete(&qfactorAsN, cf); p_Delete(&theProduct, s);
    24972509   
    … …  
    25082520    /* The following multiplication + output of the product is very slow
    25092521       in the svn/trunk SINGULAR version; see trac ticket #308.
    2510        Here, in the Spielwiese, the result is instantaneous! */
     2522       Here, in the Spielwiese, the result is instantaneous. */
    25112523    theProduct = p_Mult_q(entry, factor, s);
    25122524    p_Write(theProduct, s);
    2513     clog << "ending very special multiplication..." << endl;
     2525    clog << "...ending very special multiplication" << endl;
    25142526    n_Delete(&qfactorAsN, cf); p_Delete(&theProduct, s);
     2527       
     2528    rDelete(s); // kills 'cf' and 'r' as well
     2529  }
     2530  void test_Q_Ext_s_t_NestedFractions()
     2531  {
     2532    clog << "Start by creating Q[s, t]..." << endl;
     2533
     2534    char* n[] = {"s", "t"};
     2535    ring r = rDefault( 0, 2, n);   // Q[s, t]
     2536    TS_ASSERT_DIFFERS( r, NULLp );
     2537
     2538    PrintRing(r);
     2539
     2540    TS_ASSERT( rField_is_Domain(r) );
     2541    TS_ASSERT( rField_is_Q(r) );
     2542
     2543    TS_ASSERT( !rField_is_Zp(r) );
     2544    TS_ASSERT( !rField_is_Zp(r, 17) );
     2545
     2546    TS_ASSERT_EQUALS( rVar(r), 2);
     2547
     2548    n_coeffType type = nRegister(n_transExt, ntInitChar);
     2549    TS_ASSERT(type == n_transExt);
     2550
     2551    TransExtInfo extParam;
     2552    extParam.r = r;
     2553
     2554    clog << "Next create the rational function field Q(s, t)..." << endl;
     2555
     2556    const coeffs cf = nInitChar(type, &extParam);   // Q(s, t)
     2557
     2558    if( cf == NULL )
     2559      TS_FAIL("Could not get needed coeff. domain");
     2560
     2561    TS_ASSERT_DIFFERS( cf->cfCoeffWrite, NULLp );
     2562
     2563    if( cf->cfCoeffWrite != NULL )
     2564    {
     2565      clog << "Coeff-domain: "  << endl;
     2566      n_CoeffWrite(cf); PrintLn();
     2567    }
     2568
     2569    TS_ASSERT( !nCoeff_is_algExt(cf) );
     2570    TS_ASSERT( nCoeff_is_transExt(cf) );
     2571
     2572    clog << "Finally create the polynomial ring Q(s, t)[x, y, z]..."
     2573         << endl;
     2574
     2575    char* m[] = {"x", "y", "z"};
     2576    ring s = rDefault(cf, 3, m);   // Q(s, t)[x, y, z]
     2577    TS_ASSERT_DIFFERS(s, NULLp);
     2578
     2579    PrintRing(s);
     2580
     2581    TS_ASSERT( rField_is_Domain(s) );
     2582    TS_ASSERT( !rField_is_Q(s) );
     2583    TS_ASSERT( !rField_is_Zp(s) );
     2584    TS_ASSERT( !rField_is_Zp(s, 11) );
     2585    TS_ASSERT( !rField_is_Zp(s, 17) );
     2586    TS_ASSERT( !rField_is_GF(s) );
     2587    TS_ASSERT( rField_is_Extension(s) );
     2588    TS_ASSERT( !rField_is_GF(s, 25) );
     2589    TS_ASSERT_EQUALS(rVar(s), 3);
     2590   
     2591    /* test 1 for nested fractions, i.e. fractional coefficients: */
     2592    poly v1 = NULL;
     2593    plusTermOverQ(v1, 21, 2, 1, 1, cf->extRing);       // 21/2*s
     2594    plusTermOverQ(v1, 14, 3, 1, 0, cf->extRing);       // 21/2*s + 14/3
     2595    number v1_n = toFractionNumber(v1, cf);
     2596    PrintSized(v1_n, cf);
     2597    poly v2 = NULL;
     2598    plusTermOverQ(v2, 7, 5, 1, 1, cf->extRing);       // 7/5*s
     2599    plusTermOverQ(v2, -49, 6, 2, 1, cf->extRing);     // 7/5*s - 49/6*t
     2600    number v2_n = toFractionNumber(v2, cf);
     2601    PrintSized(v2_n, cf);
     2602    number v3_n = n_Div(v1_n, v2_n, cf);   // (45*s + 20) / (6s - 35*t)
     2603    PrintSized(v3_n, cf);
     2604    n_Delete(&v1_n, cf); n_Delete(&v2_n, cf); n_Delete(&v3_n, cf);
     2605   
     2606    /* test 2 for nested fractions, i.e. fractional coefficients: */
     2607    v1 = NULL;
     2608    plusTermOverQ(v1, 1, 2, 1, 1, cf->extRing);       // 1/2*s
     2609    plusTermOverQ(v1, 1, 1, 1, 0, cf->extRing);       // 1/2*s + 1
     2610    v2 = NULL;
     2611    plusTermOverQ(v2, 1, 1, 1, 1, cf->extRing);       // s
     2612    plusTermOverQ(v2, 2, 3, 2, 1, cf->extRing);       // s + 2/3*t
     2613    poly v3 = p_Mult_q(v1, v2, cf->extRing);   // (1/2*s + 1) * (s + 2/3*t)
     2614    number v_n = toFractionNumber(v3, cf);
     2615    PrintSized(v_n, cf);
     2616    poly w1 = NULL;
     2617    plusTermOverQ(w1, 1, 2, 1, 1, cf->extRing);       // 1/2*s
     2618    plusTermOverQ(w1, 1, 1, 1, 0, cf->extRing);       // 1/2*s + 1
     2619    poly w2 = NULL;
     2620    plusTermOverQ(w2, -7, 5, 1, 0, cf->extRing);      // -7/5
     2621    poly w3 = p_Mult_q(w1, w2, cf->extRing);   // (1/2*s + 1) * (-7/5)
     2622    number w_n = toFractionNumber(w3, cf);
     2623    PrintSized(w_n, cf);
     2624    number z_n = n_Div(v_n, w_n, cf);          // -5/7*s - 10/21*t
     2625    PrintSized(z_n, cf);
     2626    n_Delete(&v_n, cf); n_Delete(&w_n, cf); n_Delete(&z_n, cf);
    25152627   
    25162628    rDelete(s); // kills 'cf' and 'r' as well
Note: See TracChangeset for help on using the changeset viewer.