#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <limits.h>
#include <ctype.h>

#include "sysdep.h"
#include "stack.h"
#include "output.h"
#include "builtin.h"
#include "funcs.h"
#include "linear.h"
#include "polynom.h"
#include "interval.h"
#include "spread.h"
#include "express.h"
#include "udf.h"
#include "mainloop.h"
#include "edit.h"


/***************************************************************************
 *	basic math functions (+, -, *, /) for real and complex numbers
 *
 ***************************************************************************/
static void real_add (double *x, double *y, double *z)
{	*z=*x+*y;
}

static void complex_add (double *x, double *xi, double *y, double *yi,
	double *z, double *zi)
{	*z=*x+*y;
	*zi=*xi+*yi;
}

void add (header *hd, header *hd1)
/***** add
	add the values elementwise.
*****/
{   spreadf2(real_add,complex_add,interval_add,hd);
}

static void real_subtract (double *x, double *y, double *z)
{	*z=*x-*y;
}

static void complex_subtract (double *x, double *xi, double *y, double *yi,
	double *z, double *zi)
{	*z=*x-*y;
	*zi=*xi-*yi;
}

void subtract (header *hd, header *hd1)
/***** substract
	substract the values elementwise.
*****/
{   spreadf2(real_subtract,complex_subtract,interval_sub,hd);
}

static void real_multiply (double *x, double *y, double *z)
{	*z=*x*(*y);
}

void complex_multiply (double *x, double *xi, double *y, double *yi,
	double *z, double *zi)
{	double h=*x * *y - *xi * *yi;
	*zi=*x * *yi + *xi * *y;
	*z=h;
}

void dotmultiply (header *hd, header *hd1)
/***** dotmultiply
	multiply the values elementwise.
*****/
{   spreadf2(real_multiply,complex_multiply,interval_mult,hd);
}

static void real_divide (double *x, double *y, double *z)
{
#ifdef FLOAT_TEST
	if (*y==0.0) { *y=1e-10; error=1; }
#endif
	*z=*x/(*y);
}

void complex_divide (double *x, double *xi, double *y, double *yi,
	double *z, double *zi)
{	double r,h;
	r=*y * *y + *yi * *yi;
#ifdef FLOAT_TEST
	if (r==0) { r=1e-10; error=1; }
#endif
	h=(*x * *y + *xi * *yi)/r;
	*zi=(*xi * *y - *x * *yi)/r;
	*z=h;
}

void dotdivide (header *hd, header *hd1)
/***** dotdivide
	divide the values elementwise.
*****/
{   spreadf2(real_divide,complex_divide,interval_div,hd);
}

void copy_complex (double *x, double *y)
{	*x++=*y++;
	*x=*y;
}

/***************************************************************************
 *	math functions (exp, log, sqr, sin, cos, tan, asin, acos, atan)
 *
 ***************************************************************************/

static void csin (double *x, double *xi, double *z, double *zi)
{	*z=cosh(*xi)*sin(*x);
	*zi=sinh(*xi)*cos(*x);
}

void msin (header *hd)
{	spread1i(sin,csin,isin,hd);
	test_error("sin");
}

static void ccos (double *x, double *xi, double *z, double *zi)
{	*z=cosh(*xi)*cos(*x);
	*zi=-sinh(*xi)*sin(*x);
}

void mcos (header *hd)
{	spread1i(cos,ccos,icos,hd);
	test_error("cos");
}

static void ctan (double *x, double *xi, double *z, double *zi)
{	double s,si,c,ci;
	csin(x,xi,&s,&si); ccos(x,xi,&c,&ci);
	complex_divide(&s,&si,&c,&ci,z,zi);
}

#ifdef FLOAT_TEST
static double rtan (double x)
{	if (cos(x)==0.0) return 1e10;
	return tan(x);
}
#endif

void mtan (header *hd)
{	spread1i(
#ifdef FLOAT_TEST
	rtan,
#else
	tan,
#endif
	ctan,itan,hd);
	test_error("tan");
}

#ifdef FLOAT_TEST
static double ratan (double x)
{	if (x<=-M_PI && x>=M_PI) return 1e10;
	else return atan(x);
}
#endif

static void carg (double *x, double *xi, double *z)
{
#ifdef FLOAT_TEST
	if (*x==0.0 && *xi==0.0) *z=0.0;
#endif
	*z = atan2(*xi,*x);
}

#ifdef FLOAT_TEST
static double rlog (double x)
{	if (x<=0) { error=1; return 0; }
	else return log(x);
}
#endif

static void cclog (double *x, double *xi, double *z, double *zi)
{
#ifdef FLOAT_TEST
	*z=rlog(sqrt(*x * *x + *xi * *xi));
#else
	*z=log(sqrt(*x * *x + *xi * *xi));
#endif
	carg(x,xi,zi);
}


static void catan (double *x, double *xi, double *y, double *yi)
{	double h,hi,g,gi,t,ti;
	h=1-*xi; hi=*x; g=1+*xi; gi=-*x;
	complex_divide(&h,&hi,&g,&gi,&t,&ti);
	cclog(&t,&ti,&h,&hi);
	*y=hi/2; *yi=-h/2;
}

void matan (header *hd)
{	spread1i(
#ifdef FLOAT_TEST
	ratan,
#else
	atan,
#endif
	catan,iatan,hd);
	test_error("atan");
}

#ifdef FLOAT_TEST
static double rasin (double x)
{	if (x<-1 || x>1) { error=1; return 0; }
	else return asin(x);
}
#endif

static void csqrt (double *x, double *xi, double *z, double *zi)
{	double a,r;
	carg(x,xi,&a); a=a/2.0;
	r=sqrt(sqrt(*x * *x + *xi * *xi));
	*z=r*cos(a);
	*zi=r*sin(a);
}

static void casin (double *x, double *xi, double *y, double *yi)
{	double h,hi,g,gi;
	complex_multiply(x,xi,x,xi,&h,&hi);
	h=1-h; hi=-hi;
	csqrt(&h,&hi,&g,&gi);
	h=-*xi+g; hi=*x+gi;
	cclog(&h,&hi,yi,y);
	*yi=-*yi;
}

void masin (header *hd)
{	spread1(
#ifdef FLOAT_TEST
	rasin,
#else
	asin,
#endif
	casin,hd);
	test_error("asin");
}

#ifdef FLOAT_TEST
static double racos (double x)
{	if (x<-1 || x>1) { error=1; return 0; }
	else return acos(x);
}
#endif

static void cacos (double *x, double *xi, double *y, double *yi)
{	double h,hi,g,gi;
	complex_multiply(x,xi,x,xi,&h,&hi);
	h=1-h; hi=-hi;
	csqrt(&h,&hi,&g,&gi);
	hi=*xi+g; h=*x-gi;
	cclog(&h,&hi,yi,y);
	*yi=-*yi;
}

void macos (header *hd)
{	spread1(
#ifdef FLOAT_TEST
	racos,
#else
	acos,
#endif
	cacos,hd);
	test_error("acos");
}

static void cexp (double *x, double *xi, double *z, double *zi)
{	double r=exp(*x);
	*z=cos(*xi)*r;
	*zi=sin(*xi)*r;
}

#ifdef FLOAT_TEST
static void rcexp (double *x, double *xi, double *z, double *zi)
{	double r;
    if (*x>=1000) { error=1; return; }
    r=exp(*x);
	*z=cos(*xi)*r;
	*zi=sin(*xi)*r;
}
#endif

static inline double rexp (double x)
{	if (x>=1000) { error=1; return 0; }
	else return exp(x);
}

void mexp (header *hd)
{
#ifdef FLOAT_TEST
	spread1i(rexp,rcexp,iexp,hd);
#else
	spread1i(exp,cexp,iexp,hd);
#endif
	test_error("exp");
}

static double rarg (double x)
{	if (x>=0) return 0.0;
	else return M_PI;
}

void mlog (header *hd)
{
#ifdef FLOAT_TEST
	spread1i(rlog,cclog,ilog,hd);
#else
	spread1i(log,cclog,ilog,hd);
#endif
	test_error("log");
}

#ifdef FLOAT_TEST
static double rsqrt (double x)
{	if (x<0.0) { error=1; return 0; }
	else return sqrt(x);
}
#endif

void msqrt (header *hd)
{	spread1i(
#ifdef FLOAT_TEST
	rsqrt,
#else
	sqrt,
#endif
	csqrt,isqrt,hd);
	test_error("sqrt");
}

static void rrem (double *x, double *n, double *y)
{
	*y=fmod(*x,*n);
}

void mrem (header *hd)
{	spreadf2(rrem,0,0,hd);
	test_error("mod");
}

static void rmod (double *x, double *n, double *y)
{
	if (*n==0.0)
		*y=*x;
	else {
		*y=fmod(*x,*n);
		if (*y<0) *y+=*n;
	}
}

void mmod (header *hd)
{	spreadf2(rmod,0,0,hd);
	test_error("mod");
}

static void cpow (double *x, double *xi, double *y, double *yi,
	double *z, double *zi)
{	double l,li,w,wi;
	if (fabs(*x)<epsilon && fabs(*xi)<epsilon)
	{	*z=*zi=0.0; return;
	}
	cclog(x,xi,&l,&li);
	complex_multiply(y,yi,&l,&li,&w,&wi);
	cexp(&w,&wi,z,zi);
}

static void rpow (double *x, double *y, double *z)
{	int n;
	if (*x>0.0) *z=pow(*x,*y);
	else if (*x==0.0) if (*y==0.0) *z=1.0; else *z=0.0;
	else
	{	n=(int)*y;
		if (n!=*y)
		{   output("Illegal argument for ^\n"); error=1;
			return;
		}
		if (n%2) *z=-pow(-*x,n);
		else *z=pow(-*x,n);
	}
}

void mpower (header *hd)
{	spreadf2(rpow,cpow,ipow,hd);
	test_error("^");
}

/****************************************************************************
 *	degree, sign, ceil, floor functions
 *
 ****************************************************************************/
static double rdegree (double x)
{	return x/180*M_PI;
}

void mdegree (header *hd)
{	spread1(rdegree,0,hd);
	test_error("signum");
}

static double rsign (double x)
{	if (x<0) return -1;
	else if (x<=0) return 0;
	else return 1;
}

void msign (header *hd)
{	spread1(rsign,0,hd);
	test_error("sign");
}

void mceil (header *hd)
{	spread1(ceil,0,hd);
	test_error("ceil");
}

void mfloor (header *hd)
{	spread1(floor,0,hd);
	test_error("floor");
}
static double rounder;

static double rround (double x)
{	x*=rounder;
	if (x>0) x=floor(x+0.5);
	else x=-floor(-x+0.5);
	return x/rounder;
}

static void cround (double *x, double *xi, double *z, double *zi)
{	*z=rround(*x);
	*zi=rround(*xi);
}

static double frounder[]={1.0,10.0,100.0,1000.0,10000.0,100000.0,1000000.0,
10000000.0,100000000.0,1000000000.0,10000000000.0};

void mround (header *hd)
{	header *hd1;
	int n;
	hd1=next_param(hd);
	if (hd1) hd1=getvalue(hd1); if (error) return;
	if (hd1->type!=s_real) wrong_arg_in("round");
	n=(int)(*realof(hd1));
	if (n>0 && n<11) rounder=frounder[n];
	else rounder=pow(10.0,n);
	spread1(rround,cround,hd);
	test_error("round");
}

/****************************************************************************
 *	complex specific functions
 *
 ****************************************************************************/

static void cconj (double *x, double *xi, double *z, double *zi)
{	*zi=-*xi; *z=*x;
}

static double ident (double x)
{	return x;
}

void mconj (header *hd)
{	spread1(ident,cconj,hd);
	test_error("conj");
}

static void crealpart (double *x, double *xi, double *z)
{	*z=*x;
}

void mre (header *hd)
{	spread1r(ident,crealpart,hd);
	test_error("re");
}

static double zero (double x)
{	return 0.0;
}

static void cimagpart (double *x, double *xi, double *z)
{	*z=*xi;
}

void mim (header *hd)
{	spread1r(zero,cimagpart,hd);
	test_error("im");
}

void marg (header *hd)
{	spread1r(rarg,carg,hd);
	test_error("arg");
}

static void cxabs (double *x, double *xi, double *z)
{	*z=sqrt(*x * *x + *xi * *xi);
}

void mabs (header *hd)
{   spread1ir(fabs,cxabs,iabs,hd);
	test_error("abs");
}

/****************************************************************************
 *	constants
 *
 ****************************************************************************/

void mpi (header *hd)
{	new_real(M_PI,"");
}

void margn (header *hd)
{	new_real(actargn,"");
}

void mtime (header *hd)
{	new_real(myclock(),"");
}

void mfree (header *hd)
{	new_real(ramend-endlocal,"");
}

#ifndef NOSHRINK

void mshrink (header *hd)
{	header *st=hd,*result;
	unsigned long size;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_real) wrong_arg_in("shrink");
	if (*realof(hd)>LONG_MAX) wrong_arg_in("shrink");
	size=(unsigned long)*realof(hd);
	if (ramend-size<newram)
	{	output("Cannot shrink that much!\n");
		error=171; return;
	}
	if (size)
	{	if (shrink(size)) ramend-=size;
		else
		{	output("Shrink failed!\n"); error=172; return;
		}
	}
	result=new_real(ramend-ramstart,"");
	moveresult(st,result);
}

#endif

void mepsilon (header *hd)
{	new_real(epsilon,"");
}

void msetepsilon (header *hd)
{	header *stack=hd,*hd1,*result;
	hd1=getvalue(hd); if (error) return;
	if (hd1->type!=s_real) wrong_arg_in("setepsilon");
	result=new_real(epsilon,"");
	changedepsilon=epsilon=*realof(hd1);
	moveresult(stack,result);
}

void mlocalepsilon (header *hd)
{	header *stack=hd,*hd1,*result;
	hd1=getvalue(hd); if (error) return;
	if (hd1->type!=s_real) wrong_arg_in("localepsilon");
	result=new_real(epsilon,"");
	epsilon=*realof(hd1);
	moveresult(stack,result);
}

void mtype (header *hd)
{	header *st=hd,*result;
	hd=getvalue(hd); if (error) return;
	result=new_real(hd->type,"");
	moveresult(st,result);
}

void mintersects (header *hd)
{	spreadf2r(0,0,iintersects,hd);
	test_error("intersects");
}

/****************************************************************************
 *	compare operator
 *
 ****************************************************************************/

static void rgreater (double *x, double *y, double *z)
{	if (*x>*y) *z=1.0;
	else *z=0.0;
}

static void igreater (double *xa, double *xb, double *ya, double *yb,
	double *z)
{	if (*xa>*yb) *z=1.0;
	else *z=0.0;
}

void mgreater (header *hd)
{   header *st=hd,*hd1,*result,*hdv;
	hdv=getvariable(hd);
	if (hdv->type==s_string)
	{   hd=getvalue(hd);
		hd1=getvalue(nextof(st)); if (error) return;
		if (hd1->type!=s_string) wrong_arg_in("==");
		result=new_real(strcmp(stringof(hd),stringof(hd1))>0,"");
		moveresult(st,result);
	}
	else spreadf2r(rgreater,0,igreater,st);
	test_error("==");
}

static void rless (double *x, double *y, double *z)
{	if (*x<*y) *z=1.0;
	else *z=0.0;
}

static void ilesst (double *xa, double *xb, double *ya, double *yb,
	double *z)
{	if (*xb<*ya) *z=1.0;
	else *z=0.0;
}

void mless (header *hd)
{   header *st=hd,*hd1,*result,*hdv;
	hdv=getvariable(hd);
	if (hdv->type==s_string)
	{   hd=getvalue(hd);
		hd1=getvalue(nextof(st)); if (error) return;
		if (hd1->type!=s_string) wrong_arg_in("==");
		result=new_real(strcmp(stringof(hd),stringof(hd1))<0,"");
		moveresult(st,result);
	}
	else spreadf2r(rless,0,ilesst,st);
	test_error("<");
}

static void rgreatereq (double *x, double *y, double *z)
{	if (*x>=*y) *z=1.0;
	else *z=0.0;
}

static void igreatereq (double *xa, double *xb, double *ya, double *yb, double *z)
{	if (*xa>=*yb) *z=1.0;
	else *z=0.0;
}

void mgreatereq (header *hd)
{   header *st=hd,*hd1,*result,*hdv;
	hdv=getvariable(hd);
	if (hdv->type==s_string)
	{   hd=getvalue(hd);
		hd1=getvalue(nextof(st)); if (error) return;
		if (hd1->type!=s_string) wrong_arg_in("==");
		result=new_real(strcmp(stringof(hd),stringof(hd1))>=0,"");
		moveresult(st,result);
	}
	else spreadf2r(rgreatereq,0,igreatereq,st);
	test_error(">=");
}

static void rlesseq (double *x, double *y, double *z)
{	if (*x<=*y) *z=1.0;
	else *z=0.0;
}

static void ilesseq (double *xa, double *xb, double *ya, double *yb,
	double *z)
{	if (*xb<=*ya) *z=1.0;
	else *z=0.0;
}

void mlesseq (header *hd)
{   header *st=hd,*hd1,*result,*hdv;
	hdv=getvariable(hd);
	if (hdv->type==s_string)
	{   hd=getvalue(hd);
		hd1=getvalue(nextof(st)); if (error) return;
		if (hd1->type!=s_string) wrong_arg_in("==");
		result=new_real(strcmp(stringof(hd),stringof(hd1))<=0,"");
		moveresult(st,result);
	}
	else spreadf2r(rlesseq,0,ilesseq,st);
	test_error("<=");
}

static void requal (double *x, double *y, double *z)
{	if (*x==*y) *z=1.0;
	else *z=0.0;
}

static void cequal (double *x, double *xi, double *y, double *yi, double *z)
{	if (*x==*y && *xi==*yi) *z=1.0;
	else *z=0.0;
}


void mequal (header *hd)
{   header *st=hd,*hd1,*result,*hdv;
	hdv=getvariable(hd);
	if (hdv->type==s_string)
	{   hd=getvalue(hd);
		hd1=getvalue(nextof(st)); if (error) return;
		if (hd1->type!=s_string) wrong_arg_in("==");
		result=new_real(strcmp(stringof(hd),stringof(hd1))==0,"");
		moveresult(st,result);
	}
	else spreadf2r(requal,cequal,cequal,st);
	test_error("==");
}

static void runequal (double *x, double *y, double *z)
{	if (*x!=*y) *z=1.0;
	else *z=0.0;
}

static void cunequal (double *x, double *xi, double *y, double *yi, double *z)
{	if (*x!=*y || *xi!=*yi) *z=1.0;
	else *z=0.0;
}

static void iunequal (double *x, double *xi, double *y, double *yi, double *z)
{	if (*x!=*y || *xi!=*yi) *z=1.0;
	else *z=0.0;
}

void munequal (header *hd)
{   header *st=hd,*hd1,*result,*hdv;
	hdv=getvariable(hd);
	if (hdv->type==s_string)
	{   hd=getvalue(hd);
		hd1=getvalue(nextof(st)); if (error) return;
		if (hd1->type!=s_string) wrong_arg_in("==");
		result=new_real(strcmp(stringof(hd),stringof(hd1))!=0,"");
		moveresult(st,result);
	}
	else spreadf2r(runequal,cunequal,iunequal,st);
	test_error("<>");
}

static void raboutequal (double *x, double *y, double *z)
{   double rx=fabs(*x),ry=fabs(*y);
	if (rx<epsilon)
	{	if (ry<epsilon) *z=1.0;
		else *z=0.0;
	}
	else if (ry<epsilon)
	{	if (rx<epsilon) *z=1.0;
		else *z=0.0;
	}
	else if (fabs(*x-*y)/rx<epsilon) *z=1.0;
	else *z=0.0;
}

static void caboutequal
	(double *x, double *xi, double *y, double *yi, double *z)
{	raboutequal(x,y,z);
	if (*z!=0.0) raboutequal(xi,yi,z);
}

static void iaboutequal (double *x, double *xi, double *y, double *yi, double *z)
{	raboutequal(x,y,z);
	if (*z!=0.0) raboutequal(xi,yi,z);
}

void maboutequal (header *hd)
{	spreadf2r(raboutequal,caboutequal,iaboutequal,hd);
	test_error("~=");
}

/****************************************************************************
 *	logical operators
 *
 ****************************************************************************/

static void ror (double *x, double *y, double *z)
{	if (*x!=0.0 || *y!=0.0) *z=1.0;
	else *z=0.0;
}

void mor (header *hd)
{	spreadf2(ror,0,ior,hd);
	test_error("||");
}

static void rrand (double *x, double *y, double *z)
{	if (*x!=0.0 && *y!=0.0) *z=1.0;
	else *z=0.0;
}

void mand (header *hd)
{	spreadf2(rrand,0,iand,hd);
	test_error("&&");
}
static double rnot (double x)
{	if (x!=0.0) return 0.0;
	else return 1.0;
}

static void cnot (double *x, double *xi, double *r)
{	if (*x==0.0 && *xi==0.0) *r=1.0;
	else *r=0.0;
}

void mnot (header *hd)
{	spread1r(rnot,cnot,hd);
	test_error("!");
}

/****************************************************************************
 *	matrix operators
 *
 ****************************************************************************/

void mcomplex (header *hd)
{	header *st=hd,*result;
	double *m,*mr;
	long i,n;
	int c,r;
	hd=getvalue(hd);
	if (hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_cmatrix(r,c,""); if (error) return;
		n=(long)r*c;
        mr=matrixof(result)+(long)2*(n-1);
		m+=n-1;
		for (i=0; i<n; i++)
		{	*mr=*m--; *(mr+1)=0.0; mr-=2;
		}
		moveresult(st,result);
	}
	else if (hd->type==s_real)
	{	result=new_complex(*realof(hd),0.0,""); if (error) return;
		moveresult(st,result);
	}
}





/***************************************************************************
 *	events handling
 *
 ***************************************************************************/
void mwait (header *hd)
{	header *st=hd,*result;
	double now;
	int h;

	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_real) wrong_arg_in("wait");
	now=myclock();
	sys_wait(*realof(hd),&h);
	if (h==escape) { output("Interrupt\n"); error=1; return; }
	now=myclock()-now;
	if (now>*realof(hd)) now=*realof(hd);
	result=new_real(now,"");
	if (error) return;
	moveresult(st,result);
}

void mkey (header *hd)
{	int scan,key;
	key=wait_key(&scan);
	if (scan==escape) { output("Interrupt\n"); error=1; return; }
	if (key) new_real(key,"");
	else new_real(scan,"");
}

void mcode (header *hd)
{	new_real(test_code(),"");
}

void minput (header *hd)
{	header *st=hd,*result;
	char input[1024],*oldnext;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_string) wrong_arg_in("input");
retry: output(stringof(hd)); output("? ");
	nojump=1;
	edit(input);
	stringon=1;
	oldnext=next; next=input; result=scan_value(); next=oldnext;
	stringon=0;
	if (error) 
	{	output("Error in input!\n"); error=0; goto retry;
	}
	moveresult(st,result);
}

void mlineinput (header *hd)
{	header *st=hd,*result;
	char input[1024];
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_string) wrong_arg_in("lineinput");
	nojump=1;
	output(stringof(hd)); output("? ");
	edit(input);
	result=new_string(input,strlen(input),"");
	moveresult(st,result);
}

void minterpret (header *hd)
{	header *st=hd,*result;
	char *oldnext;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_string) wrong_arg_in("interpret");
	stringon=1;
	oldnext=next; next=stringof(hd); result=scan(); next=oldnext;
	stringon=0;
    if (error) { result=new_string("Error",8,""); error=0; }
	moveresult(st,result);
}

void mevaluate (header *hd)
{	header *st=hd,*result;
	char *oldnext;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_string) wrong_arg_in("evaluate");
	stringon=1;
	oldnext=next; next=stringof(hd); result=scan(); next=oldnext;
	stringon=0;
    if (error) { result=new_string("Syntax error!",16,""); }
	moveresult(st,result);
}

/***************************************************************************
 *	statistics - random numbers
 *
 ***************************************************************************/

void mfastrandom (header *hd)
{	header *st=hd,*result;
	double *m;
	int r,c;
	long k,n;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_matrix || dimsof(hd)->r!=1 || dimsof(hd)->c!=2
		|| *(m=matrixof(hd))<0 || *m>=INT_MAX
		|| *(m+1)<0 || *(m+1)>INT_MAX)
		wrong_arg_in("random");
	r=(int)*m;
	c=(int)*(m+1);
	result=new_matrix(r,c,""); if (error) return;
	m=matrixof(result);
	n=(long)c*r;
	for (k=0; k<n; k++) *m++=(double)rand()/(double)RAND_MAX;
	moveresult(st,result);
}

void mnormal (header *hd)
{	header *st=hd,*result;
	double *m,r1,r2;
	int r,c;
	long k,n;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_matrix || dimsof(hd)->r!=1 || dimsof(hd)->c!=2
		|| *(m=matrixof(hd))<0 || *m>=INT_MAX 
		|| *(m+1)<0 || *(m+1)>INT_MAX)
		wrong_arg_in("normal");
	r=(int)*m;
	c=(int)*(m+1);
	result=new_matrix(r,c,""); if (error) return;
	m=matrixof(result);
	n=(long)c*r;
	for (k=0; k<n; k++) 
	{	r1=(double)rand()/(double)RAND_MAX;
		if (r1==0.0) *m++=0.0;
		else
		{	r2=(double)rand()/(double)RAND_MAX;
			*m++=sqrt(-2*log(r1))*cos(2*M_PI*r2);
		}
	}
	moveresult(st,result);
}

static double rfak (double x)
{	int i,n;
	double res=1;
	if (x<2 || x>INT_MAX) return 1.0;
	n=(int)x;
	for (i=2; i<=n; i++) res=res*i;
	return res;
}

void mfak (header *hd)
{	spread1(rfak,0,hd);
	test_error("fak");
}

static double rlogfak (double x)
{	int i,n;
	double res=0;
	if (x<2 || x>INT_MAX) return 0.0;
	n=(int)x;
	for (i=2; i<=n; i++) res=res+log(i);
	return res;
}

void mlogfak (header *hd)
{	spread1(rlogfak,0,hd);
	test_error("logfak");
}

static void rbin (double *x, double *y, double *z)
{   long i,n,m,k;
	double res;
	n=(long)*x; m=(long)*y;
	if (m<=0) *z=1.0;
	else
	{	res=k=n-m+1;
		for (i=2; i<=m; i++) { k++; res=res*k/i; }
		*z=res;
	}
}

void mbin (header *hd)
{	spreadf2(rbin,0,0,hd);
	test_error("bin");
}

static void rlogbin (double *x, double *y, double *z)
{   long i,n,m,k;
	double res;
	n=(long)*x; m=(long)*y;
	if (m<=0) *z=0.0;
	else
	{   k=n-m+1;
		res=log(n-m+1);
		for (i=2; i<=m; i++) { k++; res+=log(k)-log(i); }
		*z=res;
	}
}

void mlogbin (header *hd)
{	spreadf2(rlogbin,0,0,hd);
	test_error("bin");
}

static void rtd (double *xa, double *yf, double *zres)
{	double t,t1,a,b,h,z,p,y,x;
	int flag=0;
	if (fabs(*xa)<epsilon) { *zres=0.5; return; }
	if (*xa<0) flag=1;
	t=*xa * *xa;
	if (t>=1) { a=1; b=*yf; t1=t; }
	else { a=*yf; b=1; t1=1/t; }
	y=2/(9*a); z=2/(9*b);
	h=pow(t1,1.0/3);
	x=fabs((1-z)*h-1+y)/sqrt(z*h*h+y);
	if (b<4) x=x*(1+0.08*x*x*x*x/(b*b*b));
	h=1+x*(0.196854+x*(0.115194+x*(0.000344+x*0.019527)));
	p=0.5/(h*h*h*h);
	if (t<0.5) *zres=p/2+0.5;
	else *zres=1-p/2;
	if (flag) *zres=1-*zres;
}

void mtd (header *hd)
{	spreadf2(rtd,0,0,hd);
	test_error("tdis");
}

static double invgauss (double a)
{	double t,c,d;
	int flag=0;
	if (a<0.5) { a=1-a; flag=1; }
	t=sqrt(-2*log(fabs(1-a)));
	c=2.515517+t*(0.802853+t*0.010328);
	d=1+t*(1.432788+t*(0.189269+t*0.001308));
	if (flag) return (c/d-t);
	else return t-c/d;
}

void minvgauss (header *hd)
{	spread1(invgauss,0,hd);
	test_error("invnormaldis");
}

static void invrtd (double *x, double *y, double *zres)
{	double z=*x,f=*y,g1,g2,g3,g4,z2;
	int flag=0;
	if (z<0.5) { flag=1; z=1-z; }
	z=invgauss(z);
	z2=z*z;
	g1=z*(1+z2)/4.0;
	g2=z*(3+z2*(16+5*z2))/96.0;
	g3=z*(-15+z2*(17+z2*(19+z2*3)))/384.0;
	g4=z*(-945+z2*(-1920+z2*(1482+z2*(776+z2*79))))/92160.0;
	*zres=(((g4/f+g3)/f+g2)/f+g1)/f+z;
	if (flag) *zres=-*zres;
}

void minvtd (header *hd)
{	spreadf2(invrtd,0,0,hd);
	test_error("invtdis");
}

static void chi (double *xa, double *yf, double *zres)
{	double ch=*xa,x,y,s,t,g,j=1;
	long i,p,f;
	f=(long)*yf;
	if (ch<epsilon) { *zres=0.0; return; }
	p=(f+1)/2;
	for (i=f; i>=2; i-=2) j=j*i;
	x=pow(ch,p)*exp(-(ch/2))/j;
	if (f%2==0) y=1;
	else y=sqrt(2/(ch*M_PI));
	s=1; t=1; g=f;
	while (t>1e-5)
	{	g=g+2;
		t=t*ch/g;
		s=s+t;
	}
	*zres=x*y*s;
}

void mchi (header *hd)
{	spreadf2(chi,0,0,hd);
	test_error("chidis");
}

static double f1,f2;

static double rfd (double F)
{	double f0,a,b,h,z,p,y,x;
	if (F<epsilon) return 0.0;
	if (F<1) { a=f2; b=f1; f0=1/F; }
	else { a=f1; b=f2; f0=F; }
	y=2/(9*a); z=2/(9*b);
	h=pow(f0,1.0/3);
	x=fabs((1-z)*h-1+y)/sqrt(z*h*h+y);
	if (b<4) x=x*(1+0.08*x*x*x*x/(b*b*b));
	h=1+x*(0.196854+x*(0.115194+x*(0.000344+x*0.019527)));
	p=0.5/(h*h*h*h);
	if (F>=1) return 1-p;
	else return p;
}

void mfdis (header *hd)
{	header *st=hd,*hd1,*hd2=0;
	hd1=next_param(st);
	if (hd1)
	{	hd2=next_param(hd1);
		hd1=getvalue(hd1);
	}
	if (hd2) hd2=getvalue(hd2);
	if (error) return;
	if (hd1->type!=s_real || hd2->type!=s_real) wrong_arg_in("fdis");
	f1=*realof(hd1);
	f2=*realof(hd2);
	spread1(rfd,0,hd);
	test_error("fdis");
}

/*****************************************************************************
 *	sorting
 *
 *****************************************************************************/
static void rmax (double *x, double *y, double *z)
{	if (*x>*y) *z=*x;
	else *z=*y;
}

void mmax (header *hd)
{   spreadf2(rmax,0,imax,hd);
	test_error("max");
}

static void rmin (double *x, double *y, double *z)
{	if (*x>*y) *z=*y;
	else *z=*x;
}

void mmin (header *hd)
{	spreadf2(rmin,0,imin,hd);
	test_error("min");
}

typedef struct { double val; int ind; } sorttyp;

static int sorttyp_compare (const sorttyp *x, const sorttyp *y)
{	if (x->val>y->val) return 1;
	else if (x->val==y->val) return 0;
	else return -1;
}

void msort (header *hd)
{	header *st=hd,*result,*result1;
	double *m,*m1;
	sorttyp *t;
	int r,c,i;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_real && hd->type!=s_matrix) wrong_arg_in("sort");
	getmatrix(hd,&r,&c,&m);
	if (c==1 || r==1) result=new_matrix(r,c,"");
	else wrong_arg_in("sort");
	if (error) return;
	result1=new_matrix(r,c,"");
	if (error) return;
	if (c==1) c=r;
	if (c==0) wrong_arg_in("sort");
	if (!freeram(c*sizeof(sorttyp)))
	{	output("Out of memory!\n"); error=600; return; 
	}
	t=(sorttyp *)newram;
	for (i=0; i<c; i++)
	{	t->val=*m++; t->ind=i; t++;
	}
	qsort(newram,c,sizeof(sorttyp),
		(int (*) (const void *, const void *))sorttyp_compare);
	m=matrixof(result); m1=matrixof(result1);
	t=(sorttyp *)newram;
	for (i=0; i<c; i++)
	{	*m++=t->val; *m1++=t->ind+1; t++;
	}
	moveresult(st,result);
	moveresult(nextof(st),result1);
}


void mstatistics (header *hd)
{	header *st=hd,*hd1,*result;
	int i,n,r,c,k;
	double *m,*mr;
	hd=getvalue(hd);
	hd1=next_param(st);
	if (hd1) hd1=getvalue(hd1); if (error) return;
	if (hd1->type!=s_real || hd->type!=s_matrix) wrong_arg_in("count");
	if (*realof(hd1)>INT_MAX || *realof(hd1)<2) wrong_arg_in("count");
	n=(int)*realof(hd1);
	getmatrix(hd,&r,&c,&m);
	if (r!=1 && c!=1) wrong_arg_in("count");
	if (c==1) c=r;
	result=new_matrix(1,n,""); if (error) return;
	mr=matrixof(result); for (i=0; i<n; i++) *mr++=0.0;
	mr=matrixof(result);
	for (i=0; i<c; i++)
	{	if (*m>=0 && *m<n)
		{	k=(int)floor(*m);
			mr[k]+=1.0;
		}
		m++;
	}
	moveresult(st,result);
}



void mmax1 (header *hd)
{	header *result,*st=hd;
	double x,*m,*mr,max,max1;
	int r,c,i,j;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_matrix(r,1,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	max=*m++;
			for (j=1; j<c; j++)
			{	x=*m++;
				if (x>max) max=x;
			}
			*mr++=max;
		}
	}
	else if (hd->type==s_interval || hd->type==s_imatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_imatrix(r,1,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	max=*m++;
			max1=*m++;
			for (j=1; j<c; j++)
			{	imax(&max,&max1,m,m+1,&max,&max1);
				m+=2;
			}
			*mr++=max;
			*mr++=max1;
		}
	}
	else wrong_arg_in("max");
	moveresult(st,result);
}

void mmin1 (header *hd)
{	header *result,*st=hd;
	double x,*m,*mr,max,max1;
	int r,c,i,j;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_matrix(r,1,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++) 
		{	max=*m; m++;
			for (j=1; j<c; j++) 
			{	x=*m++;
				if (x<max) max=x;
			}
			*mr++=max;
		}
	}
	else if (hd->type==s_interval || hd->type==s_imatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_imatrix(r,1,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	max=*m++;
			max1=*m++;
			for (j=1; j<c; j++)
			{	imin(&max,&max1,m,m+1,&max,&max1);
				m+=2;
			}
			*mr++=max;
			*mr++=max1;
		}
	}
	else wrong_arg_in("min");
	moveresult(st,result);
}



/****************************************************************************
 *	types
 *
 ****************************************************************************/
void miscomplex (header *hd)
{	header *st=hd,*result;
	hd=getvalue(hd);
	if (hd->type==s_complex || hd->type==s_cmatrix)
		result=new_real(1.0,"");
	else result=new_real(0.0,"");
	if (error) return;
	moveresult(st,result);
}

void misreal (header *hd)
{	header *st=hd,*result;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_matrix)
		result=new_real(1.0,"");
	else result=new_real(0.0,"");
	if (error) return;
	moveresult(st,result);
}

void misfunction (header *hd)
{	header *st=hd,*result;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_string
		&& (searchudf(stringof(hd))!=0 || find_builtin(stringof(hd))!=0))
			result=new_real(1.0,"");
	else result=new_real(0.0,"");
	if (error) return;
	moveresult(st,result);
}

void misvar (header *hd)
{	header *st=hd,*result;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_string
		&& searchvar(stringof(hd))!=0)
			result=new_real(1.0,"");
	else result=new_real(0.0,"");
	if (error) return;
	moveresult(st,result);
}

void misinterval (header *hd)
{	header *st=hd,*result;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_interval || hd->type==s_imatrix)
		result=new_real(1.0,"");
	else result=new_real(0.0,"");
	if (error) return;
	moveresult(st,result);
}

/****************************************************************************
 *	rounding
 *
 ****************************************************************************/

void mchar (header *hd)
{	header *st=hd,*result;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_real) wrong_arg_in("char");
	result=new_string("a",1,""); if (error) return;
	*stringof(result)=(char)*realof(hd);
	moveresult(st,result);
}

void mascii (header *hd)
{	header *st=hd,*result;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_string) wrong_arg_in("ascii");
	result=new_real(*stringof(hd),""); if (error) return;
	moveresult(st,result);
}

void merror (header *hd)
{	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_string) wrong_arg_in("error");
	output1("Error : %s\n",stringof(hd));
	error=301;
}

extern int preventoutput;

void merrlevel (header *hd)
{	header *st=hd,*res;
	char *oldnext;
	int en;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_string) wrong_arg_in("errorlevel");
	stringon=1;
	oldnext=next; next=stringof(hd);
	res=new_real(0,"");
	preventoutput=1;
	scan();
	preventoutput=0;
	next=oldnext;
	stringon=0;
	en=error; error=0;
	if (en)
	{	*realof(res)=en;
		moveresult(st,res);
	}
	else
	{	moveresult1(st,res);
	}
}

void mprintf (header *hd)
{	header *st=hd,*hd1,*result;
	char string[1024];
	hd1=next_param(hd);
	hd=getvalue(hd);
	hd1=getvalue(hd1); if (error) return;
	if (hd->type!=s_string || hd1->type!=s_real)
		wrong_arg_in("printf");
	sprintf(string,stringof(hd),*realof(hd1));
	result=new_string(string,strlen(string),""); if (error) return;
	moveresult(st,result);
}

void msetkey (header *hd)
/*****
	set a function key
*****/
{	header *st=hd,*hd1,*result;
	char *p;
	int n;
	hd=getvalue(hd); if (error) return;
	hd1=nextof(st); hd1=getvalue(hd1); if (error) return;
	if (hd->type!=s_real || hd1->type!=s_string) wrong_arg_in("setkey");
	n=(int)(*realof(hd))-1; p=stringof(hd1);
	if (n<0 || n>=10 || strlen(p)>63) wrong_arg_in("setkey");
	result=new_string(fktext[n],strlen(fktext[n]),"");
	if (error) return;
	strcpy(fktext[n],p);
	moveresult(st,result);
}


void mcd (header *hd)
{	header *st=hd,*result;
	char *path;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_string) wrong_arg_in("cd");
	path=cd(stringof(hd));
	result=new_string(path,strlen(path),"");
	moveresult(st,result);
}

static char **search_entries=NULL;
static int  search_count=0,search_current=0;

void mdir (header *hd)
{	header *st=hd,*result;
//	char *name;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_string) wrong_arg_in("dir");
	if (search_entries) {
		int i;
		for (i=0;i<search_count;i++)
			free(search_entries[i]);
		free(search_entries);
		search_entries=NULL;
		search_count=0;
		search_current=0;
	}
	scan_dir(".",stringof(hd),&search_entries,&search_count);
	if (search_entries)
		result=new_string(search_entries[0],strlen(search_entries[0]),"");
	else
		result=new_string("",0,"");
//	name=dir(stringof(hd));
//	if (name) result=new_string(name,strlen(name),"");
//	else result=new_string("",0,"");
	if (error) return;
	moveresult(st,result);
}

void mdir0 (header *hd)
{
	header *result;
	if (search_entries) {
		search_current++;
		if (search_current<search_count)
			result=new_string(search_entries[search_current],strlen(search_entries[search_current]),"");
		if (error) return;
		return;
	}
	result=new_string("",0,"");
}


void margs (header *hd)
/* return all args from realof(hd)-st argument on */
{	header *st=hd,*hd1,*result;
	int i,n;
	unsigned long size;
	hd=getvalue(hd);
	if (hd->type!=s_real) wrong_arg_in("args");
	n=(int)*realof(hd);
	if (n<1) wrong_arg_in("args");
	if (n>actargn)
	{	newram=(char *)st; return;
	}
	result=(header *)startlocal; i=1;
	while (i<n && result<(header *)endlocal)
	{	result=nextof(result); i++;
	}
	hd1=result;
	while (i<actargn+1 && hd1<(header *)endlocal)
	{	hd1=nextof(hd1); i++;
	}
	size=(char *)hd1-(char *)result;
	if (size<=0)
	{	output("Error in args!\n"); error=2021; return;
	}
	memmove((char *)st,(char *)result,size);
	newram=(char *)st+size;
}

void margs0 (header *hd)
/* return all args from realof(hd)-st argument on */
{	header *st=(header *)newram,*hd1,*result;
	int i,n;
	unsigned long size;
	n=actsp+1;
	if (actsp==0) return;
	if (n>actargn) n=actargn;
	result=(header *)startlocal; i=1;
	while (i<n && result<(header *)endlocal)
	{	result=nextof(result); i++;
	}
	hd1=result;
	while (i<actargn+1 && hd1<(header *)endlocal)
	{	hd1=nextof(hd1); i++;
	}
	size=(char *)hd1-(char *)result;
	if (size<=0)
	{	output("Error in args!\n"); error=2021; return;
	}
	memmove((char *)st,(char *)result,size);
	newram=(char *)st+size;
}

void mname (header *hd)
{	header *st=hd,*result;
	hd=getvalue(hd); if (error) return;
	result=new_string(hd->name,strlen(hd->name),"");
	moveresult(st,result);
}


#ifdef WAVES

void mplaywav (header *hd)
{	header *st=hd;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_string) wrong_arg_in("playwave");
	sys_playwav(stringof(hd));
	moveresult(st,hd);
}

#endif