source: git/kernel/numbers.cc @ 19370c

spielwiese
Last change on this file since 19370c was 45ce3a4, checked in by Hans Schönemann <hannes@…>, 17 years ago
*hannes: cleanup minpoly git-svn-id: file:///usr/local/Singular/svn/trunk@9776 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 14.0 KB
Line 
1/*****************************************
2*  Computer Algebra System SINGULAR      *
3*****************************************/
4/* $Id: numbers.cc,v 1.5 2007-01-29 16:58:32 Singular Exp $ */
5
6/*
7* ABSTRACT: interface to coefficient aritmetics
8*/
9
10#include <string.h>
11#include <stdlib.h>
12#include "mod2.h"
13#include "structs.h"
14#include "febase.h"
15#include "kstd1.h"
16#include "numbers.h"
17#include "longrat.h"
18#include "longalg.h"
19#include "modulop.h"
20#include "gnumpfl.h"
21#include "gnumpc.h"
22#include "ring.h"
23#include "ffields.h"
24#include "shortfl.h"
25#ifdef HAVE_RING2TOM
26#include "rmodulo2m.h"
27#endif
28
29//static int characteristic = 0;
30extern int IsPrime(int p);
31
32void   (*nNew)(number *a);
33number (*nInit)(int i);
34number (*nPar)(int i);
35int    (*nParDeg)(number n);
36int    (*nSize)(number n);
37int    (*nInt)(number &n);
38numberfunc nMult, nSub, nAdd, nDiv, nIntDiv, nIntMod, nExactDiv;
39number (*nNeg)(number a);
40number (*nInvers)(number a);
41void   (*nNormalize)(number &a);
42number (*nCopy)(number a);
43number (*nRePart)(number a);
44number (*nImPart)(number a);
45BOOLEAN (*nGreater)(number a,number b);
46BOOLEAN (*nEqual)(number a,number b);
47BOOLEAN (*nIsZero)(number a);
48BOOLEAN (*nIsOne)(number a);
49BOOLEAN (*nIsMOne)(number a);
50BOOLEAN (*nGreaterZero)(number a);
51void    (*nWrite)(number &a);
52char *  (*nRead)(char *s,number *a);
53void    (*nPower)(number a, int i, number * result);
54number  (*nGcd)(number a, number b, const ring r);
55number  (*nLcm)(number a, number b, const ring r);
56char * (*nName)(number n);
57void   (*n__Delete)(number *a, const ring r);
58
59/*0 implementation*/
60number nNULL; /* the 0 as constant */
61
62
63n_Procs_s *cf_root=NULL;
64
65void   nDummy1(number* d) { *d=NULL; }
66void   ndDelete(number* d, const ring r) { *d=NULL; }
67void   ndInpMult(number &a, number b, const ring r)
68{
69  number n=n_Mult(a,b,r);
70  n_Delete(&a,r);
71  a=n;
72}
73
74#ifdef LDEBUG
75void   nDBDummy1(number* d,char *f, int l) { *d=NULL; }
76#endif
77
78void   nDummy2(number& d) { }
79
80char * ndName(number n) { return NULL; }
81
82number ndPar(int i) { return nInit(0); }
83
84number ndReturn0(number n) { return nInit(0); }
85
86int    ndParDeg(number n) { return 0; }
87
88number ndGcd(number a, number b, const ring r) { return r->cf->nInit(1); }
89
90number ndIntMod(number a, number b) { return nInit(0); }
91
92number ndGetDenom(number &n, const ring r) { return n_Init(1,r); }
93
94int ndSize(number a) { return (int)nIsZero(a)==FALSE; }
95
96number ndCopy(number a) { return a; }
97number nd_Copy(number a,const ring r) { return r->cf->nCopy(a); }
98
99/*2
100* init operations for characteristic c (complete==TRUE)
101* init nDelete    for characteristic c (complete==FALSE)
102*/
103void nSetChar(ring r)
104{
105  int c=rInternalChar(r);
106
107  n__Delete= r->cf->cfDelete;
108  if (rField_is_Extension(r))
109  {
110    naSetChar(c,r);
111  }
112#ifdef HAVE_RING2TOM
113  /*----------------------ring Z / 2^m----------------*/
114  else if (rField_is_Ring_2toM(r))
115  {
116    nr2mSetExp(c, r);
117  }
118#endif 
119  else if (rField_is_Zp(r))
120  /*----------------------char. p----------------*/
121  {
122    npSetChar(c, r);
123  }
124  /* -------------- GF(p^m) -----------------------*/
125  else if (rField_is_GF(r))
126  {
127    nfSetChar(c,r->parameter);
128  }
129  /* -------------- R -----------------------*/
130  //if (c==(-1))
131  else if (rField_is_R(r))
132  {
133  }
134  /* -------------- long R -----------------------*/
135  /* -------------- long C -----------------------*/
136  else if ((rField_is_long_R(r))
137  || (rField_is_long_C(r)))
138  {
139    setGMPFloatDigits(r->float_len,r->float_len2);
140  }
141#ifdef TEST
142  /* -------------- R -----------------------*/
143  //if (c==(-1))
144  else if (!rField_is_R(r) && !rField_is_Q(r))
145  {
146    WerrorS("unknown field");
147  }
148#endif
149  nNew   = r->cf->nNew;
150  nNormalize=r->cf->nNormalize;
151  nInit  = r->cf->nInit;
152  nPar   = r->cf->nPar;
153  nParDeg= r->cf->nParDeg;
154  nInt   = r->cf->nInt;
155  nAdd   = r->cf->nAdd;
156  nSub   = r->cf->nSub;
157  nMult  = r->cf->nMult;
158  nDiv   = r->cf->nDiv;
159  nExactDiv= r->cf->nExactDiv;
160  nIntDiv= r->cf->nIntDiv;
161  nIntMod= r->cf->nIntMod;
162  nNeg   = r->cf->nNeg;
163  nInvers= r->cf->nInvers;
164  nCopy  = r->cf->nCopy;
165  nGreater = r->cf->nGreater;
166  nEqual = r->cf->nEqual;
167  nIsZero = r->cf->nIsZero;
168  nIsOne = r->cf->nIsOne;
169  nIsMOne = r->cf->nIsMOne;
170  nGreaterZero = r->cf->nGreaterZero;
171  nWrite = r->cf->nWrite;
172  nRead = r->cf->nRead;
173  nPower = r->cf->nPower;
174  nGcd  = r->cf->nGcd;
175  nLcm  = r->cf->nLcm;
176  nName= r->cf->nName;
177  nSize  = r->cf->nSize;
178  nRePart = r->cf->nRePart;
179  nImPart = r->cf->nImPart;
180  nNULL=r->cf->nNULL;
181}
182
183/*2
184* init operations for ring r
185*/
186void nInitChar(ring r)
187{
188  int c=rInternalChar(r);
189  n_coeffType t=rFieldType(r);
190
191  if (rField_is_Extension(r))
192  {
193    if (r->algring==NULL)
194    {
195      int ch=-c;
196      if (c==1) ch=0;
197      r->algring=(ring) rDefault(ch,r->P,r->parameter);
198      //r->algring->ShortOut=r->ShortOut;
199      // includes: nInitChar(r->algring);
200    }
201  }
202
203  n_Procs_s *n=cf_root;
204  while((n!=NULL)
205    && ((n->nChar!=c) || (n->type!=t)))
206      n=n->next;
207  if (n==NULL)
208  {
209    n=(n_Procs_s*)omAlloc0(sizeof(n_Procs_s));
210    n->next=cf_root;
211    n->ref=1;
212    n->nChar=c;
213    n->type=t;
214    cf_root=n;
215  }
216  else if ((n->nChar==c) && (n->type==t))
217  {
218    n->ref++;
219    r->cf=n;
220    return;
221  }
222  else
223  {
224    WerrorS("nInitChar failed");
225    return;
226  }
227  r->cf=n;
228  n->nChar = c;
229  n->nPar  = ndPar;
230  n->nParDeg=ndParDeg;
231  n->nSize = ndSize;
232  n->n_GetDenom= ndGetDenom;
233  n->nName =  ndName;
234  n->nImPart=ndReturn0;
235  n->cfDelete= ndDelete;
236  n->nNew=nDummy1;
237  n->nInpMult=ndInpMult;
238  n->cfCopy=nd_Copy;
239  n->nIntMod=ndIntMod; /* dummy !! */
240  n->nNormalize=nDummy2;
241  n->nGcd  = ndGcd;
242  n->nLcm  = ndGcd; /* tricky, isn't it ?*/
243  if (rField_is_Extension(r))
244  {
245    //naInitChar(c,TRUE,r);
246    n->cfDelete = naDelete;
247    n-> nNew       = naNew;
248    n-> nNormalize = naNormalize;
249    n->nInit       = naInit;
250    n->nPar        = naPar;
251    n->nParDeg     = naParDeg;
252    n->nInt        = naInt;
253    n->nAdd        = naAdd;
254    n->nSub        = naSub;
255    n->nMult       = naMult;
256    n->nDiv        = naDiv;
257    n->nExactDiv   = naDiv;
258    n->nIntDiv     = naIntDiv;
259    n->nNeg        = naNeg;
260    n->nInvers     = naInvers;
261    n->nCopy       = naCopy;
262    n->cfCopy      = na_Copy;
263    n->nGreater    = naGreater;
264    n->nEqual      = naEqual;
265    n->nIsZero     = naIsZero;
266    n->nIsOne      = naIsOne;
267    n->nIsMOne     = naIsMOne;
268    n->nGreaterZero= naGreaterZero;
269    n->nWrite      = naWrite;
270    n->nRead       = naRead;
271    n->nPower      = naPower;
272    n->nGcd        = naGcd;
273    n->nLcm        = naLcm;
274    n->cfSetMap    = naSetMap;
275    n->nName       = naName;
276    n->nSize       = naSize;
277    n->n_GetDenom   = naGetDenom;
278#ifdef LDEBUG
279    //n->nDBTest     = naDBTest;
280#endif
281  }
282#ifdef HAVE_RING2TOM
283  /* -------------- Z/2^m ----------------------- */
284  else if (rField_is_Ring_2toM(r))
285  {
286     nr2mInitExp(c,r);
287     n->nInit  = nr2mInit;
288     n->nCopy  = ndCopy;
289     n->nInt   = nr2mInt;
290     n->nAdd   = nr2mAdd;
291     n->nSub   = nr2mSub;
292     n->nMult  = nr2mMult;
293     n->nDiv   = nr2mDiv;
294     n->nIntDiv       = nr2mIntDiv;
295     n->nExactDiv= nr2mDiv;
296     n->nNeg   = nr2mNeg;
297     n->nInvers= nr2mInvers;
298     n->nGreater = nr2mGreater;
299     n->nEqual = nr2mEqual;
300     n->nIsZero = nr2mIsZero;
301     n->nIsOne = nr2mIsOne;
302     n->nIsMOne = nr2mIsMOne;
303     n->nGreaterZero = nr2mGreaterZero;
304     n->nWrite = nr2mWrite;
305     n->nRead = nr2mRead;
306     n->nPower = nr2mPower;
307     n->cfSetMap = nr2mSetMap;
308     n->nNormalize = nDummy2;
309     n->nLcm          = nr2mLcm;
310     n->nGcd          = nr2mGcd;
311//     n->nGetUnit = nr2mGetUnit; //TODO OLIVER
312     n->nName= ndName;
313#ifdef LDEBUG
314//     n->nDBTest=nr2mDBTest;
315#endif
316  }
317#endif
318  else if (rField_is_Q(r))
319  {
320    n->cfDelete= nlDelete;
321    n->nNew   = nlNew;
322    n->nNormalize=nlNormalize;
323    n->nInit  = nlInit;
324    n->nInt   = nlInt;
325    n->nAdd   = nlAdd;
326    n->nSub   = nlSub;
327    n->nMult  = nlMult;
328    n->nInpMult=nlInpMult;
329    n->nDiv   = nlDiv;
330    n->nExactDiv= nlExactDiv;
331    n->nIntDiv= nlIntDiv;
332    n->nIntMod= nlIntMod;
333    n->nNeg   = nlNeg;
334    n->nInvers= nlInvers;
335    n->nCopy  = nlCopy;
336    n->nGreater = nlGreater;
337    n->nEqual = nlEqual;
338    n->nIsZero = nlIsZero;
339    n->nIsOne = nlIsOne;
340    n->nIsMOne = nlIsMOne;
341    n->nGreaterZero = nlGreaterZero;
342    n->nWrite = nlWrite;
343    n->nRead = nlRead;
344    n->nPower = nlPower;
345    n->nGcd  = nlGcd;
346    n->nLcm  = nlLcm;
347    n->cfSetMap = nlSetMap;
348    n->nSize  = nlSize;
349    n->n_GetDenom = nlGetDenom;
350#ifdef LDEBUG
351    //n->nDBTest=nlDBTest;
352#endif
353  }
354  else if (rField_is_Zp(r))
355  /*----------------------char. p----------------*/
356  {
357    npInitChar(c,r);
358    n->nInit  = npInit;
359    n->nInt   = npInt;
360    n->nAdd   = npAdd;
361    n->nSub   = npSub;
362    n->nMult  = npMult;
363    n->nDiv   = npDiv;
364    n->nExactDiv= npDiv;
365    n->nNeg   = npNeg;
366    n->nInvers= npInvers;
367    n->nCopy  = ndCopy;
368    n->nGreater = npGreater;
369    n->nEqual = npEqual;
370    n->nIsZero = npIsZero;
371    n->nIsOne = npIsOne;
372    n->nIsMOne = npIsMOne;
373    n->nGreaterZero = npGreaterZero;
374    n->nWrite = npWrite;
375    n->nRead = npRead;
376    n->nPower = npPower;
377    n->cfSetMap = npSetMap;
378    /* nName= ndName; */
379    /*nSize  = ndSize;*/
380#ifdef LDEBUG
381    //n->nDBTest=npDBTest;
382#endif
383#ifdef NV_OPS
384    if (c>NV_MAX_PRIME)
385    {
386      n->nMult  = nvMult;
387      n->nDiv   = nvDiv;
388      n->nExactDiv= nvDiv;
389      n->nInvers= nvInvers;
390    }
391#endif
392  }
393  /* -------------- GF(p^m) -----------------------*/
394  else if (rField_is_GF(r))
395  {
396    //nfSetChar(c,r->parameter);
397    n->nInit  = nfInit;
398    n->nPar   = nfPar;
399    n->nParDeg= nfParDeg;
400    n->nInt   = nfInt;
401    n->nAdd   = nfAdd;
402    n->nSub   = nfSub;
403    n->nMult  = nfMult;
404    n->nDiv   = nfDiv;
405    n->nExactDiv= nfDiv;
406    n->nNeg   = nfNeg;
407    n->nInvers= nfInvers;
408    n->nCopy  = ndCopy;
409    n->nGreater = nfGreater;
410    n->nEqual = nfEqual;
411    n->nIsZero = nfIsZero;
412    n->nIsOne = nfIsOne;
413    n->nIsMOne = nfIsMOne;
414    n->nGreaterZero = nfGreaterZero;
415    n->nWrite = nfWrite;
416    n->nRead = nfRead;
417    n->nPower = nfPower;
418    n->cfSetMap = nfSetMap;
419    n->nName= nfName;
420    /*nSize  = ndSize;*/
421#ifdef LDEBUG
422    //n->nDBTest=nfDBTest;
423#endif
424  }
425  /* -------------- R -----------------------*/
426  //if (c==(-1))
427  else if (rField_is_R(r))
428  {
429    n->nInit  = nrInit;
430    n->nInt   = nrInt;
431    n->nAdd   = nrAdd;
432    n->nSub   = nrSub;
433    n->nMult  = nrMult;
434    n->nDiv   = nrDiv;
435    n->nExactDiv= nrDiv;
436    n->nNeg   = nrNeg;
437    n->nInvers= nrInvers;
438    n->nCopy  = ndCopy;
439    n->nGreater = nrGreater;
440    n->nEqual = nrEqual;
441    n->nIsZero = nrIsZero;
442    n->nIsOne = nrIsOne;
443    n->nIsMOne = nrIsMOne;
444    n->nGreaterZero = nrGreaterZero;
445    n->nWrite = nrWrite;
446    n->nRead = nrRead;
447    n->nPower = nrPower;
448    n->cfSetMap=nrSetMap;
449    /* nName= ndName; */
450    /*nSize  = ndSize;*/
451#ifdef LDEBUG
452    //n->nDBTest=nrDBTest;
453#endif
454  }
455  /* -------------- long R -----------------------*/
456  else if (rField_is_long_R(r))
457  {
458    n->cfDelete= ngfDelete;
459    n->nNew=ngfNew;
460    n->nInit  = ngfInit;
461    n->nInt   = ngfInt;
462    n->nAdd   = ngfAdd;
463    n->nSub   = ngfSub;
464    n->nMult  = ngfMult;
465    n->nDiv   = ngfDiv;
466    n->nExactDiv= ngfDiv;
467    n->nNeg   = ngfNeg;
468    n->nInvers= ngfInvers;
469    n->nCopy  = ngfCopy;
470    n->nGreater = ngfGreater;
471    n->nEqual = ngfEqual;
472    n->nIsZero = ngfIsZero;
473    n->nIsOne = ngfIsOne;
474    n->nIsMOne = ngfIsMOne;
475    n->nGreaterZero = ngfGreaterZero;
476    n->nWrite = ngfWrite;
477    n->nRead = ngfRead;
478    n->nPower = ngfPower;
479    n->cfSetMap=ngfSetMap;
480    n->nName= ndName;
481    n->nSize  = ndSize;
482#ifdef LDEBUG
483    //n->nDBTest=ngfDBTest;
484#endif
485  }
486  /* -------------- long C -----------------------*/
487  else if (rField_is_long_C(r))
488  {
489    n->cfDelete= ngcDelete;
490    n->nNew=ngcNew;
491    n->nNormalize=nDummy2;
492    n->nInit  = ngcInit;
493    n->nInt   = ngcInt;
494    n->nAdd   = ngcAdd;
495    n->nSub   = ngcSub;
496    n->nMult  = ngcMult;
497    n->nDiv   = ngcDiv;
498    n->nExactDiv= ngcDiv;
499    n->nNeg   = ngcNeg;
500    n->nInvers= ngcInvers;
501    n->nCopy  = ngcCopy;
502    n->nGreater = ngcGreater;
503    n->nEqual = ngcEqual;
504    n->nIsZero = ngcIsZero;
505    n->nIsOne = ngcIsOne;
506    n->nIsMOne = ngcIsMOne;
507    n->nGreaterZero = ngcGreaterZero;
508    n->nWrite = ngcWrite;
509    n->nRead = ngcRead;
510    n->nPower = ngcPower;
511    n->cfSetMap=ngcSetMap;
512    n->nPar=ngcPar;
513    n->nRePart=ngcRePart;
514    n->nImPart=ngcImPart;
515    /*nSize  = ndSize;*/
516#ifdef LDEBUG
517    //n->nDBTest=ngcDBTest;
518#endif
519  }
520#ifdef TEST
521  else
522  {
523    WerrorS("unknown field");
524  }
525#endif
526  if (!errorreported)
527  {
528    n->nNULL=n->nInit(0);
529    if (n->nRePart==NULL)
530      n->nRePart=n->nCopy;
531    if (n->nIntDiv==NULL)
532      n->nIntDiv=n->nDiv;
533  }
534}
535
536void nKillChar(ring r)
537{
538  if (r!=NULL)
539  {
540    if (r->cf!=NULL)
541    {
542      r->cf->ref--;
543      if (r->cf->ref<=0)
544      {
545        n_Procs_s tmp;
546        n_Procs_s* n=&tmp;
547        tmp.next=cf_root;
548        while((n->next!=NULL) && (n->next!=r->cf)) n=n->next;
549        if (n->next==r->cf)
550        {
551          n->next=n->next->next;
552          if (cf_root==r->cf) cf_root=n->next;
553          r->cf->cfDelete(&(r->cf->nNULL),r);
554          switch(r->cf->type)
555          {
556            case n_Zp:
557                 #ifdef HAVE_DIV_MOD
558                 if (r->cf->npInvTable!=NULL)
559                 omFreeSize( (ADDRESS)r->cf->npInvTable,
560                             r->cf->npPrimeM*sizeof(CARDINAL) );
561                 #else
562                 if (r->cf->npExpTable!=NULL)
563                 {
564                   omFreeSize( (ADDRESS)r->cf->npExpTable,
565                               r->cf->npPrimeM*sizeof(CARDINAL) );
566                   omFreeSize( (ADDRESS)r->cf->npLogTable,
567                               r->cf->npPrimeM*sizeof(CARDINAL) );
568                 }
569                 #endif
570                 break;
571            case n_Zp_a:
572            case n_Q_a:
573                 {
574                   number n=r->minpoly;
575                   if (n!=NULL)
576                   {
577                     r->minpoly=NULL;
578                     if (r==currRing) naMinimalPoly=NULL;
579                     naDelete(&n,r);
580                   }
581                 }
582                 break;
583
584            default:
585                 break;
586          }
587          omFreeSize((ADDRESS)r->cf, sizeof(n_Procs_s));
588          r->cf=NULL;
589        }
590        else
591        {
592          WarnS("cf_root list destroyed");
593        }
594      }
595    }
596    if (r->algring!=NULL)
597    {
598      rKill(r->algring);
599      r->algring=NULL;
600    }
601  }
602}
Note: See TracBrowser for help on using the repository browser.