source: git/factory/facIrredTest.cc @ e4fe2b

jengelh-datetimespielwiese
Last change on this file since e4fe2b was e4fe2b, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
FIX: Fixed huge BUG in cf_gmp.h CHG: starting to cleanup factory
  • 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 * @internal @version \$Id$
12 *
13 **/
14/*****************************************************************************/
15
16#include "config.h"
17#include <math.h>
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;
29  CanonicalForm buf;
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;
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
Note: See TracBrowser for help on using the repository browser.