Loading...
Searching...
No Matches
fac_multihensel.cc
Go to the documentation of this file.
/* emacs edit mode for this file is -*- C++ -*- */
/* $Id: fac_multihensel.cc 12231 2009-11-02 10:12:22Z hannes $ */
#include <config.h>
#include "assert.h"
#include "debug.h"
#include "timing.h"
#include "cf_defs.h"
#include "cf_eval.h"
#include "cf_binom.h"
#include "fac_util.h"
#include "fac_iterfor.h"
#include "cf_iter.h"
#ifndef HAVE_NTL
TIMING_DEFINE_PRINT(fac_solve);
TIMING_DEFINE_PRINT(fac_modpk);
TIMING_DEFINE_PRINT(fac_corrcoeff);
TIMING_DEFINE_PRINT(fac_extgcd);
static void
extgcdrest ( const CanonicalForm & a, const CanonicalForm & b, const CanonicalForm & s, const CanonicalForm & t, const CanonicalForm & c, CanonicalForm & S, CanonicalForm & T, const modpk & /*pk*/ )
{
CanonicalForm sigma = s * c, tau = t * c;
// divremainder( sigma, b, T, S, pk );
// T = pk( tau + T * a );
divrem( sigma, b, T, S );
T = tau + T * a;
}
static void
solveF ( const CFArray & P, const CFArray & Q, const CFArray & S, const CFArray & T, const CanonicalForm & C, const modpk & pk, int r, CFArray & a )
{
setCharacteristic( pk.getp(), pk.getk() );
CanonicalForm g, bb, b = mapinto( C );
int j;
for ( j = 1; j < r; j++ )
{
extgcdrest( mapinto( P[j] ), mapinto( Q[j] ), mapinto( S[j] ), mapinto( T[j] ), b, bb, a[j], pk );
b = bb;
}
a[r] = b;
setCharacteristic( 0 );
for ( j = 1; j <= r; j++ )
a[j] = mapinto( a[j] );
}
static CanonicalForm
evalF ( const CFArray & P, const CFArray & Q, const CFArray & A, int r )
{
CanonicalForm pprod = 1, sum = 0;
for ( int i = 1; i <= r; i++ )
{
sum += pprod * A[i] * Q[i];
pprod *= P[i];
}
return sum;
}
static CanonicalForm
derivAndEval ( const CanonicalForm & f, int n, const Variable & x, const CanonicalForm & a )
{
if ( n == 0 )
return f( a, x );
else if ( f.degree( x ) < n )
return 0;
else
{
CFIterator i;
CanonicalForm sum = 0, fact;
int min, j;
Variable v = Variable( f.level() + 1 );
for ( i = swapvar( f, x, v); i.hasTerms() && i.exp() >= n; i++ )
{
fact = 1;
min = i.exp() - n;
for ( j = i.exp(); j > min; j-- )
fact *= j;
sum += fact * i.coeff() * power( v, min );
}
return sum( a, v );
}
}
static CanonicalForm
checkCombination ( const CanonicalForm & W, const Evaluation & a, const IteratedFor & e, int k )
{
CanonicalForm dW = W;
int i, j;
for ( i = k-1; i >= 2 && ! dW.isZero(); i-- )
dW = derivAndEval( dW, e[i], Variable( i ), a[i] );
if ( ! dW.isZero() ) {
CanonicalForm fact = 1;
for ( i = 2; i < k; i++ )
for ( j = 2; j <= e[i]; j++ )
fact *= j;
dW /= fact;
}
return dW;
}
static CanonicalForm
prodCombination ( const Evaluation & a, const IteratedFor & e, int k )
{
CanonicalForm p = 1;
for ( int i = k-1; i > 1; i-- )
p *= binomialpower( Variable(i), -a[i], e[i] );
return p;
}
//static CanonicalForm check_dummy( const CFArray &a, const CFArray & P, const CFArray & Q )
//{
// int i, r = a.size();
// CanonicalForm res, prod;
// res = 0;
// prod = 1;
// for ( i = 1; i <= r; i++ )
// {
// res += prod * a[i] * Q[i];
// prod *= P[i];
// }
// return res;
//}
static bool check_e( const IteratedFor & e, int k, int m, int * n )
{
int sum = 0;
for ( int i = 2; i < k; i++ )
{
sum += e[i];
if ( e[i] > n[i] )
return false;
}
return sum == m+1;
}
static CanonicalForm modDeltak ( const CanonicalForm & f, const Evaluation & A, int k, int * n )
{
CanonicalForm result = f;
for ( int i = 2; i < k; i++ )
{
result.mod( binomialpower( Variable(i), -A[i], n[i]+1 ) );
}
return result;
}
static CFArray
findCorrCoeffs ( const CFArray & P, const CFArray & Q, const CFArray & P0, const CFArray & Q0, const CFArray & S, const CFArray & T, const CanonicalForm & C, const Evaluation & I, const modpk & pk, int r, int k, int h, int * n )
{
DEBINCLEVEL( cerr, "findCorrCoeffs" );
int i, m;
CFArray A(1,r), a(1,r);
CanonicalForm C0, Dm, g, prodcomb;
TIMING_START(fac_solve);
C0 = pk( I( C, 2, k-1 ), true );
solveF( P0, Q0, S, T, 1, pk, r, a );
TIMING_END(fac_solve);
DEBOUTLN( cerr, "trying to find correction coefficients for " << C );
DEBOUTLN( cerr, "which evaluates to " << C0 );
for ( i = 1; i <= r; i++ )
A[i] = remainder( pk( a[i] * C0 ), P0[i], pk );
DEBOUTLN( cerr, "the first approximation of the correction coefficients is " << A );
#ifdef DEBUGOUTPUT
if ( check_dummy( A, P, Q ) - C != 0 )
{
DEBOUTLN( cerr, "there is an error detected, the correction coefficients do not" );
DEBOUTLN( cerr, "fulfill equation F(A)" );
DEBOUTLN( cerr, "corresponding P " << P );
DEBOUTLN( cerr, " Q " << Q );
}
#endif
for ( m = 0; m <= h && ( m == 0 || Dm != 0 ); m++ )
{
Dm = pk( evalF( P, Q, A, r ) - C );
if ( Dm != 0 )
{
if ( k == 2 )
{
TIMING_START(fac_solve);
solveF( P0, Q0, S, T, Dm, pk, r, a );
TIMING_END(fac_solve);
for ( i = 1; i <= r; i++ )
A[i] -= a[i];
}
else
{
IteratedFor e( 2, k-1, m+1 );
while ( e.iterations_left() )
{
if ( check_e( e, k, m, n ) )
{
g = pk( checkCombination( Dm, I, e, k ) );
if ( ! g.isZero() && ! (g.mvar() > Variable(1)) )
{
prodcomb = prodCombination( I, e, k );
// Dm = Dm - g * prodcomb;
TIMING_START(fac_solve);
solveF( P0, Q0, S, T, g, pk, r, a );
TIMING_END(fac_solve);
for ( i = 1; i <= r; i++ )
{
// A[i] -= a[i] * prodcomb;
A[i] = pk( A[i] - a[i] * prodcomb );
}
}
}
e++;
}
}
}
DEBOUTLN( cerr, "the correction coefficients at step " << m );
DEBOUTLN( cerr, "are now " << A );
#ifdef DEBUGOUTPUT
if ( check_dummy( A, P, Q ) - C != 0 ) {
DEBOUTLN( cerr, "there is an error detected, the correction coefficients do not" );
DEBOUTLN( cerr, "fulfill equation F(A)" );
DEBOUTLN( cerr, "corresponding P " << P );
DEBOUTLN( cerr, " Q " << Q );
}
#endif
}
DEBDECLEVEL( cerr, "findCorrCoeffs" );
return A;
}
static bool
liftStep ( CFArray & P, int k, int r, int t, const modpk & b, const Evaluation & A, const CFArray & lcG, const CanonicalForm & Uk, int * n, int h )
{
DEBINCLEVEL( cerr, "liftStep" );
CFArray K( 1, r ), Q( 1, r ), Q0( 1, r ), P0( 1, r ), S( 1, r ), T( 1, r ), alpha( 1, r );
CanonicalForm Rm, C, D, xa = Variable(k) - A[k];
CanonicalForm xa1 = xa, xa2 = xa*xa;
CanonicalForm dummy;
int i, m;
DEBOUTLN( cerr, "we are now performing the liftstep to reach " << Variable(k) );
DEBOUTLN( cerr, "the factors so far are " << P );
DEBOUTLN( cerr, "modulus p^k= " << b.getpk() << "=" << b.getp() <<"^"<< b.getk() );
for ( i = 1; i <= r; i++ )
{
Variable vm = Variable( t + 1 );
Variable v1 = Variable(1);
K[i] = swapvar( replaceLc( swapvar( P[i], v1, vm ), swapvar( A( lcG[i], k+1, t ), v1, vm ) ), v1, vm );
P[i] = A( K[i], k, t );
}
DEBOUTLN( cerr, "lift K = " << K );
// d = degree( Uk, Variable( k ) );
TIMING_START(fac_extgcd);
Q[r] = 1;
for ( i = r; i > 1; i-- )
{
Q[i-1] = Q[i] * P[i];
P0[i] = A( P[i], 2, k-1 );
Q0[i] = A( Q[i], 2, k-1 );
extgcd( P0[i], Q0[i], S[i], T[i], b );
}
P0[1] = A( P[1], 2, k-1 );
Q0[1] = A( Q[1], 2, k-1 );
extgcd( P0[1], Q0[1], S[1], T[1], b );
TIMING_END(fac_extgcd);
for ( m = 1; m <= n[k]+1; m++ )
{
TIMING_START(fac_modpk);
Rm = modDeltak( prod( K ) - Uk, A, k, n );
TIMING_END(fac_modpk);
#ifdef DEBUGOUTPUT
if ( mod( Rm, xa1 ) != 0 )
{
DEBOUTLN( cerr, "something seems not to be ok with Rm which is " << Rm );
DEBOUTLN( cerr, "and should reduce to zero modulo " << xa1 );
}
#endif
if ( mod( Rm, xa2 ) != 0 )
{
C = derivAndEval( Rm, m, Variable( k ), A[k] );
D = 1;
for ( i = 2; i <= m; i++ ) D *= i;
C /= D;
TIMING_START(fac_corrcoeff);
alpha = findCorrCoeffs( P, Q, P0, Q0, S, T, C, A, b, r, k, h, n ); // -> h berechnen
TIMING_END(fac_corrcoeff);
// #ifdef DEBUGOUTPUT
// dummy = check_dummy( alpha, P, Q );
// if ( b(modDeltak( dummy, A, k, n )) != b(modDeltak( C, A, k, n )) )
// {
// DEBOUTLN( cerr, "lift fault" );
// DEBDECLEVEL( cerr, "liftStep" );
// return false;
// }
// #endif
for ( i = 1; i <= r; i++ )
K[i] = b(K[i] - alpha[i] * xa1);
DEBOUTLN( cerr, "the corrected K's are now " << K );
}
xa1 = xa2;
xa2 *= xa;
}
for ( i = 1; i <= r; i++ )
P[i] = K[i];
if ( prod( K ) - Uk != 0 )
{
DEBOUTLN( cerr, "the liftstep produced the wrong result" );
DEBOUTLN( cerr, "the product of the factors calculated so far is " << prod(K) );
DEBOUTLN( cerr, "and the Uk that should have been reached is " << Uk );
DEBDECLEVEL( cerr, "liftStep" );
return false;
}
DEBOUTLN( cerr, "the lift seems ok so far" );
DEBDECLEVEL( cerr, "liftStep" );
return true; // check for divisibility
}
bool
Hensel ( const CanonicalForm & U, CFArray & G, const CFArray & lcG, const Evaluation & A, const modpk & bound, const Variable & /*x*/ )
{
DEBINCLEVEL( cerr, "Hensel" );
int k, i, h, t = A.max();
bool goodeval = true;
CFArray Uk( A.min(), A.max() );
int * n = new int[t+1];
Uk[t] = U;
for ( k = t-1; k > 1; k-- )
{
Uk[k] = Uk[k+1]( A[k+1], Variable( k+1 ) );
n[k] = degree( Uk[k], Variable( k ) );
}
for ( k = A.min(); goodeval && (k <= t); k++ )
{
h = totaldegree( Uk[k], Variable(A.min()), Variable(k-1) );
for ( i = A.min(); i <= k; i++ )
n[i] = degree( Uk[k], Variable(i) );
goodeval = liftStep( G, k, G.max(), t, bound, A, lcG, Uk[k], n, h );
DEBOUTLN( cerr, "Factorization so far: " << G );
}
DEBDECLEVEL( cerr, "Hensel" );
delete[] n;
return goodeval;
}
#endif
divrem
void divrem(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r)
Definition canonicalform.cc:1032
power
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
Definition canonicalform.cc:1896
mapinto
CanonicalForm mapinto(const CanonicalForm &f)
Definition canonicalform.h:355
mod
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
degree
int degree(const CanonicalForm &f)
Definition canonicalform.h:316
CFArray
Array< CanonicalForm > CFArray
Definition canonicalform.h:397
setCharacteristic
void FACTORY_PUBLIC setCharacteristic(int c)
Definition cf_char.cc:28
swapvar
CanonicalForm FACTORY_PUBLIC swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition cf_ops.cc:168
totaldegree
int totaldegree(const CanonicalForm &f)
int totaldegree ( const CanonicalForm & f )
Definition cf_ops.cc:523
m
int m
Definition cfEzgcd.cc:128
Hensel
static int Hensel(const CanonicalForm &UU, CFArray &G, const Evaluation &AA, const CFArray &LeadCoeffs)
Definition cfEzgcd.cc:302
i
int i
Definition cfEzgcd.cc:132
k
int k
Definition cfEzgcd.cc:99
x
Variable x
Definition cfModGcd.cc:4090
p
int p
Definition cfModGcd.cc:4086
g
g
Definition cfModGcd.cc:4098
b
CanonicalForm b
Definition cfModGcd.cc:4111
tau
void tau(int **points, int sizePoints, int k)
Definition cfNewtonPolygon.cc:461
extgcd
CanonicalForm extgcd(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &a, CanonicalForm &b)
CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a,...
Definition cfUnivarGcd.cc:174
cf_binom.h
binomialpower
CanonicalForm binomialpower(const Variable &, const CanonicalForm &, int)
cf_defs.h
factory switches.
cf_eval.h
evaluate polynomials at points
cf_iter.h
Iterators for CanonicalForm's.
bound
static CanonicalForm bound(const CFMatrix &M)
Definition cf_linsys.cc:460
f
FILE * f
Definition checklibs.c:9
CFIterator
class to iterate through CanonicalForm's
Definition cf_iter.h:44
CanonicalForm
factory's main class
Definition canonicalform.h:86
CanonicalForm::isZero
CF_NO_INLINE bool isZero() const
Evaluation
class to evaluate a polynomial at points
Definition cf_eval.h:32
IteratedFor
Definition fac_iterfor.h:19
IteratedFor::iterations_left
bool iterations_left() const
Definition fac_iterfor.h:41
Variable
factory's class for variables
Definition factory.h:127
modpk
class to do operations mod p^k for int's p and k
Definition fac_util.h:23
modpk::getk
int getk() const
Definition fac_util.h:36
modpk::getp
int getp() const
Definition fac_util.h:35
debug.h
functions to print debug output
DEBINCLEVEL
#define DEBINCLEVEL(stream, msg)
Definition debug.h:44
DEBOUTLN
#define DEBOUTLN(stream, objects)
Definition debug.h:49
DEBDECLEVEL
#define DEBDECLEVEL(stream, msg)
Definition debug.h:45
alpha
Variable alpha
Definition facAbsBiFact.cc:52
result
return result
Definition facAbsBiFact.cc:76
s
const CanonicalForm int s
Definition facAbsFact.cc:51
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
j
int j
Definition facHensel.cc:110
prod
fq_nmod_poly_t prod
Definition facHensel.cc:100
fac_iterfor.h
remainder
CanonicalForm remainder(const CanonicalForm &f, const CanonicalForm &g, const modpk &pk)
Definition fac_util.cc:115
replaceLc
CanonicalForm replaceLc(const CanonicalForm &f, const CanonicalForm &c)
Definition fac_util.cc:90
fac_util.h
operations mod p^k and some other useful functions for factorization
min
static int min(int a, int b)
Definition fast_mult.cc:268
D
#define D(A)
Definition gentable.cc:128
Q0
STATIC_VAR int64 * Q0
Definition hilb.cc:68
T
STATIC_VAR jList * T
Definition janet.cc:30
G
STATIC_VAR TreeM * G
Definition janet.cc:31
h
STATIC_VAR Poly * h
Definition janet.cc:971
A
#define A
Definition sirandom.c:24
Q
#define Q
Definition sirandom.c:26
timing.h
TIMING_DEFINE_PRINT
#define TIMING_DEFINE_PRINT(t)
Definition timing.h:95
TIMING_END
#define TIMING_END(t)
Definition timing.h:93
TIMING_START
#define TIMING_START(t)
Definition timing.h:92