source: git/kernel/longalg.cc @ 3dbee61

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