source: git/kernel/longalg.cc @ 936551

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