/*
 * ordinary differential equation program copyright (C) 2011 - 2013 H.Niwa
 */

/*
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2, or (at your option)
 * any later version.

 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 * 02110-1301, USA.
 */

#include "config.h"

#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>

#include <string>
#include <complex>

#include "syserr.h"

#include "bin_node.h"
#include "gc.h"
#include "var.h"
#include "pred.h"
#include "context.h"
#include "unify.h"
#include "builtin.h"
#include "sysmodule.h"
#include "context.h"
#include "ode.h"

extern int BreakFlag;
extern int ThrowFlag;

int DoODE(Context* cx, Node* goalscar, List* module);
int Integral(Context* cx, Node* goalscar, List* module, Node* goalscar_car);
int ODEprint(Context* cx, Node* goalscar, List* module);
int ODEprintf(Context* cx, Node* goalscar, List* module);
int ODEcmd(Context* cx, Node* goalscar, List* module);

Node* GetIntegralCx(Context* cx, Node* id)
{
	Node* n = cx->integral;
	Node* np = Nil;
	
	if ((n == Nil) || (n->kind() != LIST)) {
		return Nil;
	}

	for (np = n; np->kind() == LIST; np=np->Cdr()) {
		Node* ni = np->Car();
		if (ni->kind() != LIST) {
			continue;
		}

		if (ni->Car() == id) {
			return ni;
		}
	}
	return Nil;
}

void SetIntegralCx(Context* cx, Node* id, Node* val)
{
	Node* ni = GetIntegralCx(cx, id);

	if (ni == Nil) {
		cx->integral = Cons(Cons(id, val), cx->integral);
		return;
	}

	((List*)ni)->SetCdr(val);
}

long double Simpson(long double t, long double y, long double dlt, 
			long long &flg, long double &ig, long double constant)
{
	long double rig = 0.0L;

	if (flg == 0) {
		flg = 1;
		rig = constant;
		ig = y;
		return rig;
	} else if (flg == 1) {
		flg = 2;
		rig = (ig + 2.0L*y) * dlt / 3.0L + constant;
		ig += (4.0L*y);
		return rig;
	} else {
		flg = 1;
		rig = (ig + y) * dlt / 3.0L + constant;
		ig += (2.0L*y);
		return rig;
	}
	
}

int DoODE(Context* cx, Node* goalscar, List* module)
{
	if (ListLength(goalscar) < 2) {
		syserr("usage : <ODE [VAR] (VAR INIT END STEP) PRED>\n");
		return 0;
	}

	Node*	nresult = goalscar->Cdr()->Car();
	
 	int rn;
	
	if ((rn = FuncArg(cx, nresult, goalscar, module)) <= 0) {
		syserr("for: failed in the evaluation of the argument. \n");
		return 0;
	}

	Node*	gl = Nil;
	Node*	glargs = Nil;
	if (nresult->kind() == UNDEF) {
		gl = goalscar->Cdr()->Cdr()->Cdr();
		glargs = goalscar->Cdr()->Cdr()->Car();
	} else {
		gl = goalscar->Cdr()->Cdr();
		glargs = nresult;
		nresult = Nil;
	}

	if ((rn = FuncArg(cx, glargs, goalscar, module)) <= 0) {
		syserr("for: failed in the evaluation of the argument. \n");
		return 0;
	}

	int nc = ListLength(glargs);
	if (nc != 4) {
		syserr("usage : <ODE [VAR] (VAR INIT END STEP) PRED>\n");
		return 0;
	}
		
	Node* nvar = glargs->Car()->Val();

	if (nvar->kind() != UNDEF) {
		syserr("usage : <ODE [VAR] (VAR INIT END STEP) PRED>\n");
		return 0;
	}

	glargs = glargs->Cdr();

	long double init_val = 0.0;
	long double to_val=0.0;
	long double step_val=0.0;

	Node* ninitval = glargs->Car()->Val();
	if (ninitval->kind() != ATOM) {
		syserr("usage : <ODE [VAR] (VAR INIT END STEP) PRED>\n");
		return 0;
	}
	if (!((Atom*)ninitval)->toFloat(init_val)) {
		syserr("usage : <ODE [VAR] (VAR INIT END STEP) PRED>\n");
		return 0;
	}
		
	Node* nto_val = glargs->Cdr()->Car()->Val();
	if (nto_val->kind() != ATOM) {
		syserr("usage : <ODE [VAR] (VAR INIT END STEP) PRED>\n");
		return 0;
	}
	if (!((Atom*)nto_val)->toFloat(to_val)) {
		syserr("usage : <ODE [VAR] (VAR INIT END STEP) PRED>\n");
		return 0;
	}

	if ((init_val > to_val) || (step_val < 0.0L)) {
		syserr("usage : <ODE [VAR] (VAR INIT END STEP) PRED>\n");
		return 0;
	}

	Node* nstep_val = glargs->Cdr()->Cdr()->Car()->Val();
	if (nstep_val->kind() != ATOM) {
		syserr("usage : <ODE [VAR] (VAR INIT END STEP) PRED>\n");
		return 0;
	}
	if (!((Atom*)nstep_val)->toFloat(step_val)) {
		syserr("usage : <ODE [VAR] (VAR INIT END STEP) PRED>\n");
		return 0;
	}
	
	Node* val = Nil;

	long double ldf = 0.0;
	Node* env = Nil->Cons(Nil);
	Context *cx2;
	Node* save_integral = Nil;
	int lp = 0;
	
	long double step = (to_val - init_val)/step_val;
#ifdef __CYGWIN__
        step += 0.9L;
#endif /* __CYGWIN__ */
	
	for (ldf = init_val; step >= 0.0L; ldf += step_val, step-=1.0L) {
		cx2 = new Context(module, cx->modulename);
		cx2->selfname = cx->selfname;
		cx2->ioin = cx->ioin;
		cx2->ioout = cx->ioout;
		cx2->tokenflag = cx->tokenflag;
		cx2->token = cx->token;

		cx2->ode = Cons(mka(ldf), Cons(mka(step_val), Nil));
		cx2->integral = save_integral;

		env = Nil->Cons(Nil);
		SetEnv(env, nvar);
		((Undef*)nvar)->Set(mka(ldf));

		cxpush(cx2, gl);
		cxpush(cx2, glargs);
		cxpush(cx2, env);
		cxpush(cx2, nvar);
		cxpush(cx2, goalscar);
		cxpush(cx2, val);
		cxpush(cx2, nresult);
				
		int r = Unify(cx2, gl, module);

		if (++lp == 100000) {
			GC();
			lp = 0;
		}

		if (r == 1) {
			if (cx->tokenflag) cx->token = cx2->token;

			Node*	res;
			if (gl->Car()->Car()->Eq(mka("unify"))
					|| gl->Car()->Car()->Eq(mka("obj"))) {
				res = gl->Car()->Cdr()->Cdr()->Car()->Cdr()->Car()->Val();
			} else {
				res = gl->Car()->Cdr()->Car()->Val();
			}
			
			val = Append(val, MkList(res));
		}

		cxpop(cx2);
		cxpop(cx2);
		cxpop(cx2);
		cxpop(cx2);
		cxpop(cx2);
		cxpop(cx2);
		cxpop(cx2);

		save_integral = cx2->integral;
				
		cx2->Clear();
		
		UnsetEnv(env);
		
		delete cx2;
		cx2 = 0;

		if (r != 1) {
			BreakFlag = 0;
			return r;
		}

		if (BreakFlag) {
			BreakFlag = 0;
			break;
		}
		if (ThrowFlag) {
			break;
		}
		if (ldf >= to_val) {
			break;
		}
	}		

	if (nresult->kind() == UNDEF) {
		Node* env = Nil->Cons(Nil);

		SetEnv(env, nresult);
		((Undef*)(nresult->Val()))->Set(val);

		PushStack(cx, Nil, Nil, env);
	}
			
	return 1;
}


int Integral(Context* cx, Node* goalscar, List* module)
{
	long double output;
	long double constant = 0.0;
	
	Node* g = goalscar->Cdr()->Val();
	int ll = ListLength(g);
	if ((ll != 2) && (ll != 3)) {
		syserr("usage : <integral VAR VAL [constant]> \n");
		return 0;
	}

	Node* nvar = g->Car()->Val();
	int rn;

	if ((rn = FuncArg(cx, nvar, goalscar, module)) <= 0) {
		syserr("integral: failed in the evaluation of the argument. \n");
		return 0;
	}

	if (nvar->kind() != UNDEF) {
		syserr("usage : <integral VAR VAL [constant]> \n");
		return 0;
	}

	Node* nval = g->Cdr()->Car()->Val();

	if ((rn = FuncArg(cx, nval, goalscar, module)) <= 0) {
		syserr("integral: failed in the evaluation of the argument. \n");
		return 0;
	}

	if (nval->kind() != ATOM) {
		syserr("usage : <integral VAR VAL [constant]> \n");
		return 0;
	}

	if (ll == 3) {
		Node* nconstant = g->Cdr()->Cdr()->Car()->Val();

		if ((rn = FuncArg(cx, nconstant, goalscar, module)) <= 0) {
			syserr("integral: failed in the evaluation of the argument. \n");
			return 0;
		}

		if (nconstant->kind() != ATOM) {
			syserr("usage : <integral VAR VAL [constant]> \n");
			return 0;
		}

		((Atom*)nconstant)->toFloat(constant);
	}
	
	long double t=0, val=0, delta=0, ig=0;
	
	((Atom*)nval)->toFloat(val);

	if (cx->ode != Nil) {
		int r1 = ((Atom*)(cx->ode->Car()))->toFloat(t);
		int r2 = ((Atom*)(cx->ode->Cdr()->Car()))->toFloat(delta);
	} else {
		syserr("integral is not executed without ODE. \n");
		return 0;
	}

	Node* id = goalscar->Car();
	Node* ni = GetIntegralCx(cx, id);
	long long flg=0;
	if (ni != Nil) {
		((Atom*)ni->Cdr()->Car())->toInt(flg);
		((Atom*)ni->Cdr()->Cdr()->Car())->toFloat(ig);
	}

	output = Simpson(t, val, delta, flg, ig, constant);

	SetIntegralCx(cx, id, MkList(mka(flg), mka(ig)));
	
	Node* env = Nil->Cons(Nil);
	SetEnv(env, nvar);
	((Undef*)nvar)->Set(mka(output));

	PushStack(cx, Nil, Nil, env);
	
	return 1;
}

int ODEdiff(Context* cx, Node* goalscar, List* module)
{
	long double output;
	long double constant = 0.0;
	
	Node* g = goalscar->Cdr()->Val();
	int ll = ListLength(g);
	if ((ll != 2) && (ll != 3)) {
		syserr("usage : <ODEdiff VAR VAL [constant]> \n");
		return 0;
	}

	Node* nvar = g->Car()->Val();
	int rn;

	if ((rn = FuncArg(cx, nvar, goalscar, module)) <= 0) {
		syserr("ODEdiff: failed in the evaluation of the argument. \n");
		return 0;
	}

	if (nvar->kind() != UNDEF) {
		syserr("usage : <ODEdiff VAR VAL [constant]> \n");
		return 0;
	}

	Node* nval = g->Cdr()->Car()->Val();

	if ((rn = FuncArg(cx, nval, goalscar, module)) <= 0) {
		syserr("ODEdiff: failed in the evaluation of the argument. \n");
		return 0;
	}

	if (ll == 3) {
		Node* nconstant = g->Cdr()->Cdr()->Car()->Val();

		if ((rn = FuncArg(cx, nconstant, goalscar, module)) <= 0) {
			syserr("ODEdiff: failed in the evaluation of the argument. \n");
			return 0;
		}

		if (nconstant->kind() != ATOM) {
			syserr("usage : <ODEdiff VAR VAL [constant]> \n");
			return 0;
		}

		((Atom*)nconstant)->toFloat(constant);
	}
	
	if (nval->kind() != ATOM) {
		syserr("usage : <ODEdiff VAR VAL [constant]> \n");
		return 0;
	}

	long double t=0, val=0, delta=0, ig=0;
	
	((Atom*)nval)->toFloat(val);

	if (cx->ode != Nil) {
		int r1 = ((Atom*)(cx->ode->Car()))->toFloat(t);
		int r2 = ((Atom*)(cx->ode->Cdr()->Car()))->toFloat(delta);
	} else {
		syserr("ODEdiff is not executed without ODE. \n");
		return 0;
	}

	Node* id = goalscar->Car();
	Node* ni = GetIntegralCx(cx, id);
	long long flg=0;
	if (ni != Nil) {
		((Atom*)ni->Cdr()->Car())->toInt(flg);
		((Atom*)ni->Cdr()->Cdr()->Car())->toFloat(ig);
	}

	if (flg == 0) {
		output = constant;
		flg++;
	} else {
		output = (val - ig) / delta;
	}
	ig = val;

	SetIntegralCx(cx, id, MkList(mka(flg), mka(ig)));
	
	Node* env = Nil->Cons(Nil);
	SetEnv(env, nvar);
	((Undef*)nvar)->Set(mka(output));

	PushStack(cx, Nil, Nil, env);
	
	return 1;
}

int ODEprint(Context* cx, Node* goalscar, List* module)
{
	Node* g = goalscar->Cdr()->Val();
	int ll = ListLength(g);
	if (ll < 1) {
		syserr("usage : <ODEprint count print-list> \n");
		return 0;
	}

	Node* ncount = g->Car()->Val();
	int rn;

	if ((rn = FuncArg(cx, ncount, goalscar, module)) <= 0) {
		syserr("ODEprint: failed in the evaluation of the argument. \n");
		return 0;
	}

	if (ncount->kind() != ATOM) {
		syserr("usage : <ODEprint count print-list> \n");
		return 0;
	}

	long long count = 0;
	if (((Atom*)ncount)->toInt(count) != 1) { 
		syserr("ODEprint: count is not a number. ");
		return 0;
	}
	
	if (count <= 0) {
		syserr("ODEprint: count is negative or zero. ");
		return 0;
	}
	
	Node* nlist = g->Cdr()->Val();

	if ((rn = FuncArg(cx, nlist, goalscar, module)) <= 0) {
		syserr("ODEprint: failed in the evaluation of the argument. \n");
		return 0;
	}

	if (cx->ode == Nil) {
		syserr("ODEprint is not executed without ODE. \n");
		return 0;
	}

	Node* id = goalscar->Car();
	Node* ni = GetIntegralCx(cx, id);
	long long flg=0;
	if (ni != Nil) {
		((Atom*)ni->Cdr()->Car())->toInt(flg);
	}

	if (flg == 0) {
		nlist->printcdr(cx->ioout);
		fprintf(cx->ioout, "\n");
		flg++;
	} else if (flg == count-1) {
		flg = 0;
	} else {
		flg++;
	}

	SetIntegralCx(cx, id, MkList(mka(flg)));
	
	return 1;
}


int ODEprintf(Context* cx, Node* goalscar, List* module)
{
	Node* g = goalscar->Cdr()->Val();
	int ll = ListLength(g);
	if (ll < 1) {
		syserr("usage : <ODEprintf count print-list> \n");
		return 0;
	}

	Node* ncount = g->Car()->Val();
	int rn;

	if ((rn = FuncArg(cx, ncount, goalscar, module)) <= 0) {
		syserr("ODEprintf: failed in the evaluation of the argument. \n");
		return 0;
	}

	if (ncount->kind() != ATOM) {
		syserr("usage : <ODEprintf count print-list> \n");
		return 0;
	}

	long long count = 0;
	if (((Atom*)ncount)->toInt(count) != 1) { 
		syserr("ODEprintf: count is not a number. ");
		return 0;
	}
	
	if (count <= 0) {
		syserr("ODEprintf: count is negative or zero. ");
		return 0;
	}
	
	Node* nlist = g->Cdr()->Val();

	if ((rn = FuncArg(cx, nlist, goalscar, module)) <= 0) {
		syserr("ODEprintf: failed in the evaluation of the argument. \n");
		return 0;
	}

	if (cx->ode == Nil) {
		syserr("ODEprintf is not executed without ODE. \n");
		return 0;
	}

	Node* id = goalscar->Car();
	Node* ni = GetIntegralCx(cx, id);
	long long flg=0;
	if (ni != Nil) {
		((Atom*)ni->Cdr()->Car())->toInt(flg);
	}

	if (flg == 0) {
		for (; nlist->kind() != ATOM; nlist=nlist->Cdr()) {
			nlist->Car()->print(cx->ioout);
		}
		flg++;
	} else if (flg == count-1) {
		flg = 0;
	} else {
		flg++;
	}

	SetIntegralCx(cx, id, MkList(mka(flg)));
	
	return 1;
}


int DoODEcmd(Context* cx, Node* gl, List* module)
{

	Node* env = Nil->Cons(Nil);
	Context *cx2;

	cx2 = new Context(module, cx->modulename);
	cx2->selfname = cx->selfname;
	cx2->ioin = cx->ioin;
	cx2->ioout = cx->ioout;
	cx2->tokenflag = cx->tokenflag;
	cx2->token = cx->token;

	cx2->ode = cx->ode;
	cx2->integral = cx->integral;

	cxpush(cx2, gl);
	cxpush(cx2, env);
				
	int r = Unify(cx2, gl, module);

		
	if (r == 1) {
		if (cx->tokenflag) cx->token = cx2->token;
	}

	cxpop(cx2);
	cxpop(cx2);


	cx2->Clear();
		
	UnsetEnv(env);
		
	delete cx2;
	cx2 = 0;

	if (r != 1) {
		BreakFlag = 0;
		return r;
	}

	return 1;
}


int ODEcmd(Context* cx, Node* goalscar, List* module)
{
	Node* g = goalscar->Cdr()->Val();
	int ll = ListLength(g);
	if (ll < 1) {
		syserr("usage : <ODEcmd count cmd-list> \n");
		return 0;
	}

	Node* ncount = g->Car()->Val();
	int rn;

	if ((rn = FuncArg(cx, ncount, goalscar, module)) <= 0) {
		syserr("ODEcmd: failed in the evaluation of the argument. \n");
		return 0;
	}

	if (ncount->kind() != ATOM) {
		syserr("usage : <ODEcmd count cmd-list> \n");
		return 0;
	}

	long long count = 0;
	if (((Atom*)ncount)->toInt(count) != 1) { 
		syserr("ODEcmd: count is not a number. ");
		return 0;
	}
	
	if (count <= 0) {
		syserr("ODEcmd: count is negative or zero. ");
		return 0;
	}
	
	Node* nlist = g->Cdr()->Val();

	if (cx->ode == Nil) {
		syserr("ODEcmd is not executed without ODE. \n");
		return 0;
	}

	Node* id = goalscar->Car();
	Node* ni = GetIntegralCx(cx, id);
	long long flg=0;
	if (ni != Nil) {
		((Atom*)ni->Cdr()->Car())->toInt(flg);
	}

	if (flg == 0) {
		if ((rn = DoODEcmd(cx, nlist, module)) <= 0) {
			syserr("ODEcmd: failed in the evaluation of the ODEcmd. \n");
			return 0;
		}
		flg++;
	} else if (flg == count-1) {
		flg = 0;
	} else {
		flg++;
	}

	SetIntegralCx(cx, id, MkList(mka(flg)));
	
	return 1;
}


