source: git/Singular/longalg.cc @ 32df82

fieker-DuValspielwiese
Last change on this file since 32df82 was 32df82, checked in by Hans Schönemann <hannes@…>, 27 years ago
* hannes: removed rcsid and Log: entries, added assignment module=poly corected type conversion int->module git-svn-id: file:///usr/local/Singular/svn/trunk@128 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 42.4 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: longalg.cc,v 1.3 1997-04-02 15:07:21 Singular Exp $ */
5/*
6* ABSTRACT:   algebraic numbers
7*/
8
9#include <stdio.h>
10#include <string.h>
11#include "mod2.h"
12#include "tok.h"
13#include "mmemory.h"
14#include "febase.h"
15#include "longrat.h"
16#include "modulop.h"
17#include "numbers.h"
18#include "polys.h"
19#include "ideals.h"
20#include "ipid.h"
21#include "longalg.h"
22#define  FEHLER1 1
23#define  FEHLER2 1
24#define  FEHLER3 1
25
26struct snaIdeal
27{
28  int anz;
29  alg *liste;
30};
31typedef struct snaIdeal * naIdeal;
32
33naIdeal naI=NULL;
34
35//#define RECA_SIZE (sizeof(alg)+sizeof(number))
36alg naMinimalPoly;
37int naNumbOfPar;
38char **naParNames;
39static int naIsChar0;
40static int naPrimeM;
41
42#ifdef LDEBUG
43#define naTest(a) naDBTest(a,__FILE__,__LINE__)
44BOOLEAN naDBTest(number a, char *f,int l);
45#else
46#define naTest(a)
47#endif
48
49number (*naMap)(number from);
50/* procedure variables */
51static numberfunc
52                nacMult, nacSub, nacAdd, nacDiv, nacGcd, nacLcm;
53#ifdef LDEBUG
54static void     (*nacDBDelete)(number *a,char *f,int l);
55#define         nacDelete(A) nacDBDelete(A,__FILE__,__LINE__)
56#else
57static void     (*nacDelete)(number *a);
58#endif
59       number   (*nacInit)(int i);
60static int      (*nacInt)(number &n);
61       void     (*nacNormalize)(number &a);
62static number   (*nacNeg)(number a);
63static void     (*nacWrite)(number &a);
64       number   (*nacCopy)(number a);
65static number   (*nacInvers)(number a);
66       BOOLEAN  (*nacIsZero)(number a);
67static BOOLEAN  (*nacIsOne)(number a);
68static BOOLEAN  (*nacIsMOne)(number a);
69static BOOLEAN  (*nacGreaterZero)(number a);
70static char   * (*nacRead) (char *s, number *a);
71extern void     napDelete(alg *p);
72static alg napRedp(alg q);
73static alg napTailred(alg q);
74static BOOLEAN napDivPoly(alg p, alg q);
75static int napExpi(int i, alg a, alg b);
76
77/*2
78*  sets the appropriate operators
79*/
80void naSetChar(int i, BOOLEAN complete, char ** param, int pars)
81{
82  if (naI!=NULL)
83  {
84    int j;
85    for (j=naI->anz-1; j>=0; j--)
86       napDelete (&naI->liste[j]);
87    Free((ADDRESS)naI->liste,naI->anz*sizeof(alg));
88    Free((ADDRESS)naI,sizeof(*naI));
89    naI=NULL;
90  }
91  naMap = naCopy;
92  naMinimalPoly = NULL;
93  naParNames=param;
94  naNumbOfPar=pars;
95  if (i == 1)
96  {
97    naIsChar0 = 1;
98#ifdef LDEBUG
99    nacDBDelete      = nlDBDelete;
100#else
101    nacDelete      = nlDelete;
102#endif
103    if (complete)
104    {
105      nacInit        = nlInit;
106      nacInt         = nlInt;
107      nacCopy        = nlCopy;
108      nacAdd         = nlAdd;
109      nacSub         = nlSub;
110      nacMult        = nlMult;
111      nacDiv         = nlDiv;
112      nacInvers      = nlInvers;
113      nacNormalize   = nlNormalize;
114      nacNeg         = nlNeg;
115      nacIsZero      = nlIsZero;
116      nacRead        = nlRead;
117      nacWrite       = nlWrite;
118      nacGreaterZero = nlGreaterZero;
119      nacIsOne       = nlIsOne;
120      nacIsMOne      = nlIsMOne;
121      nacGcd         = nlGcd;
122      nacLcm         = nlLcm;
123    }
124  }
125  else if (i < 0)
126  {
127    naIsChar0 = 0;
128#ifdef LDEBUG
129    nacDBDelete    = nDBDummy1;
130#else
131    nacDelete      = nDummy1;
132#endif
133    if (complete)
134    {
135      npSetChar(-i);
136      nacInit        = npInit;
137      nacInt         = npInt;
138      nacCopy        = npCopy;
139      nacAdd         = npAdd;
140      nacSub         = npSub;
141      nacMult        = npMult;
142      nacDiv         = npDiv;
143      nacInvers      = npInvers;
144      nacNormalize   = nDummy2;
145      nacNeg         = npNeg;
146      nacIsZero      = npIsZero;
147      nacRead        = npRead;
148      nacWrite       = npWrite;
149      nacGreaterZero = npGreaterZero;
150      nacIsOne       = npIsOne;
151      nacIsMOne      = npIsMOne;
152      nacGcd         = ndGcd;
153    }
154  }
155#ifdef TEST
156  else
157  {
158    Print("hier ging was schief:c=%d compl=%d param=%d\n",i,complete,param);
159  }
160#endif
161}
162
163/*============= procedure for polynomials: napXXXX =======================*/
164
165/*3
166* creates  an alg poly
167*/
168static alg napInit(int i)
169{
170  alg a = (alg)Alloc0(RECA_SIZE + naNumbOfPar * sizeof(int));
171
172  a->ko = nacInit(i);
173  if (!nacIsZero(a->ko))
174  {
175    return a;
176  }
177  else
178  {
179    Free((ADDRESS)a, RECA_SIZE + naNumbOfPar * sizeof(int));
180    return NULL;
181  }
182}
183
184/*3
185* creates  an alg poly
186*/
187alg napInitz(number z)
188{
189  alg a = (alg)Alloc0(RECA_SIZE + naNumbOfPar * sizeof(int));
190
191  a->ko = z;
192  return a;
193}
194
195/*3
196*  deletes head of  an alg poly
197*/
198static void napDelete1(alg *p)
199{
200  alg h = *p;
201
202  if (h!=NULL)
203  {
204    *p = h->ne;
205    nacDelete(&(h->ko));
206    Free((ADDRESS)h, RECA_SIZE + naNumbOfPar * sizeof(int));
207  }
208}
209
210/*3
211*  delete an alg poly
212*/
213void napDelete(alg *p)
214{
215  alg w, h = *p;
216
217  while (h!=NULL)
218  {
219    w = h;
220    h = h->ne;
221    nacDelete(&(w->ko));
222    Free((ADDRESS)w, RECA_SIZE + naNumbOfPar * sizeof(int));
223  }
224  *p = NULL;
225}
226
227
228/*3
229* copy a alg. poly
230*/
231static alg napCopy(alg p)
232{
233  alg w, a;
234
235  if (p==NULL) return NULL;
236  mmTestP(p,RECA_SIZE+naNumbOfPar*sizeof(int));
237  a = w = (alg)Alloc(RECA_SIZE + naNumbOfPar * sizeof(int));
238  memcpy(w->e, p->e, naNumbOfPar * sizeof(int));
239  w->ko = nacCopy(p->ko);
240  loop
241  {
242    p=p->ne;
243    if (p==NULL) break;
244    mmTestP(p,RECA_SIZE+naNumbOfPar*sizeof(int));
245    a->ne = (alg)Alloc(RECA_SIZE + naNumbOfPar * sizeof(int));
246    a = a->ne;
247    memcpy(a->e, p->e, naNumbOfPar * sizeof(int));
248    a->ko = nacCopy(p->ko);
249  }
250  a->ne = NULL;
251  return w;
252}
253
254/*3
255* copy a alg. poly
256*/
257static alg napCopyNeg(alg p)
258{
259  alg w, a;
260
261  if (p==NULL) return NULL;
262  mmTestP(p,RECA_SIZE+naNumbOfPar*sizeof(int));
263  a = w = (alg)Alloc(RECA_SIZE + naNumbOfPar * sizeof(int));
264  memcpy(w->e, p->e, naNumbOfPar * sizeof(int));
265  w->ko = nacNeg(nacCopy(p->ko));
266  loop
267  {
268    p=p->ne;
269    if (p==NULL) break;
270    mmTestP(p,RECA_SIZE+naNumbOfPar*sizeof(int));
271    a->ne = (alg)Alloc(RECA_SIZE + naNumbOfPar * sizeof(int));
272    a = a->ne;
273    memcpy(a->e, p->e, naNumbOfPar * sizeof(int));
274    a->ko = nacNeg(nacCopy(p->ko));
275  }
276  a->ne = NULL;
277  return w;
278}
279
280/*3
281* Compare exponents of p and q
282*/
283static int  napComp(alg p, alg q)
284{
285  int  i = 0;
286
287  mmTestP(p,RECA_SIZE+naNumbOfPar*sizeof(int));
288  mmTestP(q,RECA_SIZE+naNumbOfPar*sizeof(int));
289  while (p->e[i] == q->e[i])
290  {
291    i++;
292    if (i >= naNumbOfPar)
293      return 0;
294  }
295  if (p->e[i]  >  q->e[i])
296    return 1;
297  else
298    return - 1;
299}
300
301/*3
302*  addition of alg. polys
303*/
304alg napAdd(alg p1, alg p2)
305{
306  alg a1, p, a2, a;
307  int  c;  number t;
308
309  if (p1==NULL) return p2;
310  if (p2==NULL) return p1;
311  mmTestP(p1,RECA_SIZE+naNumbOfPar*sizeof(int));
312  mmTestP(p2,RECA_SIZE+naNumbOfPar*sizeof(int));
313  a1 = p1;
314  a2 = p2;
315  a = p  = (alg)Alloc(RECA_SIZE + naNumbOfPar * sizeof(int));
316  loop
317  {
318    c = napComp(a1, a2);
319    if (c == 1)
320    {
321      a = a->ne = a1;
322      a1 = a1->ne;
323      if (!a1)
324      {
325        a->ne= a2;
326        break;
327      }
328      mmTestP(a1,RECA_SIZE+naNumbOfPar*sizeof(int));
329    }
330    else if (c == -1)
331    {
332      a = a->ne = a2;
333      a2 = a2->ne;
334      if (!a2)
335      {
336        a->ne = a1;
337        break;
338      }
339      mmTestP(a2,RECA_SIZE+naNumbOfPar*sizeof(int));
340    }
341    else
342    {
343      t = nacAdd(a1->ko, a2->ko);
344      napDelete1(&a2);
345      if (nacIsZero(t))
346      {
347        nacDelete(&t);
348        napDelete1(&a1);
349      }
350      else
351      {
352        nacDelete(&(a1->ko));
353        a1->ko = t;
354        a = a->ne = a1;
355        a1 = a1->ne;
356      }
357      if (a1==NULL)
358      {
359        a->ne = a2;
360        break;
361      }
362      else if (a2==NULL)
363      {
364        a->ne = a1;
365        break;
366      }
367      mmTestP(a1,RECA_SIZE+naNumbOfPar*sizeof(int));
368      mmTestP(a2,RECA_SIZE+naNumbOfPar*sizeof(int));
369    }
370  }
371  a = p->ne;
372  Free((ADDRESS)p, RECA_SIZE + naNumbOfPar * sizeof(int));
373  return a;
374}
375
376
377/*3
378* multiply a alg. poly by -1
379*/
380static alg napNeg(alg a)
381{
382  alg p = a;
383
384  while (p!=NULL)
385  {
386    mmTestP(p,RECA_SIZE+naNumbOfPar*sizeof(int));
387    p->ko = nacNeg(p->ko);
388    p = p->ne;
389  }
390  return a;
391}
392
393/*3
394* returns ph * z
395*/
396static void napMultN(alg p, number z)
397{
398  number t;
399
400  while (p!=NULL)
401  {
402    mmTestP(p,RECA_SIZE+naNumbOfPar*sizeof(int));
403    t = nacMult(p->ko, z);
404    nacNormalize(t);
405    nacDelete(&p->ko);
406    p->ko = t;
407    p = p->ne;
408  }
409}
410
411/*3
412* update the polynomial a by multipying it by
413* the (number) coefficient
414* and the exponent vector (of) exp (a well initialized polynomial)
415*/
416static void napMultT(alg a, alg exp)
417{
418  int  i;
419  number t, h;
420  if (a==NULL)
421    return;
422  h = exp->ko;
423  if (nacIsOne(h))
424  {
425    do
426    {
427      for (i = naNumbOfPar - 1; i >= 0; i--)
428        a->e[i] += exp->e[i];
429      a = a->ne;
430    }
431    while (a!=NULL);
432  }
433  else
434  {
435    h = nacNeg(h);
436    if (nacIsOne(h))
437    {
438      h = nacNeg(h);
439      do
440      {
441        a->ko = nacNeg(a->ko);
442        for (i = naNumbOfPar - 1; i >= 0; i--)
443          a->e[i] += exp->e[i];
444        a = a->ne;
445      }
446      while (a!=NULL);
447    }
448    else
449    {
450      h = nacNeg(h);
451      do
452      {
453        t = nacMult(a->ko, h);
454        nacDelete(&(a->ko));
455        a->ko = t;
456        for (i = naNumbOfPar - 1; i >= 0; i--)
457          a->e[i] += exp->e[i];
458        a = a->ne;
459      }
460      while (a!=NULL);
461    }
462  }
463}
464
465/*3
466* multiplication of alg. polys
467* multiply p1 with p2, p1 and p2 are destroyed
468* do not put attention on speed: the procedure is only used in the interpreter
469*/
470static alg napMult(alg p1, alg p2)
471{
472  alg a1, a2, b, h;
473  if (p1==NULL)
474  {
475    napDelete(&p2);
476    return NULL;
477  }
478  if (p2==NULL)
479  {
480    napDelete(&p1);
481    return NULL;
482  }
483  b = NULL;
484  a1 = p1;
485  a2 = p2;
486  if (a2->ne!=NULL)
487  {
488    if (a1->ne!=NULL)
489    {
490      do
491      {
492        a1 = a1->ne;
493        a2 = a2->ne;
494      }
495      while ((a1!=NULL) && (a2!=NULL));
496      if (a1!=NULL)
497      {
498        a1 = p1;
499        a2 = p2;
500      }
501      else
502      {
503        a1 = p2;
504        a2 = p1;
505      }
506      do
507      {
508        if (a2->ne!=NULL)
509          h = napCopy(a1);
510        else
511          h = a1;
512        napMultT(h, a2);
513        b = napAdd(h, b);
514        napDelete1(&a2);
515      }
516      while (a2!=NULL);
517    }
518    else
519    {
520      napMultT(a2, a1);
521      b = napAdd(a2, b);
522      napDelete1(&a1);
523    }
524  }
525  else
526  {
527    napMultT(a1, a2);
528    b = napAdd(a1, b);
529    napDelete1(&a2);
530  }
531  return b;
532}
533
534/*3
535*  division with rest; f = g*q + r, returns r and destroy f
536*/
537static alg napRemainder(alg f, const alg  g)
538{
539  alg a, h, qq;
540
541  qq = (alg)Alloc(RECA_SIZE + naNumbOfPar * sizeof(int));
542  qq->ne = NULL;
543  a = f;
544  do
545  {
546    mmTestP(a,RECA_SIZE+naNumbOfPar*sizeof(int));
547    mmTestP(g,RECA_SIZE+naNumbOfPar*sizeof(int));
548    qq->e[0] = a->e[0] - g->e[0];
549    qq->ko = nacDiv(a->ko, g->ko);
550    qq->ko = nacNeg(qq->ko);
551    h = napCopy(g);
552    napMultT(h, qq);
553    nacDelete(&(qq->ko));
554    a = napAdd(a, h);
555  }
556  while ((a!=NULL) && (a->e[0] >= g->e[0]));
557  Free((ADDRESS)qq, RECA_SIZE + naNumbOfPar * sizeof(int));
558  return a;
559}
560
561/*3
562*  division with rest; f = g*q + r,  destroy f
563*/
564static void napDivMod(alg f, alg  g, alg *q, alg *r)
565{
566  alg a, h, b, qq;
567
568  qq = (alg)Alloc(RECA_SIZE + naNumbOfPar * sizeof(int));
569  qq->ne = b = NULL;
570  a = f;
571  do
572  {
573    mmTestP(a,RECA_SIZE+naNumbOfPar*sizeof(int));
574    mmTestP(g,RECA_SIZE+naNumbOfPar*sizeof(int));
575    qq->e[0] = a->e[0] - g->e[0];
576    qq->ko = nacDiv(a->ko, g->ko);
577    b = napAdd(b, napCopy(qq));
578    qq->ko = nacNeg(qq->ko);
579    h = napCopy(g);
580    napMultT(h, qq);
581    nacDelete(&(qq->ko));
582    a = napAdd(a, h);
583  }
584  while ((a!=NULL) && (a->e[0] >= g->e[0]));
585  Free((ADDRESS)qq, RECA_SIZE + naNumbOfPar * sizeof(int));
586  *q = b;
587  *r = a;
588}
589
590/*3
591*  returns z with z*x mod c = 1
592*/
593static alg napInvers(alg x, const alg c)
594{
595  alg y, r, qa, qn, q;
596  number t, h;
597
598  if (x->e[0] >= c->e[0])
599    x = napRemainder(x, c);
600  if (x->e[0]==0)
601  {
602    if (!nacIsOne(x->ko))
603    {
604      h = nacInit(1);
605      t = nacDiv(h, x->ko);
606      nacDelete(&(x->ko));
607      nacDelete(&h);
608      x->ko = t;
609    }
610    return x;
611  }
612  y = napCopy(c);
613  napDivMod(y, x, &qa, &r);
614  if (r->e[0]==0)
615  {
616    h = nacInit(-1);
617    t = nacDiv(h, r->ko);
618    nacNormalize(t);
619    napMultN(qa, t);
620    nacDelete(&h);
621    nacDelete(&t);
622    napDelete(&x);
623    napDelete(&r);
624    return qa;
625  }
626  y = x;
627  x = r;
628  napDivMod(y, x, &q, &r);
629  if (r->e[0]==0)
630  {
631    q = napMult(q, qa);
632    q = napAdd(q, napInit(1));
633    h = nacInit(1);
634    t = nacDiv(h, r->ko);
635    napMultN(q, t);
636    nacDelete(&h);
637    nacDelete(&t);
638    napDelete(&x);
639    napDelete(&r);
640    if (q->e[0] >= c->e[0])
641      q = napRemainder(q, c);
642    return q;
643  }
644  q = napMult(q, napCopy(qa));
645  q = napAdd(q, napInit(1));
646  qa = napNeg(qa);
647  loop
648  {
649    y = x;
650    x = r;
651    napDivMod(y, x, &qn, &r);
652    if (r->e[0]==0)
653    {
654      q = napMult(q, qn);
655      q = napNeg(q);
656      q = napAdd(q, qa);
657      h = nacInit(1);
658      t = nacDiv(h, r->ko);
659      napMultN(q, t);
660      nacDelete(&h);
661      nacDelete(&t);
662      napDelete(&x);
663      napDelete(&r);
664      if (q->e[0] >= c->e[0])
665        q = napRemainder(q, c);
666      return q;
667    }
668    y = q;
669    q = napMult(napCopy(q), qn);
670    q = napNeg(q);
671    q = napAdd(q, qa);
672    qa = y;
673  }
674}
675
676/*3
677* the degree of an alg poly (used for test of "constant" et al.)
678*/
679static int  napDeg(alg p)
680{
681  int  d = 0, i;
682
683  mmTestP(p,RECA_SIZE+naNumbOfPar*sizeof(int));
684  for (i = naNumbOfPar-1; i>=0; i--)
685    d += p->e[i];
686  return d;
687}
688
689/*3
690*writes a polynomial number
691*/
692static void napWrite(alg p)
693{
694  int  i;
695  if (p==NULL)
696    StringAppendS("0");
697  else
698  {
699    StringAppendS("(");
700    loop
701    {
702      BOOLEAN wroteCoeff=FALSE;
703      mmTestP(p,RECA_SIZE+naNumbOfPar*sizeof(int));
704      if ((napDeg(p)==0)
705      || (!nacIsOne(p->ko)))
706      {
707        nacWrite(p->ko);
708        wroteCoeff=(pShortOut==0);
709      }
710      for (i = 0; i <= naNumbOfPar - 1; i++)
711      {
712        if (p->e[i] > 0)
713        {
714          if (wroteCoeff)
715            StringAppendS("*");
716          else
717            wroteCoeff=(pShortOut==0);
718          StringAppend(naParNames[i]);
719          if (p->e[i] > 1)
720          {
721            if (pShortOut == 0)
722              StringAppendS("^");
723            StringAppend("%d", p->e[i]);
724          }
725        }
726      }
727      p = p->ne;
728      if (p==NULL)
729        break;
730      if (nacGreaterZero(p->ko))
731        StringAppendS("+");
732    }
733    StringAppendS(")");
734  }
735}
736
737
738static char *napHandleMons(char *s, int i, int *ex)
739{
740  int  j;
741  if (strncmp(s,naParNames[i],strlen(naParNames[i]))==0)
742  {
743    s+=strlen(naParNames[i]);
744    if ((*s >= '0') && (*s <= '9'))
745    {
746      s = eati(s, &j);
747      ex[i] += j;
748    }
749    else
750      ex[i]++;
751  }
752  return s;
753}
754
755/*3  reads a monomial  */
756static char  *napRead(char *s, alg *b)
757{
758  alg a;
759  int  i;
760  a = (alg)Alloc0(RECA_SIZE + naNumbOfPar * sizeof(int));
761  if ((*s >= '0') && (*s <= '9'))
762  {
763    s = nacRead(s, &(a->ko));
764    if (nacIsZero(a->ko))
765    {
766      napDelete1(&a);
767      *b = NULL;
768      return s;
769    }
770  }
771  else
772    a->ko = nacInit(1);
773  i = 0;
774  char  *olds;
775  loop
776  {
777    olds = s;
778    s = napHandleMons(s, i, a->e);
779    if (olds == s)
780      i++;
781    else
782      i = 0;
783    if ((*s == '\0') || (i >= naNumbOfPar))
784      break;
785  }
786  *b = a;
787  return s;
788}
789
790static int napExp(alg a, alg b)
791{
792  while (a->ne!=NULL) a = a->ne;
793  int m = a->e[0];
794  if (m==0) return m;
795  while (b->ne!=NULL) b = b->ne;
796  if (m > b->e[0]) m = b->e[0];
797  return m;
798}
799
800#if FEHLER2
801/*meins
802* finds the smallest i-th exponent in a and b
803* used to find it in a fraction
804*/
805static int napExpi(int i, alg a, alg b)
806{
807  if (a==NULL || b==NULL) return 0;
808  int m = a->e[i];
809  if (m==0) return 0;
810  while (a->ne != NULL)
811  {
812    a = a->ne;
813    if (m > a->e[i])
814    {
815      m = a->e[i];
816      if (m==0) return 0;
817    }
818  }
819  do
820  {
821    if (m > b->e[i])
822    {
823      m = b->e[i];
824      if (m==0) return 0;
825    }
826    b = b->ne;
827  }
828  while (b != NULL);
829  return m;
830}
831#endif
832
833static void napContent(alg ph)
834{
835  number h,d;
836  alg p;
837
838  p = ph;
839  if (nacIsOne(p->ko))
840    return;
841  h = nacCopy(p->ko);
842  p = p->ne;
843  do
844  {
845    d=nacGcd(p->ko, h);
846    if(nacIsOne(d))
847    {
848      nacDelete(&h);
849      nacDelete(&d);
850      return;
851    }
852    nacDelete(&h);
853    h = d;
854    p = p->ne;
855  }
856  while (p!=NULL);
857  h = nacInvers(d);
858  nacDelete(&d);
859  p = ph;
860  while (p!=NULL)
861  {
862    d = nacMult(p->ko, h);
863    nacDelete(&(p->ko));
864    p->ko = d;
865    p = p->ne;
866  }
867  nacDelete(&h);
868}
869
870static void napCleardenom(alg ph)
871{
872  number d, h;
873  alg p;
874
875  if (!naIsChar0)
876    return;
877  p = ph;
878  h = nacInit(1);
879  while (p!=NULL)
880  {
881    d = nacLcm(h, p->ko);
882    nacDelete(&h);
883    h = d;
884    p = p->ne;
885  }
886  if(!nacIsOne(h))
887  {
888    p = ph;
889    while (p!=NULL)
890    {
891      d=nacMult(h, p->ko);
892      nacDelete(&(p->ko));
893      p->ko = d;
894      p = p->ne;
895    }
896    nacDelete(&h);
897  }
898  napContent(ph);
899}
900
901static alg napGcd0(alg a, alg b)
902{
903  number x, y;
904  if (!naIsChar0)
905    return napInit(1);
906  x = nacCopy(a->ko);
907  if (nacIsOne(x))
908    return napInitz(x);
909  while (a->ne!=NULL)
910  {
911    a = a->ne;
912    y = nacGcd(x, a->ko);
913    nacDelete(&x);
914    x = y;
915    if (nacIsOne(x))
916      return napInitz(x);
917  }
918  do
919  {
920    y = nacGcd(x, b->ko);
921    nacDelete(&x);
922    x = y;
923    if (nacIsOne(x))
924      return napInitz(x);
925    b = b->ne;
926  }
927  while (b!=NULL);
928  return napInitz(x);
929}
930
931/*3
932* result =gcd(a,b)
933*/
934static alg napGcd(alg a, alg b)
935{
936  int i;
937  alg g, x, y, h;
938  if ((a==NULL)
939  || ((a->ne==NULL)&&(nacIsZero(a->ko))))
940  {
941    if ((b==NULL)
942    || ((b->ne==NULL)&&(nacIsZero(b->ko))))
943    {
944      return napInit(1);
945    }
946    return napCopy(b);
947  }
948  if (naMinimalPoly != NULL)
949  {
950    if (a->e[0] >= b->e[0])
951    {
952      x = a;
953      y = b;
954    }
955    else
956    {
957      x = b;
958      y = a;
959    }
960    g = napGcd0(x, y);
961    if (y->ne==NULL)
962    {
963      g->e[0] = napExp(x, y);
964      return g;
965    }
966    x = napCopy(x);
967    y = napCopy(y);
968    loop
969    {
970      h = napRemainder(x, y);
971      if (h==NULL)
972      {
973        napCleardenom(y);
974        if (!nacIsOne(g->ko)) napMultN(y, g->ko);
975        napDelete1(&g);
976        return y;
977      }
978      else if (h->ne==NULL)
979        break;
980      x = y;
981      y = h;
982    }
983    napDelete(&y);
984    napDelete1(&h);
985    g->e[0] = napExp(a, b);
986    return g;
987  }
988  x = (alg)Alloc0(RECA_SIZE + naNumbOfPar * sizeof(int));
989  g=a;
990  h=b;
991  x = napGcd0(g,h);
992  for (i=(naNumbOfPar-1); i>=0; i--)
993    x->e[i] = napExpi(i,a,b);
994  return x;
995}
996
997
998number napLcm(alg a)
999{
1000  number d, h;
1001  alg b;
1002
1003  h = nacInit(1);
1004  b = a;
1005  while (b!=NULL)
1006  {
1007    d = nacLcm(h, b->ko);
1008    nacDelete(&h);
1009    h = d;
1010    b = b->ne;
1011  }
1012  return h;
1013}
1014
1015
1016/*2
1017* meins  (for reduction in algebraic extension)
1018* checks if head of p divides head of q
1019* doesn't delete p and q
1020*/
1021BOOLEAN napDivPoly (alg p, alg q)
1022{
1023  int j=0;   /* evtl. von naNumber.. -1 abwaerts zaehlen */
1024
1025  while (p->e[j] <= q->e[j])
1026  {
1027    j++;
1028    if (j >= naNumbOfPar)
1029      return 1;
1030  }
1031  return 0;
1032}
1033
1034
1035/*2
1036* meins  (for reduction in algebraic extension)
1037* Normalform of poly with naI
1038* changes q and returns it
1039*/
1040alg napRedp (alg q)
1041{
1042  alg h = (alg)Alloc0(RECA_SIZE + naNumbOfPar * sizeof(int));
1043  int i=0,j;
1044
1045  loop
1046  {
1047    if (napDivPoly (naI->liste[i], q))
1048    {
1049      mmTestP((ADDRESS)q,RECA_SIZE + naNumbOfPar * sizeof(int));
1050      StringSetS("");
1051      napWrite(q);
1052      napWrite(naI->liste[i]);
1053      Print(StringAppendS("\n"));
1054      /* h = lt(q)/lt(naI->liste[i])*/
1055      h->ko = nacCopy(q->ko);
1056      for (j=naNumbOfPar-1; j>=0; j--)
1057        h->e[j] = q->e[j] - naI->liste[i]->e[j];
1058      h = napMult (h, napCopy(naI->liste[i]));
1059      h = napNeg (h);
1060      q = napAdd (q, napCopy(h));
1061      napDelete (&(h->ne));
1062      if (q == NULL)
1063      {
1064        napDelete(&h);
1065        return q;
1066      }
1067      /* try to reduce further */
1068      i = 0;
1069    }
1070    else
1071    {
1072      i++;
1073      if (i >= naI->anz)
1074      {
1075        napDelete(&h);
1076        return q;
1077      }
1078    }
1079  }
1080}
1081
1082
1083/*2
1084* meins  (for reduction in algebraic extension)
1085* reduces the tail of Poly q
1086* needs q != NULL
1087* changes q and returns it
1088*/
1089alg napTailred (alg q)
1090{
1091  alg h;
1092
1093  h = q->ne;
1094  while (h != NULL)
1095  {
1096     h = napRedp (h);
1097     if (h == NULL)
1098        return q;
1099     h = h->ne;
1100  }
1101  return q;
1102}
1103
1104
1105/*================ procedure for rational functions: naXXXX =================*/
1106
1107/*2
1108*  z:= i
1109*/
1110number naInit(int i)
1111{
1112  if (i!=0)
1113  {
1114    lnumber l = (lnumber)Alloc(sizeof(rnumber));
1115    l->z = napInit(i);
1116    if (l->z==NULL)
1117    {
1118      Free((ADDRESS)l, sizeof(rnumber));
1119      return NULL;
1120    }
1121    l->s = 2;
1122    l->n = NULL;
1123    return (number)l;
1124  }
1125  /*else*/
1126  return NULL;
1127}
1128
1129number  naPar(int i)
1130{
1131  lnumber l = (lnumber)Alloc(sizeof(rnumber));
1132  l->s = 2;
1133  l->z = napInit(1);
1134  l->z->e[i-1]=1;
1135  l->n = NULL;
1136  return (number)l;
1137}
1138
1139int     naParDeg(number n)     /* i := deg(n) */
1140{
1141  lnumber l = (lnumber)n;
1142  if (l==NULL) return -1;
1143  return napDeg(l->z);
1144}
1145
1146/*2
1147* convert a number to int (if possible)
1148*/
1149int naInt(number &n)
1150{
1151  lnumber l=(lnumber)n;
1152  if ((l!=NULL)&&(l->n==NULL)&&(napDeg(l->z)==0))
1153  {
1154    return nacInt(l->z->ko);
1155  }
1156  return 0;
1157}
1158
1159/*2
1160*  deletes p
1161*/
1162#ifdef LDEBUG
1163void naDBDelete(number *p,char *f, int lno)
1164#else
1165void naDelete(number *p)
1166#endif
1167{
1168  lnumber l = (lnumber) * p;
1169  if (l==NULL) return;
1170  napDelete(&(l->z));
1171  napDelete(&(l->n));
1172#if defined(MDEBUG) && defined(LDEBUG)
1173  mmDBFreeBlock((ADDRESS)l, sizeof(rnumber),f,lno);
1174#else
1175  Free((ADDRESS)l, sizeof(rnumber));
1176#endif
1177  *p = NULL;
1178}
1179
1180/*2
1181* copy p to erg
1182*/
1183number naCopy(number p)
1184{
1185  if (p==NULL) return NULL;
1186  naTest(p);
1187  lnumber erg;
1188  lnumber src = (lnumber)p;
1189  erg = (lnumber)Alloc0(sizeof(rnumber));
1190  erg->z = napCopy(src->z);
1191  erg->n = napCopy(src->n);
1192  erg->s = src->s;
1193  return (number)erg;
1194}
1195
1196/*2
1197* a dummy number: 0
1198*/
1199void naNew(number *z)
1200{
1201  *z = NULL;
1202}
1203
1204/*2
1205*  addition; lu:= la + lb
1206*/
1207number naAdd(number la, number lb)
1208{
1209  alg x, y;
1210  lnumber lu;
1211  lnumber a = (lnumber)la;
1212  lnumber b = (lnumber)lb;
1213  if (a==NULL) return naCopy(lb);
1214  if (b==NULL) return naCopy(la);
1215  mmTestP(a,sizeof(rnumber));
1216  mmTestP(b,sizeof(rnumber));
1217  lu = (lnumber)Alloc(sizeof(rnumber));
1218  if (b->n!=NULL) x = napMult(napCopy(a->z), napCopy(b->n));
1219  else            x = napCopy(a->z);
1220  if (a->n!=NULL) y = napMult(napCopy(b->z), napCopy(a->n));
1221  else            y = napCopy(b->z);
1222  lu->z = napAdd(x, y);
1223  if (lu->z==NULL)
1224  {
1225    Free((ADDRESS)lu, sizeof(rnumber));
1226    return (number)NULL;
1227  }
1228  if (a->n!=NULL)
1229  {
1230    if (b->n!=NULL) x = napMult(napCopy(a->n), napCopy(b->n));
1231    else            x = napCopy(a->n);
1232  }
1233  else
1234  {
1235    if (b->n!=NULL) x = napCopy(b->n);
1236    else            x = NULL;
1237  }
1238  if ((x!=NULL) && (napDeg(x)==0) && nacIsOne(x->ko))
1239    napDelete(&x);
1240  lu->n = x;
1241  lu->s = 0;
1242  naTest((number)lu);
1243  return (number)lu;
1244}
1245
1246/*2
1247*  subtraction; r:= la - lb
1248*/
1249number naSub(number la, number lb)
1250{
1251  lnumber lu;
1252
1253  if (lb==NULL) return naCopy(la);
1254  if (la==NULL)
1255  {
1256    //if (lb!=NULL)
1257    //{
1258      lu = (lnumber)naCopy(lb);
1259      lu->z = napNeg(lu->z);
1260      return (number)lu;
1261    //}
1262    //else
1263    //  return NULL;
1264  }
1265
1266  alg x, y;
1267  lnumber a = (lnumber)la;
1268  lnumber b = (lnumber)lb;
1269
1270  mmTestP(a,sizeof(rnumber));
1271  mmTestP(b,sizeof(rnumber));
1272  lu = (lnumber)Alloc(sizeof(rnumber));
1273  if (b->n!=NULL) x = napMult(napCopy(a->z), napCopy(b->n));
1274  else            x = napCopy(a->z);
1275  if (a->n!=NULL) y = napMult(napCopy(b->z), napCopyNeg(a->n));
1276  else            y = napCopyNeg(b->z);
1277  lu->z = napAdd(x, y);
1278  if (lu->z==NULL)
1279  {
1280    Free((ADDRESS)lu, sizeof(rnumber));
1281    return (number)NULL;
1282  }
1283  if (a->n!=NULL)
1284  {
1285    if (b->n!=NULL) x = napMult(napCopy(a->n), napCopy(b->n));
1286    else            x = napCopy(a->n);
1287  }
1288  else
1289  {
1290    if (b->n!=NULL) x = napCopy(b->n);
1291    else            x = NULL;
1292  }
1293  if ((x!=NULL)&& (napDeg(x)==0) && nacIsOne(x->ko))
1294    napDelete(&x);
1295  lu->n = x;
1296  lu->s = 0;
1297  naTest((number)lu);
1298  return (number)lu;
1299}
1300
1301/*2
1302*  multiplication; r:= la * lb
1303*/
1304number naMult(number la, number lb)
1305{
1306  if ((la==NULL) || (lb==NULL))
1307    return NULL;
1308
1309  lnumber a = (lnumber)la;
1310  lnumber b = (lnumber)lb;
1311  lnumber lo;
1312  alg x;
1313
1314  mmTestP(a,sizeof(rnumber));
1315  mmTestP(b,sizeof(rnumber));
1316  lo = (lnumber)Alloc(sizeof(rnumber));
1317  lo->z = napMult(napCopy(a->z), napCopy(b->z));
1318  if (a->n==NULL)
1319  {
1320    if (b->n==NULL)
1321      x = NULL;
1322    else
1323      x = napCopy(b->n);
1324  }
1325  else
1326  {
1327    if (b->n==NULL)
1328      x = napCopy(a->n);
1329    else
1330      x = napMult(napCopy(b->n), napCopy(a->n));
1331  }
1332  if (naMinimalPoly!=NULL)
1333  {
1334    if (lo->z->e[0] >= naMinimalPoly->e[0])
1335      lo->z = napRemainder(lo->z, naMinimalPoly);
1336    if ((x!=NULL) && (x->e[0] >= naMinimalPoly->e[0]))
1337      x = napRemainder(x, naMinimalPoly);
1338  }
1339#if FEHLER1
1340  if (naI!=NULL)
1341  {
1342    lo->z = napRedp (lo->z);
1343    if (lo->z != NULL)
1344       lo->z = napTailred (lo->z);
1345    if (x!=NULL)
1346    {
1347      x = napRedp (x);
1348      if (x!=NULL)
1349        x = napTailred (x);
1350    }
1351  }
1352#endif
1353  if ((x!=NULL) && (napDeg(x)==0) && nacIsOne(x->ko))
1354    napDelete(&x);
1355  lo->n = x;
1356  lo->s = 0;
1357  if(lo->z==NULL)
1358  {
1359    Free((ADDRESS)lo,sizeof(rnumber));
1360    lo=NULL;
1361  }
1362  naTest((number)lo);
1363  return (number)lo;
1364}
1365
1366number naIntDiv(number la, number lb)
1367{
1368  lnumber res;
1369  lnumber a = (lnumber)la;
1370  lnumber b = (lnumber)lb;
1371  if ((a==NULL) || (a->z==NULL))
1372    return NULL;
1373  if ((b==NULL) || (b->z==NULL))
1374  {
1375    WerrorS("div. by 0");
1376    return NULL;
1377  }
1378  res = (lnumber)Alloc(sizeof(rnumber));
1379  res->z = napCopy(a->z);
1380  res->n = napCopy(b->z);
1381  res->s = 0;
1382  number nres=(number)res;
1383  naNormalize(nres);
1384
1385  //napDelete(&res->n);
1386  naTest(nres);
1387  return nres;
1388}
1389
1390/*2
1391*  division; lo:= la / lb
1392*/
1393number naDiv(number la, number lb)
1394{
1395  lnumber lo;
1396  lnumber a = (lnumber)la;
1397  lnumber b = (lnumber)lb;
1398  alg x;
1399
1400  if ((a==NULL) || (a->z==NULL))
1401    return NULL;
1402
1403  if ((b==NULL) || (b->z==NULL))
1404  {
1405    WerrorS("div. by 0");
1406    return NULL;
1407  }
1408  mmTestP(a,sizeof(rnumber));
1409  mmTestP(b,sizeof(rnumber));
1410  lo = (lnumber)Alloc(sizeof(rnumber));
1411  if (b->n!=NULL)
1412    lo->z = napMult(napCopy(a->z), napCopy(b->n));
1413  else
1414    lo->z = napCopy(a->z);
1415  if (a->n!=NULL)
1416    x = napMult(napCopy(b->z), napCopy(a->n));
1417  else
1418    x = napCopy(b->z);
1419  if (naMinimalPoly!=NULL)
1420  {
1421    if (lo->z->e[0] >= naMinimalPoly->e[0])
1422      lo->z = napRemainder(lo->z, naMinimalPoly);
1423    if (x->e[0] >= naMinimalPoly->e[0])
1424      x = napRemainder(x, naMinimalPoly);
1425  }
1426#if FEHLER1
1427  if (naI!=NULL)
1428  {
1429    lo->z = napRedp (lo->z);
1430    if (lo->z != NULL)
1431       lo->z = napTailred (lo->z);
1432    if (x!=NULL)
1433    {
1434      x = napRedp (x);
1435      if (x!=NULL)
1436        x = napTailred (x);
1437    }
1438  }
1439#endif
1440  if ((napDeg(x)==0) && nacIsOne(x->ko))
1441    napDelete(&x);
1442  lo->s = 0;
1443  lo->n = x;
1444  naTest((number)lo);
1445  return (number)lo;
1446}
1447
1448/*2
1449*  za:= - za
1450*/
1451number naNeg(number za)
1452{
1453  if (za!=NULL)
1454  {
1455    lnumber e = (lnumber)za;
1456    naTest(za);
1457    e->z = napNeg(e->z);
1458  }
1459  return za;
1460}
1461
1462/*2
1463* 1/a
1464*/
1465number naInvers(number a)
1466{
1467  lnumber lo;
1468  lnumber b = (lnumber)a;
1469  alg x;
1470
1471  if (b==NULL)
1472  {
1473    WerrorS("div. by 0");
1474    return NULL;
1475  }
1476  mmTestP(b,sizeof(rnumber));
1477  lo = (lnumber)Alloc0(sizeof(rnumber));
1478  lo->s = b->s;
1479  if (b->n!=NULL)
1480    lo->z = napCopy(b->n);
1481  else
1482    lo->z = napInit(1);
1483  x = b->z;
1484  if ((napDeg(x)!=0) || !nacIsOne(x->ko))
1485    x = napCopy(x);
1486  else
1487  {
1488    lo->n = NULL;
1489    naTest((number)lo);
1490    return (number)lo;
1491  }
1492  if (naMinimalPoly!=NULL)
1493  {
1494    x = napInvers(x, naMinimalPoly);
1495    x = napMult(x, lo->z);
1496    if (x->e[0] >= naMinimalPoly->e[0])
1497      x = napRemainder(x, naMinimalPoly);
1498    lo->z = x;
1499    lo->n = NULL;
1500    lo->s = 2;
1501    while (x!=NULL)
1502    {
1503      nacNormalize(x->ko);
1504      x = x->ne;
1505    }
1506  }
1507  else
1508    lo->n = x;
1509  naTest((number)lo);
1510  return (number)lo;
1511}
1512
1513
1514BOOLEAN naIsZero(number za)
1515{
1516  lnumber zb = (lnumber)za;
1517  naTest(za);
1518#ifdef TEST
1519  if ((zb!=NULL) && (zb->z==NULL)) WerrorS("internal zero error(2)");
1520#endif
1521  return ((zb==NULL) || (zb->z==NULL));
1522}
1523
1524
1525BOOLEAN naGreaterZero(number za)
1526{
1527  lnumber zb = (lnumber)za;
1528#ifdef TEST
1529  if ((zb!=NULL) && (zb->z==NULL)) WerrorS("internal zero error(3)");
1530#endif
1531  naTest(za);
1532  return ((zb!=NULL) && (zb->z!=NULL));
1533}
1534
1535
1536/*2
1537* a = b ?
1538*/
1539BOOLEAN naEqual (number a, number b)
1540{
1541  number h;
1542  BOOLEAN bo;
1543  h = naSub(a, b);
1544  bo = naIsZero(h);
1545  naDelete(&h);
1546  return bo;
1547}
1548
1549
1550BOOLEAN naGreater (number a, number b)
1551{
1552  if (naIsZero(a))
1553    return FALSE;
1554  if (naIsZero(b))
1555    return TRUE; /* a!= 0)*/
1556  return napDeg(((lnumber)a)->n)>napDeg(((lnumber)b)->n);
1557}
1558
1559/*2
1560* reads a number
1561*/
1562char  *naRead(char *s, number *p)
1563{
1564  alg x;
1565  lnumber a;
1566  s = napRead(s, &x);
1567  if (x==NULL)
1568  {
1569    *p = NULL;
1570    return s;
1571  }
1572  *p = (number)Alloc0(sizeof(rnumber));
1573  a = (lnumber)*p;
1574  if ((naMinimalPoly!=NULL) && (x->e[0] >= naMinimalPoly->e[0]))
1575    a->z = napRemainder(x, naMinimalPoly);
1576#if FEHLER3
1577  else if (naI!=NULL)
1578  {
1579    a->z = napRedp(x);
1580    if (a->z != NULL)
1581      a->z = napTailred (a->z);
1582  }
1583#endif
1584  else
1585    a->z = x;
1586  if(a->z==NULL)
1587  {
1588    Free((ADDRESS)*p,sizeof(rnumber));
1589    *p=NULL;
1590  }
1591  else
1592  {
1593    a->n = NULL;
1594    a->s = 0;
1595    naTest(*p);
1596  }
1597  return s;
1598}
1599
1600/*2
1601* tries to convert a number to a name
1602*/
1603char * naName(number n)
1604{
1605  lnumber ph = (lnumber)n;
1606  if ((ph==NULL)||(ph->z==NULL))
1607    return NULL;
1608  int i;
1609  char *s=(char *)AllocL(4* naNumbOfPar);
1610  char *t=(char *)Alloc(8);
1611  s[0]='\0';
1612  for (i = 0; i <= naNumbOfPar - 1; i++)
1613  {
1614    if (ph->z->e[i] > 0)
1615    {
1616      if (ph->z->e[i] >1)
1617      {
1618        sprintf(t,"%s%d",naParNames[i],ph->z->e[i]);
1619        strcat(s,t);
1620      }
1621      else
1622      {
1623        strcat(s,naParNames[i]);
1624      }
1625    }
1626  }
1627  Free((ADDRESS)t,8);
1628  if (s[0]=='\0')
1629  {
1630    FreeL((ADDRESS)s);
1631    return NULL;
1632  }
1633  return s;
1634}
1635
1636/*2
1637*  writes a number
1638*/
1639void naWrite(number &phn)
1640{
1641  lnumber ph = (lnumber)phn;
1642  if ((ph==NULL)||(ph->z==NULL))
1643    StringAppendS("0");
1644  else
1645  {
1646    //phn->s = 0;
1647    naNormalize(phn);
1648    napWrite(ph->z);
1649    if (ph->n!=NULL)
1650    {
1651      StringAppendS("/");
1652      napWrite(ph->n);
1653    }
1654  }
1655}
1656
1657/*2
1658* za == 1 ?
1659*/
1660BOOLEAN naIsOne(number za)
1661{
1662  lnumber a = (lnumber)za;
1663  alg x, y;
1664  number t;
1665  if (a==NULL) return FALSE;
1666  mmTestP(a,sizeof(rnumber));
1667#ifdef TEST
1668  if (a->z==NULL) WerrorS("internal zero error(4)");
1669#endif
1670  if (a->n==NULL)
1671  {
1672    if (napDeg(a->z)==0) return nacIsOne(a->z->ko);
1673    else                 return FALSE;
1674  }
1675  x = a->z;
1676  y = a->n;
1677  do
1678  {
1679    if (napComp(x, y))
1680      return FALSE;
1681    else
1682    {
1683      t = nacSub(x->ko, y->ko);
1684      if (!nacIsZero(t))
1685      {
1686        nacDelete(&t);
1687        return FALSE;
1688      }
1689      else
1690        nacDelete(&t);
1691    }
1692    x = x->ne;
1693    y = y->ne;
1694  }
1695  while ((x!=NULL) && (y!=NULL));
1696  if ((x!=NULL) || (y!=NULL)) return FALSE;
1697  napDelete(&a->z);
1698  napDelete(&a->n);
1699  a->z = napInit(1);
1700  a->n = NULL;
1701  a->s = 2;
1702  return TRUE;
1703}
1704
1705/*2
1706* za == -1 ?
1707*/
1708BOOLEAN naIsMOne(number za)
1709{
1710  lnumber a = (lnumber)za;
1711  alg x, y;
1712  number t;
1713  if (a==NULL) return FALSE;
1714  mmTestP(a,sizeof(rnumber));
1715#ifdef TEST
1716  if (a->z==NULL) WerrorS("internal zero error(4)");
1717#endif
1718  if (a->n==NULL)
1719  {
1720    if (napDeg(a->z)==0) return nacIsMOne(a->z->ko);
1721    else                 return FALSE;
1722  }
1723  x = a->z;
1724  y = a->n;
1725  do
1726  {
1727    t = nacAdd(x->ko, y->ko);
1728    if (!nacIsZero(t))
1729    {
1730      nacDelete(&t);
1731      return FALSE;
1732    }
1733    else
1734      nacDelete(&t);
1735    x = x->ne;
1736    y = y->ne;
1737  }
1738  while ((x!=NULL) && (y!=NULL));
1739  if ((x!=NULL) || (y!=NULL)) return FALSE;
1740  napDelete(&a->z);
1741  napDelete(&a->n);
1742  a->z = napInit(-1);
1743  a->n = NULL;
1744  a->s = 2;
1745  return TRUE;
1746}
1747
1748/*2
1749* returns the i-th power of p (i>=0)
1750*/
1751void naPower(number p, int i, number *rc)
1752{
1753  number x;
1754  *rc = naInit(1);
1755  for (; i > 0; i--)
1756  {
1757    x = naMult(*rc, p);
1758    naDelete(rc);
1759    *rc = x;
1760  }
1761}
1762
1763/*2
1764* result =gcd(a,b)
1765*/
1766number naGcd(number a, number b)
1767{
1768  lnumber x, y;
1769  lnumber result = (lnumber)Alloc0(sizeof(rnumber));
1770
1771  x = (lnumber)a;
1772  y = (lnumber)b;
1773  if (naNumbOfPar == 1)
1774  {
1775    if (naMinimalPoly!=NULL)
1776    {
1777      if (x->z->ne!=NULL)
1778        result->z = napCopy(x->z);
1779      else
1780        result->z = napGcd0(x->z, y->z);
1781    }
1782    else
1783      result->z = napGcd(x->z, y->z);
1784  }
1785  else
1786    result->z = napGcd0(x->z, y->z);
1787  naTest((number)result);
1788  return (number)result;
1789}
1790
1791/*2
1792* naNumbOfPar = 1:
1793* clears denominator         algebraic case;
1794* tries to simplify ratio    transcendental case;
1795*
1796* cancels monomials
1797* occuring in denominator
1798* and enumerator  ?          naNumbOfPar != 1;
1799*/
1800void naNormalize(number &pp)
1801{
1802
1803  int  i,m;
1804  number  h, h1;
1805  alg x, y, r, x1, y1;
1806  lnumber p = (lnumber)pp;
1807
1808  if ((p==NULL) /*|| (p->s==2)*/)
1809    return;
1810  BOOLEAN normal_test=(p->s==2);
1811  x = p->z;
1812  y = p->n;
1813  p->s = 2;
1814  if ((y!=NULL) && (naMinimalPoly!=NULL))
1815  {
1816    y = napInvers(y, naMinimalPoly);
1817    x = napMult(x, y);
1818    if (x->e[0] >= naMinimalPoly->e[0])
1819    x = napRemainder(x, naMinimalPoly);
1820    p->z = x;
1821    p->n = y = NULL;
1822  }
1823  if (y==NULL)
1824  {
1825    if (naIsChar0)
1826    {
1827      r = x;
1828      while (r!=NULL)
1829      {
1830        nacNormalize(r->ko);
1831        r = r->ne;
1832      }
1833    }
1834    return;
1835  }
1836  // p->n !=NULL:
1837#if FEHLER1
1838  if (naMinimalPoly == NULL)
1839  {
1840    alg xx=x;
1841    alg yy=y;
1842    for (i=(naNumbOfPar-1); i>=0; i--)
1843    {
1844      m = napExpi(i, yy, xx);
1845      if (m != 0)          // in this case xx!=NULL!=yy
1846      {
1847        while (xx != NULL)
1848        {
1849           xx->e[i] -= m;
1850           xx = xx->ne;
1851        }
1852        while (yy != NULL)
1853        {
1854          yy->e[i] -= m;
1855          yy = yy->ne;
1856        }
1857      }
1858    }
1859  }
1860#endif
1861  if (naIsChar0)
1862  {
1863    h1 = nacInvers(y->ko);
1864    nacNormalize(h1);
1865    napMultN(x, h1);
1866    napMultN(y->ne, h1);
1867    nacDelete(&y->ko);
1868    y->ko = nacInit(1);
1869    nacDelete(&h1);
1870  }
1871  if (napDeg(y)==0)
1872  {
1873    if (nacIsOne(y->ko))
1874    {
1875      napDelete1(&y);
1876      p->n = NULL;
1877      return;
1878    }
1879    h = nacInit(1);
1880    h1 = nacDiv(h, y->ko);
1881    nacDelete(&h);
1882    nacNormalize(h1);
1883    napMultN(x, h1);
1884    nacDelete(&h1);
1885    napDelete1(&y);
1886    p->n = NULL;
1887    return;
1888  }
1889  if (naNumbOfPar != 1)
1890    return;
1891  i = napExp(x, y);
1892  if (i!=0)
1893  {
1894    r = x;
1895    while (r!=NULL)
1896    {
1897      r->e[0] -= i;
1898      r = r->ne;
1899    }
1900    r = y;
1901    while (r!=NULL)
1902    {
1903      r->e[0] -= i;
1904      r = r->ne;
1905    }
1906  }
1907  if (y->ne==NULL)
1908  {
1909    if (nacIsOne(y->ko))
1910    {
1911      if (y->e[0]==0)
1912      {
1913        napDelete1(&y);
1914        p->n = NULL;
1915      }
1916      return;
1917    }
1918    h = nacInit(1);
1919    h1 = nacDiv(h, y->ko);
1920    if (y->e[0]!=0)
1921    {
1922      nacDelete(&y->ko);
1923      y->ko = h;
1924    }
1925    else
1926    {
1927      napDelete1(&y);
1928      nacDelete(&h);
1929      p->n = NULL;
1930    }
1931    napMultN(x, h1);
1932    nacDelete(&h1);
1933    return;
1934  }
1935  if (x->e[0] >= y->e[0])
1936  {
1937    x1 = napCopy(x);
1938    y1 = napCopy(y);
1939  }
1940  else
1941  {
1942    x1 = napCopy(y);
1943    y1 = napCopy(x);
1944  }
1945  loop
1946  {
1947    r = napRemainder(x1, y1);
1948    if ((r==NULL) || (r->ne==NULL)) break;
1949    x1 = y1;
1950    y1 = r;
1951  }
1952  if (r!=NULL)
1953  {
1954    napDelete(&r);
1955    napDelete(&y1);
1956  }
1957  else
1958  {
1959    napDivMod(x, y1, &(p->z), &r);
1960    napDivMod(y, y1, &(p->n), &r);
1961    napDelete(&y1);
1962  }
1963  x = p->z;
1964  y = p->n;
1965  if (y->ne==NULL)
1966  {
1967    if (nacIsOne(y->ko))
1968    {
1969      if (y->e[0]==0)
1970      {
1971        napDelete1(&y);
1972        p->n = NULL;
1973      }
1974      return;
1975    }
1976    h = nacInit(1);
1977    h1 = nacDiv(h, y->ko);
1978    if (y->e[0]!=0)
1979    {
1980      nacDelete(&(y->ko));
1981      y->ko = h;
1982    }
1983    else
1984    {
1985      napDelete1(&y);
1986      nacDelete(&h);
1987      p->n = NULL;
1988    }
1989    napMultN(x, h1);
1990    nacDelete(&h1);
1991  }
1992}
1993
1994/*2
1995* returns in result->n 1
1996* and in     result->z the lcm(a->z,b->n)
1997*/
1998number naLcm(number la, number lb)
1999{
2000  lnumber result;
2001  lnumber a = (lnumber)la;
2002  lnumber b = (lnumber)lb;
2003  alg x;
2004  number t, bt, r;
2005  result = (lnumber)Alloc0(sizeof(rnumber));
2006  if (((naMinimalPoly==NULL) && (naI==NULL)) || !naIsChar0)
2007  {
2008    result->z = napInit(1);
2009    return (number)result;
2010  }
2011  naNormalize(lb);
2012  x = napCopy(a->z);
2013  t = napLcm(b->z);
2014  if (!nacIsOne(t))
2015  {
2016    bt = nacGcd(t, x->ko);
2017    r = nacMult(t, x->ko);
2018    nacDelete(&(x->ko));
2019    x->ko = nacDiv(r, bt);
2020    nacNormalize(x->ko);
2021    nacDelete(&bt);
2022    nacDelete(&r);
2023  }
2024  nacDelete(&t);
2025  result->z = x;
2026  naTest((number)result);
2027  return (number)result;
2028}
2029
2030/*2
2031* input: a set of constant polynomials
2032* sets the global variable naI
2033*/
2034void naSetIdeal(ideal I)
2035{
2036  int i;
2037
2038  if (idIs0(I))
2039  {
2040    for (i=naI->anz-1; i>=0; i--)
2041      napDelete(&naI->liste[i]);
2042    Free((ADDRESS)naI,sizeof(*naI));
2043    naI=NULL;
2044  }
2045  else
2046  {
2047    lnumber h;
2048    number a;
2049    alg x;
2050
2051    naI=(naIdeal)Alloc(sizeof(*naI));
2052    naI->anz=IDELEMS(I);
2053    naI->liste=(alg*)Alloc(naI->anz*sizeof(alg));
2054    for (i=IDELEMS(I)-1; i>=0; i--)
2055    {
2056      h=(lnumber)pGetCoeff(I->m[i]);
2057      /* We only need the enumerator of h, as we expect it to be a polynomial */
2058      naI->liste[i]=napCopy(h->z);
2059      /* If it isn't normalized (lc = 1) do this */
2060      if (!nacIsOne(naI->liste[i]->ko))
2061      {
2062        x=naI->liste[i];
2063        a=nacCopy(x->ko);
2064        a=nacDiv(nacInit(1),a);
2065        napMultN(x,a);
2066        nacDelete(&a);
2067      }
2068    }
2069  }
2070}
2071
2072/*2
2073* map Z/p -> Q(a)
2074*/
2075number naMapP0(number c)
2076{
2077  if (npIsZero(c)) return NULL;
2078  lnumber l=(lnumber)Alloc(sizeof(rnumber));
2079  l->s=2;
2080  l->z=(alg)Alloc0(RECA_SIZE + naNumbOfPar * sizeof(int));
2081  int i=(int)c;
2082  if (i>(naPrimeM>>2)) i-=naPrimeM;
2083  l->z->ko=nlInit(i);
2084  l->n=NULL;
2085  return (number)l;
2086}
2087
2088/*2
2089* map Q -> Q(a)
2090*/
2091number naMap00(number c)
2092{
2093  if (nlIsZero(c)) return NULL;
2094  lnumber l=(lnumber)Alloc(sizeof(rnumber));
2095  l->s=0;
2096  l->z=(alg)Alloc0(RECA_SIZE + naNumbOfPar * sizeof(int));
2097  l->z->ko=nlCopy(c);
2098  l->n=NULL;
2099  return (number)l;
2100}
2101
2102/*2
2103* map Z/p -> Z/p(a)
2104*/
2105number naMapPP(number c)
2106{
2107  if (npIsZero(c)) return NULL;
2108  lnumber l=(lnumber)Alloc(sizeof(rnumber));
2109  l->s=2;
2110  l->z=(alg)Alloc0(RECA_SIZE + naNumbOfPar * sizeof(int));
2111  l->z->ko=c; /* omit npCopy, because npCopy is a no-op */
2112  l->n=NULL;
2113  return (number)l;
2114}
2115
2116/*2
2117* map Z/p' -> Z/p(a)
2118*/
2119number naMapPP1(number c)
2120{
2121  if (npIsZero(c)) return NULL;
2122  int i=(int)c;
2123  if (i>naPrimeM) i-=naPrimeM;
2124  number n=npInit(i);
2125  if (npIsZero(n)) return NULL;
2126  lnumber l=(lnumber)Alloc(sizeof(rnumber));
2127  l->s=2;
2128  l->z=(alg)Alloc0(RECA_SIZE + naNumbOfPar * sizeof(int));
2129  l->z->ko=n;
2130  l->n=NULL;
2131  return (number)l;
2132}
2133
2134/*2
2135* map Q -> Z/p(a)
2136*/
2137number naMap0P(number c)
2138{
2139  if (nlIsZero(c)) return NULL;
2140  number n=npInit(nlInt(c));
2141  if (npIsZero(n)) return NULL;
2142  lnumber l=(lnumber)Alloc(sizeof(rnumber));
2143  l->s=2;
2144  l->z=(alg)Alloc0(RECA_SIZE + naNumbOfPar * sizeof(int));
2145  l->z->ko=n;
2146  l->n=NULL;
2147  return (number)l;
2148}
2149
2150static number (*nacMap)(number);
2151static int naParsToCopy;
2152static alg napMap(alg p)
2153{
2154  alg w, a;
2155
2156  if (p==NULL) return NULL;
2157  a = w = (alg)Alloc0(RECA_SIZE + naNumbOfPar * sizeof(int));
2158  memcpy(a->e, p->e, naParsToCopy * sizeof(int));
2159  w->ko = nacMap(p->ko);
2160  loop
2161  {
2162    p=p->ne;
2163    if (p==NULL) break;
2164    a->ne = (alg)Alloc0(RECA_SIZE + naNumbOfPar * sizeof(int));
2165    a = a->ne;
2166    memcpy(a->e, p->e, naParsToCopy * sizeof(int));
2167    a->ko = nacMap(p->ko);
2168  }
2169  a->ne = NULL;
2170  return w;
2171}
2172
2173/*2
2174* map _(a) -> _(b)
2175*/
2176number naMapQaQb(number c)
2177{
2178  if (c==NULL) return NULL;
2179  lnumber erg= (lnumber)Alloc0(sizeof(rnumber));
2180  lnumber src =(lnumber)c;
2181  erg->s=src->s;
2182  erg->z=napMap(src->z);
2183  erg->n=napMap(src->n);
2184  return (number)erg;
2185}
2186
2187BOOLEAN naSetMap(int c, char ** par, int nop, number minpol)
2188{
2189  if (currRing->ch==1) /* -> Q(a) */
2190  {
2191    if (c == 0)
2192    {
2193      nMap = naMap00;   /*Q -> Q(a)*/
2194      return TRUE;
2195    }
2196    if (c>1)
2197    {
2198      if (par==NULL)
2199      {
2200        naPrimeM = c;
2201        nMap = naMapP0;  /* Z/p -> Q(a)*/
2202        return TRUE;
2203      }
2204      else
2205      {
2206        return FALSE;   /* GF(q) -> Q(a) */
2207      }
2208    }
2209    if (c<0)
2210    {
2211      return FALSE;   /* Z/p(a) -> Q(b)*/
2212    }
2213    if (c==1)
2214    {
2215      int i;
2216      naParsToCopy=0;
2217      for(i=0;i<nop;i++)
2218      {
2219        if ((i>=currRing->P)
2220        ||(strcmp(par[i],currRing->parameter[i])!=0))
2221           return FALSE;
2222        naParsToCopy++;
2223      }
2224      nacMap=nacCopy;
2225      nMap=naMapQaQb;
2226      return TRUE;   /* Q(a) -> Q(a) */
2227    }
2228  }
2229  /*-----------------------------------------------------*/
2230  if (currRing->ch<0) /* -> Z/p(a) */
2231  {
2232    if (c == 0)
2233    {
2234      nMap = naMap0P;   /*Q -> Z/p(a)*/
2235      return TRUE;
2236    }
2237    if (c>1)
2238    {
2239      if (par==NULL)
2240      {
2241        if (c==npPrimeM)
2242        {
2243          nMap = naMapPP;  /* Z/p -> Z/p(a)*/
2244          return TRUE;
2245        }
2246        else
2247        {
2248          naPrimeM = c;
2249          nMap = naMapPP1;  /* Z/p' -> Z/p(a)*/
2250          return TRUE;
2251        }
2252      }
2253      else
2254      {
2255        return FALSE;   /* GF(q) ->Z/p(a) */
2256      }
2257    }
2258    if (c<0)
2259    {
2260      if (c==currRing->ch)
2261      {
2262        nacMap=nacCopy;
2263      }
2264      else
2265      {
2266        npMapPrime=-c;
2267        nacMap = npMapP;
2268      }
2269      int i;
2270      naParsToCopy=0;
2271      for(i=0;i<nop;i++)
2272      {
2273        if ((i>=currRing->P)
2274        ||(strcmp(par[i],currRing->parameter[i])!=0))
2275           return FALSE;
2276        naParsToCopy++;
2277      }
2278      nMap=naMapQaQb;
2279      return TRUE;   /* Z/p(a),Z/p'(a) -> Z/p(b)*/
2280    }
2281    if (c==1)
2282    {
2283      return FALSE;   /* Q(a) -> Z/p(b) */
2284    }
2285  }
2286  return FALSE;      /* default */
2287}
2288
2289/*2
2290* convert a alg number into a poly
2291*/
2292poly naPermNumber(number z, int * par_perm, int P)
2293{
2294  if (z==NULL) return NULL;
2295  poly res=NULL;
2296  poly p;
2297  napoly za=((lnumber)z)->z;
2298  do
2299  {
2300    p=pNew();
2301    pNext(p)=NULL;
2302    nNew(&pGetCoeff(p));
2303    int i;
2304    for(i=pVariables;i>=0;i--)
2305       pGetExp(p,i)=0;
2306    napoly pa=NULL;
2307    if (currRing->parameter!=NULL)
2308    {
2309      pGetCoeff(p)=(number)Alloc0(sizeof(rnumber));
2310      ((lnumber)pGetCoeff(p))->s=2;
2311      ((lnumber)pGetCoeff(p))->z=napInitz(nacCopy(napGetCoeff(za)));
2312      pa=((lnumber)pGetCoeff(p))->z;
2313    }
2314    else
2315    {
2316      pGetCoeff(p)=nCopy(napGetCoeff(za));
2317    }
2318    for(i=0;i<P;i++)
2319    {
2320      if(za->e[i]!=0)
2321      {
2322        if(par_perm[i]>0)
2323          pGetExp(p,par_perm[i])=za->e[i];
2324        else if((par_perm[i]<0)&&(pa!=NULL))
2325          pa->e[-par_perm[i]-1]=za->e[i];
2326        else
2327        {
2328          pDelete(&p);
2329          break;
2330        }
2331      }
2332    }
2333    if (p!=NULL)
2334    {
2335      pSetm(p);
2336      pTest(p);
2337      res=pAdd(res,p);
2338    }
2339    za=za->ne;
2340  }
2341  while (za!=NULL);
2342  pTest(res);
2343  return res;
2344}
2345
2346#ifdef LDEBUG
2347BOOLEAN naDBTest(number a, char *f,int l)
2348{
2349  lnumber x=(lnumber)a;
2350  if (x == NULL)
2351    return TRUE;
2352#ifdef MDEBUG
2353  if (!mmDBTestBlock(a,sizeof(rnumber),f,l))
2354    return FALSE;
2355#endif
2356  alg p = x->z;
2357  if (p==NULL)
2358  {
2359    Print("0/* in %s:%d\n",f,l);
2360    return FALSE;
2361  }
2362  while(p!=NULL)
2363  {
2364    if((naMinimalPoly!=NULL)&&(p->e[0]>naMinimalPoly->e[0])
2365    &&(p!=naMinimalPoly))
2366    {
2367      Print("deg>minpoly in %s:%d\n",f,l);
2368      return FALSE;
2369    }
2370    //if (naIsChar0 && (((int)p->ko &3) == 0) && (p->ko->s==0) && (x->s==2))
2371    //{
2372    //  Print("normalized with non-normal coeffs in %s:%d\n",f,l);
2373    //  return FALSE;
2374    //}
2375    if (naIsChar0 && nlDBTest(p->ko,f,l))
2376      return FALSE;
2377#ifdef MDEBUG
2378    if (!mmDBTestBlock(p,RECA_SIZE+naNumbOfPar*sizeof(int),f,l))
2379      return FALSE;
2380#endif
2381    p = p->ne;
2382  }
2383  p = x->n;
2384  while(p!=NULL)
2385  {
2386    if (naIsChar0 && nlDBTest(p->ko,f,l))
2387      return FALSE;
2388#ifdef MDEBUG
2389    if (!mmDBTestBlock(p,RECA_SIZE+naNumbOfPar*sizeof(int),f,l))
2390      return FALSE;
2391#endif
2392    p = p->ne;
2393  }
2394  return TRUE;
2395}
2396#endif
2397
Note: See TracBrowser for help on using the repository browser.