Changeset b809a8 in git


Ignore:
Timestamp:
Dec 7, 2007, 12:26:57 PM (16 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '2a584933abf2a2d3082034c7586d38bb6de1a30a')
Children:
e16099450a756157fc589b91903dcbe96bcb93aa
Parents:
6fd3a2f8d7423912028ac16c59ca32f3c02bc7c7
Message:
*hannes: gcd_poly


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

Legend:

Unmodified
Added
Removed

  • factory/canonicalform.h

    r6fd3a2 rb809a8  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: canonicalform.h,v 1.33 2006-05-26 11:50:37 Singular Exp $ */
     2/* $Id: canonicalform.h,v 1.34 2007-12-07 11:23:11 Singular Exp $ */
    33
    44#ifndef INCL_CANONICALFORM_H
    … …  
    217217CanonicalForm gcd ( const CanonicalForm&, const CanonicalForm& );
    218218
     219CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g );
     220
    219221CanonicalForm extgcd ( const CanonicalForm&, const CanonicalForm&, CanonicalForm&, CanonicalForm& );
    220222

  • factory/cf_gcd.cc

    r6fd3a2 rb809a8  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: cf_gcd.cc,v 1.56 2007-11-22 09:57:28 Singular Exp $ */
     2/* $Id: cf_gcd.cc,v 1.57 2007-12-07 11:26:57 Singular Exp $ */
    33
    44#include <config.h>
    … …  
    2525#endif
    2626
    27 static CanonicalForm gcd_poly( const CanonicalForm &, const CanonicalForm & );
    2827static CanonicalForm cf_content ( const CanonicalForm &, const CanonicalForm & );
    2928static bool gcd_avoid_mtaildegree ( CanonicalForm &, CanonicalForm &, CanonicalForm & );
    … …  
    450449}
    451450
    452 //{{{ static CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
     451//{{{ CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
    453452//{{{ docu
    454453//
    … …  
    465464int si_factor_reminder=1;
    466465#endif
    467 static CanonicalForm
    468 gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
     466CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
    469467{
    470468  CanonicalForm fc, gc, d1;
    … …  
    10961094  }
    10971095}
    1098 /* ------------------------------------------------------------------------ */
    1099 /* forward declarations: fin_ezgcd stuff*/
    1100 static CanonicalForm fin_ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, FFREvaluation & b, bool internal );
    1101 static bool fin_findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, FFREvaluation & b, int delta, int degF, int degG );
    1102 static CanonicalForm fin_ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, FFREvaluation & b, const modpk & bound );
    1103 
    1104 CanonicalForm fin_ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG )
    1105 {
    1106     FFREvaluation b;
    1107     return fin_ezgcd( FF, GG, b, false );
    1108 }
    1109 
    1110 static CanonicalForm fin_ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, FFREvaluation & b, bool internal )
    1111 {
    1112     CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, lcD;
    1113     CFArray DD( 1, 2 ), lcDD( 1, 2 );
    1114     int degF, degG, delta, t;
    1115     FFREvaluation bt;
    1116     bool gcdfound = false;
    1117     Variable x = Variable(1);
    1118     modpk p = modpk(getCharacteristic(), 1); // k = 1
    1119     DEBINCLEVEL( cerr, "fin_ezgcd" );
    1120     DEBOUTLN( cerr, "FF = " << FF );
    1121     DEBOUTLN( cerr, "GG = " << GG );
    1122     f = content( FF, x ); g = content( GG, x ); d = gcd( f, g );
    1123     DEBOUTLN( cerr, "f = " << f );
    1124     DEBOUTLN( cerr, "g = " << g );
    1125     F = FF / f; G = GG / g;
    1126     if ( F.isUnivariate() && G.isUnivariate() )
    1127     {
    1128         DEBDECLEVEL( cerr, "fin_ezgcd" );
    1129         if(F.mvar()==G.mvar())
    1130           d*=gcd(F,G);
    1131         return d;
    1132     }
    1133     else  if ( gcd_test_one( F, G, false ) )
    1134     {
    1135         DEBDECLEVEL( cerr, "fin_ezgcd" );
    1136         return d;
    1137     }
    1138     lcF = LC( F, x ); lcG = LC( G, x );
    1139     lcD = gcd( lcF, lcG );
    1140     delta = 0;
    1141     degF = degree( F, x ); degG = degree( G, x );
    1142     t = tmax( F.level(), G.level() );
    1143     if ( ! internal )
    1144     {
    1145         b = FFREvaluation( 2, t, FFRandom() );
    1146         b.init(); // choose random point
    1147     }
    1148     while ( ! gcdfound ) {
    1149         /// ---> A2
    1150         DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
    1151         DEBOUTLN( cerr, "F = " << F );
    1152         DEBOUTLN( cerr, "G = " << G );
    1153         if( ! fin_findeval( F, G, Fb, Gb, Db, b, delta, degF, degG ) )
    1154         {
    1155            Off(SW_USE_EZGCD_P);
    1156            d *= gcd_poly_p(F,G);
    1157            On(SW_USE_EZGCD_P);
    1158            return d;
    1159         }
    1160 
    1161         DEBOUTLN( cerr, "found evaluation b = " << b );
    1162         DEBOUTLN( cerr, "F(b) = " << Fb );
    1163         DEBOUTLN( cerr, "G(b) = " << Gb );
    1164         DEBOUTLN( cerr, "D(b) = " << Db );
    1165         delta = degree( Db );
    1166         /// ---> A3
    1167         if ( delta == 0 )
    1168         {
    1169           DEBDECLEVEL( cerr, "fin_ezgcd" );
    1170           return d;
    1171         }
    1172         /// ---> A4
    1173         //deltaold = delta;
    1174         while ( 1 ) {
    1175             bt = b;
    1176             if( ! fin_findeval( F, G, Fbt, Gbt, Dbt, bt, delta + 1, degF, degG ) )
    1177             {
    1178               Off(SW_USE_EZGCD_P);
    1179               d *= gcd_poly_p(F,G);
    1180               On(SW_USE_EZGCD_P);
    1181               return d;
    1182             }
    1183 
    1184             int dd=degree( Dbt );
    1185             if ( dd /*degree( Dbt )*/ == 0 )
    1186             {
    1187                 DEBDECLEVEL( cerr, "fin_ezgcd" );
    1188                 return d;
    1189             }
    1190             if ( dd /*degree( Dbt )*/ == delta )
    1191                 break;
    1192             else  if ( dd /*degree( Dbt )*/ < delta ) {
    1193                 delta = dd /*degree( Dbt )*/;
    1194                 b = bt;
    1195                 Db = Dbt; Fb = Fbt; Gb = Gbt;
    1196             }
    1197         }
    1198         DEBOUTLN( cerr, "now after A4, delta = " << delta );
    1199         /// ---> A5
    1200         if ( degF <= degG && delta == degF && fdivides( F, G ) )
    1201         {
    1202             DEBDECLEVEL( cerr, "fin_ezgcd" );
    1203             return d*F;
    1204         }
    1205         if ( degG < degF && delta == degG && fdivides( G, F ) )
    1206         {
    1207             DEBDECLEVEL( cerr, "fin_ezgcd" );
    1208             return d*G;
    1209         }
    1210         if ( delta != degF && delta != degG ) {
    1211             bool B_is_F;
    1212             /// ---> A6
    1213             CanonicalForm xxx;
    1214             //if ( gcd( (DD[1] = Fb / Db), Db ) == 1 ) {
    1215             DD[1] = Fb / Db;
    1216             xxx= gcd( DD[1], Db );
    1217             DEBOUTLN( cerr, "gcd((Fb/Db),Db) = " << xxx );
    1218             DEBOUTLN( cerr, "Fb/Db = " << DD[1] );
    1219             DEBOUTLN( cerr, "Db = " << Db );
    1220             if (xxx.inCoeffDomain()) {
    1221                 B = F;
    1222                 DD[2] = Db;
    1223                 lcDD[1] = lcF;
    1224                 lcDD[2] = lcF;
    1225                 B *= lcF;
    1226                 B_is_F=true;
    1227             }
    1228             //else  if ( gcd( (DD[1] = Gb / Db), Db ) == 1 ) {
    1229             else
    1230             {
    1231               DD[1] = Gb / Db;
    1232               xxx=gcd( DD[1], Db );
    1233               DEBOUTLN( cerr, "gcd((Gb/Db),Db) = " << xxx );
    1234               DEBOUTLN( cerr, "Gb/Db = " << DD[1] );
    1235               DEBOUTLN( cerr, "Db = " << Db );
    1236               if (xxx.inCoeffDomain())
    1237               {
    1238                 B = G;
    1239                 DD[2] = Db;
    1240                 lcDD[1] = lcG;
    1241                 lcDD[2] = lcG;
    1242                 B *= lcG;
    1243                 B_is_F=false;
    1244               }
    1245               else
    1246               {
    1247 #ifdef DEBUGOUTPUT
    1248                 CanonicalForm dummyres = d * fin_ezgcd_specialcase( F, G, b, p );
    1249                 DEBDECLEVEL( cerr, "fin_ezgcd" );
    1250                 return dummyres;
    1251 #else
    1252                 return d * fin_ezgcd_specialcase( F, G, b, p );
    1253 #endif
    1254               }
    1255             }
    1256             /// ---> A7
    1257             DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
    1258             DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
    1259             DEBOUTLN( cerr, "(hensel) B    = " << B );
    1260             DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
    1261             DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
    1262             DEBOUTLN( cerr, "(hensel) DD   = " << DD );
    1263             DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
    1264             gcdfound = Hensel( B, DD, lcDD, b, p, x );
    1265             DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
    1266             if (gcdfound)
    1267             {
    1268               CanonicalForm xxx=content(DD[2],Variable(1));
    1269               CanonicalForm cand=DD[2] / xxx; //content(DD[2],Variable(1));
    1270 #if 0
    1271               gcdfound= fdivides(cand,G) &&  fdivides(cand,F);
    1272 #else
    1273               if (B_is_F /*B==F*lcF*/)
    1274               {
    1275                 DEBOUTLN( cerr, "(test) G: "<<G<<" % gcd:"<<cand<<" -> " << G%cand );
    1276                 gcdfound= fdivides(cand,G);
    1277               }
    1278               else
    1279               {
    1280                 DEBOUTLN( cerr, "(test) F: "<<F<<" % gcd:"<<cand<<" -> " << F%cand);
    1281                 gcdfound= fdivides(cand,F);
    1282               }
    1283 #endif
    1284             }
    1285             /// ---> A8 (gcdfound)
    1286         }
    1287         delta++;
    1288     }
    1289     /// ---> A9
    1290     DEBDECLEVEL( cerr, "fin_ezgcd" );
    1291     return d * DD[2] / content(DD[2],Variable(1));
    1292 }
    1293 
    1294 static CanonicalForm fin_ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, FFREvaluation & b, const modpk & bound )
    1295 {
    1296     CanonicalForm d;
    1297 #if 1
    1298     Off(SW_USE_EZGCD);
    1299     //bool ntl0=isOn(SW_USE_NTL_GCD_0);
    1300     //Off(SW_USE_NTL_GCD_0);
    1301     //bool ntlp=isOn(SW_USE_NTL_GCD_P);
    1302     //Off(SW_USE_NTL_GCD_P);
    1303     d=gcd( F, G );
    1304     //if (ntl0) On(SW_USE_NTL_GCD_0);
    1305     //if (ntlp) On(SW_USE_NTL_GCD_P);
    1306     On(SW_USE_EZGCD);
    1307     return d;
    1308 #else
    1309     DEBOUTLN( cerr, "fin_ezgcd: special case" );
    1310     CanonicalForm Ft, Gt, L, LL, Fb, Gb, Db, Lb, D, Ds, Dt;
    1311     CFArray DD( 1, 2 ), lcDD( 1, 2 );
    1312     Variable x = Variable( 1 );
    1313     bool gcdfound;
    1314     Dt = 1;
    1315     /// ---> S1
    1316     DEBOUTLN( cerr, "fin_ezgcdspec: (S1)" );
    1317     Ft = fin_ezgcd( F, F.deriv( x ) );
    1318     if ( Ft.isOne() ) {
    1319         // In this case F is squarefree and we came here by bad chance
    1320         // (means: bad evaluation point).  Try again with another
    1321         // evaluation point.  Bug fix (?) by JS.  The bad example was
    1322         // gcd.debug -ocr /+USE_EZGCD/@12/CB \
    1323         //     '(16*B^8-208*B^6*C+927*B^4*C^2-1512*B^2*C^3+432*C^4)' \
    1324         //     '(4*B^7*C^2-50*B^5*C^3+208*B^3*C^4-288*B*C^5)'
    1325         if( ! b.step() ) // Error: run out of points
    1326         {
    1327           Off(SW_USE_EZGCD_P);
    1328           d = gcd_poly_p(F,G);
    1329           On(SW_USE_EZGCD_P);
    1330           return d;
    1331         }
    1332         return fin_ezgcd( F, G, b, true );
    1333     }
    1334     DEBOUTLN( cerr, "fin_ezgcdspec: (S1) done, Ft = " << Ft );
    1335     L = F / Ft;
    1336     /// ---> S2
    1337     DEBOUTLN( cerr, "fin_ezgcdspec: (S2)" );
    1338     L = fin_ezgcd( L, G, b, true );
    1339     DEBOUTLN( cerr, "fin_ezgcdspec: (S2) done, Ds = " << Ds );
    1340     Gt = G / L;
    1341     Ds = L; Dt = L;
    1342     Lb = b( L ); Gb = b( Gt ); Fb = b( Ft );
    1343     d = gcd( LC( L, x ), gcd( LC( Ft, x ), LC( Gt, x ) ) );
    1344     while ( true ) {
    1345         /// ---> S3
    1346         DEBOUTLN( cerr, "fin_ezgcdspec: (S3)" );
    1347         DEBOUTLN( cerr, "fin_ezgcdspec: Fb = " << Fb );
    1348         DEBOUTLN( cerr, "fin_ezgcdspec: Gb = " << Gb );
    1349         DD[1] = gcd( Lb, gcd( Fb, Gb ) );
    1350         /// ---> S4
    1351         DEBOUTLN( cerr, "fin_ezgcdspec: (S4)" );
    1352         if ( degree( DD[1] ) == 0 )
    1353             return Ds;
    1354         /// ---> S5
    1355         DEBOUTLN( cerr, "fin_ezgcdspec: (S5)" );
    1356         DEBOUTLN( cerr, "fin_ezgcdspec: Lb = " << Lb );
    1357         DEBOUTLN( cerr, "fin_ezgcdspec: Db = " << DD[1] );
    1358         Db = DD[1];
    1359         if ( ! (DD[2] = Lb/DD[1]).isOne() ) {
    1360             DEBOUTLN( cerr, "fin_ezgcdspec: (S55)" );
    1361             lcDD[2] = LC( L, x );
    1362             lcDD[1] = d;
    1363             DD[1] = ( DD[1] * b( d ) ) / lc( DD[1] );
    1364             DD[2] = ( DD[2] * b( lcDD[2] ) ) / lc( DD[2] );
    1365             LL = L * d;
    1366             DEBOUTLN( cerr, "fin_ezgcdspec: begin with lift." );
    1367             DEBOUTLN( cerr, "fin_ezgcdspec: B     = " << LL );
    1368             DEBOUTLN( cerr, "fin_ezgcdspec: DD    = " << DD );
    1369             DEBOUTLN( cerr, "fin_ezgcdspec: lcDD  = " << lcDD );
    1370             DEBOUTLN( cerr, "fin_ezgcdspec: b     = " << b );
    1371             DEBOUTLN( cerr, "fin_ezgcdspec: lc(B) = " << LC( LL, x ) );
    1372             DEBOUTLN( cerr, "fin_ezgcdspec: test  = " << b( LL ) - DD[1] * DD[2] );
    1373             gcdfound = Hensel( LL, DD, lcDD, b, p, x );
    1374             ASSERT( gcdfound, "fatal error in algorithm" );
    1375             Dt = pp( DD[1] );
    1376         }
    1377         DEBOUTLN( cerr, "fin_ezgcdspec: (S7)" );
    1378         Ds = Ds * Dt;
    1379         Fb = Fb / Db;
    1380         Gb = Gb / Db;
    1381     }
    1382     // this point is never reached, but to avoid compiler warnings let's return a value
    1383     return 0;
    1384 #endif
    1385 }
    1386 
    1387 static bool fin_findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, FFREvaluation & b, int delta, int degF, int degG )
    1388 {
    1389     int i;
    1390     bool ok;
    1391     DEBOUTLN( cerr, "fin_ezgcd: (fin_findeval) F = " << F  <<", G="<< G);
    1392     DEBOUTLN( cerr, "fin_ezgcd: (fin_findeval) degF = " << degF << ", degG="<<degG );
    1393     while(1) {
    1394         DEBOUTLN( cerr, "fin_ezgcd: (fin_findeval) b = " << b );
    1395         Fb = b( F );
    1396         ok = degree( Fb ) == degF;
    1397         if ( ok ) {
    1398             Gb = b( G );
    1399             ok = degree( Gb ) == degG;
    1400         }
    1401         if ( ok )
    1402         {
    1403             Db = gcd( Fb, Gb );
    1404             if ( delta > 0 )
    1405               ok = degree( Db ) < delta;
    1406         }
    1407         if ( ok )
    1408             break;
    1409 
    1410         if( b.step() ) // can choose valid point
    1411            continue;
    1412 
    1413         return false;
    1414     }
    1415     return true;
    1416 }

  • factory/ffreval.cc

    r6fd3a2 rb809a8  
    6262}
    6363
     64/* ------------------------------------------------------------------------ */
     65/* forward declarations: fin_ezgcd stuff*/
     66static CanonicalForm fin_ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, FFREvaluation & b, bool internal );
     67static bool fin_findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, FFREvaluation & b, int delta, int degF, int degG );
     68static CanonicalForm fin_ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, FFREvaluation & b, const modpk & bound );
     69
     70CanonicalForm fin_ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG )
     71{
     72    FFREvaluation b;
     73    return fin_ezgcd( FF, GG, b, false );
     74}
     75
     76static CanonicalForm fin_ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, FFREvaluation & b, bool internal )
     77{
     78    CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, lcD;
     79    CFArray DD( 1, 2 ), lcDD( 1, 2 );
     80    int degF, degG, delta, t;
     81    FFREvaluation bt;
     82    bool gcdfound = false;
     83    Variable x = Variable(1);
     84    modpk p = modpk(getCharacteristic(), 1); // k = 1
     85    DEBINCLEVEL( cerr, "fin_ezgcd" );
     86    DEBOUTLN( cerr, "FF = " << FF );
     87    DEBOUTLN( cerr, "GG = " << GG );
     88    f = content( FF, x ); g = content( GG, x ); d = gcd( f, g );
     89    DEBOUTLN( cerr, "f = " << f );
     90    DEBOUTLN( cerr, "g = " << g );
     91    F = FF / f; G = GG / g;
     92    if ( F.isUnivariate() && G.isUnivariate() )
     93    {
     94        DEBDECLEVEL( cerr, "fin_ezgcd" );
     95        if(F.mvar()==G.mvar())
     96          d*=gcd_poly(F,G);
     97        return d;
     98    }
     99    else  if ( gcd_test_one( F, G, false ) )
     100    {
     101        DEBDECLEVEL( cerr, "fin_ezgcd" );
     102        return d;
     103    }
     104    lcF = LC( F, x ); lcG = LC( G, x );
     105    lcD = gcd( lcF, lcG );
     106    delta = 0;
     107    degF = degree( F, x ); degG = degree( G, x );
     108    t = tmax( F.level(), G.level() );
     109    if ( ! internal )
     110    {
     111        b = FFREvaluation( 2, t, FFRandom() );
     112        b.init(); // choose random point
     113    }
     114    while ( ! gcdfound ) {
     115        /// ---> A2
     116        DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
     117        DEBOUTLN( cerr, "F = " << F );
     118        DEBOUTLN( cerr, "G = " << G );
     119        if( ! fin_findeval( F, G, Fb, Gb, Db, b, delta, degF, degG ) )
     120        {
     121           Off(SW_USE_EZGCD_P);
     122           d *= gcd_poly(F,G);
     123           On(SW_USE_EZGCD_P);
     124           return d;
     125        }
     126
     127        DEBOUTLN( cerr, "found evaluation b = " << b );
     128        DEBOUTLN( cerr, "F(b) = " << Fb );
     129        DEBOUTLN( cerr, "G(b) = " << Gb );
     130        DEBOUTLN( cerr, "D(b) = " << Db );
     131        delta = degree( Db );
     132        /// ---> A3
     133        if ( delta == 0 )
     134        {
     135          DEBDECLEVEL( cerr, "fin_ezgcd" );
     136          return d;
     137        }
     138        /// ---> A4
     139        //deltaold = delta;
     140        while ( 1 ) {
     141            bt = b;
     142            if( ! fin_findeval( F, G, Fbt, Gbt, Dbt, bt, delta + 1, degF, degG ) )
     143            {
     144              Off(SW_USE_EZGCD_P);
     145              d *= gcd_poly(F,G);
     146              On(SW_USE_EZGCD_P);
     147              return d;
     148            }
     149
     150            int dd=degree( Dbt );
     151            if ( dd /*degree( Dbt )*/ == 0 )
     152            {
     153                DEBDECLEVEL( cerr, "fin_ezgcd" );
     154                return d;
     155            }
     156            if ( dd /*degree( Dbt )*/ == delta )
     157                break;
     158            else  if ( dd /*degree( Dbt )*/ < delta ) {
     159                delta = dd /*degree( Dbt )*/;
     160                b = bt;
     161                Db = Dbt; Fb = Fbt; Gb = Gbt;
     162            }
     163        }
     164        DEBOUTLN( cerr, "now after A4, delta = " << delta );
     165        /// ---> A5
     166        if ( degF <= degG && delta == degF && fdivides( F, G ) )
     167        {
     168            DEBDECLEVEL( cerr, "fin_ezgcd" );
     169            return d*F;
     170        }
     171        if ( degG < degF && delta == degG && fdivides( G, F ) )
     172        {
     173            DEBDECLEVEL( cerr, "fin_ezgcd" );
     174            return d*G;
     175        }
     176        if ( delta != degF && delta != degG ) {
     177            bool B_is_F;
     178            /// ---> A6
     179            CanonicalForm xxx;
     180            //if ( gcd( (DD[1] = Fb / Db), Db ) == 1 ) {
     181            DD[1] = Fb / Db;
     182            xxx= gcd( DD[1], Db );
     183            DEBOUTLN( cerr, "gcd((Fb/Db),Db) = " << xxx );
     184            DEBOUTLN( cerr, "Fb/Db = " << DD[1] );
     185            DEBOUTLN( cerr, "Db = " << Db );
     186            if (xxx.inCoeffDomain()) {
     187                B = F;
     188                DD[2] = Db;
     189                lcDD[1] = lcF;
     190                lcDD[2] = lcF;
     191                B *= lcF;
     192                B_is_F=true;
     193            }
     194            //else  if ( gcd( (DD[1] = Gb / Db), Db ) == 1 ) {
     195            else
     196            {
     197              DD[1] = Gb / Db;
     198              xxx=gcd( DD[1], Db );
     199              DEBOUTLN( cerr, "gcd((Gb/Db),Db) = " << xxx );
     200              DEBOUTLN( cerr, "Gb/Db = " << DD[1] );
     201              DEBOUTLN( cerr, "Db = " << Db );
     202              if (xxx.inCoeffDomain())
     203              {
     204                B = G;
     205                DD[2] = Db;
     206                lcDD[1] = lcG;
     207                lcDD[2] = lcG;
     208                B *= lcG;
     209                B_is_F=false;
     210              }
     211              else
     212              {
     213#ifdef DEBUGOUTPUT
     214                CanonicalForm dummyres = d * fin_ezgcd_specialcase( F, G, b, p );
     215                DEBDECLEVEL( cerr, "fin_ezgcd" );
     216                return dummyres;
     217#else
     218                return d * fin_ezgcd_specialcase( F, G, b, p );
     219#endif
     220              }
     221            }
     222            /// ---> A7
     223            DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
     224            DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
     225            DEBOUTLN( cerr, "(hensel) B    = " << B );
     226            DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
     227            DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
     228            DEBOUTLN( cerr, "(hensel) DD   = " << DD );
     229            DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
     230            gcdfound = Hensel( B, DD, lcDD, b, p, x );
     231            DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
     232            if (gcdfound)
     233            {
     234              CanonicalForm xxx=content(DD[2],Variable(1));
     235              CanonicalForm cand=DD[2] / xxx; //content(DD[2],Variable(1));
     236#if 0
     237              gcdfound= fdivides(cand,G) &&  fdivides(cand,F);
     238#else
     239              if (B_is_F /*B==F*lcF*/)
     240              {
     241                DEBOUTLN( cerr, "(test) G: "<<G<<" % gcd:"<<cand<<" -> " << G%cand );
     242                gcdfound= fdivides(cand,G);
     243              }
     244              else
     245              {
     246                DEBOUTLN( cerr, "(test) F: "<<F<<" % gcd:"<<cand<<" -> " << F%cand);
     247                gcdfound= fdivides(cand,F);
     248              }
     249#endif
     250            }
     251            /// ---> A8 (gcdfound)
     252        }
     253        delta++;
     254    }
     255    /// ---> A9
     256    DEBDECLEVEL( cerr, "fin_ezgcd" );
     257    return d * DD[2] / content(DD[2],Variable(1));
     258}
     259
     260static CanonicalForm fin_ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, FFREvaluation & b, const modpk & bound )
     261{
     262    CanonicalForm d;
     263#if 1
     264    Off(SW_USE_EZGCD);
     265    //bool ntl0=isOn(SW_USE_NTL_GCD_0);
     266    //Off(SW_USE_NTL_GCD_0);
     267    //bool ntlp=isOn(SW_USE_NTL_GCD_P);
     268    //Off(SW_USE_NTL_GCD_P);
     269    d=gcd_poly( F, G );
     270    //if (ntl0) On(SW_USE_NTL_GCD_0);
     271    //if (ntlp) On(SW_USE_NTL_GCD_P);
     272    On(SW_USE_EZGCD);
     273    return d;
     274#else
     275    DEBOUTLN( cerr, "fin_ezgcd: special case" );
     276    CanonicalForm Ft, Gt, L, LL, Fb, Gb, Db, Lb, D, Ds, Dt;
     277    CFArray DD( 1, 2 ), lcDD( 1, 2 );
     278    Variable x = Variable( 1 );
     279    bool gcdfound;
     280    Dt = 1;
     281    /// ---> S1
     282    DEBOUTLN( cerr, "fin_ezgcdspec: (S1)" );
     283    Ft = fin_ezgcd( F, F.deriv( x ) );
     284    if ( Ft.isOne() ) {
     285        // In this case F is squarefree and we came here by bad chance
     286        // (means: bad evaluation point).  Try again with another
     287        // evaluation point.  Bug fix (?) by JS.  The bad example was
     288        // gcd.debug -ocr /+USE_EZGCD/@12/CB \
     289        //     '(16*B^8-208*B^6*C+927*B^4*C^2-1512*B^2*C^3+432*C^4)' \
     290        //     '(4*B^7*C^2-50*B^5*C^3+208*B^3*C^4-288*B*C^5)'
     291        if( ! b.step() ) // Error: run out of points
     292        {
     293          Off(SW_USE_EZGCD_P);
     294          d = gcd_poly(F,G);
     295          On(SW_USE_EZGCD_P);
     296          return d;
     297        }
     298        return fin_ezgcd( F, G, b, true );
     299    }
     300    DEBOUTLN( cerr, "fin_ezgcdspec: (S1) done, Ft = " << Ft );
     301    L = F / Ft;
     302    /// ---> S2
     303    DEBOUTLN( cerr, "fin_ezgcdspec: (S2)" );
     304    L = fin_ezgcd( L, G, b, true );
     305    DEBOUTLN( cerr, "fin_ezgcdspec: (S2) done, Ds = " << Ds );
     306    Gt = G / L;
     307    Ds = L; Dt = L;
     308    Lb = b( L ); Gb = b( Gt ); Fb = b( Ft );
     309    d = gcd( LC( L, x ), gcd( LC( Ft, x ), LC( Gt, x ) ) );
     310    while ( true ) {
     311        /// ---> S3
     312        DEBOUTLN( cerr, "fin_ezgcdspec: (S3)" );
     313        DEBOUTLN( cerr, "fin_ezgcdspec: Fb = " << Fb );
     314        DEBOUTLN( cerr, "fin_ezgcdspec: Gb = " << Gb );
     315        DD[1] = gcd( Lb, gcd( Fb, Gb ) );
     316        /// ---> S4
     317        DEBOUTLN( cerr, "fin_ezgcdspec: (S4)" );
     318        if ( degree( DD[1] ) == 0 )
     319            return Ds;
     320        /// ---> S5
     321        DEBOUTLN( cerr, "fin_ezgcdspec: (S5)" );
     322        DEBOUTLN( cerr, "fin_ezgcdspec: Lb = " << Lb );
     323        DEBOUTLN( cerr, "fin_ezgcdspec: Db = " << DD[1] );
     324        Db = DD[1];
     325        if ( ! (DD[2] = Lb/DD[1]).isOne() ) {
     326            DEBOUTLN( cerr, "fin_ezgcdspec: (S55)" );
     327            lcDD[2] = LC( L, x );
     328            lcDD[1] = d;
     329            DD[1] = ( DD[1] * b( d ) ) / lc( DD[1] );
     330            DD[2] = ( DD[2] * b( lcDD[2] ) ) / lc( DD[2] );
     331            LL = L * d;
     332            DEBOUTLN( cerr, "fin_ezgcdspec: begin with lift." );
     333            DEBOUTLN( cerr, "fin_ezgcdspec: B     = " << LL );
     334            DEBOUTLN( cerr, "fin_ezgcdspec: DD    = " << DD );
     335            DEBOUTLN( cerr, "fin_ezgcdspec: lcDD  = " << lcDD );
     336            DEBOUTLN( cerr, "fin_ezgcdspec: b     = " << b );
     337            DEBOUTLN( cerr, "fin_ezgcdspec: lc(B) = " << LC( LL, x ) );
     338            DEBOUTLN( cerr, "fin_ezgcdspec: test  = " << b( LL ) - DD[1] * DD[2] );
     339            gcdfound = Hensel( LL, DD, lcDD, b, p, x );
     340            ASSERT( gcdfound, "fatal error in algorithm" );
     341            Dt = pp( DD[1] );
     342        }
     343        DEBOUTLN( cerr, "fin_ezgcdspec: (S7)" );
     344        Ds = Ds * Dt;
     345        Fb = Fb / Db;
     346        Gb = Gb / Db;
     347    }
     348    // this point is never reached, but to avoid compiler warnings let's return a value
     349    return 0;
     350#endif
     351}
     352
     353static bool fin_findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, FFREvaluation & b, int delta, int degF, int degG )
     354{
     355    int i;
     356    bool ok;
     357    DEBOUTLN( cerr, "fin_ezgcd: (fin_findeval) F = " << F  <<", G="<< G);
     358    DEBOUTLN( cerr, "fin_ezgcd: (fin_findeval) degF = " << degF << ", degG="<<degG );
     359    while(1) {
     360        DEBOUTLN( cerr, "fin_ezgcd: (fin_findeval) b = " << b );
     361        Fb = b( F );
     362        ok = degree( Fb ) == degF;
     363        if ( ok ) {
     364            Gb = b( G );
     365            ok = degree( Gb ) == degG;
     366        }
     367        if ( ok )
     368        {
     369            Db = gcd( Fb, Gb );
     370            if ( delta > 0 )
     371              ok = degree( Db ) < delta;
     372        }
     373        if ( ok )
     374            break;
     375
     376        if( b.step() ) // can choose valid point
     377           continue;
     378
     379        return false;
     380    }
     381    return true;
     382}
Note: See TracChangeset for help on using the changeset viewer.