source: git/factory/gfops.cc @ ff3a4f

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