source: git/Singular/longalg.cc @ 49f089

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