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 | #include "config.h" |
---|
15 | #include <math.h> |
---|
16 | |
---|
17 | #include "facIrredTest.h" |
---|
18 | #include "cf_map.h" |
---|
19 | #include "cf_random.h" |
---|
20 | |
---|
21 | //returns 0 if number of zeros/k >= l |
---|
22 | double numZeros (const CanonicalForm& F, int k) |
---|
23 | { |
---|
24 | int result= 0; |
---|
25 | |
---|
26 | FFRandom FFgen; |
---|
27 | CanonicalForm buf; |
---|
28 | for (int i= 0; i < k; i++) |
---|
29 | { |
---|
30 | buf= F; |
---|
31 | for (int j= F.level(); j > 0; j++) |
---|
32 | buf= buf (FFgen.generate(), j); |
---|
33 | if (buf.isZero()) |
---|
34 | result++; |
---|
35 | } |
---|
36 | |
---|
37 | return (double) result/k; |
---|
38 | } |
---|
39 | |
---|
40 | double inverseERF (double d) |
---|
41 | { |
---|
42 | double pi= 3.141592653589793; |
---|
43 | double a= 0.140012288; |
---|
44 | |
---|
45 | double buf1; |
---|
46 | double buf= 2.0/(pi*a)+log (1.0-d*d)/2.0; |
---|
47 | buf1= buf; |
---|
48 | buf *= buf; |
---|
49 | buf -= log (1.0-d*d)/a; |
---|
50 | buf= sqrt (buf); |
---|
51 | buf -= buf1; |
---|
52 | buf= sqrt (buf); |
---|
53 | |
---|
54 | if (d < 0) |
---|
55 | buf= -buf; |
---|
56 | |
---|
57 | return buf; |
---|
58 | } |
---|
59 | |
---|
60 | //doesn't make much sense to use this if getCharacteristic() > 17 ? |
---|
61 | int probIrredTest (const CanonicalForm& F, double error) |
---|
62 | { |
---|
63 | CFMap N; |
---|
64 | CanonicalForm G= compress (F, N); |
---|
65 | int n= G.level(); |
---|
66 | int p= getCharacteristic(); |
---|
67 | |
---|
68 | double sqrtTrials= inverseERF (1-2.0*error)*sqrt (2.0); |
---|
69 | |
---|
70 | double s= sqrtTrials; |
---|
71 | |
---|
72 | double pn= pow ((double) p, (double) n); |
---|
73 | double p1= (double) 1/p; |
---|
74 | p1= p1*(1.0-p1); |
---|
75 | p1= p1/(double) pn; |
---|
76 | p1= sqrt (p1); |
---|
77 | p1 *= s; |
---|
78 | p1 += (double) 1/p; |
---|
79 | |
---|
80 | double p2= (double) (2*p-1)/(p*p); |
---|
81 | p2= p2*(1-p2); |
---|
82 | p2= p2/(double) pn; |
---|
83 | p2= sqrt (p2); |
---|
84 | p2 *= s; |
---|
85 | p2= (double) (2*p - 1)/(p*p)-p2; |
---|
86 | |
---|
87 | //no testing possible |
---|
88 | if (p2 < p1) |
---|
89 | return 0; |
---|
90 | |
---|
91 | double den= sqrt (p1*(1-p1))+sqrt (p2*(1-p2)); |
---|
92 | double num= p2-p1; |
---|
93 | |
---|
94 | sqrtTrials *= den/num; |
---|
95 | |
---|
96 | int trials= (int) floor (pow (sqrtTrials, 2.0)); |
---|
97 | |
---|
98 | double experimentalNumZeros= numZeros (G, trials); |
---|
99 | |
---|
100 | double pmiddle= sqrt (p1*p2); |
---|
101 | |
---|
102 | num= den; |
---|
103 | den= sqrt (p1*(1.0-p2))+sqrt (p2*(1.0-p1)); |
---|
104 | pmiddle *= (den/num); |
---|
105 | |
---|
106 | if (experimentalNumZeros < pmiddle) |
---|
107 | return 1; |
---|
108 | else |
---|
109 | return -1; |
---|
110 | } |
---|
111 | |
---|