source: git/factory/cf_cyclo.cc @ 82db91a

fieker-DuValspielwiese
Last change on this file since 82db91a was 82db91a, checked in by Niels Ranosch <ranosch@…>, 13 years ago
FIX: function-macros should act as functions even if "empty", esp. in if-then block.
  • Property mode set to 100644
File size: 3.6 KB
Line 
1// -*- c++ -*-
2//*****************************************************************************
3/** @file cf_cyclo.cc
4 *
5 * @author Martin Lee
6 * @date 29.01.2010
7 *
8 * Compute cyclotomic polynomials and factorize integers by brute force
9 *
10 * @par Copyright:
11 *   (c) by The SINGULAR Team, see LICENSE file
12 *
13 * @internal
14 * @version \$Id$
15 *
16**/
17//*****************************************************************************
18
19#include "config.h"
20
21#include "canonicalform.h"
22#include "cf_primes.h"
23#include "cf_util.h"
24#include "cf_assert.h"
25
26#ifdef HAVE_NTL
27#include <NTL/ZZ.h>
28#endif
29
30
31/// integer factorization using table look-ups,
32/// function may fail if integer contains primes which exceed the largest prime
33/// in our table
34int* integerFactorizer (const long integer, int& length, bool& fail)
35{
36  ASSERT (integer != 0 && integer != 1 && integer != -1,
37          "non-zero non-unit expected");
38  int* result=NULL;
39  length= 0;
40  fail= false;
41  int i= integer;
42  if (integer < 0)
43    i = -integer;
44
45  int exp= 0;
46  while ((i != 1) && (i%2 == 0))
47  {
48    i /= 2;
49    exp++;
50  }
51  if (exp != 0)
52  {
53    result= new int [exp];
54    for (int k= 0; k < exp; k++)
55      result[k]= 2;
56    length += exp;
57  }
58  if (i == 1) return result;
59
60  long j= 0;
61  exp= 0;
62  int* buf;
63  int next_prime;
64  while ((i != 1) && (j < 31937))
65  {
66    next_prime= cf_getPrime (j);
67    while ((i != 1) && (i%next_prime == 0))
68    {
69      i /= next_prime;
70      exp++;
71    }
72    if (exp != 0)
73    {
74      buf= result;
75      result= new int [length + exp];
76      for (int k= 0; k < length; k++)
77        result [k]= buf[k];
78      for (int k= 0; k < exp; k++)
79        result [k + length]= next_prime;
80      length += exp;
81    }
82    exp= 0;
83    j++;
84  }
85  if (j >= 31397)
86    fail= true;
87  ASSERT (j < 31397, "integer factorizer ran out of primes"); //sic
88  return result;
89}
90
91/// make prime factorization distinct
92static inline
93int* makeDistinct (int* factors, const int factors_length, int& length)
94{
95  length= 1;
96  int* result= new int [length];
97  int* buf;
98  result[0]= factors [0];
99  for (int i= 1; i < factors_length; i++)
100  {
101    if (factors[i - 1] != factors[i])
102    {
103      buf= result;
104      result= new int [length + 1];
105      for (int j= 0; j < length; j++)
106        result[j]= buf [j];
107      result[length]= factors[i];
108      length++;
109    }
110  }
111  return result;
112}
113
114/// compute the n-th cyclotomic polynomial,
115/// function may fail if integer_factorizer fails to factorize n
116CanonicalForm cyclotomicPoly (int n, bool& fail)
117{
118  fail= false;
119  Variable x= Variable (1);
120  CanonicalForm result= x - 1;
121  if (n == 1)
122    return result;
123  int* prime_factors;
124  int prime_factors_length;
125  int distinct_factors_length;
126  prime_factors= integerFactorizer (n, prime_factors_length, fail);
127  int* distinct_factors= makeDistinct (prime_factors, prime_factors_length,
128                                        distinct_factors_length);
129  if (fail)
130    return 1;
131  CanonicalForm buf;
132  int prod= 1;
133  for (int i= 0; i < distinct_factors_length; i++)
134  {
135    result= result (power (x, distinct_factors[i]), x)/result;
136    prod *= distinct_factors[i];
137  }
138  return result (power (x, n/prod), x);
139}
140
141#ifdef HAVE_NTL
142/// checks if alpha is a primitive element, alpha is assumed to be an algebraic
143/// variable over some finite prime field
144bool isPrimitive (const Variable& alpha, bool& fail)
145{
146  int p= getCharacteristic();
147  CanonicalForm mipo= getMipo (alpha);
148  int order= ipower(p, degree(mipo)) - 1;
149  CanonicalForm cyclo= cyclotomicPoly (order, fail);
150  if (fail)
151    return false;
152  if (mod(cyclo, mipo (Variable(1), alpha)) == 0)
153    return true;
154  else
155    return false;
156}
157
158#endif
Note: See TracBrowser for help on using the repository browser.