source: git/factory/gfops.cc @ 4dfcb1

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