Loading...
Searching...
No Matches
cf_chinese.cc File Reference

algorithms for chinese remaindering and rational reconstruction More...

`#include "config.h"`
`#include "NTLconvert.h"`
`#include "FLINTconvert.h"`
`#include "cf_assert.h"`
`#include "debug.h"`
`#include "canonicalform.h"`
`#include "cf_iter.h"`
`#include "cf_util.h"`

Go to the source code of this file.

## Functions

void | chineseRemainder (const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew)
| void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
void | chineseRemainder (const CFArray &x, const CFArray &q, CanonicalForm &xnew, CanonicalForm &qnew)
| void chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
CanonicalForm | Farey (const CanonicalForm &f, const CanonicalForm &q)
| Farey rational reconstruction.
static CanonicalForm | chin_mul_inv (const CanonicalForm a, const CanonicalForm b, int ind, CFArray &inv)
void | out_cf (const char *s1, const CanonicalForm &f, const char *s2)
void | chineseRemainderCached (const CFArray &a, const CFArray &n, CanonicalForm &xnew, CanonicalForm &prod, CFArray &inv)
void | chineseRemainderCached (const CanonicalForm &a, const CanonicalForm &q1, const CanonicalForm &b, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew, CFArray &inv)

## Detailed Description

algorithms for chinese remaindering and rational reconstruction

Used by: cf_gcd.cc, cf_linsys.cc

Header file: cf_algorithm.h

Definition in file cf_chinese.cc.

## Function Documentation

## chin_mul_inv()

CanonicalForm chin_mul_inv | ( | const CanonicalForm | a,
|  | const CanonicalForm | b,
|  | int | ind,
|  | CFArray & | inv )

inline static

Definition at line 276 of file cf_chinese.cc.

{
if (inv[ind].isZero())
{
CanonicalForm s,dummy;
(void)bextgcd( a, b, s, dummy );
inv[ind]=s;
return s;
}
else
return inv[ind];
}
bextgcd
CanonicalForm bextgcd(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &a, CanonicalForm &b)
CanonicalForm bextgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a,...
Definition canonicalform.cc:1722
b
CanonicalForm b
Definition cfModGcd.cc:4111
CanonicalForm
factory's main class
Definition canonicalform.h:86
s
const CanonicalForm int s
Definition facAbsFact.cc:51
isZero
bool isZero(const CFArray &A)
checks if entries of A are zero
Definition facSparseHensel.h:468

## chineseRemainder() [1/2]

void chineseRemainder | ( | const CanonicalForm & | x1,
|  | const CanonicalForm & | q1,
|  | const CanonicalForm & | x2,
|  | const CanonicalForm & | q2,
|  | CanonicalForm & | xnew,
|  | CanonicalForm & | qnew )

void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )

chineseRemainder - integer chinese remaindering.

Calculate xnew such that xnew=x1 (mod q1) and xnew=x2 (mod q2) and qnew = q1*q2. q1 and q2 should be positive integers, pairwise prime, x1 and x2 should be polynomials with integer coefficients. If x1 and x2 are polynomials with positive coefficients, the result is guaranteed to have positive coefficients, too.

Note: This algorithm is optimized for the case q1>>q2.

This is a standard algorithm. See, for example, Geddes/Czapor/Labahn - 'Algorithms for Computer Algebra', par. 5.6 and 5.8, or the article of M. Lauer - 'Computing by Homomorphic Images' in B. Buchberger - 'Computer Algebra - Symbolic and Algebraic Computation'.

Note: Be sure you are calculating in Z, and not in Q!

Definition at line 57 of file cf_chinese.cc.

{
DEBINCLEVEL( cerr, "chineseRemainder" );
DEBOUTLN( cerr, "log(q1) = " << q1.ilog2() );
DEBOUTLN( cerr, "log(q2) = " << q2.ilog2() );
// We calculate xnew as follows:
// xnew = v1 + v2 * q1
// where
// v1 = x1 (mod q1)
// v2 = (x2-v1)/q1 (mod q2) (*)
//
// We do one extra test to check whether x2-v1 vanishes (mod
// q2) in (*) since it is not costly and may save us
// from calculating the inverse of q1 (mod q2).
//
// u: v1 (mod q2)
// d: x2-v1 (mod q2)
// s: 1/q1 (mod q2)
//
CanonicalForm v2, v1;
CanonicalForm u, d, s, dummy;
v1 = mod( x1, q1 );
u = mod( v1, q2 );
d = mod( x2-u, q2 );
if ( d.isZero() )
{
xnew = v1;
qnew = q1 * q2;
DEBDECLEVEL( cerr, "chineseRemainder" );
return;
}
(void)bextgcd( q1, q2, s, dummy );
v2 = mod( d*s, q2 );
xnew = v1 + v2*q1;
// After all, calculate new modulus. It is important that
// this is done at the very end of the algorithm, since q1
// and qnew may refer to the same object (same is true for x1
// and xnew).
qnew = q1 * q2;
DEBDECLEVEL( cerr, "chineseRemainder" );
}
mod
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
CanonicalForm::isZero
CF_NO_INLINE bool isZero() const
CanonicalForm::ilog2
int ilog2() const
int CanonicalForm::ilog2 () const
Definition canonicalform.cc:1417
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

## chineseRemainder() [2/2]

void chineseRemainder | ( | const CFArray & | x,
|  | const CFArray & | q,
|  | CanonicalForm & | xnew,
|  | CanonicalForm & | qnew )

void chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )

chineseRemainder - integer chinese remaindering.

Calculate xnew such that xnew=x[i] (mod q[i]) and qnew is the product of all q[i]. q[i] should be positive integers, pairwise prime. x[i] should be polynomials with integer coefficients. If all coefficients of all x[i] are positive integers, the result is guaranteed to have positive coefficients, too.

This is a standard algorithm, too, except for the fact that we use a divide-and-conquer method instead of a linear approach to calculate the remainder.

Note: Be sure you are calculating in Z, and not in Q!

Definition at line 124 of file cf_chinese.cc.

{
DEBINCLEVEL( cerr, "chineseRemainder( ... CFArray ... )" );
ASSERT( x.min() == q.min() && x.size() == q.size(), "incompatible arrays" );
CFArray X(x), Q(q);
int i, j, n = x.size(), start = x.min();
DEBOUTLN( cerr, "array size = " << n );
while ( n != 1 )
{
i = j = start;
while ( i < start + n - 1 )
{
// This is a little bit dangerous: X[i] and X[j] (and
// Q[i] and Q[j]) may refer to the same object. But
// xnew and qnew in the above function are modified
// at the very end of the function, so we do not
// modify x1 and q1, resp., by accident.
chineseRemainder( X[i], Q[i], X[i+1], Q[i+1], X[j], Q[j] );
i += 2;
j++;
}
if ( n & 1 )
{
X[j] = X[i];
Q[j] = Q[i];
}
// Maybe we would get some memory back at this point if
// we would set X[j+1, ..., n] and Q[j+1, ..., n] to zero
// at this point?
n = ( n + 1) / 2;
}
xnew = X[start];
qnew = Q[q.min()];
DEBDECLEVEL( cerr, "chineseRemainder( ... CFArray ... )" );
}
CFArray
Array< CanonicalForm > CFArray
Definition canonicalform.h:397
i
int i
Definition cfEzgcd.cc:132
x
Variable x
Definition cfModGcd.cc:4090
ASSERT
#define ASSERT(expression, message)
Definition cf_assert.h:99
chineseRemainder
void chineseRemainder(const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew)
void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2,...
Definition cf_chinese.cc:57
Array::size
int size() const
Definition ftmpl_array.cc:92
Array::min
int min() const
Definition ftmpl_array.cc:98
j
int j
Definition facHensel.cc:110
Q
#define Q
Definition sirandom.c:26

## chineseRemainderCached() [1/2]

void chineseRemainderCached | ( | const CanonicalForm & | a,
|  | const CanonicalForm & | q1,
|  | const CanonicalForm & | b,
|  | const CanonicalForm & | q2,
|  | CanonicalForm & | xnew,
|  | CanonicalForm & | qnew,
|  | CFArray & | inv )

Definition at line 308 of file cf_chinese.cc.

{
CFArray A(2); A[0]=a; A[1]=b;
CFArray Q(2); Q[0]=q1; Q[1]=q2;
chineseRemainderCached(A,Q,xnew,qnew,inv);
}
chineseRemainderCached
void chineseRemainderCached(const CFArray &a, const CFArray &n, CanonicalForm &xnew, CanonicalForm &prod, CFArray &inv)
Definition cf_chinese.cc:290
A
#define A
Definition sirandom.c:24

## chineseRemainderCached() [2/2]

void chineseRemainderCached | ( | const CFArray & | a,
|  | const CFArray & | n,
|  | CanonicalForm & | xnew,
|  | CanonicalForm & | prod,
|  | CFArray & | inv )

Definition at line 290 of file cf_chinese.cc.

{
CanonicalForm p, sum=0L; prod=1L;
int i;
int len=n.size();
for (i = 0; i < len; i++) prod *= n[i];
for (i = 0; i < len; i++)
{
p = prod / n[i];
sum += a[i] * chin_mul_inv(p, n[i], i, inv) * p;
}
xnew = mod(sum , prod);
}
p
int p
Definition cfModGcd.cc:4086
chin_mul_inv
static CanonicalForm chin_mul_inv(const CanonicalForm a, const CanonicalForm b, int ind, CFArray &inv)
Definition cf_chinese.cc:276
prod
fq_nmod_poly_t prod
Definition facHensel.cc:100

## Farey()

CanonicalForm Farey | ( | const CanonicalForm & | f,
|  | const CanonicalForm & | q )

Farey rational reconstruction.

If NTL is available it uses the fast algorithm from NTL, i.e. Encarnacion, Collins.

Definition at line 202 of file cf_chinese.cc.

{
int is_rat=isOn(SW_RATIONAL);
Off(SW_RATIONAL);
Variable x = f.mvar();
CanonicalForm result = 0;
CanonicalForm c;
CFIterator i;
#ifdef HAVE_FLINT
fmpz_t FLINTq;
fmpz_init(FLINTq);
convertCF2initFmpz(FLINTq,q);
fmpz_t FLINTc;
fmpz_init(FLINTc);
fmpq_t FLINTres;
fmpq_init(FLINTres);
#elif defined(HAVE_NTL)
ZZ NTLq= convertFacCF2NTLZZ (q);
ZZ bound;
SqrRoot (bound, NTLq/2);
#else
factoryError("NTL/FLINT missing:Farey");
#endif
for ( i = f; i.hasTerms(); i++ )
{
c = i.coeff();
if ( c.inCoeffDomain())
{
#ifdef HAVE_FLINT
if (c.inZ())
{
convertCF2initFmpz(FLINTc,c);
fmpq_reconstruct_fmpz(FLINTres,FLINTc,FLINTq);
result += power (x, i.exp())*(convertFmpq2CF(FLINTres));
}
#elif defined(HAVE_NTL)
if (c.inZ())
{
ZZ NTLc= convertFacCF2NTLZZ (c);
bool lessZero= (sign (NTLc) == -1);
if (lessZero)
NTL::negate (NTLc, NTLc);
ZZ NTLnum, NTLden;
if (ReconstructRational (NTLnum, NTLden, NTLc, NTLq, bound, bound))
{
if (lessZero)
NTL::negate (NTLnum, NTLnum);
CanonicalForm num= convertNTLZZX2CF (to_ZZX (NTLnum), Variable (1));
CanonicalForm den= convertNTLZZX2CF (to_ZZX (NTLden), Variable (1));
On (SW_RATIONAL);
result += power (x, i.exp())*(num/den);
Off (SW_RATIONAL);
}
}
#else
if (c.inZ())
result += power (x, i.exp()) * Farey_n(c,q);
#endif
else
result += power( x, i.exp() ) * Farey(c,q);
}
else
result += power( x, i.exp() ) * Farey(c,q);
}
if (is_rat) On(SW_RATIONAL);
#ifdef HAVE_FLINT
fmpq_clear(FLINTres);
fmpz_clear(FLINTc);
fmpz_clear(FLINTq);
#endif
return result;
}
convertFmpq2CF
CanonicalForm convertFmpq2CF(const fmpq_t q)
conversion of a FLINT rational to CanonicalForm
Definition FLINTconvert.cc:260
convertCF2initFmpz
void convertCF2initFmpz(fmpz_t result, const CanonicalForm &f)
conversion of a factory integer to fmpz_t(init.)
Definition FLINTconvert.cc:133
convertNTLZZX2CF
CanonicalForm convertNTLZZX2CF(const ZZX &polynom, const Variable &x)
Definition NTLconvert.cc:285
convertFacCF2NTLZZ
ZZ convertFacCF2NTLZZ(const CanonicalForm &f)
NAME: convertFacCF2NTLZZX.
Definition NTLconvert.cc:670
isOn
bool isOn(int sw)
switches
Definition canonicalform.cc:1971
On
void On(int sw)
switches
Definition canonicalform.cc:1957
Off
void Off(int sw)
switches
Definition canonicalform.cc:1964
power
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
Definition canonicalform.cc:1896
num
CanonicalForm num(const CanonicalForm &f)
Definition canonicalform.h:337
den
CanonicalForm den(const CanonicalForm &f)
Definition canonicalform.h:340
Farey
CanonicalForm Farey(const CanonicalForm &f, const CanonicalForm &q)
Farey rational reconstruction.
Definition cf_chinese.cc:202
SW_RATIONAL
static const int SW_RATIONAL
set to 1 for computations over Q
Definition cf_defs.h:31
bound
static CanonicalForm bound(const CFMatrix &M)
Definition cf_linsys.cc:460
factoryError
VAR void(* factoryError)(const char *s)
Definition cf_util.cc:80
f
FILE * f
Definition checklibs.c:9
CFIterator
class to iterate through CanonicalForm's
Definition cf_iter.h:44
CanonicalForm::inCoeffDomain
bool inCoeffDomain() const
Definition canonicalform.cc:122
CanonicalForm::inZ
bool inZ() const
predicates
Definition canonicalform.cc:69
Variable
factory's class for variables
Definition factory.h:127
result
return result
Definition facAbsBiFact.cc:76
sign
static int sign(int x)
Definition ring.cc:3504

## out_cf()

void out_cf | ( | const char * | s1,
|  | const CanonicalForm & | f,
|  | const char * | s2 )

Definition at line 103 of file cf_factor.cc.

{
printf("%s",s1);
if (f.isZero()) printf("+0");
//else if (! f.inCoeffDomain() )
else if (! f.inBaseDomain() )
{
int l = f.level();
for ( CFIterator i = f; i.hasTerms(); i++ )
{
int e=i.exp();
if (i.coeff().isOne())
{
printf("+");
if (e==0) printf("1");
else
{
printf("%c",'a'+l-1);
if (e!=1) printf("^%d",e);
}
}
else
{
out_cf("+(",i.coeff(),")");
if (e!=0)
{
printf("*%c",'a'+l-1);
if (e!=1) printf("^%d",e);
}
}
}
}
else
{
if ( f.isImm() )
{
if (CFFactory::gettype()==GaloisFieldDomain)
{
long a= imm2int (f.getval());
if ( a == gf_q )
printf ("+%ld", a);
else if ( a == 0L )
printf ("+1");
else if ( a == 1L )
printf ("+%c",gf_name);
else
{
printf ("+%c",gf_name);
printf ("^%ld",a);
}
}
else
{
long l=f.intval();
if (l<0) printf("%ld",l);
else printf("+%ld",l);
}
}
else
{
#ifdef NOSTREAMIO
if (f.inZ())
{
mpz_t m;
gmp_numerator(f,m);
char * str = new char[mpz_sizeinbase( m, 10 ) + 2];
str = mpz_get_str( str, 10, m );
puts(str);
delete[] str;
mpz_clear(m);
}
else if (f.inQ())
{
mpz_t m;
gmp_numerator(f,m);
char * str = new char[mpz_sizeinbase( m, 10 ) + 2];
str = mpz_get_str( str, 10, m );
while(str[strlen(str)]<' ') { str[strlen(str)]='\0'; }
puts(str);putchar('/');
delete[] str;
mpz_clear(m);
gmp_denominator(f,m);
str = new char[mpz_sizeinbase( m, 10 ) + 2];
str = mpz_get_str( str, 10, m );
while(str[strlen(str)]<' ') { str[strlen(str)]='\0'; }
puts(str);
delete[] str;
mpz_clear(m);
}
#else
std::cout << f;
#endif
}
//if (f.inZ()) printf("(Z)");
//else if (f.inQ()) printf("(Q)");
//else if (f.inFF()) printf("(FF)");
//else if (f.inPP()) printf("(PP)");
//else if (f.inGF()) printf("(PP)");
//else
if (f.inExtension()) printf("E(%d)",f.level());
}
printf("%s",s2);
}
l
int l
Definition cfEzgcd.cc:100
m
int m
Definition cfEzgcd.cc:128
GaloisFieldDomain
#define GaloisFieldDomain
Definition cf_defs.h:18
out_cf
void out_cf(const char *s1, const CanonicalForm &f, const char *s2)
Definition cf_factor.cc:103
CFFactory::gettype
static int gettype()
Definition cf_factory.h:28
gmp_numerator
void FACTORY_PUBLIC gmp_numerator(const CanonicalForm &f, mpz_ptr result)
Definition singext.cc:20
gmp_denominator
void FACTORY_PUBLIC gmp_denominator(const CanonicalForm &f, mpz_ptr result)
Definition singext.cc:40
gf_q
VAR int gf_q
Definition gfops.cc:47
gf_name
VAR char gf_name
Definition gfops.cc:52
imm2int
static long imm2int(const InternalCF *const imm)
Definition imm.h:70
LibThread::str
char * str(leftv arg)
Definition shared.cc:699