#include #include #include #include #include #include "stack.h" #include "output.h" #include "funcs.h" #include "interval.h" #include "builtin.h" #include "express.h" #include "spread.h" #include "udf.h" #include "matrix.h" #include "linear.h" #include "command.h" #include "input.h" void scan_summe (void); extern int nosubmref,trace; int udfon=0,actargn=0,actsp=0; char *next; /**************** builtin operators ****************/ void cscalp (double *s, double *si, double *x, double *xi, double *y, double *yi) { *s += *x * *y - *xi * *yi; *si += *x * *yi + *xi * *y; } void iscalp (double *s, double *si, double *x, double *xi, double *y, double *yi) { double a,b; interval_mult(x,xi,y,yi,&a,&b); interval_add(s,si,&a,&b,s,si); } void ccopy (double *y, double *x, double *xi) { *y++=*x; *y=*xi; } void icopy (double *y, double *x, double *xi) { *y++=*x; *y=*xi; } void multiply (header *hd, header *hd1) /***** multiply matrix multiplication. *****/ { header *result,*st=hd; dims *d,*d1; double *m,*m1,*m2,*mm1,*mm2,x,y,null=0.0; int i,j,c,r,k; hd=getvalue(hd); hd1=getvalue(hd1); if (error) return; if (hd->type==s_matrix && hd1->type==s_matrix) { d=dimsof(hd); d1=dimsof(hd1); if (d->c != d1->r) { error=8; output("Cannot multiply these!\n"); return; } r=d->r; c=d1->c; result=new_matrix(r,c,""); if (error) return; m=matrixof(result); m1=matrixof(hd); m2=matrixof(hd1); for (i=0; ic,i,0); mm2=m2+j; x=0.0; for (k=0; kc; k++) { x+=(*mm1)*(*mm2); mm1++; mm2+=d1->c; } *mat(m,c,i,j)=x; } moveresult(st,result); return; } if (hd->type==s_matrix) { if (hd1->type==s_cmatrix) { d=dimsof(hd); d1=dimsof(hd1); if (d->c != d1->r) { error=8; output("Cannot multiply these!\n"); return; } r=d->r; c=d1->c; result=new_cmatrix(r,c,""); if (error) return; m=matrixof(result); m1=matrixof(hd); m2=matrixof(hd1); for (i=0; ic,i,0); mm2=m2+(long)2*j; x=0.0; y=0.0; for (k=0; kc; k++) { cscalp(&x,&y,mm1,&null,mm2,mm2+1); mm1++; mm2+=2*d1->c; } ccopy(cmat(m,c,i,j),&x,&y); } moveresult(st,result); return; } if (hd1->type==s_imatrix) { d=dimsof(hd); d1=dimsof(hd1); if (d->c != d1->r) { error=8; output("Cannot multiply these!\n"); return; } r=d->r; c=d1->c; result=new_imatrix(r,c,""); if (error) return; m=matrixof(result); m1=matrixof(hd); m2=matrixof(hd1); for (i=0; ic,i,0); mm2=m2+(long)2*j; x=0.0; y=0.0; for (k=0; kc; k++) { iscalp(&x,&y,mm1,mm1,mm2,mm2+1); mm1++; mm2+=2*d1->c; } icopy(imat(m,c,i,j),&x,&y); } moveresult(st,result); return; } } if (hd->type==s_cmatrix && hd1->type==s_matrix) { d=dimsof(hd); d1=dimsof(hd1); if (d->c != d1->r) { error=8; output("Cannot multiply these!\n"); return; } r=d->r; c=d1->c; result=new_cmatrix(r,c,""); if (error) return; m=matrixof(result); m1=matrixof(hd); m2=matrixof(hd1); for (i=0; ic,i,0); mm2=m2+j; x=0.0; y=0.0; for (k=0; kc; k++) { cscalp(&x,&y,mm1,mm1+1,mm2,&null); mm1+=2; mm2+=d1->c; } ccopy(cmat(m,c,i,j),&x,&y); } moveresult(st,result); return; } if (hd->type==s_cmatrix && hd1->type==s_cmatrix) { d=dimsof(hd); d1=dimsof(hd1); if (d->c != d1->r) { error=8; output("Cannot multiply these!\n"); return; } r=d->r; c=d1->c; result=new_cmatrix(r,c,""); if (error) return; m=matrixof(result); m1=matrixof(hd); m2=matrixof(hd1); for (i=0; ic,i,0); mm2=m2+(long)2*j; x=0.0; y=0.0; for (k=0; kc; k++) { cscalp(&x,&y,mm1,mm1+1,mm2,mm2+1); mm1+=2; mm2+=2*d1->c; } ccopy(cmat(m,c,i,j),&x,&y); } moveresult(st,result); return; } if (hd->type==s_imatrix) { if (hd1->type==s_matrix) { d=dimsof(hd); d1=dimsof(hd1); if (d->c != d1->r) { error=8; output("Cannot multiply these!\n"); return; } r=d->r; c=d1->c; result=new_imatrix(r,c,""); if (error) return; m=matrixof(result); m1=matrixof(hd); m2=matrixof(hd1); for (i=0; ic,i,0); mm2=m2+j; x=0.0; y=0.0; for (k=0; kc; k++) { iscalp(&x,&y,mm1,mm1+1,mm2,mm2); mm1+=2; mm2+=d1->c; } icopy(imat(m,c,i,j),&x,&y); } moveresult(st,result); return; } if (hd1->type==s_imatrix) { d=dimsof(hd); d1=dimsof(hd1); if (d->c != d1->r) { error=8; output("Cannot multiply these!\n"); return; } r=d->r; c=d1->c; result=new_imatrix(r,c,""); if (error) return; m=matrixof(result); m1=matrixof(hd); m2=matrixof(hd1); for (i=0; ic,i,0); mm2=m2+(long)2*j; x=0.0; y=0.0; for (k=0; kc; k++) { iscalp(&x,&y,mm1,mm1+1,mm2,mm2+1); mm1+=2; mm2+=2*d1->c; } icopy(imat(m,c,i,j),&x,&y); } moveresult(st,result); return; } } else dotmultiply(st,nextof(st)); } void divide (header *hd, header *hd1) { dotdivide(hd,hd1); } void real_invert (double *x, double *y) { *y= -*x; } void complex_invert (double *x, double *xi, double *y, double *yi) { *y= -*x; *yi= -*xi; } void invert (header *hd) /***** invert compute -matrix. *****/ { header *result,*st=hd; hd=getvalue(hd); if (error) return; if (isinterval(hd)) result=map1i(interval_invert,hd); else result=map1(real_invert,complex_invert,hd); if (!error) moveresult(st,result); } /***************** scanning ***************************/ void scan_space (void) { start: while (*next==' ' || *next==TAB) next++; if (!udfon && *next=='.' && *(next+1)=='.') { next_line(); if (error) return; goto start; } } void scan_namemax (char *name, int lmax) { int count=0; if (*next=='\"') { next++; while (*next!='\"' && *next) { *name++=*next++; count++; if (count>=lmax-1) { output("Name too long!\n"); error=11; break; } } if (*next=='\"') next++; } else if (!isalpha(*next) && *next!='_') { output1("Name expected at:\n%s\n",next); error=11; *name=0; return; } else { if (*next=='_') { *name++=*next++; count++; } while (isalpha(*next) || isdigit(*next)) { *name++=*next++; count++; if (count>=lmax-1) { output("Name too long!\n"); error=11; break; } } } *name=0; } void scan_name (char *name) { scan_namemax (name,16); } int scan_arguments (int *scount) /* look ahead for arguments */ { int count=0,olds=nosubmref,nocount=0; header *hd,*hdold; *scount=0; nosubmref=1; while (1) { scan_space(); if (*next==')' || *next==']') break; if (*next==',') hd=new_reference(0,""); else { hd=scan(); scan_space(); } if (*next=='=') { next++; hdold=(header *)newram; scan_value(); if (!error) { strcpy(hdold->name,hd->name); hdold->xor=hd->xor; moveresult(hd,hdold); nocount=1; } } else if (nocount) { output("Illegal parameter after named parameter!\n"); error=2700; } if (error) { nosubmref=olds; return 0; } while (hd<(header *)newram) { if (!nocount) count++; hd=nextof(hd); } if (count>=20) { output("To many arguments!\n"); error=56; return 0; } if (*next==';') { *scount=count; } else if (*next!=',') break; next++; } if (*next!=')' && *next!=']') { output("Error in arguments!\n"); error=19; return 0; } next++; nosubmref=olds; return count; } void scan_matrix (void) /***** scan_matrix scan a matrix from input. form: [x y z ... ; v w u ...], where x,y,z,u,v,w are sums. *****/ { header *hd,*result; dims *d; int c,r,count,i,j,complex=0,interval=0,n; unsigned long ic; unsigned long allcount; double *m,*ms,cnull[2]={0,0},*mr; hd=new_matrix(0,0,""); /* don't know how big it will be! */ if (error) return; count=0; getmatrix(hd,&r,&c,&m); ms=m; r=0; c=0; scan_space(); while (1) { scan_space(); if (*next==0) { next_line(); scan_space(); } /* matrix is allowed to pass line ends */ if (*next==';' || *next==']') /* new column starts */ { if (*next==';') next++; if (!freeramfrom(ms,count*(long)(r+1)*sizeof(double))) { output("Memory overflow!\n"); error=18; return; } if (count>c) /* this column is to long */ { if (r>0) /* expand matrix */ { for (j=count-1; j>=0; j--) { if (complex) copy_complex(cmat(ms,count,r,j), cmat(ms,c,r,j)); else if (interval) copy_interval(imat(ms,count,r,j), imat(ms,c,r,j)); else *mat(ms,count,r,j)=*mat(ms,c,r,j); } for (i=r-1; i>=0; i--) { if (i>0) for (j=c-1; j>=0; j--) { if (complex) copy_complex(cmat(ms,count,i,j), cmat(ms,c,i,j)); else if (interval) copy_interval(imat(ms,count,i,j), imat(ms,c,i,j)); else *mat(ms,count,i,j)=*mat(ms,c,i,j); } for (j=c; jtype==s_complex || result->type==s_cmatrix)) { if (interval) { output("Cannot make an interval matrix complex.\n"); error=16; return; } complex=1; /* make matrix complex up to now (oh boy!) */ allcount=((char *)m-(char *)ms)/sizeof(double); if (!freeram(allcount*sizeof(double)+result->size)) { output("Memory overflow!\n"); error=16; return; } if (allcount) { memmove(newram+allcount*sizeof(double),newram,result->size); result=(header *)((char *)result+allcount*sizeof(double)); for (ic=allcount-1; ic>0; ic--) { *(ms+(long)2*ic)=*(ms+ic); *(ms+(long)2*ic+1)=0.0; } *(ms+1)=0.0; newram=(char *)(ms+(long)2*allcount); m=(double *)newram; } hd->type=s_cmatrix; } else if (!interval && (result->type==s_interval || result->type==s_imatrix)) { if (complex) { output("Cannot make a complex matrix interval.\n"); error=16; return; } interval=1; /* make matrix complex up to now (oh boy!) */ allcount=((char *)m-(char *)ms)/sizeof(double); if (!freeram(allcount*sizeof(double)+result->size)) { output("Memory overflow!\n"); error=16; return; } if (allcount) { memmove(newram+allcount*sizeof(double),newram,result->size); result=(header *)((char *)result+allcount*sizeof(double)); for (ic=allcount-1; ic>0; ic--) { *(ms+(long)2*ic)=*(ms+ic); *(ms+(long)2*ic+1)=*(ms+ic); } *(ms+1)=*ms; newram=(char *)(ms+(long)2*allcount); m=(double *)newram; } hd->type=s_imatrix; } if (result->type==s_real) { if (interval) *m++=*aof(result); else *m++=*realof(result); count++; if (complex) *m++=0.0; else if (interval) *m++=*aof(result); } else if (result->type==s_matrix) { if (dimsof(result)->r!=1 || dimsof(result)->c<1) { output("Illegal matrix element\n"); error=16; return; } n=dimsof(result)->c; mr=matrixof(result); if (complex) { for (j=0; jc; mr=matrixof(result); for (j=0; jtype==s_cmatrix) { if (dimsof(result)->r!=1 || dimsof(result)->c<1) { output("Illegal matrix element\n"); error=16; return; } n=dimsof(result)->c; mr=matrixof(result); if (complex) { for (j=0; jtype==s_imatrix) { if (dimsof(result)->r!=1 || dimsof(result)->c<1) { output("Illegal matrix element\n"); error=16; return; } n=dimsof(result)->c; mr=matrixof(result); if (interval) { for (j=0; jtype==s_complex && complex) { *m++=*realof(result); *m++=*imagof(result); count++; } else if (result->type==s_interval && interval) { *m++=*aof(result); *m++=*bof(result); count++; } else { error=-1; output("Illegal vector!\n"); return; } if (count>=INT_MAX) { output1("Matrix has more than %d columns!\n",INT_MAX); error=17; return; } newram=(char *)m; if (!freeram(16)) { output("Memory overflow!\n"); error=16; return; } } next++; d=(dims *)(hd+1); if (c==0) r=0; d->c=c; d->r=r; if (r>=INT_MAX) { output1("Matrix has more than %d rows!\n",INT_MAX); error=18; return; } if (complex) hd->size=cmatrixsize(c,r); else if (interval) hd->size=imatrixsize(c,r); else hd->size=matrixsize(c,r); newram=(char *)hd+hd->size; } #define myisdigit(c) ((c)>='0' && (c)<='9') void scan_elementary (void) /***** scan_elemtary scan an elementary expression, like a value or variable. scans also (...). *****/ { double x; int n,nargs,hadargs=0,sp; header *hd=(header *)newram,*var,*hd1; char name[16],*s; scan_space(); if ((*next=='.' && myisdigit(*(next+1))) || myisdigit(*next)) { sscanf(next,"%lf%n",&x,&n); next+=n; if (*next=='i') /* complex number! */ { next++; new_complex(0,x,""); } else new_real(x,""); } else if (*next==2) /* a double stored in binary form */ { next++; #ifdef SPECIAL_ALIGNMENT memmove((char *)(&x),next,sizeof(double)); #else x=*((double *)next); #endif next+=sizeof(double); if (*next=='i') /* complex number! */ { next++; new_complex(0,x,""); } else new_real(x,""); } else if (*next==3) { output("Command name used as variable!\n"); error=4000; return; } else if (isalpha(*next) || *next=='_') { scan_name(name); if (error) return; scan_space(); nargs=0; if (*next=='{') { next++; scan(); if (error) return; scan_space(); if (*next!='}') { output("} missing!\n"); error=1010; return; } next++; get_element1(name,hd); goto after; } if (*next=='(' || *next=='[') /* arguments or indices */ { hadargs=(*next=='[')?2:1; next++; nargs=scan_arguments(&sp); if (error) return; } if (name[0]=='_') { if (hadargs==1 && exec_builtin(name+1,nargs,hd)); #ifdef DLL else if (hadargs==1 && exec_dll(name+1,nargs,hd)); #endif else { output1( "%s no function or variable, or wrong argument number!\n", name); error=1; return; } return; } if (hadargs==2) var=searchvar(name); else if (hadargs==1) { var=searchudf(name); if (!var) var=searchvar(name); } else var=searchvar(name); if (var && var->type==s_udf && hadargs==1) { interpret_udf(var,hd,nargs,sp); if (error) return; } else if (!var && hadargs) { if (hadargs==1 && exec_builtin(name,nargs,hd)); #ifdef DLL else if (hadargs==1 && exec_dll(name,nargs,hd)); #endif else { error=24; if (hadargs==2) output1("%s no variable!\n",name); else output1( "%s no function or variable, or wrong argument number!\n", name); return; } } else if (var && hadargs) { get_element(nargs,var,hd); } else hd=new_reference(var,name); } else if (*next=='#' && *(next+1)!='#') { next++; mindex(hd); } else if (*next=='+') { next++; scan_elementary(); } else if (*next=='-') { next++; scan_elementary(); if (!error) invert(hd); } else if (*next=='(') { next++; scan(); if (error) return; scan_space(); if (*next!=')') { output("Closing bracket ) missing!\n"); error=5; return; } newram=(char *)nextof(hd); next++; } else if (*next=='[') { next++; scan_matrix(); } else if (*next=='~' && *(next+1)!='=') { next++; scan_value(); if (error) return; newram=(char *)nextof(hd); hd1=(header *)newram; scan_space(); if (*next=='~') { next++; minterval1(hd); } else { if (*next!=',') { error=1; return; } next++; scan_value(); if (error) return; newram=(char *)nextof(hd1); scan_space(); if (*next!='~') { error=1; return; } next++; minterval(hd); } } else if (*next=='\"') { next++; s=next; while (*next!='\"' && *next!=0) next++; hd=new_string(s,next-s,""); if (*next!='\"') { output("\" missing\n"); error=1; return; } next++; } else if (*next=='\'' && *(next+1)=='\'') { next+=2; s=next; while ((*next!='\'' || *(next+1)!='\'') && *next!=0) next++; hd=new_string(s,next-s,""); if (*next!='\'' || *(next+1)!='\'') { output("\'\' missing\n"); error=1; return; } next+=2; } else error=1; after: if (error) return; /* for things, that come after an elementary expression */ if (*next=='\'') { next++; transpose(hd); return; } if (*next=='°') { next++; mdegree(hd); return; } scan_space(); if (*next=='^' || (*next=='*' && *(next+1)=='*')) { if (*next=='^') next++; else next+=2; newram=(char *)nextof(hd); scan_elementary(); if (error) return; mpower(hd); } } void scan_factor (void) { scan_elementary(); } void scan_summand (void) { header *hd=(header *)newram,*hd1; scan_space(); scan_factor(); if (error) return; while (1) { hd1=(header *)newram; scan_space(); if ((*next=='*' && *(next+1)!='*') || (*next=='.' && *(next+1)=='*')) { if (*next=='*') next++; else next+=2; scan_factor(); if (!error) dotmultiply(hd,hd1); } else if (*next=='/' || (*next=='.' && *(next+1)=='/')) { if (*next=='/') next++; else next+=2; scan_factor(); if (!error) dotdivide(hd,hd1); } else if (*next=='.') { next++; scan_factor(); if (!error) multiply(hd,hd1); } else if (*next=='\\') { next++; newram=(char *)nextof(hd); scan_factor(); if (!error) msolve(hd); } else break; if (error) break; } } void scan_summe (void) { header *hd=(header *)newram,*hd1; scan_space(); scan_summand(); if (error) return; while (1) { hd1=(header *)newram; scan_space(); if (*next=='+') { next++; scan_summand(); if (!error) add(hd,hd1); } else if (*next=='-') { next++; scan_summand(); if (!error) subtract(hd,hd1); } else break; if (error) break; } } void scan_dp (void) { header *hd=(header *)newram,*hd1,*hd2; scan_space(); if (*next==':') { next++; new_command(c_allv); return; } scan_summe(); if (*next==':') /* a vector a:b:c or a:c */ { next++; hd1=(header *)newram; scan_summe(); if (error) return; scan_space(); if (*next==':') { next++; hd2=(header *)newram; scan_summe(); if (error) return; } else { hd2=hd1; hd1=new_real(1.0,""); } if (error) return; vectorize(hd,hd1,hd2); } } void scan_compare (void) { header *hd=(header *)newram; scan_space(); if (*next=='!') { next++; scan_compare(); mnot(hd); return; } scan_dp(); if (error) return; scan_space(); if (*next=='<') { if (*(next+1)!='<') { next++; newram=(char *)nextof(hd); if (*next=='=') { next++; scan_dp(); if (error) return; mlesseq(hd); return; } else if (*next=='>') { next++; scan_dp(); if (error) return; munequal(hd); return; } scan_dp(); if (error) return; mless(hd); } else { next+=2; newram=(char *)nextof(hd); if (*next=='=') { next++; scan_dp(); if (error) return; milesseq(hd); return; } scan_dp(); if (error) return; miless(hd); } } else if (*next=='>') { next++; newram=(char *)nextof(hd); if (*next=='=') { next++; scan_dp(); if (error) return; mgreatereq(hd); return; } scan_dp(); if (error) return; mgreater(hd); } else if (*next=='=' && *(next+1)=='=') { next+=2; newram=(char *)nextof(hd); scan_dp(); if (error) return; mequal(hd); } else if (*next=='~' && *(next+1)=='=' && *(next+2)!='=') { next+=2; newram=(char *)nextof(hd); scan_dp(); if (error) return; maboutequal(hd); } else if (*next=='!' && *(next+1)=='=') { next+=2; newram=(char *)nextof(hd); scan_dp(); if (error) return; munequal(hd); } else if (*next=='_') { next++; newram=(char *)nextof(hd); scan_compare(); if (error) return; mvconcat(hd); } else if (*next=='|' && *(next+1)!='|') { next++; newram=(char *)nextof(hd); scan_compare(); if (error) return; mhconcat(hd); } } void scan_andlogical (void) { header *hd=(header *)newram; scan_compare(); if (error) return; scan_space(); if (*next=='&' && *(next+1)=='&') { next+=2; newram=(char *)nextof(hd); scan_andlogical(); if (error) return; mand(hd); } } void scan_logical (void) { header *hd=(header *)newram; scan_andlogical(); if (error) return; scan_space(); if (*next=='|' && *(next+1)=='|') { next+=2; newram=(char *)nextof(hd); scan_logical(); if (error) return; mor(hd); } } header *scan (void) { header *result=(header *)newram; scan_space(); if (*next=='{') { next++; scan_logical(); if (error) return result; while (1) { scan_space(); if (*next=='}') { next++; return result; } if (*next!=',') { output("Error in {}!\n"); error=104; return result; } next++; scan_logical(); if (error) return result; } } else { scan_logical(); } return result; } header *scan_value (void) { header *result=(header *)newram,*hd,*hd1,*marker,*nextresult, *endresults; int oldnosubmref; unsigned long size,hd1size; scan_space(); if (*next=='{') { next++; oldnosubmref=nosubmref; nosubmref=1; scan_logical(); nosubmref=oldnosubmref; hd1=getvalue(result); if (error) return result; moveresult(result,hd1); while (1) { scan_space(); if (*next=='}') { next++; return result; } if (*next!=',') { output("Error in {}!\n"); error=104; return result; } next++; hd=(header *)newram; scan_value(); if (error) return result; hd1=getvalue(hd); if (error) return result; moveresult(hd,hd1); } } else { scan_logical(); marker=result; endresults=(header *)newram; while (marker=nextresult && (char *)hd1size); newram=(char *)nextof(marker); break; } } if (hd1>nextresult && (char *)hd1size; size=hd1size-marker->size; memmove((char *)nextresult+size,(char *)nextresult, newram-(char *)nextresult); if (hd1>nextresult && (char *)hd1