[BACK]Return to isolv.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Diff for /OpenXM_contrib2/asir2000/builtin/isolv.c between version 1.5 and 1.8

version 1.5, 2005/07/14 22:46:03 version 1.8, 2019/06/04 07:11:22
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2000/builtin/isolv.c,v 1.4 2005/02/08 18:06:05 saito Exp $   * $OpenXM: OpenXM_contrib2/asir2000/builtin/isolv.c,v 1.7 2018/03/29 01:32:50 noro Exp $
  */   */
   
 #include "ca.h"  #include "ca.h"
Line 11 
Line 11 
 static void Solve(NODE, Obj *);  static void Solve(NODE, Obj *);
 static void NSolve(NODE, Obj *);  static void NSolve(NODE, Obj *);
   
   /* in builtin/vars.c */
   void Pvars();
   
   /* */
 void Solve1(P, Q, pointer *);  void Solve1(P, Q, pointer *);
 void Sturm(P, VECT *);  void Sturm(P, VECT *);
 void boundbody(P, Q *);  void boundbody(P, Q *);
 void binary(int , MAT);  void binary(int , MAT);
 void separate(Q, Q, VECT, Q, Q, int, int, MAT, int *);  void separate(Q, Q, VECT, Q, Q, int, int, MAT, int *);
 void ueval(P, Q, Q *);  void ueval(P, Q, Q *);
   int stumq(VECT, Q);
   
   
   // in engine/bf.c
   Num tobf(Num,int);
   
 struct ftab isolv_tab[] = {  struct ftab isolv_tab[] = {
         {"solve", Solve, 2},    {"solve", Solve, 2},
         {"nsolve", NSolve, 2},    {"nsolve", NSolve, 2},
         {0,0,0},    {0,0,0},
 };  };
   
 static void  static void
Line 29  Solve(arg, rp)
Line 38  Solve(arg, rp)
 NODE arg;  NODE arg;
 Obj  *rp;  Obj  *rp;
 {  {
         pointer p, Eps;    pointer p, Eps;
         pointer    root;    pointer    root;
         V       v;    V       v;
         Q       eps;    Q       eps;
   
         p = (pointer)ARG0(arg);    p = (pointer)ARG0(arg);
         if ( !p ) {    if ( !p ) {
                 *rp = 0;      *rp = 0;
                 return;      return;
         }    }
         Eps = (pointer)ARG1(arg);    Eps = (pointer)ARG1(arg);
         asir_assert(Eps, O_N, "solve");    asir_assert(Eps, O_N, "solve");
         if ( NID(Eps) != N_Q ) {    if ( NID(Eps) != N_Q ) {
                 fprintf(stderr,"solve arg 2 is required a rational number");      fprintf(stderr,"solve arg 2 is required a rational number");
                 error(" : invalid argument");      error(" : invalid argument");
                 return;      return;
         }    }
         DUPQ((Q)Eps, eps);    DUPQ((Q)Eps, eps);
         SGN(eps) = 1;    SGN(eps) = 1;
         switch (OID(p)) {    switch (OID(p)) {
                 case O_N:      case O_N:
                         *rp = 0;        *rp = 0;
                         break;        break;
                 case O_P:      case O_P:
                         Pvars(arg, &root);        Pvars(arg, &root);
                         if (NEXT(BDY((LIST)root)) != 0) {        if (NEXT(BDY((LIST)root)) != 0) {
                                 fprintf(stderr,"solve arg 1 is univariate polynormial");          fprintf(stderr,"solve arg 1 is univariate polynormial");
                                 error(" : invalid argument");          error(" : invalid argument");
                                 break;          break;
                         }        }
                         Solve1((P)p, eps, &root);        Solve1((P)p, eps, &root);
                         *rp = (Obj)root;        *rp = (Obj)root;
                         break;        break;
                 case O_LIST:      case O_LIST:
                         fprintf(stderr,"solve,");        fprintf(stderr,"solve,");
                         error(" : Sorry, not yet implement of multivars");        error(" : Sorry, not yet implement of multivars");
                         break;        break;
                 default:      default:
                         *rp = 0;        *rp = 0;
         }    }
 }  }
   
 static void  static void
 NSolve(arg, rp)  NSolve(NODE arg, Obj *rp)
 NODE arg;  
 Obj  *rp;  
 {  {
         pointer p, Eps;    pointer  p, Eps;
         pointer root;    pointer  root;
         LIST            listp;    LIST    listp;
         V                       v;    V      v;
         Q                       eps;    Q      eps;
         NODE            n, n0, m0, m, ln0;    NODE    n, n0, m0, m, ln0;
         Num             r;    Num    r;
         Itv             iv;    Itv    iv;
         BF                      breal;    BF      breal;
   
         p = (pointer)ARG0(arg);    p = (pointer)ARG0(arg);
         if ( !p ) {    if ( !p ) {
                 *rp = 0;      *rp = 0;
                 return;      return;
         }    }
         Eps = (pointer)ARG1(arg);    Eps = (pointer)ARG1(arg);
         asir_assert(Eps, O_N, "solve");    asir_assert(Eps, O_N, "solve");
         if ( NID(Eps) != N_Q ) {    if ( NID(Eps) != N_Q ) {
                 fprintf(stderr,"solve arg 2 is required a rational number");      fprintf(stderr,"solve arg 2 is required a rational number");
                 error(" : invalid argument");      error(" : invalid argument");
                 return;      return;
         }    }
         DUPQ((Q)Eps, eps);    DUPQ((Q)Eps, eps);
         SGN(eps) = 1;    SGN(eps) = 1;
         switch (OID(p)) {    switch (OID(p)) {
                 case O_N:      case O_N:
                         *rp = 0;        *rp = 0;
                         break;        break;
                 case O_P:      case O_P:
                         Pvars(arg, &root);        Pvars(arg, &root);
                         if (NEXT(BDY((LIST)root)) != 0) {        if (NEXT(BDY((LIST)root)) != 0) {
                                 fprintf(stderr,"solve arg 1 is univariate polynormial");          fprintf(stderr,"solve arg 1 is univariate polynormial");
                                 error(" : invalid argument");          error(" : invalid argument");
                                 break;          break;
                         }        }
                         Solve1((P)p, eps, &root);        Solve1((P)p, eps, &root);
                         for (m0 = BDY((LIST)root), n0 = 0; m0; m0 = NEXT(m0)) {        for (m0 = BDY((LIST)root), n0 = 0; m0; m0 = NEXT(m0)) {
                                 m = BDY((LIST)BDY(m0));          m = BDY((LIST)BDY(m0));
                                 miditvp(BDY(m), &r);          miditvp(BDY(m), &r);
                                 ToBf(r, &breal);          //ToBf(r, &breal);
                                 NEXTNODE( n0, n );          breal = (BF)tobf(r, DEFAULTPREC);
                                 MKNODE(ln0, breal, NEXT(m));          NEXTNODE( n0, n );
                                 MKLIST((LIST)listp, ln0);          MKNODE(ln0, breal, NEXT(m));
                                 BDY(n) = (pointer)listp;          MKLIST(listp, ln0);
                         }          BDY(n) = (pointer)listp;
                         NEXT(n) = 0;        }
                         MKLIST((LIST)listp,n0);        NEXT(n) = 0;
                         *rp = (pointer)listp;        MKLIST(listp,n0);
                         break;        *rp = (pointer)listp;
                 case O_LIST:        break;
                         fprintf(stderr,"solve,");      case O_LIST:
                         error(" : Sorry, not yet implement of multivars");        fprintf(stderr,"solve,");
                         break;        error(" : Sorry, not yet implement of multivars");
                 default:        break;
                         *rp = 0;      default:
         }        *rp = 0;
     }
 }  }
   
 void  void
Line 140  P    inp;
Line 148  P    inp;
 Q    eps;  Q    eps;
 pointer *rt;  pointer *rt;
 {  {
         P    p;    P    p;
         Q    up, low, a;    Q    up, low, a;
         DCP fctp, onedeg, zerodeg;    DCP fctp, onedeg, zerodeg;
         LIST listp;    LIST listp;
         VECT sseq;    VECT sseq;
         MAT  root;    MAT  root;
         int  chvu, chvl, pad, tnumb, numb, i, j;    int  chvu, chvl, pad, tnumb, numb, i, j;
         Itv  iv;    Itv  iv;
         NODE n0, n, ln0, ln;    NODE n0, n, ln0, ln;
   
         boundbody(inp, &up);    boundbody(inp, &up);
         if (!up) {    if (!up) {
                 *rt = 0;      *rt = 0;
                 return;      return;
         }    }
         Sturm(inp, &sseq);    Sturm(inp, &sseq);
         DUPQ(up,low);    DUPQ(up,low);
         SGN(low) = -1;    SGN(low) = -1;
         chvu = stumq(sseq, up);    chvu = stumq(sseq, up);
         chvl = stumq(sseq, low);    chvl = stumq(sseq, low);
         tnumb = abs(chvu - chvl);    tnumb = abs(chvu - chvl);
         MKMAT(root, tnumb, 4);    MKMAT(root, tnumb, 4);
         pad = -1;    pad = -1;
   
         fctrp(CO,inp,&fctp);    fctrp(CO,inp,&fctp);
         for (fctp = NEXT(fctp), i = 0; fctp; fctp = NEXT(fctp)) {    for (fctp = NEXT(fctp), i = 0; fctp; fctp = NEXT(fctp)) {
                 p = COEF(fctp);      p = COEF(fctp);
                 onedeg = DC(p);      onedeg = DC(p);
                 if ( !cmpq(DEG(onedeg), ONE) ) {      if ( !cmpq(DEG(onedeg), ONE) ) {
                         pad++;        pad++;
                         if ( !NEXT(onedeg) ) {        if ( !NEXT(onedeg) ) {
                                 BDY(root)[pad][0] = 0;          BDY(root)[pad][0] = 0;
                                 BDY(root)[pad][1] = 0;          BDY(root)[pad][1] = 0;
                                 BDY(root)[pad][2] = DEG(fctp);          BDY(root)[pad][2] = DEG(fctp);
                                 BDY(root)[pad][3] = p;          BDY(root)[pad][3] = p;
                         } else {        } else {
                                 divq((Q)COEF(NEXT(onedeg)),(Q)COEF(onedeg),&a);          divq((Q)COEF(NEXT(onedeg)),(Q)COEF(onedeg),&a);
                                 BDY(root)[pad][0] = a;          BDY(root)[pad][0] = a;
                                 BDY(root)[pad][1] = BDY(root)[pad][0];          BDY(root)[pad][1] = BDY(root)[pad][0];
                                 BDY(root)[pad][2] = DEG(fctp);          BDY(root)[pad][2] = DEG(fctp);
                                 BDY(root)[pad][3] = p;          BDY(root)[pad][3] = p;
                         }        }
                         continue;        continue;
                 }      }
                 boundbody(p, &up);      boundbody(p, &up);
                 Sturm(p, &sseq);      Sturm(p, &sseq);
                 DUPQ(up,low);      DUPQ(up,low);
                 SGN(low) = -1;      SGN(low) = -1;
                 chvu = stumq(sseq, up);      chvu = stumq(sseq, up);
                 chvl = stumq(sseq, low);      chvl = stumq(sseq, low);
                 numb = abs(chvu - chvl);      numb = abs(chvu - chvl);
                 separate(DEG(fctp), eps, sseq, up, low, chvu, chvl, root, &pad);      separate(DEG(fctp), eps, sseq, up, low, chvu, chvl, root, &pad);
         }    }
         for (i = 0; i < pad; i++) {    for (i = 0; i < pad; i++) {
                 for (j = i; j <= pad; j++) {      for (j = i; j <= pad; j++) {
                         if (cmpq(BDY(root)[i][0], BDY(root)[j][0]) > 0) {        if (cmpq(BDY(root)[i][0], BDY(root)[j][0]) > 0) {
                                 a = BDY(root)[i][0];          a = BDY(root)[i][0];
                                 BDY(root)[i][0] = BDY(root)[j][0];          BDY(root)[i][0] = BDY(root)[j][0];
                                 BDY(root)[j][0] = a;          BDY(root)[j][0] = a;
                                 a = BDY(root)[i][1];          a = BDY(root)[i][1];
                                 BDY(root)[i][1] = BDY(root)[j][1];          BDY(root)[i][1] = BDY(root)[j][1];
                                 BDY(root)[j][1] = a;          BDY(root)[j][1] = a;
                                 a = BDY(root)[i][2];          a = BDY(root)[i][2];
                                 BDY(root)[i][2] = BDY(root)[j][2];          BDY(root)[i][2] = BDY(root)[j][2];
                                 BDY(root)[j][2] = a;          BDY(root)[j][2] = a;
                                 a = BDY(root)[i][3];          a = BDY(root)[i][3];
                                 BDY(root)[i][3] = BDY(root)[j][3];          BDY(root)[i][3] = BDY(root)[j][3];
                                 BDY(root)[j][3] = a;          BDY(root)[j][3] = a;
                         }        }
                 }      }
         }    }
         for (i = 0; i < pad; i++) {    for (i = 0; i < pad; i++) {
                 while(cmpq(BDY(root)[i][1], BDY(root)[i+1][0]) > 0 ) {      while(cmpq(BDY(root)[i][1], BDY(root)[i+1][0]) > 0 ) {
                         binary(i, root);        binary(i, root);
                         binary(i+1, root);        binary(i+1, root);
                         if ( cmpq(BDY(root)[i][0], BDY(root)[i+1][1]) > 0 ) {        if ( cmpq(BDY(root)[i][0], BDY(root)[i+1][1]) > 0 ) {
                                 a = BDY(root)[i][0];          a = BDY(root)[i][0];
                                 BDY(root)[i][0] = BDY(root)[i+1][0];          BDY(root)[i][0] = BDY(root)[i+1][0];
                                 BDY(root)[i+1][0] = a;          BDY(root)[i+1][0] = a;
                                 a = BDY(root)[i][1];          a = BDY(root)[i][1];
                                 BDY(root)[i][1] = BDY(root)[i+1][1];          BDY(root)[i][1] = BDY(root)[i+1][1];
                                 BDY(root)[i+1][1] = a;          BDY(root)[i+1][1] = a;
                                 a = BDY(root)[i][2];          a = BDY(root)[i][2];
                                 BDY(root)[i][2] = BDY(root)[i+1][2];          BDY(root)[i][2] = BDY(root)[i+1][2];
                                 BDY(root)[i+1][2] = a;          BDY(root)[i+1][2] = a;
                                 a = BDY(root)[i][3];          a = BDY(root)[i][3];
                                 BDY(root)[i][3] = BDY(root)[i+1][3];          BDY(root)[i][3] = BDY(root)[i+1][3];
                                 BDY(root)[i+1][3] = a;          BDY(root)[i+1][3] = a;
                                 break;          break;
                         }        }
                 }      }
         }    }
         for (i = 0, n0 = 0; i <= pad; i++) {    for (i = 0, n0 = 0; i <= pad; i++) {
                 istoitv(BDY(root)[i][0], BDY(root)[i][1], &iv);      istoitv(BDY(root)[i][0], BDY(root)[i][1], &iv);
                 NEXTNODE(n0,n);      NEXTNODE(n0,n);
                 MKNODE(ln, BDY(root)[i][2], 0); MKNODE(ln0, iv, ln);      MKNODE(ln, BDY(root)[i][2], 0); MKNODE(ln0, iv, ln);
                 MKLIST(listp, ln0);BDY(n) = (pointer)listp;      MKLIST(listp, ln0);BDY(n) = (pointer)listp;
         }    }
         NEXT(n) = 0;    NEXT(n) = 0;
         MKLIST(listp,n0);    MKLIST(listp,n0);
         *rt = (pointer)listp;    *rt = (pointer)listp;
 }  }
   
 void  void
 separate(mult, eps, sseq, up, low, upn, lown, root, padp)  separate(mult, eps, sseq, up, low, upn, lown, root, padp)
 VECT sseq;  VECT sseq;
 Q               mult, eps, up, low;  Q    mult, eps, up, low;
 int     upn, lown;  int  upn, lown;
 MAT     root;  MAT  root;
 int     *padp;  int  *padp;
 {  {
         int de, midn;    int de, midn;
         Q   mid, e;    Q   mid, e;
         P   p;    P   p;
   
         de = abs(lown - upn);    de = abs(lown - upn);
         if (de == 0) return;    if (de == 0) return;
         if (de == 1) {    if (de == 1) {
                 (*padp)++;      (*padp)++;
                 BDY(root)[*padp][0] = up;      BDY(root)[*padp][0] = up;
                 BDY(root)[*padp][1] = low;      BDY(root)[*padp][1] = low;
                 BDY(root)[*padp][3] = (P *)sseq->body[0];      BDY(root)[*padp][3] = (P *)sseq->body[0];
                 subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &e );      subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &e );
                 SGN(e) = 1;      SGN(e) = 1;
                 while (cmpq(e, eps) > 0) {      while (cmpq(e, eps) > 0) {
                         binary(*padp, root);        binary(*padp, root);
                         subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &e);        subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &e);
                         SGN(e) = 1;        SGN(e) = 1;
                 }      }
                 BDY(root)[*padp][2] = mult;      BDY(root)[*padp][2] = mult;
                 return;      return;
         }    }
         addq(up, low, &mid);    addq(up, low, &mid);
         divq(mid, TWO, &mid);    divq(mid, TWO, &mid);
         midn = stumq(sseq, mid);    midn = stumq(sseq, mid);
         separate(mult, eps, sseq, low, mid, lown, midn, root, padp);    separate(mult, eps, sseq, low, mid, lown, midn, root, padp);
         separate(mult, eps, sseq, mid, up, midn, upn, root, padp);    separate(mult, eps, sseq, mid, up, midn, upn, root, padp);
 }  }
   
 void  void
Line 284  binary(indx, root)
Line 292  binary(indx, root)
 int indx;  int indx;
 MAT root;  MAT root;
 {  {
         Q       a, b, c, d, e;    Q  a, b, c, d, e;
         P       p;    P  p;
         p = (P)BDY(root)[indx][3];    p = (P)BDY(root)[indx][3];
         addq(BDY(root)[indx][0], BDY(root)[indx][1], &c);    addq(BDY(root)[indx][0], BDY(root)[indx][1], &c);
         divq(c, TWO, &d);    divq(c, TWO, &d);
         ueval(p, BDY(root)[indx][1], &a);    ueval(p, BDY(root)[indx][1], &a);
         ueval(p, d, &b);    ueval(p, d, &b);
         if (SGN(a) == SGN(b)){    if (SGN(a) == SGN(b)){
                 BDY(root)[indx][1] = d;      BDY(root)[indx][1] = d;
         } else {    } else {
                 BDY(root)[indx][0] = d;      BDY(root)[indx][0] = d;
         }    }
 }  }
   
 void  void
Line 303  Sturm(p, ret)
Line 311  Sturm(p, ret)
 P    p;  P    p;
 VECT *ret;  VECT *ret;
 {  {
         P    g1,g2,q,r,s, *t;    P    g1,g2,q,r,s, *t;
         Q    a,b,c,d,h,l,m,x;    Q    a,b,c,d,h,l,m,x;
         V    v;    V    v;
         VECT seq;    VECT seq;
         int  i,j;    int  i,j;
   
         v = VR(p);    v = VR(p);
         t = (P *)ALLOCA((deg(v,p)+1)*sizeof(P));    t = (P *)ALLOCA((deg(v,p)+1)*sizeof(P));
         g1 = t[0] = p; diffp(CO,p,v,(P *)&a); ptozp((P)a,1,&c,&g2); t[1] = g2;    g1 = t[0] = p; diffp(CO,p,v,(P *)&a); ptozp((P)a,1,&c,&g2); t[1] = g2;
         for ( i = 1, h = ONE, x = ONE; ; ) {    for ( i = 1, h = ONE, x = ONE; ; ) {
                 if ( NUM(g2) ) break;      if ( NUM(g2) ) break;
                 subq(DEG(DC(g1)),DEG(DC(g2)),&d);      subq(DEG(DC(g1)),DEG(DC(g2)),&d);
                 l = (Q)LC(g2);      l = (Q)LC(g2);
                 if ( SGN(l) < 0 ) {      if ( SGN(l) < 0 ) {
                         chsgnq(l,&a); l = a;        chsgnq(l,&a); l = a;
                 }      }
                 addq(d,ONE,&a); pwrq(l,a,&b); mulp(CO,(P)b,g1,(P *)&a);      addq(d,ONE,&a); pwrq(l,a,&b); mulp(CO,(P)b,g1,(P *)&a);
                 divsrp(CO,(P)a,g2,&q,&r);      divsrp(CO,(P)a,g2,&q,&r);
                 if ( !r ) break;      if ( !r ) break;
                 chsgnp(r,&s); r = s; i++;      chsgnp(r,&s); r = s; i++;
                 if ( NUM(r) ) {      if ( NUM(r) ) {
                         t[i] = r; break;        t[i] = r; break;
                 }      }
                 pwrq(h,d,&m); g1 = g2;      pwrq(h,d,&m); g1 = g2;
                 mulq(m,x,&a); divsp(CO,r,(P)a,&g2); t[i] = g2;      mulq(m,x,&a); divsp(CO,r,(P)a,&g2); t[i] = g2;
                 x = (Q)LC(g1);      x = (Q)LC(g1);
                 if ( SGN(x) < 0 ) {      if ( SGN(x) < 0 ) {
                         chsgnq(x,&a); x = a;        chsgnq(x,&a); x = a;
                 }      }
                 pwrq(x,d,&a); mulq(a,h,&b); divq(b,m,&h);      pwrq(x,d,&a); mulq(a,h,&b); divq(b,m,&h);
         }    }
         MKVECT(seq,i+1);    MKVECT(seq,i+1);
         for ( j = 0; j <= i; j++ ) seq->body[j] = (pointer)t[j];    for ( j = 0; j <= i; j++ ) seq->body[j] = (pointer)t[j];
         *ret = seq;    *ret = seq;
 }  }
   
 int  int
 stumq(s, val)  stumq(VECT s, Q val)
 VECT s;  
 Q val;  
 {  {
         int len, i, j, c;    int len, i, j, c;
         P   *ss;    P   *ss;
         Q   a, a0;    Q   a, a0;
   
         len = s->len;    len = s->len;
         ss = (P *)s->body;    ss = (P *)s->body;
         for ( j = 0; j < len; j++ ){    for ( j = 0; j < len; j++ ){
                 ueval(ss[j],val,&a0);      ueval(ss[j],val,&a0);
                 if (a0) break;      if (a0) break;
         }    }
         for ( i = j++, c =0; i < len; i++) {    for ( i = j++, c =0; i < len; i++) {
                 ueval( ss[i], val, &a);      ueval( ss[i], val, &a);
                 if ( a ) {      if ( a ) {
                         if( (SGN(a) > 0 && SGN(a0) < 0) || (SGN(a) < 0 && SGN(a0) > 0) ){        if( (SGN(a) > 0 && SGN(a0) < 0) || (SGN(a) < 0 && SGN(a0) > 0) ){
                                 c++;          c++;
                                 a0 = a;          a0 = a;
                         }        }
                 }      }
         }    }
         return c;    return c;
 }  }
   
 void  void
Line 371  boundbody(p, q)
Line 377  boundbody(p, q)
 P  p;  P  p;
 Q *q;  Q *q;
 {  {
         Q               t, max, tmp;    Q    t, max, tmp;
         DCP     dc;    DCP  dc;
   
         if ( !p )    if ( !p )
                 *q = 0;      *q = 0;
         else if ( p->id == O_N )    else if ( p->id == O_N )
                 *q = 0;      *q = 0;
         else {    else {
                 NEWQ(tmp);      NEWQ(tmp);
                 SGN(tmp)=1;      SGN(tmp)=1;
                 for ( dc = DC(p), max=0; dc; dc = NEXT(dc) ) {      for ( dc = DC(p), max=0; dc; dc = NEXT(dc) ) {
                         t = (Q)COEF(dc);        t = (Q)COEF(dc);
                         NM(tmp)=NM(t);        NM(tmp)=NM(t);
                         DN(tmp)=DN(t);        DN(tmp)=DN(t);
                         if ( cmpq(tmp, max) > 0 ) DUPQ(tmp, max);        if ( cmpq(tmp, max) > 0 ) DUPQ(tmp, max);
                 }      }
                 addq(ONE, max, q);      addq(ONE, max, q);
         }    }
 }  }
   
 void  void
Line 397  P p;
Line 403  P p;
 Q q;  Q q;
 Q *rp;  Q *rp;
 {  {
         Q   d, d1, a, b, t;    Q   d, d1, a, b, t;
         Q   deg, da;    Q   deg, da;
         Q   nm, dn;    Q   nm, dn;
         DCP dc;    DCP dc;
   
         if ( !p ) *rp = 0;    if ( !p ) *rp = 0;
         else if ( NUM(p) ) *rp = (Q)p;    else if ( NUM(p) ) *rp = (Q)p;
         else {    else {
                 if ( q ) {      if ( q ) {
                         NTOQ( DN(q), 1, dn );        NTOQ( DN(q), 1, dn );
                         NTOQ( NM(q), SGN(q), nm );        NTOQ( NM(q), SGN(q), nm );
                 } else {      } else {
                         dn = 0;        dn = 0;
                         nm = 0;        nm = 0;
                 }      }
                 if ( !dn ) {      if ( !dn ) {
                         dc = DC(p); t = (Q)COEF(dc);        dc = DC(p); t = (Q)COEF(dc);
                         for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc= NEXT(dc) ) {        for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc= NEXT(dc) ) {
                                 subq(d, DEG(dc), &d1); pwrq(nm, d1, &a);          subq(d, DEG(dc), &d1); pwrq(nm, d1, &a);
                                 mulq(t,a,&b); addq(b,(Q)COEF(dc),&t);          mulq(t,a,&b); addq(b,(Q)COEF(dc),&t);
                         }        }
                         if ( d ) {        if ( d ) {
                                 pwrq(nm,d,&a); mulq(t,a,&b); t = b;          pwrq(nm,d,&a); mulq(t,a,&b); t = b;
                         }        }
                         *rp = t;        *rp = t;
                 } else {      } else {
                         dc = DC(p); t = (Q)COEF(dc);        dc = DC(p); t = (Q)COEF(dc);
                         for ( d=deg= DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc= NEXT(dc) ) {        for ( d=deg= DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc= NEXT(dc) ) {
                                 subq(d, DEG(dc), &d1); pwrq(nm, d1, &a);          subq(d, DEG(dc), &d1); pwrq(nm, d1, &a);
                                 mulq(t,a,&b);          mulq(t,a,&b);
                                 subq(deg, DEG(dc), &d1); pwrq(dn, d1, &a);          subq(deg, DEG(dc), &d1); pwrq(dn, d1, &a);
                                 mulq(a, (Q)COEF(dc), &da);          mulq(a, (Q)COEF(dc), &da);
                                 addq(b,da,&t);          addq(b,da,&t);
                         }        }
                         if ( d ) {        if ( d ) {
                                 pwrq(nm,d,&a); mulq(t,a,&b); t = b;          pwrq(nm,d,&a); mulq(t,a,&b); t = b;
                         }        }
                         *rp = t;        *rp = t;
                 }      }
         }    }
 }  }
 #endif  #endif
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.8
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>