Return to isolv.c CVS log

Up to [local] / OpenXM_contrib2 / asir2018 / builtin

fix interval arithmetic

/*
 * $OpenXM: OpenXM_contrib2/asir2018/builtin/isolv.c,v 1.4 2019/12/24 10:26:38 kondoh Exp $
 */

#include "ca.h"
#include "parse.h"
#include "version.h"

#if defined(INTERVAL) 

/* in Q.c */
void dupz(Z, Z*);

static void Solve(NODE, Obj *);
static void NSolve(NODE, Obj *);

/* in builtin/vars.c */
void Pvars();

/* */
void Solve1(P, Q, pointer *);
void Sturm(P, VECT *);
void boundbody(P, Q *);
void binary(int , MAT);
void separate(Q, Q, VECT, Q, Q, int, int, MAT, int *);
void ueval(P, Q, Q *);
int stumq(VECT, Q);


// in engine/bf.c
Num tobf(Num,int);

struct ftab isolv_tab[] = {
  {"solve", Solve, 2},
  {"nsolve", NSolve, 2},
  {0,0,0},
};

static void
Solve(arg, rp)
NODE arg;
Obj  *rp;
{
  //pointer p, Eps;
  pointer	root;
  V			v;
  Q			eps;
  Q			Eps;
  P			p;


  p = (P)ARG0(arg);
  if ( !p ) {
    *rp = 0;
    return;
  }
  Eps = (Q)ARG1(arg);
  asir_assert(Eps, O_N, "solve");
  if ( NID(Eps) != N_Q ) {
    fprintf(stderr,"solve arg 2 is required a rational number");
    error(" : invalid argument");
    return;
  }
  //DUPQ((Q)Eps, eps);
  //SGN(eps) = 1;
  NEWQ(eps); mpq_abs(BDY(eps),BDY(Eps));
 
  switch (OID(p)) {
    case O_N:
      *rp = 0;
      break;
    case O_P:
      Pvars(arg, &root);
      if (NEXT(BDY((LIST)root)) != 0) {
        fprintf(stderr,"solve arg 1 is univariate polynormial");
        error(" : invalid argument");
        break;
      }
      Solve1(p, eps, &root);
      *rp = (Obj)root;
      break;
    case O_LIST:
      fprintf(stderr,"solve,");
      error(" : Sorry, not yet implement of multivars");
      break;
    default:
      *rp = 0;
  }
}

static void
NSolve(NODE arg, Obj *rp)
{
  P			p;
  Q			Eps;
  Q			eps;
  pointer  root;
  LIST    listp;
  V      v;
  NODE    n, n0, m0, m, ln0;
  Num    r;
  Itv    iv;
  BF      breal;

  p = (P)ARG0(arg);
  if ( !p ) {
    *rp = 0;
    return;
  }
  Eps = (Q)ARG1(arg);
  asir_assert(Eps, O_N, "solve");
  if ( NID(Eps) != N_Q ) {
    fprintf(stderr,"solve arg 2 is required a rational number");
    error(" : invalid argument");
    return;
  }
  //DUPQ((Q)Eps, eps);
  //SGN(eps) = 1;
  NEWQ(eps); mpq_abs(BDY(eps),BDY(Eps));

  switch (OID(p)) {
    case O_N:
      *rp = 0;
      break;
    case O_P:
      Pvars(arg, &root);
      if (NEXT(BDY((LIST)root)) != 0) {
        fprintf(stderr,"solve arg 1 is univariate polynormial");
        error(" : invalid argument");
        break;
      }
      Solve1(p, eps, &root);
      for (m0 = BDY((LIST)root), n0 = 0; m0; m0 = NEXT(m0)) {
        m = BDY((LIST)BDY(m0));
        miditvp(BDY(m), &r);
        //ToBf(r, &breal);
        breal = (BF)tobf(r, DEFAULTPREC);
        NEXTNODE( n0, n );
        MKNODE(ln0, breal, NEXT(m));
        MKLIST(listp, ln0);
        BDY(n) = (pointer)listp;
      }
      NEXT(n) = 0;
      MKLIST(listp,n0);
      *rp = (pointer)listp;
      break;
    case O_LIST:
      fprintf(stderr,"solve,");
      error(" : Sorry, not yet implement of multivars");
      break;
    default:
      *rp = 0;
  }
}

void
Solve1(inp, eps, rt)
P    inp;
Q    eps;
pointer *rt;
{
  P    p;
  Q    up, low, a;
  DCP fctp, onedeg, zerodeg;
  LIST listp;
  VECT sseq;
  MAT  root;
  int  chvu, chvl, pad, tnumb, numb, i, j;
  Itv  iv;
  NODE n0, n, ln0, ln;

  boundbody(inp, &up);
  if (!up) {
    *rt = 0;
    return;
  }
  Sturm(inp, &sseq);
  //DUPQ(up,low);
  //SGN(low) = -1;
  NEWQ(low); mpq_neg(BDY(low),BDY(up));
  chvu = stumq(sseq, up);
  chvl = stumq(sseq, low);
  tnumb = abs(chvu - chvl);
  MKMAT(root, tnumb, 4);
  pad = -1;

  fctrp(CO,inp,&fctp);
  for (fctp = NEXT(fctp), i = 0; fctp; fctp = NEXT(fctp)) {
    p = COEF(fctp);
    onedeg = DC(p);
    if ( !cmpz(DEG(onedeg), ONE) ) {
      pad++;
      if ( !NEXT(onedeg) ) {
        BDY(root)[pad][0] = 0;
        BDY(root)[pad][1] = 0;
        BDY(root)[pad][2] = DEG(fctp);
        BDY(root)[pad][3] = p;
      } else {
        divq((Q)COEF(NEXT(onedeg)),(Q)COEF(onedeg),&a);
        BDY(root)[pad][0] = a;
        BDY(root)[pad][1] = BDY(root)[pad][0];
        BDY(root)[pad][2] = DEG(fctp);
        BDY(root)[pad][3] = p;
      }
      continue;
    }
    boundbody(p, &up);
    Sturm(p, &sseq);
    //DUPQ(up,low);
    //SGN(low) = -1;
    NEWQ(low); mpq_neg(BDY(low),BDY(up));
    chvu = stumq(sseq, up);
    chvl = stumq(sseq, low);
    numb = abs(chvu - chvl);
    separate((Q)DEG(fctp), eps, sseq, up, low, chvu, chvl, root, &pad);
  }
  for (i = 0; i < pad; i++) {
    for (j = i; j <= pad; j++) {
      if (cmpq(BDY(root)[i][0], BDY(root)[j][0]) > 0) {
        a = BDY(root)[i][0];
        BDY(root)[i][0] = BDY(root)[j][0];
        BDY(root)[j][0] = a;
        a = BDY(root)[i][1];
        BDY(root)[i][1] = BDY(root)[j][1];
        BDY(root)[j][1] = a;
        a = BDY(root)[i][2];
        BDY(root)[i][2] = BDY(root)[j][2];
        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++) {
    istoitv(BDY(root)[i][0], BDY(root)[i][1], &iv);
    NEXTNODE(n0,n);
    MKNODE(ln, BDY(root)[i][2], 0); MKNODE(ln0, iv, ln);
    MKLIST(listp, ln0);BDY(n) = (pointer)listp;
  }
  NEXT(n) = 0;
  MKLIST(listp,n0);
  *rt = (pointer)listp;
}

void
separate(mult, eps, sseq, up, low, upn, lown, root, padp)
VECT sseq;
Q    mult, eps, up, low;
int  upn, lown;
MAT  root;
int  *padp;
{
  int de, midn;
  Q   mid, e;
  P   p;
  Q tmp_e;

  de = abs(lown - upn);
  if (de == 0) return;
  if (de == 1) {
    (*padp)++;
    BDY(root)[*padp][0] = up;
    BDY(root)[*padp][1] = low;
    BDY(root)[*padp][3] = (P *)sseq->body[0];
    subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &tmp_e );
    //SGN(e) = 1;
    absq(tmp_e, &e);
    while (cmpq(e, eps) > 0) {
      binary(*padp, root);
      subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &tmp_e);
      //SGN(e) = 1;
      absq(tmp_e, &e);
    }
    BDY(root)[*padp][2] = mult;
    return;
  }
  addq(up, low, &mid);
  divq(mid, (Q)TWO, &mid);
  midn = stumq(sseq, mid);
  separate(mult, eps, sseq, low, mid, lown, midn, 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, (Q)TWO, &d);
  ueval(p, BDY(root)[indx][1], &a);
  ueval(p, d, &b);
  if (sgnq(a) == sgnq(b)){
    BDY(root)[indx][1] = d;
  } else {
    BDY(root)[indx][0] = d;
  }
}

void
Sturm(p, ret)
P    p;
VECT *ret;
{
  P    g1,g2,q,r,s, *t;
  Q    a,b,c,d,h,l,m,x;
  V    v;
  VECT seq;
  int  i,j;

  v = VR(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;
  NEWQ(h);
  NEWQ(x);
  dupz(ONE, (Z *)&h);
  dupz(ONE, (Z *)&x);
  for ( i = 1; ; ) {
    if ( NUM(g2) ) break;
    subz(DEG(DC(g1)),DEG(DC(g2)),(Z*)&d);
    l = (Q)LC(g2);
    if ( sgnq(l) < 0 ) {
      chsgnq(l,&a); l = a;
    }
    addq(d,(Q)ONE,&a); pwrq(l,a,&b); mulp(CO,(P)b,g1,(P *)&a);
    divsrp(CO,(P)a,g2,&q,&r);
    if ( !r ) break;
    chsgnp(r,&s); r = s; i++;
    if ( NUM(r) ) {
      t[i] = r; break;
    }
    pwrq(h,d,&m); g1 = g2;
    mulq(m,x,&a); divsp(CO,r,(P)a,&g2); t[i] = g2;
    x = (Q)LC(g1);
    if ( sgnq(x) < 0 ) {
      chsgnq(x,&a); x = a;
    }
    pwrq(x,d,&a); mulq(a,h,&b); divq(b,m,&h);
  }
  MKVECT(seq,i+1);
  for ( j = 0; j <= i; j++ ) seq->body[j] = (pointer)t[j];
  *ret = seq;
}

int
stumq(VECT s, Q val)
{
  int len, i, j, c;
  P   *ss;
  Q   a, a0;

  len = s->len;
  ss = (P *)s->body;
  for ( j = 0; j < len; j++ ){
    ueval(ss[j],val,&a0);
    if (a0) break;
  }
  for ( i = j++, c =0; i < len; i++) {
    ueval( ss[i], val, &a);
    if ( a ) {
      if( (sgnq(a) > 0 && sgnq(a0) < 0) || (sgnq(a) < 0 && sgnq(a0) > 0) ){
        c++;
        a0 = a;
      }
    }
  }
  return c;
}

void
boundbody(p, q)
P  p;
Q *q;
{
  Q    t, max, tmp;
  DCP  dc;

  if ( !p )
    *q = 0;
  else if ( p->id == O_N )
    *q = 0;
  else {
    //NEWQ(tmp);
    //SGN(tmp)=1;
    NEWQ(max);
    for ( dc = DC(p); dc; dc = NEXT(dc) ) {
      t = (Q)COEF(dc);
      if ( cmpq(t, max) > 0 ) { //DUPQ(tmp, max);
        mpq_abs(BDY(max),BDY(t));
      }
    }
    addq((Q)ONE, max, q);
    //mpq_clear(max);
  }
}

void
ueval(p, q, rp)
P p;
Q q;
Q *rp;
{
  Z   d, d1;
  Z   deg;
  Q   da;
  Q   a, b, t;
  Z   nm, dn;
  DCP dc;

  if ( !p ) *rp = 0;
  else if ( NUM(p) ) *rp = (Q)p;
  else {
    if ( q ) {
      //NTOQ( DN(q), 1, dn );
      dnq(q, &dn);
      //NTOQ( NM(q), sgnq(q), nm );
      nmq(q, &nm);
    } else {
      dn = 0;
      nm = 0;
    }
    if ( !dn ) {
      dc = DC(p); t = (Q)COEF(dc);
      for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc= NEXT(dc) ) {
        subz(d, DEG(dc), &d1); pwrq((Q)nm, (Q)d1, &a);
        mulq(t,a,&b); addq(b,(Q)COEF(dc),&t);
      }
      if ( d ) {
        pwrq((Q)nm,(Q)d,&a); mulq(t,a,&b); t = b;
      }
      *rp = t;
    } else {
      dc = DC(p); t = (Q)COEF(dc);
      for ( d=deg= DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc= NEXT(dc) ) {
        subz(d, DEG(dc), &d1); pwrq((Q)nm, (Q)d1, &a);
        mulq(t,a,&b);
        subz(deg, DEG(dc), &d1); pwrq((Q)dn, (Q)d1, &a);
        mulq(a, (Q)COEF(dc), &da);
        addq(b,da,&t);
      }
      if ( d ) {
        pwrz(nm,d,(Z*)&a); mulq(t,a,&b); t = b;
      }
      *rp = t;
    }
  }
}
#endif