[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.2 and 1.3

version 1.2, 2003/10/20 07:32:19 version 1.3, 2003/10/23 01:32:59
Line 1 
Line 1 
 /*  /*
  * $OpenXM$   * $OpenXM: OpenXM_contrib2/asir2000/builtin/isolv.c,v 1.2 2003/10/20 07:32:19 saito Exp $
  */   */
   
 #include "ca.h"  #include "ca.h"
Line 14  static void NSolve(NODE, Obj *);
Line 14  static void NSolve(NODE, Obj *);
 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 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 *);
   
Line 75  NSolve(arg, rp)
Line 76  NSolve(arg, rp)
 NODE arg;  NODE arg;
 Obj  *rp;  Obj  *rp;
 {  {
         pointer p, Eps;          pointer p, Eps;
         LIST    root, listp;          pointer root, 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 ) {
Line 116  Obj  *rp;
Line 117  Obj  *rp;
                                 ToBf(r, &breal);                                  ToBf(r, &breal);
                                 NEXTNODE( n0, n );                                  NEXTNODE( n0, n );
                                 MKNODE(ln0, breal, NEXT(m));                                  MKNODE(ln0, breal, NEXT(m));
                                 MKLIST(listp, ln0);                                  MKLIST((LIST)listp, ln0);
                                 BDY(n) = (pointer)listp;                                  BDY(n) = (pointer)listp;
                         }                          }
                         NEXT(n) = 0;                          NEXT(n) = 0;
                         MKLIST(listp,n0);                          MKLIST((LIST)listp,n0);
                         *rp = (pointer)listp;                          *rp = (pointer)listp;
                         break;                          break;
                 case O_LIST:                  case O_LIST:
Line 159  pointer *rt;
Line 160  pointer *rt;
         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, 3);          MKMAT(root, tnumb, 4);
         pad = -1;          pad = -1;
   
         fctrp(CO,inp,&fctp);          fctrp(CO,inp,&fctp);
Line 172  pointer *rt;
Line 173  pointer *rt;
                                 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;
                         } 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;
                         }                          }
                         continue;                          continue;
                 }                  }
Line 201  pointer *rt;
Line 204  pointer *rt;
                                 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];
                                   BDY(root)[i][3] = BDY(root)[j][3];
                                   BDY(root)[j][3] = a;
                         }                          }
                 }                  }
         }          }
           for (i = 0; i < pad; i++) {
                   while(cmpq(BDY(root)[i][1], BDY(root)[i+1][0]) > 0 ) {
                           binary(i, root);
                           binary(i+1, root);
                           if ( cmpq(BDY(root)[i][0], BDY(root)[i+1][1]) > 0 ) {
                                   a = BDY(root)[i][0];
                                   BDY(root)[i][0] = BDY(root)[i+1][0];
                                   BDY(root)[i+1][0] = a;
                                   a = BDY(root)[i][1];
                                   BDY(root)[i][1] = BDY(root)[i+1][1];
                                   BDY(root)[i+1][1] = a;
                                   a = BDY(root)[i][2];
                                   BDY(root)[i][2] = BDY(root)[i+1][2];
                                   BDY(root)[i+1][2] = a;
                                   a = BDY(root)[i][3];
                                   BDY(root)[i][3] = BDY(root)[i+1][3];
                                   BDY(root)[i+1][3] = a;
                                   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);
Line 224  MAT root;
Line 251  MAT root;
 int     *padp;  int     *padp;
 {  {
         int de, midn;          int de, midn;
         Q   mid, a, b, c, d, e;          Q   mid, e;
         P   p;          P   p;
   
         de = abs(lown - upn);          de = abs(lown - upn);
Line 233  int *padp;
Line 260  int *padp;
                 (*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];
                 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;
                 p = (P *)sseq->body[0];  
                 while (cmpq(e, eps) > 0) {                  while (cmpq(e, eps) > 0) {
                         addq(BDY(root)[*padp][0], BDY(root)[*padp][1], &c);                          binary(*padp, root);
                         divq(c, TWO, &d);  
                         ueval(p, BDY(root)[*padp][1], &a);  
                         ueval(p, d, &b);  
                         if (SGN(a) == SGN(b))  
                                 BDY(root)[*padp][1] = d;  
                         else BDY(root)[*padp][0] = d;  
                         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;
                 }                  }
Line 255  int *padp;
Line 276  int *padp;
         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
   binary(indx, root)
   int indx;
   MAT root;
   {
           Q       a, b, c, d, e;
           P       p;
           p = (P*)BDY(root)[indx][3];
           addq(BDY(root)[indx][0], BDY(root)[indx][1], &c);
           divq(c, TWO, &d);
           ueval(p, BDY(root)[indx][1], &a);
           ueval(p, d, &b);
           if (SGN(a) == SGN(b)){
                   BDY(root)[indx][1] = d;
           } else {
                   BDY(root)[indx][0] = d;
           }
 }  }
   
 void  void
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>