source: git/coeffs/ffields.cc @ d0a51ee

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