source: git/factory/facIrredTest.cc @ 4a5e77

fieker-DuValspielwiese
Last change on this file since 4a5e77 was 4a5e77, checked in by Martin Lee <martinlee84@…>, 13 years ago
added probabilistic irreducibility test, available via system ("probIrredTest",poly,"error") git-svn-id: file:///usr/local/Singular/svn/trunk@14042 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100755
File size: 2.3 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 p= getCharacteristic();
27  int result= 0;
28
29  FFRandom FFgen;
30  CanonicalForm buf;
31  for (int i= 0; i < k; i++)
32  {
33    buf= F;
34    for (int j= F.level(); j > 0; j++) 
35      buf= buf (FFgen.generate(), j);
36    if (buf.isZero())
37      result++;
38  }
39
40  return (double) result/k;
41}
42
43double inverseERF (double d)
44{
45  double pi= 3.141592653589793;
46  double a= 0.140012288;
47
48  double buf1;
49  double buf= 2.0/(pi*a)+log (1.0-d*d)/2.0;
50  buf1= buf;
51  buf *= buf;
52  buf -= log (1.0-d*d)/a;
53  buf= sqrt (buf);
54  buf -= buf1;
55  buf= sqrt (buf);
56
57  if (d < 0)
58    buf= -buf;
59
60  return buf;
61}
62
63//doesn't make much sense to use this if getCharacteristic() > 17 ?
64int probIrredTest (const CanonicalForm& F, double error)
65{
66  CFMap N;
67  CanonicalForm G= compress (F, N);
68  int n= G.level();
69  int p= getCharacteristic();
70
71  double sqrtTrials= inverseERF (1-2.0*error)*sqrt (2.0);
72
73  double s= sqrtTrials;
74
75  double pn= pow ((double) p, (double) n);
76  double p1= (double) 1/p;
77  p1= p1*(1.0-p1);
78  p1= p1/(double) pn;
79  p1= sqrt (p1);
80  p1 *= s;
81  p1 += (double) 1/p;
82
83  double p2= (double) (2*p-1)/(p*p);
84  p2= p2*(1-p2);
85  p2= p2/(double) pn;
86  p2= sqrt (p2);
87  p2 *= s;
88  p2= (double) (2*p - 1)/(p*p)-p2;
89
90  //no testing possible
91  if (p2 < p1)
92    return 0;
93
94  double den= sqrt (p1*(1-p1))+sqrt (p2*(1-p2));
95  double num= p2-p1;
96
97  sqrtTrials *= den/num;
98
99  int trials= (int) floor (pow (sqrtTrials, 2.0));
100
101  double l= 2.0;
102  double experimentalNumZeros= numZeros (G, trials);
103
104  double pmiddle= sqrt (p1*p2);
105
106  num= den;
107  den= sqrt (p1*(1.0-p2))+sqrt (p2*(1.0-p1));
108  pmiddle *= (den/num);
109
110  if (experimentalNumZeros < pmiddle)
111    return 1;
112  else
113    return -1;
114}
115
Note: See TracBrowser for help on using the repository browser.