Loading...
Searching...
No Matches
facIrredTest.cc
Go to the documentation of this file.
/*****************************************************************************\
* Computer Algebra System SINGULAR
\*****************************************************************************/
/** @file facIrredTest.cc
*
* This file implements a probabilistic irreducibility test for polynomials over
* Z/p.
*
* @author Martin Lee
*
**/
/*****************************************************************************/
#include "config.h"
#include <cmath>
#include "facIrredTest.h"
#include "cf_map.h"
#include "cf_random.h"
//returns 0 if number of zeros/k >= l
double numZeros (const CanonicalForm& F, int k)
{
int result= 0;
FFRandom FFgen;
CanonicalForm buf;
for (int i= 0; i < k; i++)
{
buf= F;
for (int j= F.level(); j > 0; j++)
buf= buf (FFgen.generate(), j);
if (buf.isZero())
result++;
}
return (double) result/k;
}
double inverseERF (double d)
{
double pi= 3.141592653589793;
double a= 0.140012288;
double buf1;
double buf= 2.0/(pi*a)+log (1.0-d*d)/2.0;
buf1= buf;
buf *= buf;
buf -= log (1.0-d*d)/a;
buf= sqrt (buf);
buf -= buf1;
buf= sqrt (buf);
if (d < 0)
buf= -buf;
return buf;
}
//doesn't make much sense to use this if getCharacteristic() > 17 ?
int probIrredTest (const CanonicalForm& F, double error)
{
CFMap N;
CanonicalForm G= compress (F, N);
int n= G.level();
int p= getCharacteristic();
double sqrtTrials= inverseERF (1-2.0*error)*sqrt (2.0);
double s= sqrtTrials;
double pn= pow ((double) p, (double) n);
double p1= (double) 1/p;
p1= p1*(1.0-p1);
p1= p1/(double) pn;
p1= sqrt (p1);
p1 *= s;
p1 += (double) 1/p;
double p2= (double) (2*p-1)/(p*p);
p2= p2*(1-p2);
p2= p2/(double) pn;
p2= sqrt (p2);
p2 *= s;
p2= (double) (2*p - 1)/(p*p)-p2;
//no testing possible
if (p2 < p1)
return 0;
double den= sqrt (p1*(1-p1))+sqrt (p2*(1-p2));
double num= p2-p1;
sqrtTrials *= den/num;
int trials= (int) floor (pow (sqrtTrials, 2.0));
double experimentalNumZeros= numZeros (G, trials);
double pmiddle= sqrt (p1*p2);
num= den;
den= sqrt (p1*(1.0-p2))+sqrt (p2*(1.0-p1));
pmiddle *= (den/num);
if (experimentalNumZeros < pmiddle)
return 1;
else
return -1;
}
pow
Rational pow(const Rational &a, int e)
Definition GMPrat.cc:411
num
CanonicalForm num(const CanonicalForm &f)
Definition canonicalform.h:337
den
CanonicalForm den(const CanonicalForm &f)
Definition canonicalform.h:340
getCharacteristic
int FACTORY_PUBLIC getCharacteristic()
Definition cf_char.cc:70
N
const CanonicalForm CFMap CFMap & N
Definition cfEzgcd.cc:56
i
int i
Definition cfEzgcd.cc:132
k
int k
Definition cfEzgcd.cc:99
p
int p
Definition cfModGcd.cc:4086
compress
CanonicalForm compress(const CanonicalForm &f, CFMap &m)
CanonicalForm compress ( const CanonicalForm & f, CFMap & m ).
Definition cf_map.cc:210
cf_map.h
map polynomials
cf_random.h
generate random integers, random elements of finite fields
CFMap
class CFMap
Definition cf_map.h:85
CanonicalForm
factory's main class
Definition canonicalform.h:86
CanonicalForm::level
int level() const
level() returns the level of CO.
Definition canonicalform.cc:576
FFRandom
generate random elements in F_p
Definition cf_random.h:44
FFRandom::generate
CanonicalForm generate() const
Definition cf_random.cc:65
result
return result
Definition facAbsBiFact.cc:76
s
const CanonicalForm int s
Definition facAbsFact.cc:51
buf1
CanonicalForm buf1
Definition facFqBivar.cc:76
j
int j
Definition facHensel.cc:110
numZeros
double numZeros(const CanonicalForm &F, int k)
evaluate F at k random points in Z/p^n and count the number of zeros that occur
Definition facIrredTest.cc:24
probIrredTest
int probIrredTest(const CanonicalForm &F, double error)
given some error probIrredTest detects irreducibility or reducibility of F with confidence level 1-er...
Definition facIrredTest.cc:63
inverseERF
double inverseERF(double d)
Definition facIrredTest.cc:42
facIrredTest.h
This file provides a probabilistic irreducibility test for polynomials over Z/p.
G
STATIC_VAR TreeM * G
Definition janet.cc:31
pi
#define pi
Definition libparse.cc:1145
sqrt
gmp_float sqrt(const gmp_float &a)
Definition mpr_complex.cc:319
log
gmp_float log(const gmp_float &a)
Definition mpr_complex.cc:335
error
#define error(a)
Definition mpr_numeric.cc:966
buf
int status int void * buf
Definition si_signals.h:69