source: git/factory/cf_cyclo.cc @ 72bfc8

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