Changeset a7ec94 in git


Ignore:
Timestamp:
Apr 3, 2006, 10:10:32 AM (18 years ago)
Author:
Wilfred Pohl <pohl@…>
Branches:
(u'spielwiese', '17f1d200f27c5bd38f5dfc6e8a0879242279d1d8')
Children:
b66fdff25ea59b6eccb8449f7f51cd88ee9c8abd
Parents:
4c1175b6e3158e332e93d0a68e192139454c5616
Message:
new version


git-svn-id: file:///usr/local/Singular/svn/trunk@9052 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • factory/cf_gcd.cc

    r4c1175 ra7ec94  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: cf_gcd.cc,v 1.42 2006-02-20 15:47:16 Singular Exp $ */
     2/* $Id: cf_gcd.cc,v 1.43 2006-04-03 08:10:32 pohl Exp $ */
    33
    44#include <config.h>
    … …  
    1313#include "cf_primes.h"
    1414#include "cf_algorithm.h"
    15 #include "cf_map.h"
    16 #include "sm_sparsemod.h"
    1715#include "fac_util.h"
    1816#include "ftmpl_functions.h"
     17#ifdef HAVE_NTL
     18#undef HAVE_NTL
     19#endif
    1920
    2021#ifdef HAVE_NTL
    2122#include <NTL/ZZX.h>
    2223#include "NTLconvert.h"
    23 bool isPurePoly(const CanonicalForm & f);
     24bool isPurePoly(const CanonicalForm & );
     25static CanonicalForm gcd_univar_ntl0( const CanonicalForm &, const CanonicalForm & );
     26static CanonicalForm gcd_univar_ntlp( const CanonicalForm &, const CanonicalForm & );
    2427#endif
    2528
    26 static CanonicalForm gcd_poly( const CanonicalForm & f, const CanonicalForm& g, const bool modularflag );
    27 static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g );
     29static CanonicalForm gcd_poly( const CanonicalForm &, const CanonicalForm & );
     30static CanonicalForm cf_content ( const CanonicalForm &, const CanonicalForm & );
     31static bool gcd_avoid_mtaildegree ( CanonicalForm &, CanonicalForm &, CanonicalForm & );
     32static void cf_prepgcd( const CanonicalForm &, const CanonicalForm &, int &, int &, int & );
    2833
    2934bool
    … …  
    6570    return gcd( e( F ), e( G ) ).degree() < 1;
    6671}
    67 
    68 //{{{ static CanonicalForm balance ( const CanonicalForm & f, const CanonicalForm & q )
    69 //{{{ docu
    70 //
    71 // balance() - map f from positive to symmetric representation
    72 //   mod q.
    73 //
    74 // This makes sense for univariate polynomials over Z only.
    75 // q should be an integer.
    76 //
    77 // Used by gcd_poly_univar0().
    78 //
    79 //}}}
    80 static CanonicalForm
    81 balance ( const CanonicalForm & f, const CanonicalForm & q )
    82 {
    83     Variable x = f.mvar();
    84     CanonicalForm result = 0, qh = q / 2;
    85     CanonicalForm c;
    86     CFIterator i;
    87     for ( i = f; i.hasTerms(); i++ ) {
    88         c = mod( i.coeff(), q );
    89         if ( c > qh )
    90             result += power( x, i.exp() ) * (c - q);
    91         else
    92             result += power( x, i.exp() ) * c;
    93     }
    94     return result;
    95 }
    96 //}}}
    9772
    9873//{{{ static CanonicalForm icontent ( const CanonicalForm & f, const CanonicalForm & c )
    … …  
    193168//}}}
    194169
     170//{{{ static CanonicalForm balance ( const CanonicalForm & f, const CanonicalForm & q )
     171//{{{ docu
     172//
     173// balance() - map f from positive to symmetric representation
     174//   mod q.
     175//
     176// This makes sense for univariate polynomials over Z only.
     177// q should be an integer.
     178//
     179// Used by gcd_poly_univar0().
     180//
     181//}}}
     182static CanonicalForm
     183balance ( const CanonicalForm & f, const CanonicalForm & q )
     184{
     185    Variable x = f.mvar();
     186    CanonicalForm result = 0, qh = q / 2;
     187    CanonicalForm c;
     188    CFIterator i;
     189    for ( i = f; i.hasTerms(); i++ ) {
     190        c = mod( i.coeff(), q );
     191        if ( c > qh )
     192            result += power( x, i.exp() ) * (c - q);
     193        else
     194            result += power( x, i.exp() ) * c;
     195    }
     196    return result;
     197}
     198//}}}
     199
    195200static CanonicalForm
    196201gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
    197202{
    198 #ifdef HAVE_NTL
    199   if (isOn(SW_USE_NTL_GCD_P) && isPurePoly(F) && isPurePoly(G))
    200   {
    201     if ( getCharacteristic() > 0 )
    202     {
    203       //CanonicalForm cf=F.lc();
    204       //CanonicalForm f=F / cf;
    205       //CanonicalForm cg=G.lc();
    206       //CanonicalForm g= G / cg;
    207       zz_pContext ccc(getCharacteristic());
    208       ccc.restore();
    209       zz_p::init(getCharacteristic());
    210       zz_pX F1=convertFacCF2NTLzzpX(F);
    211       zz_pX G1=convertFacCF2NTLzzpX(G);
    212       zz_pX R=GCD(F1,G1);
    213       return  convertNTLzzpX2CF(R,F.mvar());
    214     }
    215     else
    216     {
    217       CanonicalForm f=F ;
    218       CanonicalForm g=G ;
    219       bool rat=isOn( SW_RATIONAL );
    220       if ( isOn( SW_RATIONAL ) )
    221       {
    222          DEBOUTLN( cerr, "NTL_gcd: ..." );
    223          CanonicalForm cdF = bCommonDen( F );
    224          CanonicalForm cdG = bCommonDen( G );
    225          Off( SW_RATIONAL );
    226          CanonicalForm l = lcm( cdF, cdG );
    227          On( SW_RATIONAL );
    228          f *= l, g *= l;
    229       }
    230       DEBOUTLN( cerr, "NTL_gcd: f=" << f );
    231       DEBOUTLN( cerr, "NTL_gcd: g=" << g );
    232       ZZX F1=convertFacCF2NTLZZX(f);
    233       ZZX G1=convertFacCF2NTLZZX(g);
    234       ZZX R=GCD(F1,G1);
    235       CanonicalForm r=convertNTLZZX2CF(R,F.mvar());
    236       DEBOUTLN( cerr, "NTL_gcd: -> " << r );
    237       if (rat)
    238       {
    239         r /= r.lc();
    240         DEBOUTLN( cerr, "NTL_gcd2: -> " << r );
    241       }
    242       return r;
    243     }
    244   }
    245 #endif
    246203  CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
    247204  int p, i, n;
    … …  
    322279    }
    323280    else
    324       return gcd_poly( F, G, false );
     281      return gcd_poly( F, G );
    325282    DEBOUTLN( cerr, "another try ..." );
    326283  }
    327284}
    328285
    329 CanonicalForm
    330 gcd_poly1( const CanonicalForm & f, const CanonicalForm & g, const bool modularflag )
     286static CanonicalForm
     287gcd_poly_p( const CanonicalForm & f, const CanonicalForm & g )
    331288{
    332289    CanonicalForm pi, pi1;
    333     Variable v = f.mvar();
    334 
    335     if ( f.degree( v ) >= g.degree( v ) )
     290    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
     291    int delta = degree( f ) - degree( g );
     292
     293    if ( delta >= 0 )
    336294    {
    337295        pi = f; pi1 = g;
    … …  
    339297    else
    340298    {
    341         pi = g; pi1 = f;
    342     }
    343     CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
    344     int delta;
     299        pi = g; pi1 = f; delta = -delta;
     300    }
    345301    Ci = content( pi ); Ci1 = content( pi1 );
    346302    pi1 = pi1 / Ci1; pi = pi / Ci;
    347303    C = gcd( Ci, Ci1 );
    348     if ( pi.isUnivariate() && pi1.isUnivariate() )
    349     {
     304    if ( !( pi.isUnivariate() && pi1.isUnivariate() ) )
     305    {
     306      if ( gcd_test_one( pi1, pi, true ) )
     307        return C;
     308    }
    350309#ifdef HAVE_NTL
    351       if (( modularflag) ||
    352       ((isOn(SW_USE_NTL_GCD_P)||isOn(SW_USE_NTL_GCD_0))
    353        && isPurePoly(pi) && isPurePoly(pi1)))
    354          return gcd_poly_univar0(pi, pi1, true) * C;
    355 #else
    356       if ( modularflag)
    357         return gcd_poly_univar0( pi, pi1, true ) * C;
     310    else
     311    {
     312      if ( isOn(SW_USE_NTL_GCD_P) && isPurePoly(pi) && isPurePoly(pi1) )
     313        return gcd_univar_ntlp(pi, pi1 ) * C;
     314    }
    358315#endif
    359     }
    360     else if ( gcd_test_one( pi1, pi, true ) )
    361       return C;
    362     delta = degree( pi, v ) - degree( pi1, v );
     316    Variable v = f.mvar();
    363317    Hi = power( LC( pi1, v ), delta );
    364318    if ( (delta+1) % 2 )
    … …  
    385339}
    386340
    387 //{{{ static CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g, bool modularflag )
     341static CanonicalForm
     342gcd_poly_0( const CanonicalForm & f, const CanonicalForm & g )
     343{
     344    CanonicalForm pi, pi1;
     345    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
     346    int delta = degree( f ) - degree( g );
     347
     348    if ( delta >= 0 )
     349    {
     350        pi = f; pi1 = g;
     351    }
     352    else
     353    {
     354        pi = g; pi1 = f; delta = -delta;
     355    }
     356    Ci = content( pi ); Ci1 = content( pi1 );
     357    pi1 = pi1 / Ci1; pi = pi / Ci;
     358    C = gcd( Ci, Ci1 );
     359    if ( pi.isUnivariate() && pi1.isUnivariate() )
     360    {
     361#ifdef HAVE_NTL
     362        if ( isOn(SW_USE_NTL_GCD_0) && isPurePoly(pi) && isPurePoly(pi1) )
     363            return gcd_univar_ntl0(pi, pi1 ) * C;
     364#endif
     365        return gcd_poly_univar0( pi, pi1, true ) * C;
     366    }
     367    else if ( gcd_test_one( pi1, pi, true ) )
     368      return C;
     369    Variable v = f.mvar();
     370    Hi = power( LC( pi1, v ), delta );
     371    if ( (delta+1) % 2 )
     372        bi = 1;
     373    else
     374        bi = -1;
     375    while ( degree( pi1, v ) > 0 ) {
     376        pi2 = psr( pi, pi1, v );
     377        pi2 = pi2 / bi;
     378        pi = pi1; pi1 = pi2;
     379        if ( degree( pi1, v ) > 0 ) {
     380            delta = degree( pi, v ) - degree( pi1, v );
     381            if ( (delta+1) % 2 )
     382                bi = LC( pi, v ) * power( Hi, delta );
     383            else
     384                bi = -LC( pi, v ) * power( Hi, delta );
     385            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
     386        }
     387    }
     388    if ( degree( pi1, v ) == 0 )
     389        return C;
     390    else
     391        return C * pp( pi );
     392}
     393
     394//{{{ static CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
    388395//{{{ docu
    389396//
    … …  
    394401// characteristic and settings of SW_USE_EZGCD and SW_USE_SPARSEMOD, resp.
    395402//
    396 // modularflag is reached down to gcd_poly1() without change in case of
    397 // zero characteristic.
    398 //
    399403// Used by gcd() and gcd_poly_univar0().
    400404//
    … …  
    404408#endif
    405409static CanonicalForm
    406 gcd_poly ( const CanonicalForm & f, const CanonicalForm & g, const bool modularflag )
    407 {
     410gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
     411{
     412    CanonicalForm fc, gc, d1;
     413    int mp, cc, p1, pe;
     414    mp = f.level()+1;
     415    cf_prepgcd( f, g, cc, p1, pe);
     416    if ( cc != 0 )
     417    {
     418        if ( cc > 0 )
     419        {
     420            fc = replacevar( f, Variable(cc), Variable(mp) );
     421            gc = g;
     422        }
     423        else
     424        {
     425            fc = replacevar( g, Variable(-cc), Variable(mp) );
     426            gc = f;
     427        }
     428        return cf_content( fc, gc );
     429    }
     430// now each appearing variable is in f and g
     431    fc = f;
     432    gc = g;
     433    if( gcd_avoid_mtaildegree ( fc, gc, d1 ) )
     434        return d1;
    408435    if ( getCharacteristic() != 0 )
    409436    {
    410       if ( f.isUnivariate() && g.isUnivariate() )
    411       {
    412 #ifdef HAVE_NTL
    413         if ( (isOn(SW_USE_NTL_GCD_P))
    414         && isPurePoly(f) && isPurePoly(g))
    415            return gcd_poly_univar0(f, g, true);
    416 #else
    417         return gcd_poly_univar0( f, g, true );
    418 #endif
    419       }
    420       // now: f or g is not univariate
    421       if (f.level() != g.level())
    422       {
    423         CFMap M, N;
    424         compress( f, g, M, N );
    425         CanonicalForm fM = M(f);
    426         CanonicalForm gM = M(g);
    427         if ( fM.mvar() != gM.mvar() )
    428         {
    429           if ( fM.mvar() > gM.mvar() )
    430             return N( cf_content( fM, gM ) );
    431           else
    432             return N( cf_content( gM, fM ) );
    433         }
     437        if ( p1 == fc.level() )
     438            fc = gcd_poly_p( fc, gc );
    434439        else
    435           return N( gcd_poly1( fM, gM, false ) );
    436       }
    437       else
    438         return gcd_poly1( f, g, false );
    439     }
    440     else if ( isOn( SW_USE_EZGCD ) && ! ( f.isUnivariate() && g.isUnivariate() ) )
    441     {
    442         CFMap M, N;
    443         compress( f, g, M, N );
    444 #if 0
    445         CanonicalForm r=N( ezgcd( M(f), M(g) ) );
    446         if ((f%r!=0) || (g % r !=0))
    447         {
    448            if (si_factor_reminder)
    449            printf("ezgcd failed, trying gcd_poly1\n");
    450            return gcd_poly1( f, g, modularflag);
     440        {
     441            fc = replacevar( fc, Variable(p1), Variable(mp) );
     442            gc = replacevar( gc, Variable(p1), Variable(mp) );
     443            fc = replacevar( gcd_poly_p( fc, gc ), Variable(mp), Variable(p1) );
     444        }
     445    }
     446    else if ( isOn( SW_USE_EZGCD ) && !f.isUnivariate() )
     447    {
     448        if ( pe == 1 )
     449            fc = ezgcd( fc, gc );
     450        else if ( pe > 0 )// no variable at position 1
     451        {
     452            fc = replacevar( fc, Variable(pe), Variable(1) );
     453            gc = replacevar( gc, Variable(pe), Variable(1) );
     454            fc = replacevar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
    451455        }
    452456        else
    453            return r;
    454 #else
    455          return N( ezgcd( M(f), M(g) ) );
    456 #endif
    457     }
    458     else if ( isOn( SW_USE_SPARSEMOD )
    459     && ! ( f.isUnivariate() && g.isUnivariate() ) )
    460     {
    461 #if 0
    462         CanonicalForm r=sparsemod( f, g );
    463         if ((f%r!=0) || (g % r !=0))
    464         {
    465            if (si_factor_reminder)
    466            printf("sparsemod failed, trying gcd_poly1\n");
    467            return r;
    468            //return gcd_poly1( f, g, modularflag);
    469         }
    470         else
    471           return r;
    472 #else
    473         return sparsemod( f, g );
    474 #endif
    475     }
    476     else
    477     {
    478         return gcd_poly1( f, g, modularflag );
    479     }
     457        {
     458            pe = -pe;
     459            fc = swapvar( fc, Variable(pe), Variable(1) );
     460            gc = swapvar( gc, Variable(pe), Variable(1) );
     461            fc = swapvar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
     462        }
     463    }
     464    else
     465    {
     466        fc = gcd_poly_0( fc, gc );
     467    }
     468    if ( d1.degree() > 0 )
     469        fc *= d1;
     470    return fc;
    480471}
    481472//}}}
    … …  
    498489        CanonicalForm result = g;
    499490        while ( i.hasTerms() && ! result.isOne() ) {
    500             result = gcd( result, i.coeff() );
     491            result = gcd( i.coeff(), result );
    501492            i++;
    502493        }
    … …  
    504495    }
    505496    else
    506         if ( f.sign() < 0 )
    507             return -f;
    508         else
    509             return f;
     497        return abs( f );
    510498}
    511499//}}}
    … …  
    522510content ( const CanonicalForm & f )
    523511{
    524     return cf_content( f, 0 );
     512    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) ) {
     513        CFIterator i = f;
     514        CanonicalForm result = abs( i.coeff() );
     515        i++;
     516        while ( i.hasTerms() && ! result.isOne() ) {
     517            result = gcd( i.coeff(), result );
     518            i++;
     519        }
     520        return result;
     521    }
     522    else
     523        return abs( f );
    525524}
    526525//}}}
    … …  
    594593//}}}
    595594
    596 CanonicalForm
     595static CanonicalForm
    597596gcd ( const CanonicalForm & f, const CanonicalForm & g )
    598597{
    599     if ( f.isZero() )
    600         if ( g.lc().sign() < 0 )
    601             return -g;
     598    bool b = f.isZero();
     599    if ( b || g.isZero() )
     600    {
     601        if ( b )
     602            return abs( g );
    602603        else
    603             return g;
    604     else  if ( g.isZero() )
    605         if ( f.lc().sign() < 0 )
    606             return -f;
    607         else
    608             return f;
    609     else  if ( f.inBaseDomain() )
    610         if ( g.inBaseDomain() )
    611             return bgcd( f, g );
    612         else
    613             return cf_content( g, f );
    614     else  if ( g.inBaseDomain() )
    615         return cf_content( f, g );
    616     else  if ( f.mvar() == g.mvar() )
     604            return abs( f );
     605    }
     606    if ( f.inPolyDomain() || g.inPolyDomain() )
     607    {
     608        if ( f.mvar() != g.mvar() )
     609        {
     610            if ( f.mvar() > g.mvar() )
     611                return cf_content( f, g );
     612            else
     613                return cf_content( g, f );
     614        }
    617615        if ( f.inExtension() && getReduce( f.mvar() ) )
    618616            return 1;
    619         else {
     617        else
     618        {
    620619            if ( divides( f, g ) )
    621                 if ( f.lc().sign() < 0 )
    622                     return -f;
    623                 else
    624                     return f;
     620                return abs( f );
    625621            else  if ( divides( g, f ) )
    626                 if ( g.lc().sign() < 0 )
    627                     return -g;
    628                 else
    629                     return g;
    630             if ( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) {
     622                return abs( g );
     623            if ( !( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) )
     624            {
     625                CanonicalForm d;
     626                do{ d = gcd_poly( f, g ); }
     627                while ((!divides(d,f)) || (!divides(d,g)));
     628                return abs( d );
     629            }
     630            else
     631            {
    631632                CanonicalForm cdF = bCommonDen( f );
    632633                CanonicalForm cdG = bCommonDen( g );
    … …  
    636637                CanonicalForm F = f * l, G = g * l;
    637638                Off( SW_RATIONAL );
    638                 do { l = gcd_poly( F, G, true ); }
     639                do { l = gcd_poly( F, G ); }
    639640                while ((!divides(l,F)) || (!divides(l,G)));
    640641                On( SW_RATIONAL );
    641                 if ( l.lc().sign() < 0 )
    642                     return -l;
    643                 else
    644                     return l;
    645             }
    646             else {
    647                 CanonicalForm d;
    648                 do{ d = gcd_poly( f, g, getCharacteristic()==0 ); }
    649                 while ((!divides(d,f)) || (!divides(d,g)));
    650                 if ( d.lc().sign() < 0 )
    651                     return -d;
    652                 else
    653                     return d;
    654             }
    655         }
    656     else  if ( f.mvar() > g.mvar() )
    657         return cf_content( f, g );
    658     else
    659         return cf_content( g, f );
     642                return abs( l );
     643            }
     644        }
     645    }
     646    if ( f.inBaseDomain() && g.inBaseDomain() )
     647        return bgcd( f, g );
     648    else
     649        return 1;
    660650}
    661651
    … …  
    674664{
    675665    if ( f.isZero() || g.isZero() )
    676         return f;
     666        return 0;
    677667    else
    678668        return ( f / gcd( f, g ) ) * g;
    679669}
    680670//}}}
     671
     672#ifdef HAVE_NTL
     673
     674static CanonicalForm
     675gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G )
     676{
     677    ZZX F1=convertFacCF2NTLZZX(F);
     678    ZZX G1=convertFacCF2NTLZZX(G);
     679    ZZX R=GCD(F1,G1);
     680    return convertNTLZZX2CF(R,F.mvar());
     681}
     682
     683static CanonicalForm
     684gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
     685{
     686    zz_pContext ccc(getCharacteristic());
     687    ccc.restore();
     688    zz_p::init(getCharacteristic());
     689    zz_pX F1=convertFacCF2NTLzzpX(F);
     690    zz_pX G1=convertFacCF2NTLzzpX(G);
     691    zz_pX R=GCD(F1,G1);
     692    return  convertNTLzzpX2CF(R,F.mvar());
     693}
     694
     695#endif
     696
     697static bool
     698gcd_avoid_mtaildegree ( CanonicalForm & f1, CanonicalForm & g1, CanonicalForm & d1 )
     699{
     700    bool rdy = true;
     701    int df = f1.taildegree();
     702    int dg = g1.taildegree();
     703
     704    d1 = d1.genOne();
     705    if ( dg == 0 )
     706    {
     707        if ( df == 0 )
     708            return false;
     709        else
     710        {
     711            if ( f1.degree() == df )
     712                d1 = cf_content( g1, LC( f1 ) );
     713            else
     714            {
     715                f1 /= power( f1.mvar(), df );
     716                rdy = false;
     717            }
     718        }
     719    }
     720    else
     721    {
     722        if ( df == 0)
     723        {
     724            if ( g1.degree() == dg )
     725                d1 = cf_content( f1, LC( g1 ) );
     726            else
     727            {
     728                g1 /= power( g1.mvar(), dg );
     729                rdy = false;
     730            }
     731        }
     732        else
     733        {
     734            if ( df > dg )
     735                d1 = power( f1.mvar(), dg );
     736            else
     737                d1 = power( f1.mvar(), df );
     738            if ( f1.degree() == df )
     739            {
     740                if (g1.degree() == dg)
     741                    d1 *= gcd( LC( f1 ), LC( g1 ) );
     742                else
     743                {
     744                    g1 /= power( g1.mvar(), dg);
     745                    d1 *= cf_content( g1, LC( f1 ) );
     746                }
     747            }
     748            else
     749            {
     750                f1 /= power( f1.mvar(), df );
     751                if ( g1.degree() == dg )
     752                    d1 *= cf_content( f1, LC( g1 ) );
     753                else
     754                {
     755                    g1 /= power( g1.mvar(), dg );
     756                    rdy = false;
     757                }
     758            }
     759        }
     760    }
     761    return rdy;
     762}
     763
     764/*
     765*  compute positions p1 and pe of optimal variables:
     766*    pe is used in "ezgcd" and
     767*    p1 in "gcd_poly1"
     768*/
     769static
     770void optvalues ( const int * df, const int * dg, const int n, int & p1, int &pe )
     771{
     772    int i, o1, oe;
     773    if ( df[n] > dg[n] )
     774    {
     775        o1 = df[n]; oe = dg[n];
     776    }
     777    else
     778    {
     779        o1 = dg[n]; oe = df[n];
     780    }
     781    i = n-1;
     782    while ( i > 0 )
     783    {
     784        if ( df[i] != 0 )
     785        {
     786            if ( df[i] > dg[i] )
     787            {
     788                if ( o1 > df[i]) { o1 = df[i]; p1 = i; }
     789                if ( oe <= dg[i]) { oe = dg[i]; pe = i; }
     790            }
     791            else
     792            {
     793                if ( o1 > dg[i]) { o1 = dg[i]; p1 = i; }
     794                if ( oe <= df[i]) { oe = df[i]; pe = i; }
     795            }
     796        }
     797        i--;
     798    }
     799}
     800
     801/*
     802*  make some changes of variables, see optvalues
     803*/
     804static void
     805cf_prepgcd( const CanonicalForm & f, const CanonicalForm & g, int & cc, int & p1, int &pe )
     806{
     807    int i, k, n;
     808    n = f.level();
     809    cc = 0;
     810    p1 = pe = n;
     811    if ( n == 1 )
     812        return;
     813    int * degsf = new int[n+1];
     814    int * degsg = new int[n+1];
     815    for ( i = n; i > 0; i-- )
     816    {
     817        degsf[i] = degsg[i] = 0;
     818    }
     819    degsf = degrees( f, degsf );
     820    degsg = degrees( g, degsg );
     821
     822    k = 0;
     823    for ( i = n-1; i > 0; i-- )
     824    {
     825        if ( degsf[i] == 0 )
     826        {
     827            if ( degsg[i] != 0 )
     828            {
     829                cc = -i;
     830                break;
     831            }
     832        }
     833        else
     834        {
     835            if ( degsg[i] == 0 )
     836            {
     837                cc = i;
     838                break;
     839            }
     840            else k++;
     841        }
     842    }
     843
     844    if ( ( cc == 0 ) && ( k != 0 ) )
     845        optvalues( degsf, degsg, n, p1, pe );
     846    if ( ( pe != 1 ) && ( degsf[1] != 0 ) )
     847        pe = -pe;
     848   
     849    delete [] degsf;
     850    delete [] degsg;
     851}
Note: See TracChangeset for help on using the changeset viewer.