source: git/kernel/longalg.cc @ b506079

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