Changeset f11d7b in git


Ignore:
Timestamp:
Jan 19, 2004, 12:26:20 PM (20 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '648d28f488f6ff08f5607ff229b9ad9e4a5b93c2')
Children:
b44d16e5c775914df07cc382ca039bb32caae77a
Parents:
bc9f64c8ededf224f4b2323f83f51d42ba00716e
Message:
*hannes: gcd via NTL


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

Legend:

Unmodified
Added
Removed
  • factory/NTLconvert.cc

    rbc9f64c rf11d7b  
    1 /* $Id: NTLconvert.cc,v 1.12 2003-10-15 17:19:38 Singular Exp $ */
     1/* $Id: NTLconvert.cc,v 1.13 2004-01-19 11:26:19 Singular Exp $ */
    22#include <config.h>
    33
     
    2020#include <NTL/ZZXFactoring.h>
    2121#include <NTL/ZZ_pXFactoring.h>
     22#include <NTL/lzz_pXFactoring.h>
    2223#include <NTL/GF2XFactoring.h>
    2324#include <NTL/ZZ_pEXFactoring.h>
     25#include <NTL/lzz_pEXFactoring.h>
    2426#include <NTL/GF2EXFactoring.h>
    2527#include <NTL/tools.h>
     
    146148  return ntl_poly;
    147149}
     150zz_pX convertFacCF2NTLzzpX(CanonicalForm f)
     151{
     152  zz_pX ntl_poly;
     153
     154  CFIterator i;
     155  i=f;
     156
     157  int j=0;
     158  int NTLcurrentExp=i.exp();
     159  int largestExp=i.exp();
     160  int k;
     161
     162  // we now build up the NTL-polynomial
     163  ntl_poly.SetMaxLength(largestExp+1);
     164
     165  for (;i.hasTerms();i++)
     166  {
     167    for (k=NTLcurrentExp;k>i.exp();k--)
     168    {
     169      SetCoeff(ntl_poly,k,0);
     170    }
     171    NTLcurrentExp=i.exp();
     172
     173    CanonicalForm c=i.coeff();
     174    if (!c.isImm()) c.mapinto(); //c%= getCharacteristic();
     175    if (!c.isImm())
     176    {  //This case will never happen if the characteristic is in fact a prime
     177       // number, since all coefficients are represented as immediates
     178       #ifndef NOSTREAMIO
     179       cout<<"convertFacCF2NTLzz_pX: coefficient not immediate! : "<<f<<"\n";
     180       #else
     181       printf("convertFacCF2NTLzz_pX: coefficient not immediate!, char=%d\n",
     182              getCharacteristic());
     183       #endif
     184       exit(1);
     185    }
     186    else
     187    {
     188      SetCoeff(ntl_poly,NTLcurrentExp,c.intval());
     189    }
     190    NTLcurrentExp--;
     191  }
     192
     193  //Set the remaining coefficients of ntl_poly to zero.
     194  // This is necessary, because NTL internally
     195  // also stores powers with zero coefficient,
     196  // whereas factory stores tuples of degree and coefficient
     197  //leaving out tuples if the coefficient equals zero
     198  for (k=NTLcurrentExp;k>=0;k--)
     199  {
     200    SetCoeff(ntl_poly,k,0);
     201  }
     202
     203  //normalize the polynomial and return it
     204  ntl_poly.normalize();
     205
     206  return ntl_poly;
     207}
    148208
    149209////////////////////////////////////////////////////////////////////////////////
     
    261321}
    262322
     323CanonicalForm convertNTLzzpX2CF(zz_pX poly,Variable x)
     324{
     325  //printf("convertNTLzzpX2CF\n");
     326  CanonicalForm bigone;
     327
     328
     329  if (deg(poly)>0)
     330  {
     331    // poly is non-constant
     332    bigone=0;
     333    bigone.mapinto();
     334    // Compute the canonicalform coefficient by coefficient,
     335    // bigone summarizes the result.
     336    for (int j=0;j<deg(poly)+1;j++)
     337    {
     338      if (coeff(poly,j)!=0)
     339      {
     340        bigone+=(power(x,j)*CanonicalForm(to_long(rep(coeff(poly,j)))));
     341      }
     342    }
     343  }
     344  else
     345  {
     346    // poly is immediate
     347    bigone=CanonicalForm(to_long(rep(coeff(poly,0))));
     348    bigone.mapinto();
     349  }
     350  return bigone;
     351}
     352
     353CanonicalForm convertNTLZZX2CF(ZZX polynom,Variable x)
     354{
     355  //printf("convertNTLZZX2CF\n");
     356  CanonicalForm bigone;
     357
     358  // Go through the vector e and build up the CFFList
     359  // As usual bigone summarizes the result
     360  bigone=0;
     361  ZZ coefficient;
     362
     363  for (int j=0;j<=deg(polynom);j++)
     364  {
     365    coefficient=coeff(polynom,j);
     366    if (!IsZero(coefficient))
     367    {
     368      bigone += (power(x,j)*convertZZ2CF(coefficient));
     369    }
     370  }
     371  return bigone;
     372}
    263373
    264374////////////////////////////////////////////////////////////////////////////////
     
    357467  {
    358468    rueckgabe.append(CFFactor(convertNTLZZpX2CF(e[i].a,x),e[i].b));
     469  }
     470  if(isOn(SW_USE_NTL_SORT)) rueckgabe.sort(NTLcmpCF);
     471  // the multiplicity at pos 1
     472  if (!IsOne(multi))
     473    rueckgabe.insert(CFFactor(CanonicalForm(to_long(rep(multi))),1));
     474  return rueckgabe;
     475}
     476CFFList convertNTLvec_pair_zzpX_long2FacCFFList
     477                                  (vec_pair_zz_pX_long e,zz_p multi,Variable x)
     478{
     479  //printf("convertNTLvec_pair_zzpX_long2FacCFFList\n");
     480  CFFList rueckgabe;
     481  zz_pX polynom;
     482  long exponent;
     483  CanonicalForm bigone;
     484
     485  // Maybe, e may additionally be sorted with respect to increasing degree of x
     486  // but this is not
     487  //important for the factorization, but nevertheless would take computing time,
     488  // so it is omitted
     489
     490
     491  // Go through the vector e and compute the CFFList
     492  // again bigone summarizes the result
     493  for (int i=e.length()-1;i>=0;i--)
     494  {
     495    rueckgabe.append(CFFactor(convertNTLzzpX2CF(e[i].a,x),e[i].b));
    359496  }
    360497  if(isOn(SW_USE_NTL_SORT)) rueckgabe.sort(NTLcmpCF);
     
    629766  CanonicalForm bigone;
    630767
    631   // Maybe, e may additionally be sorted with respect to increasing degree of x, but this is not
    632   //important for the factorization, but nevertheless would take computing time, so it is omitted
    633 
    634 
    635768  // Go through the vector e and build up the CFFList
    636769  // As usual bigone summarizes the result
    637770  for (int i=e.length()-1;i>=0;i--)
    638771  {
    639     bigone=0;
    640772    ZZ coefficient;
    641773    polynom=e[i].a;
    642774    exponent=e[i].b;
    643     long coeff_long;
    644 
    645     for (int j=0;j<deg(polynom)+1;j++)
    646     {
    647       coefficient=coeff(polynom,j);
    648       if (!IsZero(coefficient))
    649       {
    650         bigone += (power(x,j)*convertZZ2CF(coefficient));
    651       }
    652     }
    653 
     775    bigone=convertNTLZZX2CF(polynom,x);
    654776    //append the converted polynomial to the list
    655777    rueckgabe.append(CFFactor(bigone,exponent));
     
    684806{
    685807  return convertNTLZZpX2CF(rep(coefficient),x);
     808}
     809CanonicalForm convertNTLzzpE2CF(zz_pE coefficient,Variable x)
     810{
     811  return convertNTLzzpX2CF(rep(coefficient),x);
    686812}
    687813
     
    751877  return rueckgabe;
    752878}
     879CFFList convertNTLvec_pair_zzpEX_long2FacCFFList(vec_pair_zz_pEX_long e,zz_pE multi,Variable x,Variable alpha)
     880{
     881  CFFList rueckgabe;
     882  zz_pEX polynom;
     883  long exponent;
     884  CanonicalForm bigone;
     885
     886  // Maybe, e may additionally be sorted with respect to increasing degree of x, but this is not
     887  //important for the factorization, but nevertheless would take computing time, so it is omitted
     888
     889  // Go through the vector e and build up the CFFList
     890  // As usual bigone summarizes the result during every loop
     891  for (int i=e.length()-1;i>=0;i--)
     892  {
     893    bigone=0;
     894
     895    polynom=e[i].a;
     896    exponent=e[i].b;
     897
     898    for (int j=0;j<deg(polynom)+1;j++)
     899    {
     900      if (IsOne(coeff(polynom,j)))
     901      {
     902        bigone+=power(x,j);
     903      }
     904      else
     905      {
     906        CanonicalForm coefficient=convertNTLzzpE2CF(coeff(polynom,j),alpha);
     907        if (coeff(polynom,j)!=0)
     908        {
     909          bigone += (power(x,j)*coefficient);
     910        }
     911      }
     912    }
     913    //append the computed polynomials together with its exponent to the CFFList
     914    rueckgabe.append(CFFactor(bigone,exponent));
     915  }
     916  if(isOn(SW_USE_NTL_SORT)) rueckgabe.sort(NTLcmpCF);
     917  // Start by appending the multiplicity
     918  if (!IsOne(multi))
     919    rueckgabe.insert(CFFactor(convertNTLzzpE2CF(multi,alpha),1));
     920
     921  //return the computed CFFList
     922  return rueckgabe;
     923}
    753924
    754925////////////////////////////////////////////////////////////////////////////////
     
    9031074  return result;
    9041075}
     1076zz_pEX convertFacCF2NTLzz_pEX(CanonicalForm f, zz_pX mipo)
     1077{
     1078  zz_pE::init(mipo);
     1079  zz_pEX result;
     1080  CFIterator i;
     1081  i=f;
     1082
     1083  int j=0;
     1084  int NTLcurrentExp=i.exp();
     1085  int largestExp=i.exp();
     1086  int k;
     1087
     1088  result.SetMaxLength(largestExp+1);
     1089  for(;i.hasTerms();i++)
     1090  {
     1091    for(k=NTLcurrentExp;k>i.exp();k--) SetCoeff(result,k,0);
     1092    NTLcurrentExp=i.exp();
     1093    CanonicalForm c=i.coeff();
     1094    zz_pX cc=convertFacCF2NTLzzpX(c);
     1095    //ZZ_pE ccc;
     1096    //conv(ccc,cc);
     1097    SetCoeff(result,NTLcurrentExp,to_zz_pE(cc));
     1098    NTLcurrentExp--;
     1099  }
     1100  for(k=NTLcurrentExp;k>=0;k--) SetCoeff(result,k,0);
     1101  result.normalize();
     1102  return result;
     1103}
    9051104#endif
  • factory/NTLconvert.h

    rbc9f64c rf11d7b  
    1 /* $Id: NTLconvert.h,v 1.6 2003-08-28 14:15:18 Singular Exp $ */
     1/* $Id: NTLconvert.h,v 1.7 2004-01-19 11:26:20 Singular Exp $ */
    22#ifndef INCL_NTLCONVERT_H
    33#define INCL_NTLCONVERT_H
     
    2323#include <NTL/ZZXFactoring.h>
    2424#include <NTL/ZZ_pXFactoring.h>
     25#include <NTL/lzz_pXFactoring.h>
    2526#include <NTL/GF2XFactoring.h>
    2627#include "int_int.h"
     
    3435
    3536ZZ_pX convertFacCF2NTLZZpX(CanonicalForm f);
     37zz_pX convertFacCF2NTLzzpX(CanonicalForm f);
    3638GF2X convertFacCF2NTLGF2X(CanonicalForm f);
    3739CanonicalForm convertNTLZZpX2CF(ZZ_pX poly,Variable x);
     40CanonicalForm convertNTLzzpX2CF(zz_pX poly,Variable x);
    3841CanonicalForm convertNTLGF2X2CF(GF2X poly,Variable x);
     42CanonicalForm convertNTLZZX2CF(ZZX polynom,Variable x);
    3943CFFList convertNTLvec_pair_ZZpX_long2FacCFFList(vec_pair_ZZ_pX_long e,ZZ_p multi,Variable x);
     44CFFList convertNTLvec_pair_zzpX_long2FacCFFList(vec_pair_zz_pX_long e,zz_p multi,Variable x);
    4045
    4146CFFList convertNTLvec_pair_GF2X_long2FacCFFList(vec_pair_GF2X_long e,GF2 multi,Variable x);
  • factory/cf_gcd.cc

    rbc9f64c rf11d7b  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: cf_gcd.cc,v 1.21 2001-01-30 13:41:12 Singular Exp $ */
     2/* $Id: cf_gcd.cc,v 1.22 2004-01-19 11:26:20 Singular Exp $ */
    33
    44#include <config.h>
     
    1818#include "ftmpl_functions.h"
    1919
     20#ifdef HAVE_NTL
     21#include "NTLconvert.h"
     22#endif
     23
    2024static CanonicalForm gcd_poly( const CanonicalForm & f, const CanonicalForm& g, bool modularflag );
    2125
     
    170174gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
    171175{
    172     CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
    173     int p, i, n;
    174 
    175     if ( primitive ) {
    176         f = F;
    177         g = G;
    178         c = 1;
    179     }
    180     else {
    181         CanonicalForm cF = content( F ), cG = content( G );
    182         f = F / cF;
    183         g = G / cG;
    184         c = bgcd( cF, cG );
    185     }
    186     cg = gcd( f.lc(), g.lc() );
    187     cl = ( f.lc() / cg ) * g.lc();
     176#ifdef HAVE_NTL
     177  if (isOn(SW_USE_NTL_GCD))
     178  {
     179    if ( getCharacteristic() > 0 )
     180    {
     181      zz_pContext ccc(getCharacteristic());
     182      ccc.restore();
     183      zz_p::init(getCharacteristic());
     184      zz_pX F1=convertFacCF2NTLzzpX(F);
     185      zz_pX G1=convertFacCF2NTLzzpX(G);
     186      zz_pX R=GCD(F1,G1);
     187      return convertNTLzzpX2CF(R,F.mvar());
     188    }
     189    else
     190    {
     191      ZZX F1=convertFacCF2NTLZZX(F);
     192      ZZX G1=convertFacCF2NTLZZX(G);
     193      ZZX R=GCD(F1,G1);
     194      return convertNTLZZX2CF(R,F.mvar());
     195    }
     196  }
     197#endif
     198  CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
     199  int p, i, n;
     200
     201  if ( primitive )
     202  {
     203    f = F;
     204    g = G;
     205    c = 1;
     206  }
     207  else
     208  {
     209    CanonicalForm cF = content( F ), cG = content( G );
     210    f = F / cF;
     211    g = G / cG;
     212    c = bgcd( cF, cG );
     213  }
     214  cg = gcd( f.lc(), g.lc() );
     215  cl = ( f.lc() / cg ) * g.lc();
    188216//     B = 2 * cg * tmin(
    189217//         maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
    190218//         maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
    191219//         )+1;
    192     M = tmin( maxNorm(f), maxNorm(g) );
    193     BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
    194     q = 0;
    195     i = cf_getNumSmallPrimes() - 1;
    196     while ( true ) {
    197         B = BB;
    198         while ( i >= 0 && q < B ) {
    199             p = cf_getSmallPrime( i );
    200             i--;
    201             while ( i >= 0 && mod( cl, p ) == 0 ) {
    202                 p = cf_getSmallPrime( i );
    203                 i--;
    204             }
    205             setCharacteristic( p );
    206             Dp = gcd( mapinto( f ), mapinto( g ) );
    207             Dp = ( Dp / Dp.lc() ) * mapinto( cg );
    208             setCharacteristic( 0 );
    209             if ( Dp.degree() == 0 ) return c;
    210             if ( q.isZero() ) {
    211                 D = mapinto( Dp );
    212                 q = p;
    213                 B = power(CanonicalForm(2),D.degree())*M+1;
    214             }
    215             else {
    216                 if ( Dp.degree() == D.degree() ) {
    217                     chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
    218                     q = newq;
    219                     D = newD;
    220                 }
    221                 else if ( Dp.degree() < D.degree() ) {
    222                     // all previous p's are bad primes
    223                     q = p;
    224                     D = mapinto( Dp );
    225                     B = power(CanonicalForm(2),D.degree())*M+1;
    226                 }
    227                 // else p is a bad prime
    228             }
     220  M = tmin( maxNorm(f), maxNorm(g) );
     221  BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
     222  q = 0;
     223  i = cf_getNumSmallPrimes() - 1;
     224  while ( true )
     225  {
     226    B = BB;
     227    while ( i >= 0 && q < B )
     228    {
     229      p = cf_getSmallPrime( i );
     230      i--;
     231      while ( i >= 0 && mod( cl, p ) == 0 )
     232      {
     233        p = cf_getSmallPrime( i );
     234        i--;
     235      }
     236      setCharacteristic( p );
     237      Dp = gcd( mapinto( f ), mapinto( g ) );
     238      Dp = ( Dp / Dp.lc() ) * mapinto( cg );
     239      setCharacteristic( 0 );
     240      if ( Dp.degree() == 0 )
     241        return c;
     242      if ( q.isZero() )
     243      {
     244        D = mapinto( Dp );
     245        q = p;
     246        B = power(CanonicalForm(2),D.degree())*M+1;
     247      }
     248      else
     249      {
     250        if ( Dp.degree() == D.degree() )
     251        {
     252          chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
     253          q = newq;
     254          D = newD;
    229255        }
    230         if ( i >= 0 ) {
    231             // now balance D mod q
    232             D = pp( balance( D, q ) );
    233             if ( divides( D, f ) && divides( D, g ) )
    234                 return D * c;
    235             else
    236                 q = 0;
     256        else if ( Dp.degree() < D.degree() )
     257        {
     258          // all previous p's are bad primes
     259          q = p;
     260          D = mapinto( Dp );
     261          B = power(CanonicalForm(2),D.degree())*M+1;
    237262        }
    238         else {
    239             return gcd_poly( F, G, false );
    240         }
    241         DEBOUTLN( cerr, "another try ..." );
    242     }
     263        // else p is a bad prime
     264      }
     265    }
     266    if ( i >= 0 )
     267    {
     268      // now balance D mod q
     269      D = pp( balance( D, q ) );
     270      if ( divides( D, f ) && divides( D, g ) )
     271        return D * c;
     272      else
     273        q = 0;
     274    }
     275    else
     276      return gcd_poly( F, G, false );
     277    DEBOUTLN( cerr, "another try ..." );
     278  }
    243279}
    244280
Note: See TracChangeset for help on using the changeset viewer.