source: git/kernel/ffields.cc @ 341696

spielwiese
Last change on this file since 341696 was 341696, checked in by Hans Schönemann <hannes@…>, 14 years ago
Adding Id property to all files git-svn-id: file:///usr/local/Singular/svn/trunk@12231 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 13.3 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6* ABSTRACT: finite fields with a none-prime number of elements (via tables)
7*/
8
9#include <string.h>
10#include "mod2.h"
11#include <mylimits.h>
12#include "febase.h"
13#include "omalloc.h"
14#include "numbers.h"
15#include "ring.h"
16#include "ffields.h"
17
18int nfCharQ=0;  /* the number of elemts: q*/
19int nfM1;       /*representation of -1*/
20int nfCharP=0;  /* the characteristic: p*/
21static int nfCharQ1=0; /* q-1 */
22CARDINAL *nfPlus1Table=NULL; /* the table i=log(z^i) -> log(z^i+1) */
23char * nfParameter;          /*  the name of the primitive element */
24/* the q's from the table 'fftable' */
25unsigned short fftable[]={
26    4,  8, 16, 32, 64, 128, 256, 512,1024,2048,4096,8192,16384, 32768,
27/*2^2 2^3 2^4 2^5 2^6  2^7  2^8  2^9 2^10 2^11 2^12 2^13  2^14  2^15*/
28    9, 27, 81,243,729,2187, 6561,19683,59049,
29/*3^2 3^3 3^4 3^5 3^6  3^7  3^8   3^9  3^10*/
30   25,125,625,3125,15625,
31/*5^2 5^3 5^4 5^5  5^6*/
32   49,343,2401,16807,
33/*7^2 7^3  7^4 7^5*/
34   121,1331, 14641,
35/*11^2 11^3  11^4*/
36  169, 2197, 28561,
37/*13^2 13^3  13^4*/
38  289, 4913,
39/*17^2 17^3*/
40  361, 6859,
41/*19^2 19^3*/
42  529, 12167,
43/*23^2 23^3*/
44  841, 24389,
45/*29^2 29^3*/
46  961, 29791,
47/*31^2 31^3*/
48  1369, 50653,
49/*37^2  37^3*/
50  1681, /*41^2*/
51  1849, /*43^2*/
52  2209, /*47^2*/
53  2809, /*53^2*/
54  3481, /*59^2*/
55  3721, /*61^2*/
56  4489, /*67^2*/
57  5041, /*71^2*/
58  5329, /*73^2*/
59  6241, /*79^2*/
60  6889, /*83^2*/
61  7921, /*89^2*/
62  9409, /*97^2*/
63  10201, /*101^2*/
64  10609, /*103^2*/
65  11449, /*107^2*/
66  11881, /*109^2*/
67  12769, /*113^2*/
68  16129, /*127^2*/
69  17161, /*131^2*/
70  18769, /*137^2*/
71  19321, /*139^2*/
72  22201, /*149^2*/
73  22801, /*151^2*/
74  24649, /*157^2*/
75  26569, /*163^2*/
76  27889, /*167^2*/
77  29929, /*173^2*/
78  32041, /*179^2*/
79  32761, /*181^2*/
80  36481, /*191^2*/
81  37249, /*193^2*/
82  38809, /*197^2*/
83  39601, /*199^2*/
84  49729, /*223^2*/
85  44521, /*211^2*/
86  51529, /*227^2*/
87  52441, /*229^2*/
88  54289, /*233^2*/
89  57121, /*239^2*/
90  58081, /*241^2*/
91  63001, /*251^2*/
92  0 };
93
94/*1
95* numbers in GF(p^n):
96* let nfCharQ=q=nfCharP^n=p^n
97* GF(q)\{0} will be generated by powers of an element Z
98* Z^i will be represented by the int i, 1 by the int 0, 0 by the int q=nfChar
99*/
100
101#ifdef LDEBUG
102/*2
103* debugging: is a a valid representation of a number ?
104*/
105BOOLEAN nfDBTest (number a, const char *f, const int l)
106{
107  if (((long)a<0L) || ((long)a>(long)nfCharQ))
108  {
109    Print("wrong %d in %s:%d\n",(int)((long)a),f,l);
110    return FALSE;
111  }
112  int i=0;
113  do
114  {
115    if (nfPlus1Table[i]>nfCharQ)
116    {
117      Print("wrong table %d=%d in %s:%d\n",i,nfPlus1Table[i],f,l);
118      return FALSE;
119    }
120    i++;
121  } while (i<nfCharQ);
122  return TRUE;
123}
124#define nfTest(N) nfDBTest(N,__FILE__,__LINE__)
125#endif
126
127/*2
128* k >= 0 ?
129*/
130BOOLEAN nfGreaterZero (number k)
131{
132#ifdef LDEBUG
133  nfTest(k);
134#endif
135  return !nfIsZero(k) && !nfIsMOne(k);
136}
137
138/*2
139* a*b
140*/
141number nfMult (number a,number b)
142{
143#ifdef LDEBUG
144  nfTest(a);
145  nfTest(b);
146#endif
147  if (((long)a == (long)nfCharQ) || ((long)b == (long)nfCharQ))
148    return (number)(long)nfCharQ;
149  /*else*/
150  int i=(int)((long)a+(long)b);
151  if (i>=nfCharQ1) i-=nfCharQ1;
152#ifdef LDEBUG
153  nfTest((number)(long)i);
154#endif
155  return (number)(long)i;
156}
157
158/*2
159* int -> number
160*/
161number nfInit (int i, const ring r)
162{
163  // Hmm .. this is just to prevent initialization
164  // from nfInitChar to go into an infinite loop
165  if (i==0) return (number)(long)nfCharQ;
166  while (i <  0)    i += nfCharP;
167  while (i >= nfCharP) i -= nfCharP;
168  if (i==0) return (number)(long)nfCharQ;
169  CARDINAL c=0;
170  while (i>1)
171  {
172    c=nfPlus1Table[c];
173    i--;
174  }
175#ifdef LDEBUG
176  nfTest((number)(long)c);
177#endif
178  return (number)(long)c;
179}
180
181/*
182* the generating element `z`
183*/
184number nfPar (int i)
185{
186  return (number)1;
187}
188
189/*2
190* the degree of the "alg. number"
191*/
192int nfParDeg(number n)
193{
194#ifdef LDEBUG
195  nfTest(n);
196#endif
197  if((long)nfCharQ == (long)n) return -1;
198  return (int)((long)n);
199}
200
201/*2
202* number -> int
203*/
204int nfInt (number &n, const ring r)
205{
206  return 0;
207}
208
209/*2
210* a + b
211*/
212number nfAdd (number a, number b)
213{
214/*4 z^a+z^b=z^b*(z^(a-b)+1), if a>=b; *
215*          =z^a*(z^(b-a)+1)  if a<b  */
216#ifdef LDEBUG
217  nfTest(a);
218  nfTest(b);
219#endif
220  if ((long)nfCharQ == (long)a) return b;
221  if ((long)nfCharQ == (long)b) return a;
222  long zb,zab,r;
223  if ((long)a >= (long)b)
224  {
225    zb = (long)b;
226    zab = (long)a-(long)b;
227  }
228  else
229  {
230    zb = (long)a;
231    zab = (long)b-(long)a;
232  }
233#ifdef LDEBUG
234  nfTest((number)zab);
235#endif
236  if (nfPlus1Table[zab]==nfCharQ) r=(long)nfCharQ; /*if z^(a-b)+1 =0*/
237  else
238  {
239    r= zb+(long)nfPlus1Table[zab];
240    if(r>=(long)nfCharQ1) r-=(long)nfCharQ1;
241  }
242#ifdef LDEBUG
243  nfTest((number)r);
244#endif
245  return (number)r;
246}
247
248/*2
249* a - b
250*/
251number nfSub (number a, number b)
252{
253  number mb = nfNeg(b);
254  return nfAdd(a,mb);
255}
256
257/*2
258* a == 0 ?
259*/
260BOOLEAN nfIsZero (number  a)
261{
262#ifdef LDEBUG
263  nfTest(a);
264#endif
265  return (long)nfCharQ == (long)a;
266}
267
268/*2
269* a == 1 ?
270*/
271BOOLEAN nfIsOne (number a)
272{
273#ifdef LDEBUG
274  nfTest(a);
275#endif
276  return 0L == (long)a;
277}
278
279/*2
280* a == -1 ?
281*/
282BOOLEAN nfIsMOne (number a)
283{
284#ifdef LDEBUG
285  nfTest(a);
286#endif
287  if (0L == (long)a) return FALSE; /* special handling of char 2*/
288  return (long)nfM1 == (long)a;
289}
290
291/*2
292* a / b
293*/
294number nfDiv (number a,number b)
295{
296#ifdef LDEBUG
297  nfTest(b);
298#endif
299  if ((long)b==(long)nfCharQ)
300  {
301    WerrorS(nDivBy0);
302    return (number)((long)nfCharQ);
303  }
304#ifdef LDEBUG
305  nfTest(a);
306#endif
307  if ((long)a==(long)nfCharQ)
308    return (number)((long)nfCharQ);
309  /*else*/
310  long s = (long)a - (long)b;
311  if (s < 0L)
312    s += (long)nfCharQ1;
313#ifdef LDEBUG
314  nfTest((number)s);
315#endif
316  return (number)s;
317}
318
319/*2
320* 1 / c
321*/
322number  nfInvers (number c)
323{
324#ifdef LDEBUG
325  nfTest(c);
326#endif
327  if ((long)c==(long)nfCharQ)
328  {
329    WerrorS(nDivBy0);
330    return (number)((long)nfCharQ);
331  }
332#ifdef LDEBUG
333  nfTest(((number)((long)nfCharQ1-(long)c)));
334#endif
335  return (number)((long)nfCharQ1-(long)c);
336}
337
338/*2
339* -c
340*/
341number nfNeg (number c)
342{
343/*4 -z^c=z^c*(-1)=z^c*nfM1*/
344#ifdef LDEBUG
345  nfTest(c);
346#endif
347  if ((long)nfCharQ == (long)c) return c;
348  long i=(long)c+(long)nfM1;
349  if (i>=(long)nfCharQ1) i-=(long)nfCharQ1;
350#ifdef LDEBUG
351  nfTest((number)i);
352#endif
353  return (number)i;
354}
355
356/*2
357* a > b ?
358*/
359BOOLEAN nfGreater (number a,number b)
360{
361#ifdef LDEBUG
362  nfTest(a);
363  nfTest(b);
364#endif
365  return (long)a != (long)b;
366}
367
368/*2
369* a == b ?
370*/
371BOOLEAN nfEqual (number a,number b)
372{
373#ifdef LDEBUG
374  nfTest(a);
375  nfTest(b);
376#endif
377  return (long)a == (long)b;
378}
379
380/*2
381* write via StringAppend
382*/
383void nfWrite (number &a)
384{
385#ifdef LDEBUG
386  nfTest(a);
387#endif
388  if ((long)a==(long)nfCharQ)  StringAppendS("0");
389  else if ((long)a==0L)   StringAppendS("1");
390  else if (nfIsMOne(a))   StringAppendS("-1");
391  else
392  {
393    StringAppendS(nfParameter);
394    if ((long)a!=1L)
395    {
396      if(currRing->ShortOut==0)  StringAppendS("^");
397      StringAppend("%d",(int)((long)a));
398    }
399  }
400}
401
402/*2
403*
404*/
405char * nfName(number a)
406{
407#ifdef LDEBUG
408  nfTest(a);
409#endif
410  char *s;
411  if (((long)a==(long)nfCharQ) || ((long)a==0L)) return NULL;
412  else if ((long)a==1L)
413  {
414    return omStrDup(nfParameter);
415  }
416  else
417  {
418    s=(char *)omAlloc(4+strlen(nfParameter));
419    sprintf(s,"%s%d",nfParameter,(int)((long)a));
420  }
421  return s;
422}
423/*2
424* c ^ i with i>=0
425*/
426void nfPower (number a, int i, number * result)
427{
428#ifdef LDEBUG
429  nfTest(a);
430#endif
431  if (i==0)
432  {
433    //*result=nfInit(1);
434    *result = (number)0L;
435  }
436  else if (i==1)
437  {
438    *result = a;
439  }
440  else
441  {
442    nfPower(a,i-1,result);
443    *result = nfMult(a,*result);
444  }
445#ifdef LDEBUG
446  nfTest(*result);
447#endif
448}
449
450/*4
451* read an integer (with reduction mod p)
452*/
453static const char* nfEati(const char *s, int *i)
454{
455  if (*s >= '0' && *s <= '9')
456  {
457    *i = 0;
458    do
459    {
460      *i *= 10;
461      *i += *s++ - '0';
462      if (*i > (MAX_INT_VAL / 10)) *i = *i % nfCharP;
463    }
464    while (*s >= '0' && *s <= '9');
465    if (*i >= nfCharP) *i = *i % nfCharP;
466  }
467  else *i = 1;
468  return s;
469}
470
471/*2
472* read a number
473*/
474const char * nfRead (const char *s, number *a)
475{
476  int i;
477  number z;
478  number n;
479
480  s = nfEati(s, &i);
481  z=nfInit(i, currRing);
482  *a=z;
483  if (*s == '/')
484  {
485    s++;
486    s = nfEati(s, &i);
487    n=nfInit(i, currRing);
488    *a = nfDiv(z,n);
489  }
490  if (strncmp(s,nfParameter,strlen(nfParameter))==0)
491  {
492    s+=strlen(nfParameter);
493    if ((*s >= '0') && (*s <= '9'))
494    {
495      s=eati(s,&i);
496      while (i>=nfCharQ1) i-=nfCharQ1;
497    }
498    else
499      i=1;
500    z=(number)(long)i;
501    *a=nfMult(*a,z);
502  }
503#ifdef LDEBUG
504  nfTest(*a);
505#endif
506  return s;
507}
508
509#ifdef HAVE_FACTORY
510int gf_tab_numdigits62 ( int q );
511int convertback62 ( char * p, int n );
512#else
513static int gf_tab_numdigits62 ( int q )
514{
515    if ( q < 62 )          return 1;
516    else  if ( q < 62*62 ) return 2;
517    /*else*/               return 3;
518}
519
520static int convback62 ( char c )
521{
522    if ( c >= '0' && c <= '9' )
523        return int(c) - int('0');
524    else  if ( c >= 'A' && c <= 'Z' )
525        return int(c) - int('A') + 10;
526    else
527        return int(c) - int('a') + 36;
528}
529
530static int convertback62 ( char * p, int n )
531{
532    int r = 0;
533    for ( int j = 0; j < n; j++ )
534        r = r * 62 + convback62( p[j] );
535    return r;
536}
537#endif
538
539int nfMinPoly[16];
540
541void nfShowMipo()
542{
543  int i=nfMinPoly[0];
544  int j=0;
545  loop
546  {
547    j++;
548    if (nfMinPoly[j]!=0)
549      StringAppend("%d*%s^%d",nfMinPoly[j],nfParameter,i);
550    i--;
551    if(i<0) break;
552    if (nfMinPoly[j]!=0)
553      StringAppendS("+");
554  }
555}
556
557static void nfReadMipo(char *s)
558{
559  const char *l=strchr(s,';')+1;
560  char *n;
561  int i=strtol(l,&n,10);
562  l=n;
563  int j=1;
564  nfMinPoly[0]=i;
565  while(i>=0)
566  {
567    nfMinPoly[j]=strtol(l,&n,10);
568    if (l==n) break;
569    l=n;
570    j++;
571    i--;
572  }
573  if (i>=0)
574  {
575    WerrorS("error in reading minpoly from gftables");
576  }
577}
578
579/*2
580* init global variables from files 'gftables/%d'
581*/
582void nfSetChar(int c, char **param)
583{
584  //Print("GF(%d)\n",c);
585  nfParameter=param[0];
586  if ((c==nfCharQ)||(c==-nfCharQ))
587    /*this field is already set*/  return;
588  int i=0;
589  while ((fftable[i]!=c) && (fftable[i]!=0)) i++;
590  if (fftable[i]==0)
591    return;
592  if (nfCharQ > 1)
593  {
594    omFreeSize( (ADDRESS)nfPlus1Table,nfCharQ*sizeof(CARDINAL) );
595    nfPlus1Table=NULL;
596  }
597  if ((c>1) || (c<0))
598  {
599    if (c>1) nfCharQ = c;
600    else     nfCharQ = -c;
601    char buf[100];
602    sprintf(buf,"gftables/%d",nfCharQ);
603    FILE * fp = feFopen(buf,"r",NULL,TRUE);
604    if (fp==NULL)
605    {
606      return;
607    }
608    if(!fgets( buf, sizeof(buf), fp)) return;
609    if(strcmp(buf,"@@ factory GF(q) table @@\n")!=0)
610    {
611      goto err;
612    }
613    if(!fgets( buf, sizeof(buf), fp))
614    {
615      goto err;
616    }
617    int q;
618    sscanf(buf,"%d %d",&nfCharP,&q);
619    nfReadMipo(buf);
620    nfCharQ1=nfCharQ-1;
621    //Print("nfCharQ=%d,nfCharQ1=%d,mipo=>>%s<<\n",nfCharQ,nfCharQ1,buf);
622    nfPlus1Table= (CARDINAL *)omAlloc( (nfCharQ)*sizeof(CARDINAL) );
623    int digs = gf_tab_numdigits62( nfCharQ );
624    char * bufptr;
625    int i = 1;
626    int k;
627    while ( i < nfCharQ )
628    {
629      fgets( buf, sizeof(buf), fp);
630      //( strlen( buffer ) == (size_t)digs * 30, "illegal table" );
631      bufptr = buf;
632      k = 0;
633      while ( (i < nfCharQ) && (k < 30) )
634      {
635        nfPlus1Table[i] = convertback62( bufptr, digs );
636        if(nfPlus1Table[i]>nfCharQ)
637        {
638          Print("wrong entry %d: %d(%c%c%c)\n",i,nfPlus1Table[i],bufptr[0],bufptr[1],bufptr[2]);
639        }
640        bufptr += digs;
641        if (nfPlus1Table[i]==nfCharQ)
642        {
643          if(i==nfCharQ1)
644          {
645            nfM1=0;
646          }
647          else
648          {
649            nfM1=i;
650          }
651        }
652        i++; k++;
653      }
654    }
655    nfPlus1Table[0]=nfPlus1Table[nfCharQ1];
656  }
657  else
658    nfCharQ=0;
659#ifdef LDEBUG
660  nfTest((number)0);
661#endif
662  return;
663err:
664  Werror("illegal GF-table %d",nfCharQ);
665}
666
667/*2
668* map Z/p -> GF(p,n)
669*/
670number nfMapP(number c)
671{
672  return nfInit((int)((long)c), currRing);
673}
674
675/*2
676* map GF(p,n1) -> GF(p,n2), n1 < n2, n1 | n2
677*/
678int nfMapGG_factor;
679number nfMapGG(number c)
680{
681  int i=(long)c;
682  i*= nfMapGG_factor;
683  while (i >nfCharQ1) i-=nfCharQ1;
684  return (number)((long)i);
685}
686/*2
687* map GF(p,n1) -> GF(p,n2), n1 > n2, n2 | n1
688*/
689number nfMapGGrev(number c)
690{
691  int ex=(int)((long)c);
692  if ((ex % nfMapGG_factor)==0)
693    return (number)(((long)ex) / ((long)nfMapGG_factor));
694  else
695    return (number)(long)nfCharQ; /* 0 */
696}
697
698/*2
699* set map function nMap ... -> GF(p,n)
700*/
701nMapFunc nfSetMap(ring src, ring dst)
702{
703  if (rField_is_GF(src,nfCharQ))
704  {
705    return ndCopy;   /* GF(p,n) -> GF(p,n) */
706  }
707  if (rField_is_GF(src))
708  {
709    int q=src->ch;
710    if ((nfCharQ % q)==0) /* GF(p,n1) -> GF(p,n2), n2 > n1 */
711    {
712      // check if n2 is a multiple of n1
713      int n1=1;
714      int qq=nfCharP;
715      while(qq!=q) { qq *= nfCharP; n1++; }
716      int n2=1;
717      qq=nfCharP;
718      while(qq!=nfCharQ) { qq *= nfCharP; n2++; }
719      Print("map %d^%d -> %d^%d\n",nfCharP,n1,nfCharP,n2);
720      if ((n2 % n1)==0)
721      {
722        int save_ch=currRing->ch;
723        char **save_par=currRing->parameter;
724        nfSetChar(src->ch,src->parameter);
725        int nn=nfPlus1Table[0];
726        nfSetChar(save_ch,save_par);
727        nfMapGG_factor= nfPlus1Table[0] / nn;
728        Print("nfMapGG_factor=%d (%d / %d)\n",nfMapGG_factor, nfPlus1Table[0], nn);
729        return nfMapGG;
730      }
731      else if ((n1 % n2)==0)
732      {
733        nfMapGG_factor= (n1/n2);
734        return nfMapGGrev;
735      }
736      else
737        return NULL;
738    }
739  }
740  if (rField_is_Zp(src,nfCharP))
741  {
742    return nfMapP;    /* Z/p -> GF(p,n) */
743  }
744  return NULL;     /* default */
745}
Note: See TracBrowser for help on using the repository browser.