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

Diff for /OpenXM_contrib2/asir2000/plot/calc.c between version 1.1 and 1.15

version 1.1, 1999/12/03 07:39:12 version 1.15, 2020/10/04 03:14:09
Line 1 
Line 1 
 /* $OpenXM: OpenXM/src/asir99/plot/calc.c,v 1.1.1.1 1999/11/10 08:12:34 noro Exp $ */  /*
    * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
    * All rights reserved.
    *
    * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
    * non-exclusive and royalty-free license to use, copy, modify and
    * redistribute, solely for non-commercial and non-profit purposes, the
    * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
    * conditions of this Agreement. For the avoidance of doubt, you acquire
    * only a limited right to use the SOFTWARE hereunder, and FLL or any
    * third party developer retains all rights, including but not limited to
    * copyrights, in and to the SOFTWARE.
    *
    * (1) FLL does not grant you a license in any way for commercial
    * purposes. You may use the SOFTWARE only for non-commercial and
    * non-profit purposes only, such as academic, research and internal
    * business use.
    * (2) The SOFTWARE is protected by the Copyright Law of Japan and
    * international copyright treaties. If you make copies of the SOFTWARE,
    * with or without modification, as permitted hereunder, you shall affix
    * to all such copies of the SOFTWARE the above copyright notice.
    * (3) An explicit reference to this SOFTWARE and its copyright owner
    * shall be made on your publication or presentation in any form of the
    * results obtained by use of the SOFTWARE.
    * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
    * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
    * for such modification or the source code of the modified part of the
    * SOFTWARE.
    *
    * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
    * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
    * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
    * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
    * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
    * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
    * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
    * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
    * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
    * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
    * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
    * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
    * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
    * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
    * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
    * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
    * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
    *
    * $OpenXM: OpenXM_contrib2/asir2000/plot/calc.c,v 1.14 2018/03/29 01:32:55 noro Exp $
   */
 #include "ca.h"  #include "ca.h"
   #include "parse.h"
 #include "ifplot.h"  #include "ifplot.h"
 #include <math.h>  #include <math.h>
 #if PARI  #if defined(PARI)
 #include "genpari.h"  #include "genpari.h"
 #endif  #endif
   
 double usubstrp(P,double);  #ifndef MAXSHORT
   #define MAXSHORT ((short)0x7fff)
   #endif
   
 void calc(tab,can)  void todouble(Obj,Obj *);
 double **tab;  void Psetprec();
 struct canvas *can;  
 {  
         double x,y,xmin,ymin,xstep,ystep;  
         int ix,iy;  
         Real r,rx,ry;  
         Obj fr,g;  
         int w,h;  
         V vx,vy;  
         Obj t,s;  
   
         initmarker(can,"Evaluating...");  void calc(double **tab,struct canvas *can,int nox){
         MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);    //memory_plot,IFPLOTD,INEQND,INEQNANDD,INEQNORD
         vx = can->vx;    //INEQNXORD,conplotmainD
         vy = can->vy;    double x,y,xstep,ystep;
         w = can->width; h = can->height;    int ix,iy;
         xmin = can->xmin; xstep = (can->xmax-can->xmin)/w;    Real r,rx,ry;
         ymin = can->ymin; ystep = (can->ymax-can->ymin)/h;    Obj fr,g,t,s;
         MKReal(1.0,rx); MKReal(1.0,ry); /* dummy real */  
         for( ix = 0, x = xmin; ix < w ; ix++, x += xstep ) {    if(!nox)initmarker(can,"Evaluating...");
 #if 0    todouble((Obj)can->formula,&fr);
                 MKReal(x,r); substp(CO,fr,vx,(P)r,&g);    xstep=(can->xmax-can->xmin)/can->width;
                 marker(can,DIR_X,ix);    ystep=(can->ymax-can->ymin)/can->height;
                 for( iy = 0, y = ymin; iy < h ; iy++, y += ystep )    MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
                         tab[ix][iy] = usubstrp(g,y);    BDY(rx)=can->xmin;
 #endif    substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&t); devalr(CO,t,&g);
                 BDY(rx) = x; substr(CO,0,fr,vx,x?(P)rx:0,&t); devalr(CO,t,&g);    BDY(ry)=can->ymin;
                 marker(can,DIR_X,ix);    substr(CO,0,g,can->vy,can->ymin?(Obj)ry:0,&t); devalr(CO,t,&s);
                 for( iy = 0, y = ymin; iy < h ; iy++, y += ystep ) {    can->vmax=can->vmin=ToReal(s);
                         BDY(ry) = y;    for(ix=0,x=can->xmin; ix<can->width; ix++,x+=xstep){
                         substr(CO,0,g,vy,y?(P)ry:0,&t); devalr(CO,t,&s);      BDY(rx)=x; substr(CO,0,fr,can->vx,x?(Obj)rx:0,&t);
                         tab[ix][iy] = ToReal(s);      devalr(CO,t,&g);
                 }      if(!nox)marker(can,DIR_X,ix);
         }      for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
         BDY(ry)=y;
         substr(CO,0,g,can->vy,y?(Obj)ry:0,&t);
         devalr(CO,t,&s);
         tab[ix][iy]=ToReal(s);
         if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
         if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
       }
     }
 }  }
   
 double usubstrp(p,r)  void calcq(double **tab,struct canvas *can,int nox){
 P p;    //IFPLOTQ,INEQNQ,INEQNANDQ,INEQNORQ,INEQNXORQ
 double r;    //plotoverD
 {    Q dx,dy,xstep,ystep,q1,w,h,c;
         double t;    P g,g1,f1,f2,x,y;
         DCP dc;    int ix,iy;
         int d;    Obj fr,gm,t,s;
         double pwrreal0();    Real r,rx,ry;
   
         if ( !p )    todouble((Obj)can->formula,&fr);
                 t = 0.0;    MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
         else if ( NUM(p) )    BDY(rx)=can->xmin;
                 t = BDY((Real)p);    substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&t); devalr(CO,t,&gm);
         else {    BDY(ry)=can->ymin;
                 dc = DC(p); t = BDY((Real)COEF(dc));    substr(CO,0,gm,can->vy,can->ymin?(Obj)ry:0,&t); devalr(CO,t,&s);
                 for ( d = QTOS(DEG(dc)), dc = NEXT(dc); dc;    can->vmax=can->vmin=ToReal(s);
                         d = QTOS(DEG(dc)), dc = NEXT(dc) ) {  
                         t = t*pwrreal0(r,(d-QTOS(DEG(dc))))+BDY((Real)COEF(dc));    subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
                 }    subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
                 if ( d )    MKV(can->vx,x); mulp(CO,(P)xstep,x,(P *)&t);
                         t *= pwrreal0(r,d);    addp(CO,(P)can->qxmin,(P)t,(P *)&s); substp(CO,can->formula,can->vx,(P)s,&f1);
         }    MKV(can->vy,y); mulp(CO,(P)ystep,y,(P *)&t);
         return t;    addp(CO,(P)can->qymin,(P)t,(P *)&s); substp(CO,f1,can->vy,(P)s,&f2);
     ptozp(f2,1,&c,&g);
     if(!nox) initmarker(can,"Evaluating...");
     for(iy=0;iy<can->height;iy++){
       marker(can,DIR_Y,iy);
       STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,(P *)&t); ptozp((P)t,1,&c,&g1);
       for(ix=0;ix<can->width;ix++){
         STOQ(ix,q1);substp(CO,g1,can->vx,(P)q1,(P *)&t);
         devalr(CO,t,&s);
         tab[ix][iy]=ToReal(s);
         if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
         if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
       }
     }
 }  }
   
 void qcalc(tab,can)  void calcb(double **tab,struct canvas *can,int nox){
 char **tab;    //IFPLOTB,INEQNB,INEQNANDB,INEQNORB,INEQNXORB
 struct canvas *can;    Q dx,dy,xstep,ystep,q1,w,h,c;
 {    P g,g1,f1,f2,x,y,t,s;
         Q dx,dy,w,h,xstep,ystep,c,q1,q2;    int ix,iy,*a,*pa;
         P g,g1,f1,f2,x,y,t,s;    VECT ss;
         DCP dc;    Obj fr,gm,tm,sm;
         int ix,iy;    Real r,rx,ry;
         int *a,*pa;  
         char *px;  
         VECT ss;  
   
     todouble((Obj)can->formula,&fr);
         subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);    MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
         subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);    BDY(rx)=can->xmin;
         MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);    substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&tm); devalr(CO,tm,&gm);
         addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);    BDY(ry)=can->ymin;
         MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);    substr(CO,0,gm,can->vy,can->ymin?(Obj)ry:0,&tm); devalr(CO,tm,&sm);
         addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);    can->vmax=can->vmin=ToReal(sm);
         ptozp(f2,1,&c,&g);  
         a = (int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));    for(iy=0;iy<can->height;iy++)for(ix=0;ix<can->width;ix++)tab[ix][iy]=1.0;
         initmarker(can,"Horizontal scan...");    subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
         for( ix = 0; ix < can->width; ix++ ) {    subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
                 marker(can,DIR_X,ix);    MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
                 STOQ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);    addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
                 if ( !g1 )    MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
                         for ( iy = 0; iy < can->height; iy++ )    addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
                                 tab[ix][iy] = 1;    ptozp(f2,1,&c,&g);
                 else if ( !NUM(g1) ) {    a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
                         strum(CO,g1,&ss); seproot(ss,0,can->height,a);    for(iy=0;iy<can->height;iy++)for(ix=0;ix<can->width;ix++)tab[ix][iy]=1.0;
                         for ( iy = 0, pa = a; iy < can->height; iy++, pa++ )    subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
                                 if ( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )    subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
                                         tab[ix][iy] = 1;    MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
                 }    addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
         }    MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
         initmarker(can,"Vertical scan...");    addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
         for( iy = 0; iy < can->height; iy++ ) {    ptozp(f2,1,&c,&g);
                 marker(can,DIR_Y,iy);    a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
                 STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);    for(ix=0;ix<can->width;ix++){
                 if ( !g1 )      STOQ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);
                         for ( ix = 0; ix < can->width; ix++ )      if(!g1)for(iy=0;iy<can->height;iy++)tab[ix][iy]=0.0;
                                 tab[ix][iy] = 1;      else if(!NUM(g1)){
                 else if ( !NUM(g1) ) {        sturmseq(CO,g1,&ss);
                         strum(CO,g1,&ss); seproot(ss,0,can->width,a);        seproot(ss,0,can->width,a);
                         for ( ix = 0, pa = a; ix < can->width; ix++, pa++ )        for(iy=0,pa=a;iy<can->height;iy++,pa++){
                                 if ( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )          if(*pa<0||(*(pa+1)>=0&&(*pa>*(pa+1))))tab[ix][iy]=0.0;
                                         tab[ix][iy] = 1;          else {
                 }            STOQ(iy,q1);substp(CO,g1,can->vy,(P)q1,&t);
         }            devalr(CO,(Obj)t,(Obj *)&s);
             tab[ix][iy]=ToReal(s);
             if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
             if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
           }
         }
       }
     }
     for(iy=0;iy<can->height;iy++){
       STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);
       if(!g1) for(ix=0;ix<can->width;ix++)tab[ix][iy]=0.0;
       else if(!NUM(g1)){
         sturmseq(CO,g1,&ss);
         seproot(ss,0,can->height,a);
         for(ix=0,pa=a;ix<can->width;ix++,pa++){
           if(tab[ix][iy]!=0.0){
             if(*pa<0||(*(pa+1)>=0&&(*pa>*(pa+1))))tab[ix][iy]=0.0;
             else {
               STOQ(ix,q1);substp(CO,g1,can->vx,(P)q1,&t);
               devalr(CO,(Obj)t,(Obj *)&s);
               tab[ix][iy]=ToReal(s);
               if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
               if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
             }
           }
         }
       }
     }
 }  }
   
 strum(vl,p,rp)  double usubstrp(P p,double r){
 VL vl;    DCP dc;
 P p;    int d;
 VECT *rp;    double t,pwrreal0();
 {  
         P g1,g2,q,r,s;  
         P *t;  
         V v;  
         VECT ret;  
         int i,j;  
         Q a,b,c,d,h,l,m,x;  
   
         v = VR(p); t = (P *)ALLOCA((deg(v,p)+1)*sizeof(P));    if(!p) t=0.0;
         g1 = t[0] = p; diffp(vl,p,v,(P *)&a); ptozp((P)a,1,&c,&g2); t[1] = g2;    else if(NUM(p))t=BDY((Real)p);
         for ( i = 1, h = ONE, x = ONE; ; ) {    else {
                 if ( NUM(g2) )      dc=DC(p); t=BDY((Real)COEF(dc));
                         break;      for(d=QTOS(DEG(dc)),dc=NEXT(dc);dc;d=QTOS(DEG(dc)),dc=NEXT(dc)){
                 subq(DEG(DC(g1)),DEG(DC(g2)),&d);        t=t*pwrreal0(r,(d-QTOS(DEG(dc))))+BDY((Real)COEF(dc));
                 l = (Q)LC(g2);      }
                 if ( SGN(l) < 0 ) {      if(d)t*=pwrreal0(r,d);
                         chsgnq(l,&a); l = a;    }
                 }    return t;
                 addq(d,ONE,&a); pwrq(l,a,&b); mulp(vl,(P)b,g1,(P *)&a);  
                 divsrp(vl,(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(vl,r,(P)a,&g2); t[i] = g2;  
                 x = (Q)LC(g1);  
                 if ( SGN(x) < 0 ) {  
                         chsgnq(x,&a); x = a;  
                 }  
                 pwrq(x,d,&a); mulq(a,h,&b); divq(b,m,&h);  
         }  
         MKVECT(ret,i+1);  
         for ( j = 0; j <= i; j++ )  
                 ret->body[j] = (pointer)t[j];  
         *rp = ret;  
 }  }
   
 seproot(s,min,max,ar)  void qcalc(char **tab,struct canvas *can){
 VECT s;    //qifplotmain(Old type)
 int min,max;    Q dx,dy,w,h,xstep,ystep,c,q1;
 int *ar;    P g,g1,f1,f2,x,y,t,s;
 {    int ix,iy;
         P f;    int *a,*pa;
         P *ss;    VECT ss;
         Q q,t;  
         int i,j,k;  
   
         ss = (P *)s->body; f = ss[0];    subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
         for ( i = min; i <= max; i++ ) {    subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
                 STOQ(i,q); usubstqp(f,q,&t);    MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
                 if ( !t )    addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
                         ar[i] = -1;    MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
                 else {    addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
                         ar[i] = numch(s,q,t); break;    ptozp(f2,1,&c,&g);
                 }    a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
         }    initmarker(can,"Horizontal scan...");
         if ( i > max )    for( ix=0; ix < can->width; ix++ ){
                 return;      marker(can,DIR_X,ix);
         for ( j = max; j >= min; j-- ) {      STOQ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);
                 STOQ(j,q); usubstqp(f,q,&t);      if( !g1 )
                 if ( !t )        for(iy=0; iy < can->height; iy++ )
                         ar[j] = -1;          tab[ix][iy]=1;
                 else {      else if( !NUM(g1) ){
                         if ( i != j )        sturmseq(CO,g1,&ss); seproot(ss,0,can->height,a);
                                 ar[j] = numch(s,q,t);        for(iy=0, pa=a; iy < can->height; iy++, pa++ )
                         break;          if( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )
                 }            tab[ix][iy]=1;
         }      }
         if ( j <= i+1 )    }
                 return;    initmarker(can,"Vertical scan...");
         if ( ar[i] == ar[j] ) {    for( iy=0; iy < can->height; iy++ ){
                 for ( k = i+1; k < j; k++ )      marker(can,DIR_Y,iy);
                         ar[k] = ar[i];      STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);
                 return;      if( !g1 )
         }        for(ix=0; ix < can->width; ix++ )
         k = (i+j)/2;          tab[ix][iy]=1;
         seproot(s,i,k,ar);      else if( !NUM(g1) ){
         seproot(s,k,j,ar);        sturmseq(CO,g1,&ss); seproot(ss,0,can->width,a);
         for(ix=0, pa=a; ix < can->width; ix++, pa++ )
           if( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )
             tab[ix][iy]=1;
       }
     }
 }  }
   
 numch(s,n,a0)  void sturmseq(VL vl,P p,VECT *rp){
 VECT s;    P g1,g2,q,r,s,*t;
 Q n,a0;    V v;
 {    VECT ret;
         int len,i,c;    int i,j;
         Q a;    Q a,b,c,d,h,l,m,x;
         P *ss;  
   
         len = s->len; ss = (P *)s->body;    v=VR(p);t=(P *)ALLOCA((deg(v,p)+1)*sizeof(P));
         for ( i = 1, c = 0; i < len; i++ ) {    g1=t[0]=p;diffp(vl,p,v,(P *)&a);ptozp((P)a,1,&c,&g2);t[1]=g2;
                 usubstqp(ss[i],n,&a);    for(i=1,h=ONE,x=ONE;;){
                 if ( a ) {      if(NUM(g2)) break;
                         if ( (SGN(a)>0 && SGN(a0)<0) || (SGN(a)<0 && SGN(a0)>0) )      subq(DEG(DC(g1)),DEG(DC(g2)),&d);
                                 c++;      l=(Q)LC(g2);
                         a0 = a;      if(SGN(l)<0){
                 }        chsgnq(l,&a);l=a;
         }      }
         return c;      addq(d,ONE,&a);pwrq(l,a,&b);mulp(vl,(P)b,g1,(P *)&a);
       divsrp(vl,(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(vl,r,(P)a,&g2);t[i]=g2;
       x=(Q)LC(g1);
       if(SGN(x)<0){
         chsgnq(x,&a);x=a;
       }
       pwrq(x,d,&a);mulq(a,h,&b);divq(b,m,&h);
     }
     MKVECT(ret,i+1);
     for(j=0;j<=i;j++)
       ret->body[j]=(pointer)t[j];
     *rp=ret;
 }  }
   
 usubstqp(p,r,v)  void seproot(VECT s,int min,int max,int *ar){
 P p;    P f,*ss;
 Q r;    Q q,t;
 Q *v;    int i,j,k;
 {  
         Q d,d1,a,b,t;  
         DCP dc;  
   
         if ( !p )    ss=(P *)s->body;f=ss[0];
                 *v = 0;    for(i=min;i<=max;i++){
         else if ( NUM(p) )      STOQ(i,q);usubstqp(f,q,&t);
                 *v = (Q)p;      if(!t)ar[i]=-1;
         else {      else {
                 dc = DC(p); t = (Q)COEF(dc);        ar[i]=numch(s,q,t);break;
                 for ( d = DEG(dc), dc = NEXT(dc); dc;      }
                         d = DEG(dc), dc = NEXT(dc) ) {    }
                         subq(d,DEG(dc),&d1); pwrq(r,d1,&a);    if(i>max) return;
                         mulq(t,a,&b); addq(b,(Q)COEF(dc),&t);    for(j=max;j>= min;j--){
                 }      STOQ(j,q); usubstqp(f,q,&t);
                 if ( d ) {      if(!t)ar[j]=-1;
                         pwrq(r,d,&a); mulq(t,a,&b); t = b;      else {
                 }        if(i!=j)ar[j]=numch(s,q,t);
                 *v = t;        break;
         }      }
     }
     if(j<=i+1) return;
     if(ar[i]==ar[j]){
       for(k=i+1;k<j;k++)ar[k]=ar[i];
       return;
     }
     k=(i+j)/2;
     seproot(s,i,k,ar);
     seproot(s,k,j,ar);
 }  }
   
 void plotcalc(can)  int numch(VECT s,Q n,Q a0){
 struct canvas *can;    int len,i,c;
 {    Q a;
         double x,xmin,xstep,ymax,ymin,dy;    P *ss;
         int ix;  
         Real r;  
         Obj fr;  
         double usubstrp();  
         int w,h;  
         double *tab;  
         POINT *pa;  
         Real rx;  
         Obj t,s;  
   
         MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);    len=s->len;ss=(P *)s->body;
         w = can->width; h = can->height;    for(i=1,c=0;i<len;i++){
         xmin = can->xmin; xstep = (can->xmax-can->xmin)/w;      usubstqp(ss[i],n,&a);
         tab = (double *)ALLOCA(w*sizeof(double));      if(a){
         MKReal(1,rx); /* dummy real number */        if((SGN(a)>0 && SGN(a0)<0)||(SGN(a)<0&&SGN(a0)>0))c++;
         for( ix = 0, x = xmin; ix < w ; ix++, x += xstep ) {        a0=a;
                 /* full substitution */      }
                 BDY(rx) = x;    }
                 substr(CO,0,fr,can->vx,x?(P)rx:0,&s); devalr(CO,(Obj)s,&t);    return c;
                 if ( t && (OID(t)!=O_N || NID((Num)t)!=N_R) )  }
                         error("plotcalc : invalid evaluation");  
                 tab[ix] = ToReal((Num)t);  
 #if 0  
                 tab[ix] = usubstrp(fr,x);  
 #endif  
         }  
         if ( can->ymax == can->ymin ) {  
                 for ( ymax = ymin = tab[0], ix = 1; ix < w; ix++ ) {  
                         if ( tab[ix] > ymax )  
                                 ymax = tab[ix];  
                         if ( tab[ix] < ymin )  
                                 ymin = tab[ix];  
                 }  
                 can->ymax = ymax; can->ymin = ymin;  
         } else {  
                 ymax = can->ymax; ymin = can->ymin;  
         }  
         dy = ymax-ymin;  
         can->pa = (struct pa *)MALLOC(sizeof(struct pa));  
         can->pa[0].length = w;  
         can->pa[0].pos = pa = (POINT *)MALLOC(w*sizeof(POINT));  
         for ( ix = 0; ix < w; ix++ ) {  
 #ifndef MAXSHORT  
 #define MAXSHORT ((short)0x7fff)  
 #endif  
                 double t;  
   
                 XC(pa[ix]) = ix;  void usubstqp(P p,Q r,Q *v){
                 t = (h - 1)*(ymax - tab[ix])/dy;    Q d,d1,a,b,t;
                 if ( t > MAXSHORT )    DCP dc;
                         YC(pa[ix]) = MAXSHORT;  
                 else if ( t < -MAXSHORT )    if(!p)
                         YC(pa[ix]) = -MAXSHORT;      *v=0;
                 else    else if(NUM(p))*v=(Q)p;
                         YC(pa[ix]) = t;    else {
         }      dc=DC(p);t=(Q)COEF(dc);
       for(d=DEG(dc),dc=NEXT(dc);dc;d=DEG(dc),dc=NEXT(dc)){
         subq(d,DEG(dc),&d1);pwrq(r,d1,&a);
         mulq(t,a,&b);addq(b,(Q)COEF(dc),&t);
       }
       if(d){
         pwrq(r,d,&a);mulq(t,a,&b);t=b;
       }
       *v=t;
     }
 }  }
   
   Num tobf(Num,int);
   
   void plotcalcbf(struct canvas *can){
     Obj fr,s,t;
     Num xmin,xmax,ymin,ymax,xstep;
     Num u,v,ha,dx,dy,x;
     Num *tab;
     Real r;
     Q w,h1;
     int ix;
     POINT *pa;
     double rr;
     Q prec;
     NODE arg;
   
     STOQ(can->prec,prec); arg = mknode(1,prec); Psetprec(arg,&t);
     evalr(CO,(Obj)can->formula,can->prec,&fr);
     MKReal(can->xmin,r); xmin = tobf((Num)r,can->prec);
     MKReal(can->xmax,r); xmax = tobf((Num)r,can->prec);
     MKReal(can->ymin,r); ymin = tobf((Num)r,can->prec);
     MKReal(can->ymax,r); ymax = tobf((Num)r,can->prec);
     STOQ(can->width,w);
     subbf(xmax,xmin,&dx); divbf(dx,(Num)w,&xstep);
     tab=(Num *)MALLOC(can->width*sizeof(Num));
     for(ix=0,x=xmin;ix<can->width;ix++){
       substr(CO,0,fr,can->vx,(Obj)x,(Obj *)&s);
       evalr(CO,(Obj)s,can->prec,&t);
       if(t&&(OID(t)!=O_N))
         error("plotcalcbf : invalid evaluation");
       tab[ix]=(Num)t;
       addbf(x,xstep,&u); x = u;
     }
     if(!cmpbf(ymax,ymin)){
       for(ymax=ymin=tab[0],ix=1;ix<can->width;ix++){
         if(cmpbf(tab[ix],ymax)>0)ymax=tab[ix];
         if(cmpbf(tab[ix],ymin)<0)ymin=tab[ix];
       }
       can->ymax=ToReal(ymax);can->ymin=ToReal(ymin);
     }
     subbf(ymax,ymin,&dy);
     can->pa=(struct pa *)MALLOC(sizeof(struct pa));
     can->pa[0].length=can->width;
     can->pa[0].pos=pa=(POINT *)MALLOC(can->width*sizeof(POINT));
     STOQ(can->height-1,h1);
     for(ix=0;ix<can->width;ix++){
       XC(pa[ix])=ix;
       subbf(ymax,tab[ix],&u); divbf(u,dy,&v); mulbf(v,(Num)h1,&u);
       rr = ToReal(u);
       if(rr>MAXSHORT)YC(pa[ix])=MAXSHORT;
       else if(rr<-MAXSHORT)YC(pa[ix])=-MAXSHORT;
       else YC(pa[ix])=(long)rr;
     }
   }
   
   void plotcalc(struct canvas *can){
     //plot,memory_plot,plotover,plot_resize
     double x,xmin,xstep,ymax,ymin,dy,*tab,usubstrp();
     int ix,w,h;
     Real r,rx;
     Obj fr,t,s;
     POINT *pa;
   
     if ( can->prec ) {
       plotcalcbf(can);
       return;
     }
     todouble((Obj)can->formula,&fr);
     w=can->width;h=can->height;
     xmin=can->xmin;xstep=(can->xmax-can->xmin)/w;
     tab=(double *)ALLOCA(w*sizeof(double));
     MKReal(1,rx); // dummy real number
     for(ix=0,x=xmin;ix<w;ix++,x+=xstep){
       // full substitution
       BDY(rx)=x;
       substr(CO,0,fr,can->vx,x?(Obj)rx:0,&s);
       devalr(CO,(Obj)s,&t);
       if(t&&(OID(t)!=O_N||NID((Num)t)!=N_R))
         error("plotcalc : invalid evaluation");
       tab[ix]=ToReal((Num)t);
     }
     if(can->ymax==can->ymin){
       for(ymax=ymin=tab[0],ix=1;ix<w;ix++){
         if(tab[ix]>ymax)ymax=tab[ix];
         if(tab[ix]<ymin)ymin=tab[ix];
       }
       can->ymax=ymax;can->ymin=ymin;
     } else {
       ymax=can->ymax;ymin=can->ymin;
     }
     dy=ymax-ymin;
     can->pa=(struct pa *)MALLOC(sizeof(struct pa));
     can->pa[0].length=w;
     can->pa[0].pos=pa=(POINT *)MALLOC(w*sizeof(POINT));
     for(ix=0;ix<w;ix++){
       double t;
       XC(pa[ix])=ix;
       t=(h-1)*(ymax-tab[ix])/dy;
       if(t>MAXSHORT)YC(pa[ix])=MAXSHORT;
       else if(t<-MAXSHORT)YC(pa[ix])=-MAXSHORT;
       else YC(pa[ix])=(long)t;
     }
   }
   
   void polarcalc(struct canvas *can){
     double xmax,xmin,ymax,ymin,dx,dy,pmin,pstep,tr,p,*tabx,*taby;
     double usubstrp();
     int i,nstep,w,h;
     POINT *pa;
     Real r;
     Obj fr,t,s;
   
     todouble((Obj)can->formula,&fr);
     w=can->width; h=can->height; nstep=can->nzstep;
     pmin=can->zmin; pstep=(can->zmax-can->zmin)/nstep;
     tabx=(double *)ALLOCA(nstep*sizeof(double));
     taby=(double *)ALLOCA(nstep*sizeof(double));
     MKReal(1,r); // dummy real number
   
     for(i=0,p=pmin;i<nstep;i++,p+= pstep){
       // full substitution
       BDY(r)=p;
       substr(CO,0,fr,can->vx,p?(Obj)r:0,&s);
       devalr(CO,(Obj)s,&t);
       if(t&&(OID(t)!=O_N||NID((Num)t)!=N_R))
         error("polarcalc : invalid evaluation");
       tr=ToReal((Num)t);
       tabx[i]=tr*cos(p);
       taby[i]=tr*sin(p);
     }
     xmax=xmin=tabx[0];
     ymax=ymin=taby[0];
     for(i=1;i<nstep;i++){
       if(tabx[i]>xmax)xmax=tabx[i];
       if(tabx[i]<xmin)xmin=tabx[i];
       if(taby[i]>ymax)ymax=taby[i];
       if(taby[i]<ymin)ymin=taby[i];
     }
     can->xmax=xmax;can->xmin=xmin;
     can->ymax=ymax;can->ymin=ymin;
     dx=xmax-xmin;
     dy=ymax-ymin;
     can->pa=(struct pa *)MALLOC(sizeof(struct pa));
     can->pa[0].length=nstep;
     can->pa[0].pos=pa=(POINT *)MALLOC(w*sizeof(POINT));
     for(i=0;i<nstep;i++){
       XC(pa[i])=(w-1)*(tabx[i]-xmin)/dx;
       YC(pa[i])=(h-1)*(ymax-taby[i])/dy;
     }
   }
   
   void polarcalcNG(struct canvas *can){
     //polarplotNG
     double xmax,xmin,ymax,ymin,dx,dy,pmin,pstep,tr,p, *tabx,*taby;
     double usubstrp();
     int i,ix,iy,nstep,w,h;
     POINT *pa;
     Real r;
     Obj fr,t,s;
   
     todouble((Obj)can->formula,&fr);
     w=can->width; h=can->height; nstep=can->nzstep;
     pmin=can->zmin; pstep=(can->zmax-can->zmin)/nstep;
     tabx=(double *)ALLOCA(nstep*sizeof(double));
     taby=(double *)ALLOCA(nstep*sizeof(double));
     MKReal(1,r); // dummy real number
   
     for(i=0,p=pmin;i<nstep;i++,p+= pstep){
       // full substitution
       BDY(r)=p;
       substr(CO,0,fr,can->vx,p?(Obj)r:0,&s);
       devalr(CO,(Obj)s,&t);
       if(t&&(OID(t)!=O_N||NID((Num)t)!=N_R))
         error("polarcalc : invalid evaluation");
       tr=ToReal((Num)t);
       tabx[i]=tr*cos(p);
       taby[i]=tr*sin(p);
       if(i==0){
         xmax=xmin=tabx[0];
         ymax=ymin=taby[0];
       } else {
         if(tabx[i]>xmax)xmax=tabx[i];
         if(tabx[i]<xmin)xmin=tabx[i];
         if(taby[i]>ymax)ymax=taby[i];
         if(taby[i]<ymin)ymin=taby[i];
       }
     }
     can->xmax=xmax;can->xmin=xmin;
     can->ymax=ymax;can->ymin=ymin;
     dx=xmax-xmin;
     dy=ymax-ymin;
     can->pa=(struct pa *)MALLOC(sizeof(struct pa));
     can->pa[0].length=nstep;
     can->pa[0].pos=pa=(POINT *)MALLOC(w*sizeof(POINT));
     for(i=0;i<nstep;i++){
       XC(pa[i])=(w-1)*(tabx[i]-xmin)/dx;
       YC(pa[i])=(h-1)*(ymax-taby[i])/dy;
     }
   }
   
   /*
   void ineqncalc(double **tab,struct canvas *can,int nox){
     double x,y,xmin,ymin,xstep,ystep;
     int ix,iy,w,h;
     Real r,rx,ry;
     Obj fr,g,t,s;
     V vx,vy;
   
     if(!nox) initmarker(can,"Evaluating...");
     todouble((Obj)can->formula,&fr);
     vx=can->vx;vy=can->vy;
     w=can->width;h=can->height;
     xmin=can->xmin;xstep=(can->xmax-can->xmin)/w;
     ymin=can->ymin;ystep=(can->ymin-can->ymin)/h;
     MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
   
     for(ix=0,x=xmin;ix<=w;ix++,x+=xstep){
       BDY(rx)=x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t);
       devalr(CO,t,&g);
       if(!nox) marker(can,DIR_X,ix);
       for(iy=0,y=ymin;iy<=h;iy++,y+=ystep){
         BDY(ry)=y;
         substr(CO,0,g,vy,y?(Obj)ry:0,&t);
         devalr(CO,t,&s);
         tab[ix][iy]=ToReal(s);
       }
     }
   }
   */
   
   #if defined(INTERVAL)
   void itvcalc(double **mask, struct canvas *can, int nox){
     //ITVIFPLOT
     double x,y,xstep,ystep,dx,dy,wx,wy;
     int idv,ix,iy,idx,idy;
     Itv ity,itx,ddx,ddy;
     Real r,rx,ry,rx1,ry1,rdx,rdy,rdx1,rdy1;
     V vx,vy;
     Obj fr,g,t,s;
   
     idv=can->division;
     todouble((Obj)can->formula,&fr);
     vx=can->vx; vy=can->vy;
     xstep=(can->xmax-can->xmin)/can->width;
     ystep=(can->ymax-can->ymin)/can->height;
     if(idv!=0){
       wx=xstep/can->division;
       wy=ystep/can->division;
     }
     MKReal(can->ymin,ry1);
     for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
       ry=ry1;
       MKReal(y+ystep,ry1);
       istoitv((Num)(ry1),(Num)ry,&ity);
       substr(CO,0,(Obj)fr,vy,(Obj)ity,&t);
       MKReal(can->xmin,rx1);
       for(ix=0,x=can->xmin; ix<can->width; ix++,x+=xstep){
         rx=rx1;
         MKReal(x+xstep,rx1);
         istoitv((Num)(rx1),(Num)rx,&itx);
         substr(CO,0,(Obj)fr,vx,(Obj)itx,&t);
         MKReal(can->ymin,ry1);
         for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
           ry=ry1;
           MKReal(y+ystep,ry1);
           istoitv((Num)ry,(Num)ry1,&ity);
           substr(CO,0,(Obj)t,vy,(Obj)ity,&g);
           if(compnum(0,0,(Num)g))mask[ix][iy]=-1;
           else {
             mask[ix][iy]=0;
   /*
             if(idv==0) mask[ix][iy]=0;
             else {
               MKReal(y,rdy1);
               for(idy=0,dy=y;idy<idv;dy+=wy,idy++){
                 rdy=rdy1;
                 MKReal(dy+wy,rdy1);
                 istoitv((Num)rdy,(Num)rdy1,&ddy);
                 substr(CO,0,(Obj)fr,vy,(Obj)ddy,&t);
                 MKReal(x,rdx1);
                 for(idx=0,dx=x;idx<idx;dx+=wx,idx++){
                   rdx=rdx1;
                   MKReal(dx+wx,rdx1);
                   istoitv((Num)rdx,(Num)rdx1,&ddx);
                   substr(CO,0,(Obj)t,vx,(Obj)ddx,&g);
                   if(!compnum(0,0,(Num)g)){
                     mask[ix][iy]=0;
                     break;
                   }
                 }
                 if(mask[ix][iy]==0)break;
               }
             }
   */
           }
         }
       }
     }
   }
   #endif
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.15
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>