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

Diff for /OpenXM_contrib2/asir2018/builtin/isolv.c between version 1.3 and 1.4

version 1.3, 2019/10/17 03:03:12 version 1.4, 2019/12/24 10:26:38
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2018/builtin/isolv.c,v 1.2 2019/06/04 07:11:23 kondoh Exp $   * $OpenXM: OpenXM_contrib2/asir2018/builtin/isolv.c,v 1.3 2019/10/17 03:03:12 kondoh Exp $
  */   */
   
 #include "ca.h"  #include "ca.h"
Line 8 
Line 8 
   
 #if defined(INTERVAL)  #if defined(INTERVAL)
   
   /* in Q.c */
   void dupz(Z, Z*);
   
 static void Solve(NODE, Obj *);  static void Solve(NODE, Obj *);
 static void NSolve(NODE, Obj *);  static void NSolve(NODE, Obj *);
   
Line 38  Solve(arg, rp)
Line 41  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;
     Q                     Eps;
     P                     p;
   
   
       *rp = 0;    p = (P)ARG0(arg);
 #if 0  
   p = (pointer)ARG0(arg);  
   if ( !p ) {    if ( !p ) {
     *rp = 0;      *rp = 0;
     return;      return;
   }    }
   Eps = (pointer)ARG1(arg);    Eps = (Q)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;
     NEWQ(eps); mpq_abs(BDY(eps),BDY(Eps));
   
   switch (OID(p)) {    switch (OID(p)) {
     case O_N:      case O_N:
       *rp = 0;        *rp = 0;
Line 71  Obj  *rp;
Line 76  Obj  *rp;
         error(" : invalid argument");          error(" : invalid argument");
         break;          break;
       }        }
       Solve1((P)p, eps, &root);        Solve1(p, eps, &root);
       *rp = (Obj)root;        *rp = (Obj)root;
       break;        break;
     case O_LIST:      case O_LIST:
Line 81  Obj  *rp;
Line 86  Obj  *rp;
     default:      default:
       *rp = 0;        *rp = 0;
   }    }
 #endif  
 }  }
   
 static void  static void
 NSolve(NODE arg, Obj *rp)  NSolve(NODE arg, Obj *rp)
 {  {
   pointer  p, Eps;    P                     p;
     Q                     Eps;
     Q                     eps;
   pointer  root;    pointer  root;
   LIST    listp;    LIST    listp;
   V      v;    V      v;
   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;
   
       *rp = 0;    p = (P)ARG0(arg);
 #if 0  
   
   p = (pointer)ARG0(arg);  
   if ( !p ) {    if ( !p ) {
     *rp = 0;      *rp = 0;
     return;      return;
   }    }
   Eps = (pointer)ARG1(arg);    Eps = (Q)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;
     NEWQ(eps); mpq_abs(BDY(eps),BDY(Eps));
   
   switch (OID(p)) {    switch (OID(p)) {
     case O_N:      case O_N:
       *rp = 0;        *rp = 0;
Line 125  NSolve(NODE arg, Obj *rp)
Line 129  NSolve(NODE arg, Obj *rp)
         error(" : invalid argument");          error(" : invalid argument");
         break;          break;
       }        }
       Solve1((P)p, eps, &root);        Solve1(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);
Line 171  pointer *rt;
Line 175  pointer *rt;
     return;      return;
   }    }
   Sturm(inp, &sseq);    Sturm(inp, &sseq);
   DUPQ(up,low);    //DUPQ(up,low);
   SGN(low) = -1;    //SGN(low) = -1;
     NEWQ(low); mpq_neg(BDY(low),BDY(up));
   chvu = stumq(sseq, up);    chvu = stumq(sseq, up);
   chvl = stumq(sseq, low);    chvl = stumq(sseq, low);
   tnumb = abs(chvu - chvl);    tnumb = abs(chvu - chvl);
Line 183  pointer *rt;
Line 188  pointer *rt;
   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 ( !cmpz(DEG(onedeg), ONE) ) {
       pad++;        pad++;
       if ( !NEXT(onedeg) ) {        if ( !NEXT(onedeg) ) {
         BDY(root)[pad][0] = 0;          BDY(root)[pad][0] = 0;
Line 201  pointer *rt;
Line 206  pointer *rt;
     }      }
     boundbody(p, &up);      boundbody(p, &up);
     Sturm(p, &sseq);      Sturm(p, &sseq);
     DUPQ(up,low);      //DUPQ(up,low);
     SGN(low) = -1;      //SGN(low) = -1;
       NEWQ(low); mpq_neg(BDY(low),BDY(up));
     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((Q)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++) {
Line 269  int  *padp;
Line 275  int  *padp;
   int de, midn;    int de, midn;
   Q   mid, e;    Q   mid, e;
   P   p;    P   p;
     Q tmp_e;
   
   de = abs(lown - upn);    de = abs(lown - upn);
   if (de == 0) return;    if (de == 0) return;
Line 277  int  *padp;
Line 284  int  *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], &tmp_e );
     SGN(e) = 1;      //SGN(e) = 1;
       absq(tmp_e, &e);
     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], &tmp_e);
       SGN(e) = 1;        //SGN(e) = 1;
         absq(tmp_e, &e);
     }      }
     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, (Q)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);
Line 303  MAT root;
Line 312  MAT root;
   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, (Q)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 (sgnq(a) == sgnq(b)){
     BDY(root)[indx][1] = d;      BDY(root)[indx][1] = d;
   } else {    } else {
     BDY(root)[indx][0] = d;      BDY(root)[indx][0] = d;
Line 327  VECT *ret;
Line 336  VECT *ret;
   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; ; ) {    NEWQ(h);
     NEWQ(x);
     dupz(ONE, (Z *)&h);
     dupz(ONE, (Z *)&x);
     for ( i = 1; ; ) {
     if ( NUM(g2) ) break;      if ( NUM(g2) ) break;
     subq(DEG(DC(g1)),DEG(DC(g2)),&d);      subz(DEG(DC(g1)),DEG(DC(g2)),(Z*)&d);
     l = (Q)LC(g2);      l = (Q)LC(g2);
     if ( SGN(l) < 0 ) {      if ( sgnq(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,(Q)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++;
Line 344  VECT *ret;
Line 357  VECT *ret;
     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 ( sgnq(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);
Line 370  stumq(VECT s, Q val)
Line 383  stumq(VECT s, Q val)
   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( (sgnq(a) > 0 && sgnq(a0) < 0) || (sgnq(a) < 0 && sgnq(a0) > 0) ){
         c++;          c++;
         a0 = a;          a0 = a;
       }        }
Line 392  Q *q;
Line 405  Q *q;
   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) ) {      NEWQ(max);
       for ( dc = DC(p); dc; dc = NEXT(dc) ) {
       t = (Q)COEF(dc);        t = (Q)COEF(dc);
       NM(tmp)=NM(t);        if ( cmpq(t, max) > 0 ) { //DUPQ(tmp, max);
       DN(tmp)=DN(t);          mpq_abs(BDY(max),BDY(t));
       if ( cmpq(tmp, max) > 0 ) DUPQ(tmp, max);        }
     }      }
     addq(ONE, max, q);      addq((Q)ONE, max, q);
       //mpq_clear(max);
   }    }
 }  }
   
Line 410  P p;
Line 425  P p;
 Q q;  Q q;
 Q *rp;  Q *rp;
 {  {
   Q   d, d1, a, b, t;    Z   d, d1;
   Q   deg, da;    Z   deg;
   Q   nm, dn;    Q   da;
     Q   a, b, t;
     Z   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 );        dnq(q, &dn);
         //NTOQ( NM(q), sgnq(q), nm );
         nmq(q, &nm);
     } else {      } else {
       dn = 0;        dn = 0;
       nm = 0;        nm = 0;
Line 428  Q *rp;
Line 447  Q *rp;
     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);          subz(d, DEG(dc), &d1); pwrq((Q)nm, (Q)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((Q)nm,(Q)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);          subz(d, DEG(dc), &d1); pwrq((Q)nm, (Q)d1, &a);
         mulq(t,a,&b);          mulq(t,a,&b);
         subq(deg, DEG(dc), &d1); pwrq(dn, d1, &a);          subz(deg, DEG(dc), &d1); pwrq((Q)dn, (Q)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;          pwrz(nm,d,(Z*)&a); mulq(t,a,&b); t = b;
       }        }
       *rp = t;        *rp = t;
     }      }
   }    }
 #endif  
 }  }
 #endif  #endif
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>