#include <stdlib.h>
/* lots of changes - good luck */
/* NOTE!!! I am changing the parser a great deal here
I am making parameters 200-300
kernels are 1000
variables are 1100 - 2100 for now
shift/delay adds 3000 to them
currently 600 is the maximum to get it up to eg 700
delete all .o files change xpplim.h by adding 100 to MAXODE,MAXODE1
add 50 to MAXPRIMEVAR
better make worksize bigger - I am not sure how much cant hurt to
be too big - 300000 = 2.4MB memory which is not too much these days
change is_uvar - change 17 to 18
is_var - " "
The following are the currently used indices
0-99 functions of 1 variable like sin, cos, etc
100-199 fixed functions of 2 variables like + -, etc
200-399 parameters
600-699 networks
700-799 lookup tables
800-899 dummy arguments for functions
900-999 special stuff
1000-1099 volterra kernels
10000-11000 variables
2400-2499 user defined functions
3200-3300 shifted constants
20000-21000 shifted variables
NOTE - I will shortly make user-defined functions 2400 instead of 400
so that I have the block from 200-599 for parameter
*/
#ifndef WCTYPE
#include <ctype.h>
#else
#include <wctype.h>
#endif
#include <math.h>
/* #include <malloc.h> */
#include <stdio.h>
#include <string.h>
#include "parserslow.h"
#include "xpplim.h"
#include "getvar.h"
#define MAXEXPLEN 1024
#define THOUS 10000
#define DOUB_EPS 2.23E-15
#define POP stack[--stack_pointer]
double zippy;
#define PUSH(a) zippy=(a); stack[stack_pointer++]=zippy;
#define DFNORMAL 1
#define DFFP 2
#define DFSTAB 3
#define COM(a) my_symb[toklist[(a)]].com
int ERROUT;
extern int DelayFlag;
int NDELAYS=0;
/*double pow2(); */
double get_delay();
double delay_stab_eval();
double lookup(),network_value();
double atof(),poidev();
double ndrand48();
double ker_val();
double hom_bcs();
double BoxMuller;
int BoxMullerFlag=0;
int RandSeed=12345678;
extern int newseed;
extern int del_stab_flag;
double CurrentIndex=0;
int SumIndex=1;
/*************************
RPN COMPILER *
**************************/
/*****************************
* PARSER.C *
*
*
* parses any algebraic expression
* and converts to an integer array
* to be interpreted by the rpe_val
* function.
*
* the main data structure is a contiguous
* list of symbols with their priorities
* and their symbol value
*
* on the first pass, the expression is converted to
* a list of integers without any checking except for
* valid symbols and for numbers
*
* next this list of integers is converted to an RPN expression
* for evaluation.
*
*
* 6/95 stuff added to add names to namelist without compilation
*************************************************************/
/* INIT_RPN */
/*
int matherr(e)
struct exception *e;
{
char *c;
switch (e->type) {
case DOMAIN: c = "domain error"; break;
case SING : c = "argument singularity"; break;
case OVERFLOW: c = "overflow range"; break;
case UNDERFLOW: c = "underflow range"; break;
default: c = "(unknown error"; break;
}
fprintf(stderr, "math exception : %d\n", e->type);
fprintf(stderr, " name : %s\n", e->name);
fprintf(stderr, " arg 1: %e\n", e->arg1);
fprintf(stderr, " arg 2: %e\n", e->arg2);
fprintf(stderr, " ret : %e\n", e->retval);
return 1;
}
*/
init_rpn()
{
int i;
ERROUT = 1;
NCON = 0;
NFUN = 0;
NVAR = 0;
NKernel=0;
MaxPoints=4000;
NSYM = STDSYM;
two_args();
one_arg();
add_con("PI", M_PI);
add_con("I'",0.0);
/* add_con("c___1",0.0);
add_con("c___2",0.0);
add_con("c___3",0.0); */
SumIndex=NCON-1;
init_table();
if (newseed==1) RandSeed=time(0);
nsrand48(RandSeed);
}
/* FREE_UFUNS */
free_ufuns()
{
int i;
for(i=0;i<NFUN;i++)
{
free(ufun[i]);
free(ufun_def[i]);
}
}
duplicate_name(junk)
char *junk;
{
int i;
find_name(junk,&i);
if(i>=0){
if(ERROUT)printf("%s is a duplicate name\n",junk);
return(1);
}
return(0);
}
/* ADD_CONSTANT */
add_constant(junk)
char *junk;
{
int len;
char string[100];
if(duplicate_name(junk)==1)return(1);
if(NCON>=397)
{
if(ERROUT)printf("too many constants !!\n");
return(1);
}
convert(junk,string);
len=strlen(string);
if(len<1){
printf("Empty parameter - remove spaces\n");
return 1;
}
if(len>MXLEN)len=MXLEN;
strncpy(my_symb[NSYM].name,string,len);
my_symb[NSYM].name[len]='\0';
my_symb[NSYM].len=len;
my_symb[NSYM].pri=10;
my_symb[NSYM].arg=0;
my_symb[NSYM].com=200+NCON-1;
NSYM++;
return(0);
}
get_var_index(name)
char *name;
{
int in=-1;
int type,com;
find_name(name,&type);
if(type<0)return -1;
com=my_symb[type].com;
if(is_uvar(com))
{
return(com-10000);
}
return(-1);
}
/* GET_TYPE */
get_type(index)
int index;
{
return(my_symb[index].com);
}
/* ADD_CON */
add_con(name,value)
char *name;
double value;
{
if(NCON>=397)
{
if(ERROUT)printf("too many constants !!\n");
return(1);
}
constants[NCON]=value;
NCON++;
return(add_constant(name));
}
add_kernel(name,mu,expr)
char *name,*expr;
double mu;
{
char string[100];
char bexpr[MAXEXPLEN],kexpr[MAXEXPLEN];
int len,i,in=-1;
if(duplicate_name(name)==1)return(1);
if(NKernel==MAXKER){
printf("Too many kernels..\n");
return(1);
}
if(mu<0||mu>=1.0){
printf(" mu must lie in [0,1.0) \n");
return(1);
}
convert(name,string);
len=strlen(string);
if(len>MXLEN)len=MXLEN;
strncpy(my_symb[NSYM].name,string,len);
my_symb[NSYM].name[len]='\0';
my_symb[NSYM].len=len;
my_symb[NSYM].pri=10;
my_symb[NSYM].arg=0;
my_symb[NSYM].com=1000+NKernel;
kernel[NKernel].k_n1=0.0;
kernel[NKernel].mu=mu;
kernel[NKernel].k_n=0.0;
kernel[NKernel].k_n1=0.0;
kernel[NKernel].flag=0;
for(i=0;i<strlen(expr);i++)
if(expr[i]=='#')in=i;
if(in==0||in==(strlen(expr)-1)){
printf("Illegal use of convolution...\n");
return(1);
}
if(in>0){
kernel[NKernel].flag=CONV;
kernel[NKernel].expr=(char *)malloc(strlen(expr)+2-in);
kernel[NKernel].kerexpr=(char *)malloc(in+1);
for(i=0;i<in;i++)kernel[NKernel].kerexpr[i]=expr[i];
kernel[NKernel].kerexpr[in]=0;
for(i=in+1;i<strlen(expr);i++)kernel[NKernel].expr[i-in-1]=expr[i];
kernel[NKernel].expr[strlen(expr)-in-1]=0;
printf("Convolving %s with %s\n",
kernel[NKernel].kerexpr,kernel[NKernel].expr);
}
else {
kernel[NKernel].expr=(char *)malloc(strlen(expr)+2);
strcpy(kernel[NKernel].expr,expr);
}
NSYM++;
NKernel++;
return(0);
}
/* ADD_VAR */
add_var(junk, value)
char *junk;
double value;
{
char string[100];
int len;
/* printf(" variable - %s \n",junk); */
if(duplicate_name(junk)==1)return(1);
if(NVAR>=MAXODE1)
{
if(ERROUT)printf("too many variables !!\n");
return(1);
}
convert(junk,string);
len=strlen(string);
if(len>MXLEN)len=MXLEN;
strncpy(my_symb[NSYM].name,string,len);
my_symb[NSYM].name[len]='\0';
my_symb[NSYM].len=len;
my_symb[NSYM].pri=10;
my_symb[NSYM].arg=0;
my_symb[NSYM].com=10000+NVAR;
NSYM++;
variables[NVAR]=value;
NVAR++;
return(0);
}
/* ADD_EXPR */
add_expr(expr,command, length )
char *expr;
int *command, *length;
{
char dest[1024];
int my_token[1024];
int err,i;
convert(expr,dest);
/* printf(" Making token ...\n"); */
err=make_toks(dest,my_token);
/* i=0;
while(1){
printf(" %d %d \n",i,my_token[i]);
if(my_token[i]==ENDTOK)break;
i++;
} */
if(err!=0)return(1);
err = alg_to_rpn(my_token,command);
if(err!=0)return(1);
i=0;
while(command[i]!=ENDEXP)i++;
*length=i+1;
/* for(i=0;i<*length;i++)printf("%d \n",command[i]); */
return(0);
}
add_net_name(index,name)
int index;
char *name;
{
char string[50];
int len=strlen(name);
printf(" Adding net %s %d \n",name,index);
if(duplicate_name(name)==1)return(1);
convert(name,string);
if(len>MXLEN)len=MXLEN;
strncpy(my_symb[NSYM].name,string,len);
my_symb[NSYM].name[len]='\0';
my_symb[NSYM].len=len;
my_symb[NSYM].pri=10;
my_symb[NSYM].arg=1;
my_symb[NSYM].com=600+index;
NSYM++;
return(0);
}
add_vect_name(index,name)
int index;
char *name;
{
char string[50];
int len=strlen(name);
printf(" Adding vector %s %d \n",name,index);
if(duplicate_name(name)==1)return(1);
convert(name,string);
if(len>MXLEN)len=MXLEN;
strncpy(my_symb[NSYM].name,string,len);
my_symb[NSYM].name[len]='\0';
my_symb[NSYM].len=len;
my_symb[NSYM].pri=10;
my_symb[NSYM].arg=1;
my_symb[NSYM].com=VECT_ROOT+index;
NSYM++;
return(0);
}
/* ADD LOOKUP TABLE */
add_2d_table(name,file)
char *name,*file;
{
printf(" TWO D NOT HERE YET \n");
return(1);
}
add_file_table(index,file)
char *file;
int index;
{
if(load_table(file,index)==0)
{
if(ERROUT)printf("Problem with creating table !!\n");
return(1);
}
return(0);
}
add_table_name(index,name)
char *name;
int index;
{
char string[50];
int len=strlen(name);
if(duplicate_name(name)==1)return(1);
convert(name,string);
if(len>MXLEN)len=MXLEN;
strncpy(my_symb[NSYM].name,string,len);
my_symb[NSYM].name[len]='\0';
my_symb[NSYM].len=len;
my_symb[NSYM].pri=10;
my_symb[NSYM].arg=1;
my_symb[NSYM].com=700+index;
set_table_name(name,index);
NSYM++;
return(0);
}
/* ADD LOOKUP TABLE */
add_form_table(index,nn,xlo,xhi,formula)
char *formula;
double xlo,xhi;
int nn;
int index;
{
if(create_fun_table(nn,xlo,xhi,formula,index)==0)
{
if(ERROUT)printf("Problem with creating table !!\n");
return(1);
}
return(0);
}
set_old_arg_names(narg)
int narg;
{
int i;
for(i=0;i<narg;i++){
sprintf(my_symb[FIRST_ARG+i].name,"ARG%d",i+1);
my_symb[FIRST_ARG+i].len=4;
}
}
set_new_arg_names(narg,args)
char args[10][11];
int narg;
{
int i;
for(i=0;i<narg;i++){
strcpy(my_symb[FIRST_ARG+i].name,args[i]);
my_symb[FIRST_ARG+i].len=strlen(args[i]);
}
}
/* NEW ADD_FUN for new form_ode code */
add_ufun_name(name,index,narg)
char *name;
int index,narg;
{
char string[50];
int len=strlen(name);
if(duplicate_name(name)==1)return(1);
if(index>=MAXUFUN)
{
if(ERROUT)printf("too many functions !!\n");
return(1);
}
printf(" Added user fun %s \n",name);
convert(name,string);
if(len>MXLEN)len=MXLEN;
strncpy(my_symb[NSYM].name,string,len);
my_symb[NSYM].name[len]='\0';
my_symb[NSYM].len=len;
my_symb[NSYM].pri=10;
my_symb[NSYM].arg=narg;
my_symb[NSYM].com=2400+index;
NSYM++;
strcpy(ufun_names[index],name);
return (0);
}
fixup_endfun(u,l,narg)
int *u;
int l,narg;
{
u[l-1]=ENDFUN;
u[l]=narg;
u[l+1]=ENDEXP;
}
add_ufun_new(index,narg,rhs,args)
char *rhs,args[MAXARG][11];
int narg,index;
{
int i,l;
int end;
if(narg>MAXARG){
printf("Maximal arguments exceeded \n");
return(1);
}
if((ufun[index]=(int *)malloc(1024))==NULL)
{
if(ERROUT)printf("not enough memory!!\n");
return(1);
}
if((ufun_def[index]=(char *)malloc(MAXEXPLEN))==NULL)
{
if(ERROUT)printf("not enough memory!!\n");
return(1);
}
ufun_arg[index].narg=narg;
for(i=0;i<narg;i++)
strcpy(ufun_arg[index].args[i],args[i]);
set_new_arg_names(narg,args);
if(add_expr(rhs,ufun[index],&end)==0)
{
ufun[index][end-1]=ENDFUN;
ufun[index][end]=narg;
ufun[index][end+1]=ENDEXP;
strcpy(ufun_def[index],rhs);
l=strlen(ufun_def[index]);
ufun_def[index][l]=0;
narg_fun[index]=narg;
set_old_arg_names(narg);
return(0);
}
set_old_arg_names(narg);
if(ERROUT)printf(" ERROR IN FUNCTION DEFINITION\n");
return(1);
}
/* ADD_UFUN */
add_ufun(junk,expr,narg)
char *junk, *expr;
int narg;
{
char string[50];
int i,l;
int end;
int len=strlen(junk);
if(duplicate_name(junk)==1)return(1);
if(NFUN>=MAXUFUN)
{
if(ERROUT)printf("too many functions !!\n");
return(1);
}
if((ufun[NFUN]=(int *)malloc(1024))==NULL)
{
if(ERROUT)printf("not enough memory!!\n");
return(1);
}
if((ufun_def[NFUN]=(char *)malloc(MAXEXPLEN))==NULL)
{
if(ERROUT)printf("not enough memory!!\n");
return(1);
}
convert(junk,string);
if(add_expr(expr,ufun[NFUN],&end)==0)
{
if(len>MXLEN)len=MXLEN;
strncpy(my_symb[NSYM].name,string,len);
my_symb[NSYM].name[len]='\0';
my_symb[NSYM].len=len;
my_symb[NSYM].pri=10;
my_symb[NSYM].arg=narg;
my_symb[NSYM].com=2400+NFUN;
NSYM++;
ufun[NFUN][end-1]=ENDFUN;
ufun[NFUN][end]=narg;
ufun[NFUN][end+1]=ENDEXP;
strcpy(ufun_def[NFUN],expr);
l=strlen(ufun_def[NFUN]);
ufun_def[NFUN][l-1]=0;
strcpy(ufun_names[NFUN],junk);
narg_fun[NFUN]=narg;
for(i=0;i<narg;i++){
sprintf(ufun_arg[NFUN].args[i],"ARG%d",i+1);
}
NFUN++;
return(0);
}
if(ERROUT)printf(" ERROR IN FUNCTION DEFINITION\n");
return(1);
}
check_num(tok,value)
int *tok;
double value;
{
int bob,in,m,i;
for(i=0;i<NSYM;i++){
if(strncmp(my_symb[i].name,"NUM##",5)==0){
bob=my_symb[i].com;
in=bob%100;
m=bob/100;
if(m==10)in+=100;
if(constants[in]==value){
*tok=i;
return(1);
}
}
}
return(0);
}
/* is_ufun */
is_ufun( x)
int x;
{
if((x/100)==UFUN)return(1);
else return(0);
}
/* IS_UCON */
is_ucon( x)
int x;
{
int xx=x/100;
if(xx>1&&xx<6)return(1);
else return(0);
}
/* IS_UVAR */
is_uvar(x)
int x;
{
int y=x/100;
if((y>=100)&&(y<200))return(1);
else return(0);
}
isvar(y)
int y;
{
if((y>=100)&&(y<200))return(1);
return 0;
}
iscnst(y)
int y;
{
if(y>1&&y<6)return(1);
return 0;
}
isker(y)
int y;
{
if(y==10)return(1);
return 0;
}
is_kernel(x)
int x;
{
if((x/100)==10)return(1);
else return(0);
}
is_lookup(x)
int x;
{
if((x/100)==7)return(1);
else return(0);
}
find_lookup(name)
char *name;
{
int index,com;
find_name(name,&index);
if(index==-1)return(-1);
com=my_symb[index].com;
if(is_lookup(com))return(com%100);
return(-1);
}
/* FIND_NAME */
find_name(string, index)
char *string;
int *index;
{
char junk[100];
int i,len;
convert(string,junk);
len=strlen(junk);
for(i=0;i<NSYM;i++)
{
if(len==my_symb[i].len)
if(strncmp(my_symb[i].name,junk,len)==0)break;
}
if(i<NSYM)
*index=i;
else *index=-1;
}
get_param_index(name)
char *name;
{
int type,com;
find_name(name,&type);
if(type<0)return(-1);
com=my_symb[type].com;
if(is_ucon(com))
{
return(com-200);
}
return(-1);
}
/* GET_VAL */
get_val(name,value)
char *name;
double *value;
{
int type,com;
*value=0.0;
find_name(name,&type);
if(type<0)return(0);
com=my_symb[type].com;
if(is_ucon(com))
{
*value=constants[com-200];
return(1);
}
if(is_uvar(com))
{
*value=variables[com-10000];
return(1);
}
return(0);
}
/* SET_VAL */
set_val(name, value)
char *name;
double value;
{
int type,com;
find_name(name,&type);
if(type<0)return(0);
com=my_symb[type].com;
if(is_ucon(com))
{
constants[com-200]=value;
return(1);
}
if(is_uvar(com))
{
variables[com-10000]=value;
return(1);
}
return(0);
}
set_ivar(i, value)
int i;
double value;
{
SETVAR(i,value);
}
double get_ivar(i)
int i;
{ return(GETVAR(i));
}
alg_to_rpn(toklist,command)
int *toklist,*command;
{
int tokstak[500],comptr=0,tokptr=0,lstptr=0,temp;
int ncomma=0,delflag=0,deltok;
int loopstk[100];
int lptr=0;
int nif=0,nthen=0,nelse=0;
int newtok,oldtok,zip;
int i,my_com,my_arg,jmp;
char ch;
int jj;
tokstak[0]=STARTTOK;
tokptr=1;
oldtok=STARTTOK;
while(1)
{
getnew:
newtok=toklist[lstptr++];
/* for(zip=0;zip<tokptr;zip++)
printf("%d %d\n",zip,tokstak[zip]); */
/* check for delay symbol */
if(newtok==DELSYM)
{
temp=my_symb[toklist[lstptr+1]].com;
/* !! */ if(is_uvar(temp))
{
my_symb[LASTTOK].com=temp+10000; /* create a temporary sybol */
NDELAYS++;
toklist[lstptr+1]=LASTTOK;
my_symb[LASTTOK].pri=10;
}
else
{
printf("Illegal use of DELAY \n");
return(1);
}
}
/* check for delshft symbol */
if(newtok==DELSHFTSYM)
{
temp=my_symb[toklist[lstptr+1]].com;
/* !! */ if(is_uvar(temp))
{
my_symb[LASTTOK].com=temp+10000; /* create a temporary sybol */
NDELAYS++;
toklist[lstptr+1]=LASTTOK;
my_symb[LASTTOK].pri=10;
}
else
{
printf("Illegal use of DELAY Shift \n");
return(1);
}
}
/* check for shift */
if(newtok==SHIFTSYM||newtok==ISHIFTSYM)
{
temp=my_symb[toklist[lstptr+1]].com;
/* !! */ if(is_uvar(temp) || is_ucon(temp))
{
if(is_uvar(temp))my_symb[LASTTOK].com=temp+10000;
if(is_ucon(temp))my_symb[LASTTOK].com=temp+3000;
/* create a temporary sybol */
toklist[lstptr+1]=LASTTOK;
my_symb[LASTTOK].pri=10;
}
else
{
printf("Illegal use of SHIFT \n");
return(1);
}
}
next:
if((newtok==ENDTOK)&&(oldtok==STARTTOK))break;
if(newtok==LPAREN)
{
tokstak[tokptr]=LPAREN;
tokptr++;
oldtok=LPAREN;
goto getnew;
}
if(newtok==RPAREN)
{
switch(oldtok)
{
case LPAREN:
tokptr--;
oldtok=tokstak[tokptr-1];
goto getnew;
case COMMA:
tokptr--;
ncomma++;
oldtok=tokstak[tokptr-1];
goto next;
}
}
if((newtok==COMMA)&&(oldtok==COMMA))
{
tokstak[tokptr]=COMMA;
tokptr++;
goto getnew;
}
if(my_symb[oldtok%THOUS].pri>=my_symb[newtok%THOUS].pri)
{
command[comptr]=my_symb[oldtok%THOUS].com;
/* printf("com(%d)=%d\n",comptr,command[comptr]); */
if((my_symb[oldtok%THOUS].arg==2)&&
(my_symb[oldtok%THOUS].com/100==1))
ncomma--;
my_com=command[comptr];
comptr++;
/* New code 3/95 */
if(my_com==NUMSYM){
/* printf("tp=%d ",tokptr); */
tokptr--;
/* printf(" ts[%d]=%d ",tokptr,tokstak[tokptr]); */
command[comptr]=tokstak[tokptr-1];
/* printf("xcom(%d)=%d\n",comptr,command[comptr]); */
comptr++;
tokptr--;
command[comptr]=tokstak[tokptr-1];
/* printf("xcom(%d)=%d\n",comptr,command[comptr]); */
comptr++;
}
/* end new code 3/95 */
if(my_com==SUMSYM){
loopstk[lptr]=comptr;
comptr++;
lptr++;
ncomma-=1;
}
if(my_com==ENDSUM){
lptr--;
jmp=comptr-loopstk[lptr]-1;
command[loopstk[lptr]]=jmp;
}
if(my_com==MYIF){
loopstk[lptr]=comptr; /* add some space for jump */
comptr++;
lptr++;
nif++;
}
if(my_com==MYTHEN){
/* First resolve the if jump */
lptr--;
jmp=comptr-loopstk[lptr]; /* -1 is old */
command[loopstk[lptr]]=jmp;
/* Then set up for the then jump */
loopstk[lptr]=comptr;
lptr++;
comptr++;
nthen++;
}
if(my_com==MYELSE){
lptr--;
jmp=comptr-loopstk[lptr]-1;
command[loopstk[lptr]]=jmp;
nelse++;
}
if(my_com==ENDDELAY||my_com==ENDSHIFT||my_com==ENDISHIFT){
ncomma-=1;
}
if(my_com==ENDDELSHFT)
ncomma-=2;
/* if(my_com==CONV||my_com==DCONV){
ncomma-=1;
} */
/* CHECK FOR USER FUNCTION */
if(is_ufun(my_com))
{
my_arg=my_symb[oldtok].arg;
command[comptr]=my_arg;
comptr++;
ncomma=ncomma+1-my_arg;
}
/* USER FUNCTION OKAY */
tokptr--;
oldtok=tokstak[tokptr-1];
goto next;
}
/* NEW code 3/95 */
if(newtok==NUMTOK){
tokstak[tokptr++]=toklist[lstptr++];
tokstak[tokptr++]=toklist[lstptr++];
}
/* end 3/95 */
tokstak[tokptr]=newtok;
oldtok=newtok;
tokptr++;
goto getnew;
}
if(ncomma!=0){
printf("Illegal number of arguments\n");
return(1);
}
if((nif!=nelse)||(nif!=nthen)){
printf("If statement missing ELSE or THEN \n");
return(1);
}
command[comptr]=my_symb[ENDTOK].com;
/* pr_command(command); */
return(0);
}
pr_command(command)
int *command;
{
int i=0;
int token;
while(1){
token=command[i];
printf("%d %d \n",i,token);
if(token==ENDEXP)return;
i++;
}
}
show_where(string,index)
char *string;
int index;
{
char junk[MAXEXPLEN];
int i;
for(i=0;i<index;i++)junk[i]=' ';
junk[index]='^';
junk[index+1]=0;
printf("%s\n%s\n",string,junk);
}
function_sym(token) /* functions should have ( after them */
{
int com=my_symb[token].com;
int i1=com/100;
if(i1==0&&!unary_sym(token))return(1); /* single variable functions */
if(i1==1&&!binary_sym(token))return(1); /* two-variable function */
if(i1==UFUN||i1==7||i1==6||i1==5)return(1);
if(token==DELSHFTSYM||token==DELSYM||token==SHIFTSYM||token==ISHIFTSYM||com==MYIF||com==MYTHEN||com==MYELSE
||com==SUMSYM||com==ENDSUM)return(1);
return(0);
}
unary_sym(token)
{
if(token==9||token==55)return(1);
return(0);
}
binary_sym(token)
int token;
{
if(token>2&&token<9)return(1);
if(token>43&&token<51)return(1);
if(token==54)return(1);
return(0);
}
pure_number(token)
int token;
{
int com=my_symb[token].com;
int i1=com/100;
/* !! */ if(token==NUMTOK||isvar(i1)||iscnst(i1)||isker(i1)||i1==8||token==INDX)
return(1);
return(0);
}
gives_number(token)
int token;
{
int com=my_symb[token].com;
int i1=com/100;
if(token==INDX)return(1);
if(token==NUMTOK)return(1);
if(i1==0&&!unary_sym(token))return(1); /* single variable functions */
if(i1==1&&!binary_sym(token))return(1); /* two-variable function */
/* !! */ if(i1==8||isvar(i1)||iscnst(i1)||i1==7||i1==6||i1==5||isker(i1)||i1==UFUN)return(1);
if(com==MYIF||token==DELSHFTSYM||token==DELSYM||token==SHIFTSYM||token==ISHIFTSYM||com==SUMSYM)return(1);
return(0);
}
check_syntax(oldtoken,newtoken) /* 1 is BAD! */
int oldtoken,newtoken;
{
int com1=my_symb[oldtoken].com,com2=my_symb[newtoken].com;
/* if the first symbol or ( or binary symbol then must be unary symbol or
something that returns a number or another (
*/
if(unary_sym(oldtoken)||oldtoken==COMMA||oldtoken==STARTTOK
||oldtoken==LPAREN||binary_sym(oldtoken))
{
if(unary_sym(newtoken)||gives_number(newtoken)||newtoken==LPAREN)return(0);
return(1);
}
/* if this is a regular function, then better have (
*/
if(function_sym(oldtoken)){
if(newtoken==LPAREN)return(0);
return(1);
}
/* if we have a constant or variable or ) or kernel then better
have binary symbol or "then" or "else" as next symbol
*/
if(pure_number(oldtoken)){
if(binary_sym(newtoken)||newtoken==RPAREN
||newtoken==COMMA||newtoken==ENDTOK)
return(0);
return(1);
}
if(oldtoken==RPAREN){
if(binary_sym(newtoken)||newtoken==RPAREN
||newtoken==COMMA||newtoken==ENDTOK)return(0);
if(com2==MYELSE||com2==MYTHEN||com2==ENDSUM)return(0);
return(1);
}
printf("Bad token %d \n",oldtoken);
return(1);
}
/******************************
* PARSER *
******************************/
make_toks(dest,my_token)
char *dest;
int *my_token;
{
char num[40],junk[10];
double value;
int old_tok=STARTTOK,tok_in=0;
int index=0,err,token,nparen=0,lastindex=0;
union /* WARNING -- ASSUMES 32 bit int and 64 bit double */
{
struct {
int int1;
int int2;
} pieces;
struct {
double z;
} num;
} encoder;
while(dest[index]!='\0')
{
lastindex=index;
find_tok(dest,&index,&token);
if((token==MINUS)&&
((old_tok==STARTTOK)||(old_tok==COMMA)||(old_tok==LPAREN)))
token=NEGATE;
if(token==LPAREN)++nparen;
if(token==RPAREN)--nparen;
if(token==NSYM)
{
if(do_num(dest,num,&value,&index)){
show_where(dest,index);
return(1);
}
/* new code 3/95 */
encoder.num.z=value;
my_token[tok_in++]=NUMTOK;
my_token[tok_in++]=encoder.pieces.int1;
my_token[tok_in++]=encoder.pieces.int2;
if(check_syntax(old_tok,NUMTOK)==1){
printf("Illegal syntax \n");
show_where(dest,lastindex);
return(1);
}
old_tok=NUMTOK;
}
else
{
my_token[tok_in++]=token;
if(check_syntax(old_tok,token)==1){
printf("Illegal syntax (Ref:%d %d) \n",old_tok,token);
show_where(dest,lastindex);
tokeninfo(old_tok);
tokeninfo(token);
return(1);
}
old_tok=token;
}
}
my_token[tok_in++]=ENDTOK;
if(check_syntax(old_tok,ENDTOK)==1){
printf("Premature end of expression \n");
show_where(dest,lastindex);
return(1);
}
if(nparen!=0)
{
if(ERROUT)printf(" parentheses don't match\n");
return(1);
}
return(0);
}
tokeninfo(tok)
int tok;
{
printf(" %s %d %d %d %d \n",
my_symb[tok].name,my_symb[tok].len,my_symb[tok].com,
my_symb[tok].arg,my_symb[tok].pri);
}
do_num(source,num,value,ind)
char *source,*num;
double *value;
int *ind;
{
int j=0,i=*ind,error=0;
int ndec=0,nexp=0,ndig=0;
char ch,oldch;
oldch='\0';
*value=0.0;
while(1)
{
ch=source[i];
if(((ch=='+')||(ch=='-'))&&(oldch!='E'))break;
if((ch=='*')||(ch=='^')||(ch=='/')||(ch==',')||(ch==')')||(ch=='\0')
|| (ch=='|') || (ch=='>') || (ch=='<') || (ch=='&')
|| (ch=='='))break;
if((ch=='E')||(ch=='.')||(ch=='+')||(ch=='-')||isdigit(ch))
{
if(isdigit(ch))ndig++;
switch(ch)
{
case 'E':
nexp++;
if((nexp==2)||(ndig==0))goto err;break;
case '.':
ndec++;
if((ndec==2)||(nexp==1))goto err;break;
}
num[j]=ch;
j++;
i++;
oldch=ch;
}
else
{
err:
num[j]=ch;
j++;
error=1;
break;
}
}
num[j]='\0';
if(error==0)*value=atof(num);
else
if(ERROUT)printf(" illegal expression: %s\n",num);
*ind=i;
return(error);
}
convert(source,dest)
char *source,*dest;
{
char ch;
int i=0,j=0;
while (1)
{
ch=source[i];
if(!isspace(ch))dest[j++]=ch;
i++;
if(ch=='\0')break;
}
strupr(dest);
}
find_tok(source,index,tok)
char *source;
int *index,*tok;
{
int i=*index,maxlen=0,symlen;
int k,j,my_tok,match;
my_tok=NSYM;
for(k=0;k<NSYM;k++)
{
symlen=my_symb[k].len;
if(symlen<=maxlen)continue;
match=1;
for(j=0;j<symlen;j++)
{
if(source[i+j]!=my_symb[k].name[j])
{
match=0;
break;
}
}
if(match!=0)
{
my_tok=k;
maxlen=symlen;
}
}
*index=*index+maxlen;
*tok=my_tok;
}
double pmod(x,y)
double x,y;
{
double z=fmod(x,y);
if(z<0)z+=y;
return(z);
}
two_args()
{
fun2[4]=atan2;
fun2[5]=pow;
fun2[6]=max;
fun2[7]=min;
/* fun2[8]=fmod; */
fun2[8]=pmod; /* This always gives an answer in [0,y) for mod(x,y) */
fun2[9]=dand;
fun2[10]=dor;
fun2[11]=dgt;
fun2[12]=dlt;
fun2[13]=deq;
fun2[14]=dge;
fun2[15]=dle;
fun2[16]=dne;
fun2[17]=normal;
fun2[18]=bessel_j;
fun2[19]=bessel_y;
}
/* These are the Bessel Functions; if you dont have them then
return some sort of dummy value or else write a program
to compute them
*/
double bessel_j(x,y)
double x,y;
{
int n=(int)x;
return(jn(n,y));
}
double bessel_y(x,y)
double x,y;
{
int n=(int)x;
return(yn(n,y));
}
/*
double pow2(z,w)
double z,w;
{
return(pow(z,w));
double sign=1.0;
if(floor(w)==w){
if(z<0.0)sign=-1.0;
if(fabs(fmod(w,2.0))<1.e-10)sign=1.0;
return(sign*pow(fabs(z),w));
}
else
return(pow(fabs(z),w));
}
*/
/*********************************************
FANCY DELAY HERE *-------------------------<<<
*********************************************/
char *com_name(com)
int com;
{
int i;
for( i=0;i<NSYM;i++)
if( my_symb[i].com == com ) break;
if( i < NSYM )
return my_symb[i].name;
else
return "";
}
double do_shift(shift,variable)
double shift,variable;
{
int it, in;
int i=(int)(variable),ish=(int)shift;
/* printf( "shifting %d (%s) by %d to %d (%s)\n",
* (int)variable, com_name((int)variable), (int)shift, i, com_name(i) );
*/
if(i<0) return(0.0);
it=i/100;
switch(it){
case 2:
case 3:
case 4:
case 5:
in=i-200;
in+=ish;
if(in>NCON)
return 0.0;
else
return constants[in];
break;
case 100:
case 101:
case 102:
case 103:
case 104:
case 105:
case 106:
case 107:
case 108:
case 109:
case 110:
case 111:
case 112:
case 113:
case 114:
case 115:
case 116:
case 117:
case 118:
case 119:
in=i-10000;
in+=ish;
if(in>MAXODE)
return 0.0;
else
return variables[in];
default:
printf("This can't happen: Invalid symbol index for SHIFT\n");
return 0.0;
}
}
double do_ishift(shift,variable)
double shift,variable;
{
int it, in;
int i=(int)(variable),ish=(int)shift;
/* printf( "shifting %d (%s) by %d to %d (%s)\n",
* (int)variable, com_name((int)variable), (int)shift, i, com_name(i) );
*/
return variable+shift;
}
double do_delay_shift(delay,shift,variable)
double delay,shift,variable;
{
int it, in;
int i=(int)(variable),ish=(int)shift;
if(i<0) return(0.0);
in=i-10000+ish;
if(in>MAXODE)
return 0.0;
if(del_stab_flag>0){
if(DelayFlag&&delay>0.0)
return(get_delay(in-1,delay));
return(variables[in]);
}
return(delay_stab_eval(delay,in));
}
double do_delay(delay,i)
double delay,i;
{
double z;
int variable, it, in;
variable = i-10000;
if(del_stab_flag>0){
if(DelayFlag&&delay>0.0)
return(get_delay(variable-1,delay));
return(variables[variable]);
}
return(delay_stab_eval(delay,(int)variable));
}
/*
double Exp(z)
double z;
{
if(z>700)return(1.01423e+304);
return(exp(z));
}
double Ln(z)
double z;
{
if(z<1e-320)return(-736.82724);
return(log(z));
}
double Log10(z)
double z;
{
if(z<1e-320)return(-320.);
return(log10(z));
}
*/
one_arg()
{
fun1[0]=sin;
fun1[1]=cos;
fun1[2]=tan;
fun1[3]=asin;
fun1[4]=acos;
fun1[5]=atan;
fun1[6]=sinh;
fun1[7]=tanh;
fun1[8]=cosh;
fun1[9]=fabs;
fun1[10]=exp;
fun1[11]=log;
fun1[12]=log10;
fun1[13]=sqrt;
fun1[14]=neg;
fun1[15]=recip;
fun1[16]=heaviside;
fun1[17]=signum;
fun1[18]=floor;
fun1[19]=rndom;
fun1[20]=dnot;
fun1[21]=erf;
fun1[22]=erfc;
fun1[23]=hom_bcs;
fun1[24]=poidev;
}
double normal(mean,std)
double mean,std;
{
double fac,r,v1,v2;
if(BoxMullerFlag==0){
do {
v1=2.0*ndrand48()-1.0;
v2=2.0*ndrand48()-1.0;
r=v1*v1+v2*v2;
} while(r>=1.0);
fac=sqrt(-2.0*log(r)/r);
BoxMuller=v1*fac;
BoxMullerFlag=1;
return(v2*fac*std+mean);
}
else {
BoxMullerFlag=0;
return(BoxMuller*std+mean);
}
}
double max(x,y)
double x,y;
{
return(((x>y)?x:y));
}
double min(x,y)
double x,y;
{
return(((x<y)?x:y));
}
double neg(z)
double z;
{
return(-z);
}
double recip(z)
double z;
{
return(1.00/z);
}
double heaviside(z)
double z;
{
float w=1.0;
if(z<0)w=0.0;
return(w);
}
double rndom( z)
double z;
{
/* return (z*(double)rand()/32767.00); */
return(z*ndrand48());
}
double signum(z)
double z;
{
if(z<0.0)return(-1.0);
if(z>0.0)return(1.0);
return(0.0);
}
/* logical stuff */
double dnot(x)
double x;
{
return((double)(x==0.0));
}
double dand(x,y)
double x,y;
{
return((double)(x&&y));
}
double dor(x,y)
double x,y;
{
return((double)(x||y));
}
double dge(x,y)
double x,y;
{
return((double)(x>=y));
}
double dle(x,y)
double x,y;
{
return((double)(x<=y));
}
double deq(x,y)
double x,y;
{
return((double)(x==y));
}
double dne(x,y)
double x,y;
{
return((double)(x!=y));
}
double dgt(x,y)
double x,y;
{
return((double)(x>y));
}
double dlt(x,y)
double x,y;
{
return((double)(x<y));
}
/* end of logical stuff */
double evaluate(equat)
int *equat;
{
uptr=0;
stack_pointer=0;
return(eval_rpn(equat));
}
double eval_rpn(equat)
int *equat;
{
int i,it,in,j,*tmpeq;
int sumcount,ii,new[500],is;
int isum,aind;
int low,high,ind,iloop,ijmp;
double temp,temx,temy,temz;
double sum;
union /* WARNING -- ASSUMES 32 bit int and 64 bit double */
{
struct {
int int1;
int int2;
} pieces;
struct {
double z;
} num;
} encoder;
while((i=*equat++)!=ENDEXP)
{
switch(i)
{
case NUMSYM:
encoder.pieces.int2=*equat++;
encoder.pieces.int1=*equat++;
PUSH(encoder.num.z);
break;
case ENDFUN:
i=*equat++;
uptr-=i;
break;
case MYIF:
temx=POP;
ijmp=*equat++;
if(temx==0.0)equat+=ijmp;
break;
case MYTHEN:
ijmp=*equat++;
equat+=ijmp;
break;
case MYELSE:
break;
case ENDDELSHFT:
temx=POP;
temy=POP;
temz=POP;
PUSH(do_delay_shift(temx,temy,temz));
break;
case ENDDELAY:
temx=POP;
temy=POP;
PUSH(do_delay(temx,temy));
break;
case ENDSHIFT:
temx=POP;
temy=POP;
PUSH(do_shift(temx,temy));
break;
case ENDISHIFT:
temx=POP;
temy=POP;
PUSH(do_ishift(temx,temy));
break;
case SUMSYM:
temx=POP;
high=(int)temx;
temx=POP;
low=(int)temx;
ijmp=*equat++;
sum=0.0;
if(low<=high){
for(is=low;is<=high;is++){
tmpeq=equat;
constants[SumIndex]=(double)is;
sum+=eval_rpn(tmpeq);
}
}
equat+=ijmp;
PUSH(sum);
break;
case ENDSUM:
return(POP);
case INDXCOM:
PUSH(CurrentIndex);
break;
default:
{
it=i/100;
in=i%100;
switch(it)
{
case 0: PUSH(fun1[in](POP));
break;
case 1:
{
if(in==0){temx=POP;temy=POP;PUSH(temx+temy);goto bye;}
if(in==2){temx=POP;temy=POP;PUSH(temx*temy);goto bye;}
if(in==1){temx=POP;temy=POP;PUSH(temy-temx);goto bye;}
if(in==3){temx=POP;if(temx==0.0)temx=DOUB_EPS;
temy=POP;PUSH(temy/temx);goto bye;}
temx=POP;
temy=POP;
PUSH(fun2[in](temy,temx));break;
}
case 2:
case 3:
case 4:
case 5:
PUSH(constants[i-200]);break;
case 6: PUSH(network_value(POP,in));break;
case 7: PUSH(lookup(POP,in));break;
case 8:
PUSH(ustack[uptr-1-in]); break;
case 10: PUSH(ker_val(in));break;
case 100:
case 101:
case 102:
case 103:
case 104:
case 105:
case 106:
case 107:
case 108:
case 109:
case 110:
case 111:
case 112:
case 113:
case 114:
case 115:
case 116:
case 117:
case 118:
case 119:
PUSH(variables[i-10000]); break;
/* indexes for shift and delay operators... */
case 32:
case 33:
case 34:
case 35:
PUSH((double)(i-3000)); break;
case 200:
case 201:
case 202:
case 203:
case 204:
case 205:
case 206:
case 207:
case 208:
case 209:
case 210:
case 211:
case 212:
case 213:
case 214:
case 215:
case 216:
case 217:
case 218:
case 219:
PUSH((double)(i-10000)); break;
case UFUN: i=*equat++;
for(j=0;j<i;j++)
{
ustack[uptr]=POP;
uptr++;
}
PUSH(eval_rpn(ufun[in]));
break;
}
bye: j=0;
}
}
}
return(POP);
}
/* STRING STUFF */
#ifndef STRUPR
strupr(s)
char *s;
{
int i=0;
while(s[i])
{
if(islower(s[i]))s[i]-=32;
i++;
}
}
strlwr(s)
char *s;
{
int i=0;
while(s[i])
{
if(isupper(s[i]))s[i]+=32;
i++;
}
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1