My Project
Loading...
Searching...
No Matches
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
24double 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
42double 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 ?
63int probIrredTest (const CanonicalForm& F, double error)
64{
65 CFMap 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
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:411
CanonicalForm num(const CanonicalForm &f)
CanonicalForm den(const CanonicalForm &f)
int FACTORY_PUBLIC getCharacteristic()
Definition: cf_char.cc:70
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
int p
Definition: cfModGcd.cc:4078
CanonicalForm compress(const CanonicalForm &f, CFMap &m)
CanonicalForm compress ( const CanonicalForm & f, CFMap & m )
Definition: cf_map.cc:210
map polynomials
generate random integers, random elements of finite fields
class CFMap
Definition: cf_map.h:85
factory's main class
Definition: canonicalform.h:86
int level() const
level() returns the level of CO.
generate random elements in F_p
Definition: cf_random.h:44
CanonicalForm generate() const
Definition: cf_random.cc:68
return result
Definition: facAbsBiFact.cc:75
const CanonicalForm int s
Definition: facAbsFact.cc:51
CanonicalForm buf1
Definition: facFqBivar.cc:73
int j
Definition: facHensel.cc:110
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
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
double inverseERF(double d)
Definition: facIrredTest.cc:42
This file provides a probabilistic irreducibility test for polynomials over Z/p.
STATIC_VAR TreeM * G
Definition: janet.cc:31
#define pi
Definition: libparse.cc:1145
gmp_float sqrt(const gmp_float &a)
Definition: mpr_complex.cc:327
gmp_float log(const gmp_float &a)
Definition: mpr_complex.cc:343
#define error(a)
Definition: mpr_numeric.cc:966
int status int void * buf
Definition: si_signals.h:59