source: git/Singular/longalg.cc @ 25ea20d

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