source: git/Singular/longalg.cc @ 35aab3

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