File: polys.cc

package info (click to toggle)
singular 1%3A4.1.1-p2%2Bds-3
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 35,860 kB
  • sloc: cpp: 288,280; ansic: 17,387; lisp: 4,242; yacc: 1,654; python: 1,608; makefile: 1,424; lex: 1,387; perl: 632; sh: 567; xml: 182
file content (155 lines) | stat: -rw-r--r-- 4,277 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
#include "kernel/mod2.h"

#include "omalloc/omalloc.h"
#include "misc/options.h"

#include "polys.h"
#include "kernel/ideals.h"
#include "kernel/ideals.h"
#include "polys/clapsing.h"

/// Widely used global variable which specifies the current polynomial ring for Singular interpreter and legacy implementatins.
/// @Note: one should avoid using it in newer designs, for example due to possible problems in parallelization with threads.
ring  currRing = NULL;

void rChangeCurrRing(ring r)
{
  //------------ set global ring vars --------------------------------
  currRing = r;
  if( r != NULL )
  {
    rTest(r);
    //------------ global variables related to coefficients ------------
    assume( r->cf!= NULL );
    nSetChar(r->cf);
    //------------ global variables related to polys
    p_SetGlobals(r); // also setting TEST_RINGDEP_OPTS
    //------------ global variables related to factory -----------------
  }
}

poly p_Divide(poly p, poly q, const ring r)
{
  assume(q!=NULL);
  if (q==NULL)
  {
    WerrorS("div. by 0");
    return NULL;
  }
  if (p==NULL)
  {
    p_Delete(&q,r);
    return NULL;
  }
  if (pNext(q)!=NULL)
  { /* This means that q != 0 consists of at least two terms*/
    if(p_GetComp(p,r)==0)
    {
      if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
      &&(!rField_is_Ring(r)))
      {
        poly res=singclap_pdivide(p, q, r);
        p_Delete(&p,r);
        p_Delete(&q,r);
        return res;
      }
      else
      {
        ideal vi=idInit(1,1); vi->m[0]=q;
        ideal ui=idInit(1,1); ui->m[0]=p;
        ideal R; matrix U;
        ring save_ring=currRing;
        if (r!=currRing) rChangeCurrRing(r);
        int save_opt;
        SI_SAVE_OPT1(save_opt);
        si_opt_1 &= ~(Sy_bit(OPT_PROT));
        ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
        SI_RESTORE_OPT1(save_opt);
        if (r!=save_ring) rChangeCurrRing(save_ring);
        matrix T = id_Module2formatedMatrix(m,1,1,r);
        p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
        id_Delete((ideal *)&T,r);
        id_Delete((ideal *)&U,r);
        id_Delete(&R,r);
        //vi->m[0]=NULL; ui->m[0]=NULL;
        id_Delete(&vi,r);
        id_Delete(&ui,r);
        return p;
      }
    }
    else
    {
      int comps=p_MaxComp(p,r);
      ideal I=idInit(comps,1);
      poly h;
      int i;
      // conversion to a list of polys:
      while (p!=NULL)
      {
        i=p_GetComp(p,r)-1;
        h=pNext(p);
        pNext(p)=NULL;
        p_SetComp(p,0,r);
        I->m[i]=p_Add_q(I->m[i],p,r);
        p=h;
      }
      // division and conversion to vector:
      h=NULL;
      p=NULL;
      for(i=comps-1;i>=0;i--)
      {
        if (I->m[i]!=NULL)
        {
          if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
          &&(!rField_is_Ring(r)))
            h=singclap_pdivide(I->m[i],q,r);
          else
          {
            ideal vi=idInit(1,1); vi->m[0]=q;
            ideal ui=idInit(1,1); ui->m[0]=I->m[i];
            ideal R; matrix U;
            ring save_ring=currRing;
            if (r!=currRing) rChangeCurrRing(r);
            int save_opt;
            SI_SAVE_OPT1(save_opt);
            si_opt_1 &= ~(Sy_bit(OPT_PROT));
            ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
            SI_RESTORE_OPT1(save_opt);
            if (r!=save_ring) rChangeCurrRing(save_ring);
            matrix T = id_Module2formatedMatrix(m,1,1,r);
            h=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
            id_Delete((ideal*)&T,r);
            id_Delete((ideal*)&U,r);
            id_Delete(&R,r);
            vi->m[0]=NULL; ui->m[0]=NULL;
            id_Delete(&vi,r);
            id_Delete(&ui,r);
          }
          p_SetCompP(h,i+1,r);
          p=p_Add_q(p,h,r);
        }
      }
      id_Delete(&I,r);
      p_Delete(&q,r);
      return p;
    }
  }
  else
  { /* This means that q != 0 consists of just one term,
       or that currRing is over a coefficient ring. */
#ifdef HAVE_RINGS
    if (!rField_is_Domain(currRing))
    {
      WerrorS("division only defined over coefficient domains");
      return NULL;
    }
    if (pNext(q)!=NULL)
    {
      WerrorS("division over a coefficient domain only implemented for terms");
      return NULL;
    }
#endif
    return p_DivideM(p,q,r);
  }
  return FALSE;
}