source: git/coeffs/numbers.cc @ 3aae0e

spielwiese
Last change on this file since 3aae0e was 3aae0e, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
Finished shortfl. Notice: CopyMap change
  • Property mode set to 100644
File size: 15.4 KB
Line 
1/*****************************************
2*  Computer Algebra System SINGULAR      *
3*****************************************/
4/* $Id$ */
5
6/*
7* ABSTRACT: interface to coefficient aritmetics
8*/
9
10#include <string.h>
11#include <stdlib.h>
12#include "coeffs.h"
13#include "numbers.h"
14#include "longrat.h"
15#include "longalg.h"
16#include "modulop.h"
17#include "gnumpfl.h"
18#include "gnumpc.h"
19#include "ffields.h"
20#include "shortfl.h"
21#include "longtrans.h"
22#ifdef HAVE_RINGS
23#include <kernel/rmodulo2m.h>
24#include <kernel/rmodulon.h>
25#include <kernel/rintegers.h>
26
27extern omBin gmp_nrz_bin;
28#endif
29
30#ifndef assume
31#  define assume(a) if(!(a)){ Werror( "Assumption: is wrong: %s\n", #a ); };
32#endif
33
34
35//static int characteristic = 0;
36extern int IsPrime(int p);
37
38/*0 implementation*/
39number nNULL; /* the 0 as constant */
40
41n_Procs_s *cf_root=NULL;
42
43void   nNew(number* d) { *d=NULL; }
44void   ndDelete(number* d, const coeffs r) { *d=NULL; }
45void   ndInpMult(number &a, number b, const coeffs r)
46{
47  number n=n_Mult(a,b,r);
48  n_Delete(&a,r);
49  a=n;
50}
51number ndInpAdd(number &a, number b, const coeffs r)
52{
53  number n=n_Add(a,b,r);
54  n_Delete(&a,r);
55  a=n;
56  return a;
57}
58
59#ifdef LDEBUG
60void   nDBDummy1(number* d,char *f, int l) { *d=NULL; }
61BOOLEAN ndDBTest(number a, const char *f, const int l)
62{
63  return TRUE;
64}
65#endif
66
67void   nDummy2(number& d, const coeffs r) { }
68
69char * ndName(number n, const coeffs r) { return NULL; }
70
71number ndPar(int i, const coeffs r) { return n_Init(0,r); }
72
73number ndReturn0(number n, const coeffs r) { return n_Init(0,r); }
74
75int    ndParDeg(number n, const coeffs r) { return 0; }
76
77number ndGcd(number a, number b, const coeffs r) { return n_Init(1,r); }
78
79number ndIntMod(number a, number b, const coeffs r) { return n_Init(0,r); }
80
81number ndGetDenom(number &n, const coeffs r) { return n_Init(1,r); }
82number ndGetNumerator(number &a,const coeffs r) { return n_Copy(a,r); }
83
84int ndSize(number a, const coeffs r) { return (int)n_IsZero(a,r)==FALSE; }
85
86number ndCopy(number a, const coeffs) { return a; }
87number ndCopyMap(number a, const coeffs r, const coeffs aRing)
88{
89  assume( getCoeffType(r) == getCoeffType(aRing) );
90  assume( nField_has_simple_Alloc(r) && nField_has_simple_Alloc(aRing) );
91 
92  return a;
93}
94
95number nd_Copy(number a, const coeffs r) { return n_Copy(a, r); }
96
97#ifdef HAVE_RINGS
98BOOLEAN ndDivBy(number a, number b) { return TRUE; } // assume a,b !=0
99int ndDivComp(number a, number b) { return 2; }
100BOOLEAN ndIsUnit(number a) { return !nIsZero(a); }
101number  ndExtGcd (number a, number b, number *s, number *t) { return nInit(1); }
102#endif
103
104/*2
105* init operations for characteristic c (complete==TRUE)
106* init nDelete    for characteristic c (complete==FALSE)
107*/
108void nSetChar(coeffs r)
109{
110  int c=rInternalChar(r);
111
112  /*--------------------- Q -----------------*/
113  if (nField_is_Q(r))
114  {
115    r->nInit_bigint=nl_Copy;
116  }
117  /*--------------------- Q_a/ Zp_a -----------------*/
118  else if (nField_is_Extension(r))
119  {
120    if (r->minpoly != NULL)
121    {
122      naSetChar(c,r);
123      if (rField_is_Q_a(r)) r->nInit_bigint=naMap00;
124      if (rField_is_Zp_a(r)) r->nInit_bigint=naMap0P;
125    }
126    else
127    {
128      ntSetChar(c,r);
129      if (rField_is_Q_a(r)) r->nInit_bigint=ntMap00;
130      if (rField_is_Zp_a(r)) r->nInit_bigint=ntMap0P;
131    }
132  }
133#ifdef HAVE_RINGS
134  /*----------------------ring Z / 2^m----------------*/
135  else if (nField_is_Ring_2toM(r))
136  {
137    nr2mSetExp(c, r);
138    r->nInit_bigint=nr2mMapQ;
139  }
140  /*----------------------ring Z ----------------*/
141  else if (nField_is_Ring_Z(r))
142  {
143    nrzSetExp(c, r);
144    r->nInit_bigint=nrzMapQ;
145  }
146  /*----------------------ring Z / n----------------*/
147  else if (nField_is_Ring_ModN(r))
148  {
149    nrnSetExp(c, r);
150    r->nInit_bigint=nrnMapQ;
151  }
152  /*----------------------ring Z / n----------------*/
153  else if (nField_is_Ring_PtoM(r))
154  {
155    nrnSetExp(c, r);
156    r->nInit_bigint=nrnMapQ;
157  }
158#endif
159  else if (nField_is_Zp(r))
160  /*----------------------char. p----------------*/
161  {
162    npSetChar(c, r);
163    r->nInit_bigint=npMap0;
164  }
165  /* -------------- GF(p^m) -----------------------*/
166  else if (nField_is_GF(r))
167  {
168    nfSetChar(c,r->parameter);
169    r->nInit_bigint=ndReturn0; // not impl.
170  }
171  /* -------------- R -----------------------*/
172  //if (c==(-1))
173  else if (nField_is_R(r))
174  {
175    r->nInit_bigint=nrMapQ;
176  }
177  /* -------------- long R -----------------------*/
178  /* -------------- long C -----------------------*/
179  else if ((nField_is_long_R(r))
180  || (nField_is_long_C(r)))
181  {
182    setGMPFloatDigits(r->float_len,r->float_len2);
183    if (nField_is_long_R(r)) r->nInit_bigint=ngfMapQ;
184    else                     r->nInit_bigint=ngcMapQ;
185  }
186#ifdef TEST
187  /* -------------- R -----------------------*/
188  //if (c==(-1))
189  else if (!nField_is_R(r) && !nField_is_Q(r))
190  {
191    WerrorS("unknown field");
192  }
193#endif
194}
195
196/*2
197* init operations for coeffs r
198*/
199void nInitChar(coeffs r)
200{
201  int c=nInternalChar(r);
202  n_coeffType t=nFieldType(r);
203
204  if (nField_is_Extension(r))
205  {
206    if (r->algring==NULL)
207    {
208      int ch=-c;
209      if (c==1) ch=0;
210      r->algring=(ring) rDefault(ch,r->P,r->parameter);
211      //r->algring->ShortOut=r->ShortOut;
212      // includes: nInitChar(r->algring);
213    }
214  }
215
216  n_coeffType t=rFieldType(r);
217  if ((r->cf!=NULL) && (r->cf->nChar==c) && (r->cf->type==t))
218  { r->cf->ref++; return; }
219
220  n_Procs_s *n=cf_root;
221  while((n!=NULL)
222    && ((n->nChar!=c) || (n->type!=t)))
223      n=n->next;
224  if (n==NULL)
225  {
226    n=(n_Procs_s*)omAlloc0(sizeof(n_Procs_s));
227    n->next=cf_root;
228    n->ref=1;
229    n->nChar=c;
230    n->type=t;
231    cf_root=n;
232  }
233  else if ((n->nChar==c) && (n->type==t))
234  {
235    n->ref++;
236    r=n;
237    return;
238  }
239  else
240  {
241    WerrorS("nInitChar failed");
242    return;
243  }
244  r=n;
245  n->nChar = c;
246  n->nPar  = ndPar;
247  n->nParDeg=ndParDeg;
248  n->nSize = ndSize;
249  n->cfGetDenom= ndGetDenom;
250  n->cfGetNumerator= ndGetNumerator;
251  n->nName =  ndName;
252  n->nImPart=ndReturn0;
253  n->cfDelete= ndDelete;
254  n->nInpMult=ndInpMult;
255  n->cfCopy=nd_Copy;
256  n->nIntMod=ndIntMod; /* dummy !! */
257  n->nNormalize=nDummy2;
258  n->nGcd  = ndGcd;
259  n->nLcm  = ndGcd; /* tricky, isn't it ?*/
260#ifdef HAVE_RINGS
261  n->nDivComp = ndDivComp;
262  n->nDivBy = ndDivBy;
263  n->nIsUnit = ndIsUnit;
264  n->nExtGcd = ndExtGcd;
265  n->nGetUnit = (nMapFunc)NULL;
266#endif
267  if (nField_is_Extension(r))
268  {
269    //ntInitChar(c,TRUE,r);
270    n->cfDelete       = ntDelete;
271    n->nNormalize     = ntNormalize;
272    n->cfInit         = ntInit;
273    n->nPar           = ntPar;
274    n->nParDeg        = ntParDeg;
275    n->n_Int          = ntInt;
276    n->nAdd           = ntAdd;
277    n->nSub           = ntSub;
278    n->nMult          = ntMult;
279    n->nDiv           = ntDiv;
280    n->nExactDiv      = ntDiv;
281    n->nIntDiv        = ntIntDiv;
282    n->nNeg           = ntNeg;
283    n->nInvers        = ntInvers;
284    //n->nCopy          = ntCopy;
285    n->cfCopy         = nt_Copy;
286    n->nGreater       = ntGreater;
287    n->nEqual         = ntEqual;
288    n->nIsZero        = ntIsZero;
289    n->nIsOne         = ntIsOne;
290    n->nIsMOne        = ntIsMOne;
291    n->nGreaterZero   = ntGreaterZero;
292    n->cfWrite        = ntWrite;
293    n->nRead          = ntRead;
294    n->nPower         = ntPower;
295    n->nGcd           = ntGcd;
296    n->nLcm           = ntLcm;
297    n->cfSetMap       = ntSetMap;
298    n->nName          = ntName;
299    n->nSize          = ntSize;
300    n->cfGetDenom     = napGetDenom;
301    n->cfGetNumerator = napGetNumerator;
302#ifdef LDEBUG
303    n->nDBTest        = naDBTest;
304#endif
305  }
306#ifdef HAVE_RINGS
307  /* -------------- Z/2^m ----------------------- */
308  else if (nField_is_Ring_2toM(r))
309  {
310    nr2mCoeffInit(n, c, r->cf);
311  }
312  /* -------------- Z/n ----------------------- */
313  else if (nField_is_Ring_ModN(r) || nField_is_Ring_PtoM(r))
314  {
315    nrnCoeffInit(n, c, r->cf);
316  }
317  /* -------------- Z ----------------------- */
318  else if (nField_is_Ring_Z(r))
319  {
320     n->cfInit  = nrzInit;
321     n->cfDelete= nrzDelete;
322     n->nCopy  = nrzCopy;
323     n->cfCopy = cfrzCopy;
324     n->nSize  = nrzSize;
325     n->n_Int  = nrzInt;
326     n->nAdd   = nrzAdd;
327     n->nSub   = nrzSub;
328     n->nMult  = nrzMult;
329     n->nDiv   = nrzDiv;
330     n->nIntDiv = nrzIntDiv;
331     n->nIntMod = nrzIntMod;
332     n->nExactDiv= nrzDiv;
333     n->nNeg   = nrzNeg;
334     n->nInvers= nrzInvers;
335     n->nDivBy = nrzDivBy;
336     n->nDivComp = nrzDivComp;
337     n->nGreater = nrzGreater;
338     n->nEqual = nrzEqual;
339     n->nIsZero = nrzIsZero;
340     n->nIsOne = nrzIsOne;
341     n->nIsMOne = nrzIsMOne;
342     n->nGreaterZero = nrzGreaterZero;
343     n->cfWrite = nrzWrite;
344     n->nRead = nrzRead;
345     n->nPower = nrzPower;
346     n->cfSetMap = nrzSetMap;
347     n->nNormalize = nDummy2;
348     n->nLcm          = nrzLcm;
349     n->nGcd          = nrzGcd;
350     n->nIsUnit = nrzIsUnit;
351     n->nGetUnit = nrzGetUnit;
352     n->nExtGcd = nrzExtGcd;
353     n->nName= ndName;
354#ifdef LDEBUG
355     n->nDBTest=ndDBTest; // not yet implemented: nrzDBTest;
356#endif
357  }
358#endif
359  else if (nField_is_Q(r))
360  {
361    n->cfDelete= nlDelete;
362    n->nNormalize=nlNormalize;
363    n->cfInit = nlInit;
364    n->n_Int  = nlInt;
365    n->nAdd   = nlAdd;
366    n->nSub   = nlSub;
367    n->nMult  = nlMult;
368    n->nInpMult=nlInpMult;
369    n->nDiv   = nlDiv;
370    n->nExactDiv= nlExactDiv;
371    n->nIntDiv= nlIntDiv;
372    n->nIntMod= nlIntMod;
373    n->nNeg   = nlNeg;
374    n->nInvers= nlInvers;
375    n->nCopy  = nlCopy;
376    n->nGreater = nlGreater;
377    n->nEqual = nlEqual;
378    n->nIsZero = nlIsZero;
379    n->nIsOne = nlIsOne;
380    n->nIsMOne = nlIsMOne;
381    n->nGreaterZero = nlGreaterZero;
382    n->cfWrite = nlWrite;
383    n->nRead = nlRead;
384    n->nPower = nlPower;
385    n->nGcd  = nlGcd;
386    n->nLcm  = nlLcm;
387    n->cfSetMap = nlSetMap;
388    n->nSize  = nlSize;
389    n->cfGetDenom = nlGetDenom;
390    n->cfGetNumerator = nlGetNumerator;
391#ifdef LDEBUG
392    n->nDBTest=nlDBTest;
393#endif
394  }
395  else if (nField_is_Zp(r))
396  /*----------------------char. p----------------*/
397  {
398    npInitChar(c,r);
399    n->cfInit = npInit;
400    n->n_Int  = npInt;
401    n->nAdd   = npAdd;
402    n->nSub   = npSub;
403    n->nMult  = npMult;
404    n->nDiv   = npDiv;
405    n->nExactDiv= npDiv;
406    n->nNeg   = npNeg;
407    n->nInvers= npInvers;
408    n->nCopy  = ndCopy;
409    n->nGreater = npGreater;
410    n->nEqual = npEqual;
411    n->nIsZero = npIsZero;
412    n->nIsOne = npIsOne;
413    n->nIsMOne = npIsMOne;
414    n->nGreaterZero = npGreaterZero;
415    n->cfWrite = npWrite;
416    n->nRead = npRead;
417    n->nPower = npPower;
418    n->cfSetMap = npSetMap;
419    /* nName= ndName; */
420    /*nSize  = ndSize;*/
421#ifdef LDEBUG
422    n->nDBTest=npDBTest;
423#endif
424#ifdef NV_OPS
425    if (c>NV_MAX_PRIME)
426    {
427      n->nMult  = nvMult;
428      n->nDiv   = nvDiv;
429      n->nExactDiv= nvDiv;
430      n->nInvers= nvInvers;
431      n->nPower= nvPower;
432      n->nInpMult= nvInpMult;
433    }
434#endif
435  }
436  /* -------------- GF(p^m) -----------------------*/
437  else if (nField_is_GF(r))
438  {
439    //nfSetChar(c,r->parameter);
440    n->cfInit = nfInit;
441    n->nPar   = nfPar;
442    n->nParDeg= nfParDeg;
443    n->n_Int  = nfInt;
444    n->nAdd   = nfAdd;
445    n->nSub   = nfSub;
446    n->nMult  = nfMult;
447    n->nDiv   = nfDiv;
448    n->nExactDiv= nfDiv;
449    n->nNeg   = nfNeg;
450    n->nInvers= nfInvers;
451    n->nCopy  = ndCopy;
452    n->nGreater = nfGreater;
453    n->nEqual = nfEqual;
454    n->nIsZero = nfIsZero;
455    n->nIsOne = nfIsOne;
456    n->nIsMOne = nfIsMOne;
457    n->nGreaterZero = nfGreaterZero;
458    n->cfWrite = nfWrite;
459    n->nRead = nfRead;
460    n->nPower = nfPower;
461    n->cfSetMap = nfSetMap;
462    n->nName= nfName;
463    /*nSize  = ndSize;*/
464#ifdef LDEBUG
465    n->nDBTest=nfDBTest;
466#endif
467  }
468  /* -------------- R -----------------------*/
469  //if (c==(-1))
470  else if (nField_is_R(r))
471  {
472    n->cfInit = nrInit;
473    n->n_Int  = nrInt;
474    n->nAdd   = nrAdd;
475    n->nSub   = nrSub;
476    n->nMult  = nrMult;
477    n->nDiv   = nrDiv;
478    n->nExactDiv= nrDiv;
479    n->nNeg   = nrNeg;
480    n->nInvers= nrInvers;
481    n->nCopy  = ndCopy;
482    n->nGreater = nrGreater;
483    n->nEqual = nrEqual;
484    n->nIsZero = nrIsZero;
485    n->nIsOne = nrIsOne;
486    n->nIsMOne = nrIsMOne;
487    n->nGreaterZero = nrGreaterZero;
488    n->cfWrite = nrWrite;
489    n->nRead = nrRead;
490    n->nPower = nrPower;
491    n->cfSetMap=nrSetMap;
492    /* nName= ndName; */
493    n->nSize = nrSize;
494#ifdef LDEBUG
495    n->nDBTest=ndDBTest; // not yet implemented: nrDBTest;
496#endif
497  }
498  /* -------------- long R -----------------------*/
499  else if (nField_is_long_R(r))
500  {
501    n->cfDelete= ngfDelete;
502    n->cfInit = ngfInit;
503    n->n_Int  = ngfInt;
504    n->nAdd   = ngfAdd;
505    n->nSub   = ngfSub;
506    n->nMult  = ngfMult;
507    n->nDiv   = ngfDiv;
508    n->nExactDiv= ngfDiv;
509    n->nNeg   = ngfNeg;
510    n->nInvers= ngfInvers;
511    n->nCopy  = ngfCopy;
512    n->nGreater = ngfGreater;
513    n->nEqual = ngfEqual;
514    n->nIsZero = ngfIsZero;
515    n->nIsOne = ngfIsOne;
516    n->nIsMOne = ngfIsMOne;
517    n->nGreaterZero = ngfGreaterZero;
518    n->cfWrite = ngfWrite;
519    n->nRead = ngfRead;
520    n->nPower = ngfPower;
521    n->cfSetMap=ngfSetMap;
522    n->nName= ndName;
523    n->nSize  = ngfSize;
524#ifdef LDEBUG
525    n->nDBTest=ndDBTest; // not yet implemented: ngfDBTest
526#endif
527  }
528  /* -------------- long C -----------------------*/
529  else if (nField_is_long_C(r))
530  {
531    n->cfDelete= ngcDelete;
532    n->nNormalize=nDummy2;
533    n->cfInit = ngcInit;
534    n->n_Int  = ngcInt;
535    n->nAdd   = ngcAdd;
536    n->nSub   = ngcSub;
537    n->nMult  = ngcMult;
538    n->nDiv   = ngcDiv;
539    n->nExactDiv= ngcDiv;
540    n->nNeg   = ngcNeg;
541    n->nInvers= ngcInvers;
542    n->nCopy  = ngcCopy;
543    n->nGreater = ngcGreater;
544    n->nEqual = ngcEqual;
545    n->nIsZero = ngcIsZero;
546    n->nIsOne = ngcIsOne;
547    n->nIsMOne = ngcIsMOne;
548    n->nGreaterZero = ngcGreaterZero;
549    n->cfWrite = ngcWrite;
550    n->nRead = ngcRead;
551    n->nPower = ngcPower;
552    n->cfSetMap=ngcSetMap;
553    n->nPar=ngcPar;
554    n->nRePart=ngcRePart;
555    n->nImPart=ngcImPart;
556    n->nSize = ngcSize;
557#ifdef LDEBUG
558    n->nDBTest=ndDBTest; // not yet implemented: ngcDBTest
559#endif
560  }
561#ifdef TEST
562  else
563  {
564    WerrorS("unknown field");
565  }
566#endif
567#ifdef HAVE_RINGS
568  if (n->nGetUnit==(nMapFunc)NULL) n->nGetUnit=n->nCopy;
569#endif
570  if (!errorreported)
571  {
572    n->nNULL=n->cfInit(0,r);
573    if (n->nRePart==NULL)
574      n->nRePart=n->cfCopy;
575    if (n->nIntDiv==NULL)
576      n->nIntDiv=n->nDiv;
577  }
578}
579
580void nKillChar(coeffs r)
581{
582  if (r!=NULL)
583  {
584    r->ref--;
585    if (r->ref<=0)
586    {
587      n_Procs_s tmp;
588      n_Procs_s* n=&tmp;
589      tmp.next=cf_root;
590      while((n->next!=NULL) && (n->next!=r)) n=n->next;
591      if (n->next==r)
592      {
593        n->next=n->next->next;
594        if (cf_root==r) cf_root=n->next;
595        r->cfDelete(&(r->nNULL),r);
596        switch(r->type)
597        {
598          case n_Zp:
599               #ifdef HAVE_DIV_MOD
600               if (r->npInvTable!=NULL)
601               omFreeSize( (void *)r->npInvTable,
602                           r->npPrimeM*sizeof(unsigned short) );
603               #else
604               if (r->npExpTable!=NULL)
605               {
606                 omFreeSize( (void *)r->npExpTable,
607                             r->npPrimeM*sizeof(unsigned short) );
608                 omFreeSize( (void *)r->npLogTable,
609                             r->npPrimeM*sizeof(unsigned short) );
610               }
611               #endif
612               break;
613          case n_Zp_a:
614          case n_Q_a:
615               {
616                 number n=r->minpoly;
617                 if (n!=NULL)
618                 {
619                   r->minpoly=NULL;
620                   naDelete(&n,r);
621                 }
622               }
623               break;
624
625            default:
626                 break;
627          }
628          omFreeSize((void *)r, sizeof(n_Procs_s));
629          r=NULL;
630        }
631        else
632        {
633          WarnS("cf_root list destroyed");
634        }
635      }
636      r->cf=NULL;
637    }
638    //if (r->algring!=NULL)
639    //{
640    //  rKill(r->algring);
641    //  r->algring=NULL;
642    //}
643  }
644}
Note: See TracBrowser for help on using the repository browser.