source: git/Singular/longalg.cc @ 18dd47

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