source: git/factory/gfops.cc @ 386b3d

spielwiese
Last change on this file since 386b3d was 8710ff0, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: 64bit integers in factory by Adi Popescu
  • Property mode set to 100644
File size: 6.7 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2
3#include "config.h"
4
5#ifdef HAVE_CSTDIO
6#include <cstdio>
7#include <cstdlib>
8#else
9#include <stdio.h>
10#include <stdlib.h>
11#endif
12#include <string.h>
13
14#include "cf_assert.h"
15
16#include "cf_defs.h"
17#include "gf_tabutil.h"
18#include "cf_util.h"
19#include "canonicalform.h"
20#include "variable.h"
21#include "gfops.h"
22
23#ifdef SINGULAR
24#include "singext.h"
25#endif
26
27
28const int gf_maxtable = 63001;
29const int gf_maxbuffer = 200;
30
31const int gf_primes_len = 42;
32#ifndef NOASSERT
33static unsigned short gf_primes [] =
34{
35      2,   3,   5,   7,  11,  13,  17,  19,
36     23,  29,  31,  37,  41,  43,  47,  53,
37     59,  61,  67,  71,  73,  79,  83,  89,
38     97, 101, 103, 107, 109, 113, 127, 131,
39    137, 139, 149, 151, 157, 163, 167, 173,
40    179, 181, 191, 193, 197, 199, 223, 211,
41    227, 229, 233, 239, 241, 251
42};
43#endif
44
45int gf_q = 0;
46int gf_p = 0;
47int gf_n = 0;
48int gf_q1 = 0;
49int gf_m1 = 0;
50char gf_name = 'Z';
51
52unsigned short * gf_table = 0;
53
54CanonicalForm gf_mipo(0);
55
56static CanonicalForm intVec2CF ( int degree, int * coeffs, int level )
57{
58    int i;
59    CanonicalForm result;
60    for ( i = 0; i <= degree; i++ )
61    {
62        result += CanonicalForm( coeffs[ i ] ) * power( Variable( level ), degree - i );
63    }
64    return result;
65}
66
67static char *gftable_dir;
68extern "C" {
69  void set_gftable_dir(char *d){
70    gftable_dir = d;
71  }
72}
73
74static void gf_get_table ( int p, int n )
75{
76    char buffer[gf_maxbuffer];
77    int q = ipower( p, n );
78
79    // do not read the table a second time
80    if ( gf_q == q )
81    {
82        return;
83    }
84
85    if ( gf_table == 0 )
86        gf_table = new unsigned short[gf_maxtable];
87
88/*#ifdef SINGULAR
89    // just copy the table if Singular already read it
90    //printf("init_gf(gf_get_table) q=%d, nfCharQ=%d\n",q,nfCharQ);
91    if ( q == nfCharQ )
92    {
93        gf_p = p; gf_n = n;
94        gf_q = q; gf_q1 = q - 1;
95        gf_m1 = nfM1;
96        gf_mipo = intVec2CF( nfMinPoly[0], nfMinPoly + 1, 1 );
97        (void)memcpy( gf_table, nfPlus1Table, gf_q * sizeof( unsigned short ) );
98        gf_table[gf_q] = 0;
99        return;
100    }
101#endif*/
102
103    // try to open file
104    char *gffilename;
105    FILE * inputfile;
106    if (gftable_dir)
107    {
108      sprintf( buffer, "/gftable.%d.%d", p, n );
109      gffilename = (char *)malloc(strlen(gftable_dir) + strlen(buffer) + 1);
110      STICKYASSERT(gffilename,"out of memory");
111      strcpy(gffilename,gftable_dir);
112      strcat(gffilename,buffer);
113      inputfile = fopen( gffilename, "r" );
114    }
115    else
116    {
117#ifndef SINGULAR
118      sprintf( buffer, GFTABLEDIR "/gftable.%d.%d", p, n );
119      gffilename = buffer;
120      inputfile = fopen( buffer, "r" );
121#else
122      sprintf( buffer, "gftables/%d", q );
123      gffilename = buffer;
124      inputfile = feFopen( buffer, "r" );
125#endif
126    }
127    if (!inputfile)
128    {
129      fprintf(stderr,"can not open GF(q) addition table: %s\n",gffilename);
130      STICKYASSERT(inputfile, "can not open GF(q) table");
131    }
132
133    // read ID
134    char * bufptr;
135    char * success;
136    success = fgets( buffer, gf_maxbuffer, inputfile );
137    STICKYASSERT( success, "illegal table (reading ID)" );
138    STICKYASSERT( strcmp( buffer, "@@ factory GF(q) table @@\n" ) == 0, "illegal table" );
139    // read p and n from file
140    int pFile, nFile;
141    success = fgets( buffer, gf_maxbuffer, inputfile );
142    STICKYASSERT( success, "illegal table (reading p and n)" );
143    sscanf( buffer, "%d %d", &pFile, &nFile );
144    STICKYASSERT( p == pFile && n == nFile, "illegal table" );
145    // skip (sic!) factory-representation of mipo
146    // and terminating "; "
147    bufptr = (char *)strchr( buffer, ';' ) + 2;
148    // read simple representation of mipo
149    int i, degree;
150    sscanf( bufptr, "%d", &degree );
151    bufptr = (char *)strchr( bufptr, ' ' ) + 1;
152    int * mipo = new int[degree + 1];
153    for ( i = 0; i <= degree; i++ )
154    {
155        sscanf( bufptr, "%d", mipo + i );
156        bufptr = (char *)strchr( bufptr, ' ' ) + 1;
157    }
158
159    gf_p = p; gf_n = n;
160    gf_q = q; gf_q1 = q-1;
161    gf_mipo = intVec2CF( degree, mipo, 1 );
162    delete [] mipo;
163
164    // now for the table
165    int k, digs = gf_tab_numdigits62( gf_q );
166    i = 1;
167    while ( i < gf_q )
168    {
169        success = fgets( buffer, gf_maxbuffer, inputfile );
170        STICKYASSERT( strlen( buffer ) - 1 == (size_t)digs * 30, "illegal table" );
171        bufptr = buffer;
172        k = 0;
173        while ( i < gf_q && k < 30 )
174        {
175            gf_table[i] = convertback62( bufptr, digs );
176            bufptr += digs;
177            if ( gf_table[i] == gf_q )
178            {
179                if ( i == gf_q1 )
180                    gf_m1 = 0;
181                else
182                    gf_m1 = i;
183            }
184            i++; k++;
185        }
186    }
187    gf_table[0] = gf_table[gf_q1];
188    gf_table[gf_q] = 0;
189
190    (void)fclose( inputfile );
191}
192
193#ifndef NOASSERT
194static bool gf_valid_combination ( int p, int n )
195{
196    int i = 0;
197    while ( i < gf_primes_len && gf_primes[i] != p ) i++;
198    if ( i == gf_primes_len )
199        return false;
200    else
201    {
202        i = n;
203        int a = 1;
204        while ( a < gf_maxtable && i > 0 )
205        {
206            a *= p;
207            i--;
208        }
209        if ( i > 0 || a > gf_maxtable )
210            return false;
211        else
212            return true;
213    }
214}
215#endif
216
217void gf_setcharacteristic ( int p, int n, char name )
218{
219    ASSERT( gf_valid_combination( p, n ), "illegal immediate GF(q)" );
220    gf_name = name;
221    gf_get_table( p, n );
222}
223
224long gf_gf2ff ( long a )
225{
226    if ( gf_iszero( a ) )
227        return 0;
228    else
229    {
230        // starting from z^0=1, step through the table
231        // counting the steps until we hit z^a or z^0
232        // again.  since we are working in char(p), the
233        // latter is guaranteed to be fulfilled.
234        long i = 0, ff = 1;
235        do
236        {
237            if ( i == a )
238                return ff;
239            ff++;
240            i = gf_table[i];
241        } while ( i != 0 );
242        return -1;
243    }
244}
245
246int gf_gf2ff ( int a )
247{
248    if ( gf_iszero( a ) )
249        return 0;
250    else
251    {
252        // starting from z^0=1, step through the table
253        // counting the steps until we hit z^a or z^0
254        // again.  since we are working in char(p), the
255        // latter is guaranteed to be fulfilled.
256        int i = 0, ff = 1;
257        do
258        {
259            if ( i == a )
260                return ff;
261            ff++;
262            i = gf_table[i];
263        } while ( i != 0 );
264        return -1;
265    }
266}
267
268bool gf_isff ( long a )
269{
270    if ( gf_iszero( a ) )
271        return true;
272    else
273    {
274        // z^a in GF(p) iff (z^a)^p-1=1
275        return gf_isone( gf_power( a, gf_p - 1 ) );
276    }
277}
278
279bool gf_isff ( int a )
280{
281    if ( gf_iszero( a ) )
282        return true;
283    else
284    {
285        // z^a in GF(p) iff (z^a)^p-1=1
286        return gf_isone( gf_power( a, gf_p - 1 ) );
287    }
288}
Note: See TracBrowser for help on using the repository browser.