My Project  debian-1:4.1.1-p2+ds-4
facIrredTest.cc
Go to the documentation of this file.
1 /*****************************************************************************\
2  * Computer Algebra System SINGULAR
3 \*****************************************************************************/
4 /** @file facIrredTest.cc
5  *
6  * This file implements a probabilistic irreducibility test for polynomials over
7  * Z/p.
8  *
9  * @author Martin Lee
10  *
11  **/
12 /*****************************************************************************/
13 
14 
15 #include "config.h"
16 
17 #include <cmath>
18 
19 #include "facIrredTest.h"
20 #include "cf_map.h"
21 #include "cf_random.h"
22 
23 //returns 0 if number of zeros/k >= l
24 double numZeros (const CanonicalForm& F, int k)
25 {
26  int result= 0;
27 
28  FFRandom FFgen;
30  for (int i= 0; i < k; i++)
31  {
32  buf= F;
33  for (int j= F.level(); j > 0; j++)
34  buf= buf (FFgen.generate(), j);
35  if (buf.isZero())
36  result++;
37  }
38 
39  return (double) result/k;
40 }
41 
42 double inverseERF (double d)
43 {
44  double pi= 3.141592653589793;
45  double a= 0.140012288;
46 
47  double buf1;
48  double buf= 2.0/(pi*a)+log (1.0-d*d)/2.0;
49  buf1= buf;
50  buf *= buf;
51  buf -= log (1.0-d*d)/a;
52  buf= sqrt (buf);
53  buf -= buf1;
54  buf= sqrt (buf);
55 
56  if (d < 0)
57  buf= -buf;
58 
59  return buf;
60 }
61 
62 //doesn't make much sense to use this if getCharacteristic() > 17 ?
63 int probIrredTest (const CanonicalForm& F, double error)
64 {
65  CFMap N;
66  CanonicalForm G= compress (F, N);
67  int n= G.level();
68  int p= getCharacteristic();
69 
70  double sqrtTrials= inverseERF (1-2.0*error)*sqrt (2.0);
71 
72  double s= sqrtTrials;
73 
74  double pn= pow ((double) p, (double) n);
75  double p1= (double) 1/p;
76  p1= p1*(1.0-p1);
77  p1= p1/(double) pn;
78  p1= sqrt (p1);
79  p1 *= s;
80  p1 += (double) 1/p;
81 
82  double p2= (double) (2*p-1)/(p*p);
83  p2= p2*(1-p2);
84  p2= p2/(double) pn;
85  p2= sqrt (p2);
86  p2 *= s;
87  p2= (double) (2*p - 1)/(p*p)-p2;
88 
89  //no testing possible
90  if (p2 < p1)
91  return 0;
92 
93  double den= sqrt (p1*(1-p1))+sqrt (p2*(1-p2));
94  double num= p2-p1;
95 
96  sqrtTrials *= den/num;
97 
98  int trials= (int) floor (pow (sqrtTrials, 2.0));
99 
100  double experimentalNumZeros= numZeros (G, trials);
101 
102  double pmiddle= sqrt (p1*p2);
103 
104  num= den;
105  den= sqrt (p1*(1.0-p2))+sqrt (p2*(1.0-p1));
106  pmiddle *= (den/num);
107 
108  if (experimentalNumZeros < pmiddle)
109  return 1;
110  else
111  return -1;
112 }
113 
error
void error(const char *fmt,...)
Definition: emacs.cc:54
FFRandom
generate random elements in F_p
Definition: cf_random.h:43
j
int j
Definition: facHensel.cc:105
k
int k
Definition: cfEzgcd.cc:92
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
result
return result
Definition: facAbsBiFact.cc:76
num
CanonicalForm num(const CanonicalForm &f)
Definition: canonicalform.h:330
CFMap
class CFMap
Definition: cf_map.h:84
sqrt
gmp_float sqrt(const gmp_float &a)
Definition: mpr_complex.cc:327
amp::floor
const signed long floor(const ampf< Precision > &x)
Definition: amp.h:774
N
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:48
cf_random.h
getCharacteristic
int getCharacteristic()
Definition: cf_char.cc:51
CanonicalForm
factory's main class
Definition: canonicalform.h:77
pi
#define pi
Definition: libparse.cc:1142
i
int i
Definition: cfEzgcd.cc:125
buf
int status int void * buf
Definition: si_signals.h:58
compress
CanonicalForm compress(const CanonicalForm &f, CFMap &m)
CanonicalForm compress ( const CanonicalForm & f, CFMap & m )
Definition: cf_map.cc:210
buf1
CanonicalForm buf1
Definition: facFqBivar.cc:71
den
CanonicalForm den(const CanonicalForm &f)
Definition: canonicalform.h:333
inverseERF
double inverseERF(double d)
Definition: facIrredTest.cc:42
cf_map.h
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
log
gmp_float log(const gmp_float &a)
Definition: mpr_complex.cc:343
pow
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:414
FFRandom::generate
CanonicalForm generate() const
Definition: cf_random.cc:56
p
int p
Definition: cfModGcd.cc:4019
s
const CanonicalForm int s
Definition: facAbsFact.cc:55
G
static TreeM * G
Definition: janet.cc:32
CanonicalForm::level
int level() const
level() returns the level of CO.
Definition: canonicalform.cc:543
facIrredTest.h