source: git/factory/facIrredTest.cc @ 8a30b1

spielwiese
Last change on this file since 8a30b1 was 72bfc8, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: deleted @internal
  • Property mode set to 100755
File size: 2.2 KB
Line 
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
22double 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
40double 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 ?
61int 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
Note: See TracBrowser for help on using the repository browser.