/*******************************************************************
Copyright (C) 2016 AMPL Optimization, Inc.; written by David M. Gay.

Permission to use, copy, modify, and distribute this software and its
documentation for any purpose and without fee is hereby granted,
provided that the above copyright notice appear in all copies and that
both that the copyright notice and this permission notice and warranty
disclaimer appear in supporting documentation.

The author and AMPL Optimization, Inc. disclaim all warranties with
regard to this software, including all implied warranties of
merchantability and fitness.  In no event shall the author be liable
for any special, indirect or consequential damages or any damages
whatsoever resulting from loss of use, data or profits, whether in an
action of contract, negligence or other tortious action, arising out
of or in connection with the use or performance of this software.
*******************************************************************/

#include "jac2dim.h"
#include "opnos.hd"

#ifdef __cplusplus
extern "C" {
#endif

extern real con2ival_nomap_ASL(ASL*, int, real*, fint*);
extern void con2grd_nomap_ASL(ASL*, int, real *, real*, fint*);

#define Egulp 400

 static int com11, n_com1, ncom_togo, nvar0, nvinc, nvref;
 static int *vrefnext, *vrefx;
 static expr_if *if2list_end;
 static expr_va *varg2list_end;

#define NFHASH 23

 static derp *last_d;
 static expr *last_e;

#define efunc efunc2

 static relo *relolist, *relo2list;
 static expr_if *iflist, *if2list;
 static expr_va *varglist, *varg2list;
 static real one = 1.;
 static int nderp;
#undef nzc
 static int amax1, imap_len, k_seen, lasta, lasta0, lasta00, lastj,
	max_var, nocopy, nv0, nv01, nv0b, nv0c, nv1, nzc, nzclim;
 static int co_first = 1;
 static int *imap, *zc, *zci;

 static expr *(*holread) ANSI((EdRead*));
 extern real f_OPCPOW ANSI((expr* A_ASL));

 static ASL_fgh *asl;

 static void
ed_reset(ASL *a)
{

	a->i.memLast = a->i.memNext = 0;
	k_seen = 0;
	nvref = 0;
	vrefnext = vrefx = 0;
	relolist = relo2list = 0;
	last_d = 0;
	iflist = if2list = if2list_end = 0;
	varglist = varg2list = varg2list_end = 0;
	imap = 0;
	amax1 = imap_len = lastj = nocopy = 0;
	com11 = lasta = lasta0 = lasta00 = n_com1 = nderp = 0;
	last_e = 0;
	co_first = 1;
	}

 static void
fscream(EdRead *R, const char *name, int nargs, const char *kind)
{
	scream(R, ASL_readerr_argerr,
		"line %ld: attempt to call %s with %d %sargs\n",
		R->Line, name, nargs, kind);
	}

 static void
sorry_CLP(EdRead *R, const char *what)
{
	fprintf(Stderr,
		"Sorry, %s cannot handle %s.\n",
		progname ? progname : "", what);
	exit_ASL(R,ASL_readerr_CLP);
	}

#ifdef Double_Align
#define memadj(x) x
#else
#define memadj(x) (((x) + (sizeof(long)-1)) & ~(sizeof(long)-1))
#endif

 extern char* f_OPHOL ANSI((expr* A_ASL));
 extern efunc f_OPPLTERM, f_OPVARVAL, f_OPFUNCALL;
 extern sfunc f_OPIFSYM;

 static void
new_derp(int a, int b, real *c)
{
	derp *d;
	if (a == nv1)
		return;
	nderp++;
	d = (derp *)mem(sizeof(derp));
	d->next = last_d;
	last_d = d;
	d->a.i = a;
	d->b.i = b;
	d->c.rp = c;
	}

 static derp *
new_relo(expr *e, derp *Dnext, int *ap)
{
	relo *r;
	derp *d;

	if (last_d != Dnext) {
		*ap = e->a;
		for(d = last_d; d->next != Dnext; d = d->next);
		d->next = 0;
		}
	else {
		last_d = 0;
		new_derp(e->a, *ap = lasta++, &one);
		}
	if (!last_d)
		return 0;
	r = (relo *)mem(sizeof(relo));
	r->next = relolist;
	r->next2 = relo2list;
	relo2list = relolist = r;
	r->D = r->Dcond = last_d;
	r->Dnext = Dnext;
	return r->D;
	}

 static relo *
new_relo1(derp *Dnext)
{
	relo *r;

	r = (relo *)mem(sizeof(relo));
	r->next = relolist;
	relolist = r;
	r->D = 0;
	r->Dnext = Dnext;
	return r;
	}

 static expr *
new_expr(int opcode, expr *L, expr *R, int deriv)
{
	expr *rv;
	extern efunc f_OP1POW, f_OP2POW, f_OPPOW;
	efunc *o;
	int L1, R1, i;

	o = r_ops[opcode];
	if (o == f_OPPOW) {
		if (R->op == f_OPNUM)
			if (((expr_n *)R)->v == 2.) {
				o = f_OP2POW;
				R = 0;
				}
			else
				o = f_OP1POW;
		else if (L->op == f_OPNUM)
			o = f_OPCPOW;
		}
	rv = (expr *)mem(sizeof(expr));
	rv->op = o;
	rv->L.e = L;
	rv->R.e = R;
	if (deriv) {
		L1 = L && L->op != f_OPNUM;
		R1 = R && R->op != f_OPNUM;
		if (L1 | R1) {
			rv->a = lasta++;
			if (L1)
				new_derp(L->a, rv->a, &rv->dL);
			if (R1)
				new_derp(R->a, rv->a, &rv->dR);
			rv->bak = last_e;
			last_e = rv;
			if (R)
				rv->dL2 = rv->dLR = rv->dR2 = 0;
			else if (o == f_OP2POW)
				rv->dL2 = 2;
			else
				rv->dL2 = 0;
			if (L1)
			    if (R1)
				switch(opcode) {
					case PLUS:	i = Hv_plusLR;	break;
					case MINUS:	i = Hv_minusLR;	break;
					case MULT:	i = Hv_timesLR;	break;
					case DIV:	i = Hv_divLR;	break;
					default:	i = Hv_binaryLR;
					}
			    else switch(opcode) {
					case PLUS:
					case MINUS:	i = Hv_plusL;	break;
					case MULT:	i = Hv_timesL;	break;
					case UMINUS:	i = Hv_negate;	break;
					default:	i = Hv_unary;
					}
			else
				switch(opcode) {
					case PLUS:	i = Hv_plusR;	break;
					case MINUS:	i = Hv_minusR;	break;
					case MULT:	i = Hv_timesR;	break;
					default:	i = Hv_binaryR;
					}
			rv->dO.i = i;
			}
		}
	return rv;
	}

 static char op_type[] = {
#include "op_type.hd"
	};

 static expr *
eread(EdRead *R, int deriv)
{
	arglist *al;
	argpair *ap, *da, *sap;
	char **sa;
	char *dig;
	de *d;
	derp *dsave;
	efunc *op;
	expr *L, *arg, **args, **args1, **argse, *rv;
	expr_f *rvf;
	expr_if *rvif;
	expr_n *rvn;
	expr_va *rva;
	func_info *fi;
	int *at, *nn, *nn0;
	int a0, a1, i, i1, j, j1, k, kd, kd2, ks, numargs, symargs;
	int (*Xscanf)(EdRead*, const char*, ...);
	long L1;
	plterm *p;
	real *b, **fh, *hes, r, *ra;
	unsigned int Ls;
	static real dvalue[] = {
#include "dvalue.hd"
		};

	Xscanf = xscanf;
	switch(edag_peek(R)) {
		case 'f':
			if (Xscanf(R, "%d %d", &i, &j) != 2)
				badline(R);
			fi = funcs[i];
			if (fi->nargs >= 0) {
				if (fi->nargs != j) {
 bad_nargs:
					fscream(R, fi->name, j, "");
					}
				}
			else if (-(1+fi->nargs) > j)
				goto bad_nargs;
			rvf = (expr_f *)mem(sizeof(expr_f)
					+ (j-1)*sizeof(expr *));
			rvf->op = f_OPFUNCALL;
			rvf->fi = fi;
			rvf->dO.i = Hv_func;
			args = rvf->args;
			argse = args + j;
			k = ks = symargs = numargs = 0;
			while(args < argse) {
				arg = *args++ = eread(R, deriv);
				if ((op = arg->op) == (efunc*)f_OPHOL)
					symargs++;
				else if (op == (efunc*)f_OPIFSYM)
					ks++;
				else {
					numargs++;
					if (op != f_OPNUM)
						k++;
					}
				}
			symargs += ks;
			if (symargs && !(fi->ftype & 1))
				fscream(R, fi->name, symargs, "symbolic ");
			if (deriv) {
				kd2 = (kd = k) ? numargs*(numargs+1) >> 1 : 0;
				rvf->a = lasta++;
				}
			else {
				kd = kd2 = 0;
				rvf->a = nv1;
				}
			Ls = sizeof(arglist)
					+ (k + kd + ks)*sizeof(argpair)
					+ kd*kd*sizeof(real *)
					+ (numargs+kd+kd2)*sizeof(real)
					+ symargs*sizeof(char *)
					+ j*sizeof(int);
			if (kd)
				Ls += numargs*sizeof(real);
			ra = (real *)mem(Ls);
			dig = 0;
			if (kd < numargs && kd)
				dig = (char*)mem(numargs);
			b = kd ? ra + numargs : ra;
			hes = b + numargs;
			al = rvf->al = (arglist *)(hes + kd2);
			al->n = numargs + symargs;
			al->nr = numargs;
			al->ra = ra;
			if (kd) {
				al->derivs = b;
				al->hes = hes;
				memset(b, 0, (numargs + kd2)*sizeof(real));
				}
			else
				al->derivs = al->hes = 0;
			al->dig = dig;
			al->funcinfo = fi->funcinfo;
			al->AE = asl->i.ae;
			al->sa = (Const char**)(sa = (char **)(al + 1));
			ap = rvf->ap = (argpair *)(sa + symargs);
			rvf->da = da = ap + k;
			sap = rvf->sap = da + kd;
			rvf->fh = fh = (real **)(sap + ks);
			at = al->at = (int *)(fh + kd*kd);
			symargs = numargs = i = 0;
			nn = nn0 = (int *)b;
			for(args = rvf->args; args < argse; at++) {
				arg = *args++;
				if ((op = arg->op) == (efunc*)f_OPHOL) {
					*at = --symargs;
					*sa++ = ((expr_h *)arg)->sym;
					}
				else if (op == (efunc*)f_OPIFSYM) {
					*at = --symargs;
					sap->e = arg;
					(sap++)->u.s = sa++;
					}
				else {
					*at = numargs++;
					j = 1;
					if (op == f_OPNUM)
						*ra = ((expr_n *)arg)->v;
					else  {
						ap->e = arg;
						(ap++)->u.v = ra;
						if (kd) {
							j = 0;
							new_derp(arg->a,
								rvf->a, b);
							*b = 0;
							da->e = arg;
							(da++)->u.v = b;
							*nn++ = i;
							}
						}
					if (dig)
						*dig++ = j;
					b++;
					ra++;
					i++;
					}
				}
			rvf->ape = ap;
			rvf->sape = sap;
			rvf->dae = da;
			if (kd) {
				rvf->bak = last_e;
				last_e = (expr  *)rvf;
				rvf->dO.i = Hv_func;
				for(i1 = 0; i1 < kd; i1++) {
					i = nn0[i1];
					for(j1 = 0; j1 < kd; j1++) {
						j = nn0[j1];
						*fh++ = &hes[i >= j ? (i*(i+1)>>1)+j
								    : (j*(j+1)>>1)+i];
						}
					}
				}
			return (expr *)rvf;

		case 'h':
			return holread(R);

		case 's':
			if (Xscanf(R, "%hd", &L1) != 1)
				badline(R);
			r = L1;
			goto have_r;

		case 'l':
			if (Xscanf(R, "%ld", &L1) != 1)
				badline(R);
			r = L1;
			goto have_r;

		case 'n':
			if (Xscanf(R, "%lf", &r) != 1)
				badline(R);
 have_r:
			rvn = (expr_n *)mem(sizeof(expr_n));
			rvn->op = f_OPNUM_ASL;
			rvn->v = r;
			return (expr *)rvn;

		case 'o':
			break;

		case 'v':
			if (Xscanf(R,"%d",&k) != 1 || k < 0)
				badline(R);
			if (k >= nvar0)
				k += nvinc;
			if (k > max_var)
				badline(R);
			if (k < nv01 && deriv && !zc[k]++)
				zci[nzc++] = k;
			return (expr *)(var_e + k);

		default:
			badline(R);
		}

	if (Xscanf(R, asl->i.opfmt, &k) != 1 || k < 0 || k >= sizeof(op_type))
		badline(R);
	switch(op_type[k]) {

		case 1:	/* unary */
			rv = new_expr(k, eread(R, deriv), 0, deriv);
			rv->dL = dvalue[k];	/* for UMINUS, FLOOR, CEIL */
			return rv;

		case 2:	/* binary */
			if (dvalue[k] == 11)
				deriv = 0;
			L = eread(R, deriv);
			rv = new_expr(k, L, eread(R, deriv), deriv);
			rv->dL = 1.;
			rv->dR = dvalue[k];	/* for PLUS, MINUS, REM */
			return rv;

		case 3:	/* vararg (min, max) */
			i = -1;
			Xscanf(R, "%d", &i);
			if (i <= 0)
				badline(R);
			rva = (expr_va *)mem(sizeof(expr_va));
			rva->op = r_ops[k];
			rva->L.d = d =
				(de *)mem(i*sizeof(de) + sizeof(expr *));
			rva->next = varglist;
			varglist = varg2list = rva;
			if (!last_d) {
				new_derp(lasta, lasta, &one);
				lasta++;
				}
			rva->d0 = dsave = last_d;
			rva->bak = last_e;
			a0 = a1 = lasta;
			for(j = 0; i > 0; i--, d++) {
				last_d = dsave;
				last_e = 0;
				d->e = L = eread(R, deriv);
				d->ee = last_e;
				if (L->op == f_OPNUM || L->a == nv1) {
					d->d = dsave;
					d->dv.i = nv1;
					}
				else {
					if (deriv)
						d->d = new_relo(L, dsave,
								&d->dv.i);
					if (a1 < lasta)
						a1 = lasta;
					lasta = a0;
					}
				}
			d->e = 0;	/* sentinnel expr * */
			rva->a = lasta = a1;
			last_d = dsave;
			if (deriv) {
				new_derp(0, lasta++, &one);
				/* f_MINLIST or f_MAXLIST will replace the 0 */
				rva->R.D = last_d;
				nocopy = 1;
				last_e = (expr *)rva;
				rva->dO.i = Hv_vararg;
				}
			else {
				rva->R.D = 0;
				last_e = rva->bak;
				}
			return (expr *)rva;

		case 4: /* piece-wise linear */
			i = -1;
			Xscanf(R, "%d", &i);
			if (i <= 1)
				badline(R);
			j = 2*i - 1;
			p = (plterm *)mem(sizeof(plterm) + (j-1)*sizeof(real));
			p->n = i;
			b = p->bs;
			do {
				switch(edag_peek(R)) {
					case 's':
						if (Xscanf(R,"%hd",&L1) != 1)
							badline(R);
						r = L1;
						break;
					case 'l':
						if (Xscanf(R,"%ld",&L1) != 1)
							badline(R);
						r = L1;
						break;
					case 'n':
						if (Xscanf(R,"%lf",&r) == 1)
							break;
					default:
						badline(R);
					}
				*b++ = r;
				}
				while(--j > 0);
			if (b[-2] <= 0.)
				p->z = 2*i - 2;
			else {
				b = p->bs + 1;
				while(*b <= 0.)
					b += 2;
				p->z = (b - p->bs) - 1;
				}
			if (edag_peek(R) != 'v'
			 || Xscanf(R, "%d", &k) != 1
			 || k < 0 || k >= max_var)
				badline(R);
			if (k >= nvar0)
				k += nvinc;
			rv = (expr *)mem(sizeof(expr));
			rv->op = f_OPPLTERM;
			rv->L.p = p;
			rv->R.e = (expr *)(var_e + k);
			if (deriv) {
				new_derp(k, rv->a = lasta++, &rv->dL);
				rv->bak = last_e;
				last_e = rv;
				rv->dO.i = Hv_plterm;
				}
			return rv;

		case 5: /* if */
			rvif = (expr_if *)mem(sizeof(expr_if));
			rvif->op = r_ops[k];
			rvif->next = iflist;
			iflist = if2list = rvif;
			if (!last_d) {
				new_derp(lasta, lasta, &one);
				lasta++;
				}
			rvif->d0 = dsave = last_d;
			rvif->bak = last_e;
			rvif->e = eread(R, 0);
			last_e = 0;
			a0 = lasta;
			rvif->T = L = eread(R, deriv);
			rvif->Te = last_e;
			j = 0;
			if (L->op == f_OPNUM) {
				rvif->dT = dsave;
				rvif->Tv.i = nv1;
				}
			else if ((j = deriv))
				rvif->dT = new_relo(L, dsave, &rvif->Tv.i);
			a1 = lasta;
			lasta = a0;
			last_d = dsave;
			last_e = 0;
			rvif->F = L = eread(R, deriv);
			rvif->Fe = last_e;
			if (L->op == f_OPNUM) {
				rvif->dF = dsave;
				rvif->Fv.i = nv1;
				}
			else if (j)
				rvif->dF = new_relo(L, dsave, &rvif->Fv.i);
			if (lasta < a1)
				lasta = a1;
			last_d = dsave;
			if (j) {
				new_derp(0, rvif->a = lasta++, &one);
				rvif->D = last_d;
				nocopy = 1;
				last_e = (expr *)rvif;
				rvif->dO.i = Hv_if;
				}
			else {
				rvif->a = nv1;
				rvif->D = 0;
				last_e = rvif->bak;
				}
			return (expr *)rvif;

		case 11: /* OPCOUNT */
			deriv = 0;
			/* no break */
		case 6: /* sumlist */
			i = 0;
			Xscanf(R, "%d", &i);
			if (i <= 2 && (op_type[k] == 6 || i < 1))
				badline(R);
			rv = (expr *)mem(sizeof(expr) - sizeof(real)
					+ (deriv ? i+i+1 : i)*sizeof(expr *));
			rv->op = r_ops[k];
			rv->a = deriv ? lasta++ : nv1;
			rv->L.ep = args = (expr **)&rv->dR;
			if (deriv) {
				j = 0;
				do {
					*args++ = L = eread(R, deriv);
					if (L->op != f_OPNUM) {
						new_derp(L->a, rv->a, &one);
						j++;
						}
					}
					while(--i > 0);
				if (j) {
					rv->bak = last_e;
					last_e = rv;
					rv->dO.i = Hv_sumlist;
					argse = args;
					args1 = rv->L.ep;
					while(args1 < args) {
						arg = *args1++;
						if (arg->op != f_OPNUM)
							*argse++ = arg;
						}
					*argse = 0;
					}
				}
			else do
				*args++ = eread(R, deriv);
				while(--i > 0);
			rv->R.ep = args;
			return rv;
			}
	badline(R);
	return 0;
	}

 static list *
new_list(list *nxt)
{
	list *rv = (list *)mem(sizeof(list));
	rv->next = nxt;
	return rv;
	}

 static list *
crefs(void)
{
	int i;
	list *rv = 0;

	while(nzc > 0) {
		if ((i = zci[--nzc]) >= nv0) {
			rv = new_list(rv);
			rv->item.i = i;
			}
		zc[i] = 0;
		}
	return rv;
	}

 static funnel *
funnelfix(funnel *f)
{
	cexp *ce;
	funnel *fnext, *fprev;

	for(fprev = 0; f; f = fnext) {
		fnext = f->next;
		f->next = fprev;
		fprev = f;
		ce = f->ce;
		ce->z.i = ce->d->b.i;
		}
	return fprev;
	}

 static derp *
derpadjust(derp *d0, int a, derp *dnext)
{
	derp *d, *d1;
	int *r, *re;
	relo *rl;
	expr_if *il, *ile;
	expr_va *vl, *vle;
	de *de1;

	if (!(d = d0))
		return dnext;
	r = imap + lasta0;
	re = imap + lasta;
	while(r < re)
		*r++ = a++;
	if (amax < a)
		amax = a;
	r = imap;
	for(;; d = d1) {
		d->a.i = r[d->a.i];
		d->b.i = r[d->b.i];
		if (!(d1 = d->next))
			break;
		}
	d->next = dnext;
	if ((rl = relo2list)) {
		relo2list = 0;
		do {
			d = rl->Dcond;
			do {
				d->a.i = r[d->a.i];
				d->b.i = r[d->b.i];
				}
				while((d = d->next));
			}
			while((rl = rl->next2));
		}
	if (if2list != if2list_end) {
		ile = if2list_end;
		if2list_end = il = if2list;
		do {
			il->Tv.i = r[il->Tv.i];
			il->Fv.i = r[il->Fv.i];
			}
			while((il = il->next) != ile);
		}
	if (varg2list != varg2list_end) {
		vle = varg2list_end;
		varg2list_end = vl = varg2list;
		do {
			for(de1 = vl->L.d; de1->e; de1++)
				de1->dv.i = r[de1->dv.i];
			}
			while((vl = vl->next) != vle);
		}
	return d0;
	}

 static derp *
derpcopy(cexp *ce, derp *dnext)
{
	derp	*d, *dprev;
	int	*map;
	derp		d00;

	if (!(d = ce->d))
		return dnext;
	map = imap;
	for(dprev = &d00; d; d = d->next) {
		new_derp(map[d->a.i], map[d->b.i], d->c.rp);
		dprev = dprev->next = last_d;
		}
	dprev->next = dnext;
	return d00.next;
	}

 static void
imap_alloc(void)
{
	int i, *r, *re;

	if (imap) {
		imap_len += lasta;
		imap = (int *)Realloc(imap, imap_len*Sizeof(int));
		return;
		}
	imap_len = amax1 > lasta ? amax1 : lasta;
	imap_len += 100;
	r = imap = (int *)Malloc(imap_len*Sizeof(int));
	for(i = 0, re = r + nv1+1; r < re;)
		*r++ = i++;
	}

 static int
compar(const void *a, const void *b, void *v)
{
	Not_Used(v);
	return *(int*)a - *(int*)b;
	}

 static void
comsubs(int alen, cde *d)
{
	list *L;
	int a, i, j, k;
	int *r, *re;
	cexp *ce;
	derp *D, *dnext;
	relo *R;

	D = last_d;
	a = lasta00;
	dnext = 0;
	R = 0;
	for(i = j = 0; i < nzc; i++)
		if ((k = zci[i]) >= nv0)
			zci[j++] = k;
		else
			zc[k] = 0;
	if ((nzc = j)) {
		for(i = 0; i < nzc; i++)
			for(L = cexps[zci[i]-nv0].cref; L; L = L->next)
				if (!zc[L->item.i]++)
					zci[nzc++] = L->item.i;
		if (nzc > 1) {
			if (nzc < nzclim)
				qsortv(zci, nzc, sizeof(int), compar, NULL);
			else for(i = nv0, j = 0; i < max_var; i++)
				if (zc[i])
					zci[j++] = i;
			}
		if (nzc > 0) {
			R = new_relo1(dnext);
			i = 0;
			do {
				j = zci[i];
				zc[j] = 0;
				ce = &cexps[j - nv0];
				if (ce->funneled)
					imap[var_e[j].a] = a++;
				else {
					r = imap + ce->z.i;
					re = r + ce->zlen;
					while(r < re)
						*r++ = a++;
					}
				dnext = R->D = derpcopy(ce, R->D);
				}
				while(++i < nzc);
			nzc = 0;
			}
		}
	if (D || R) {
		if (!R)
			R = new_relo1(dnext);
		D = R->D = derpadjust(D, a, R->D);
		if (d->e->op != f_OPVARVAL)
			d->e->a = imap[d->e->a];
		}
	d->d = D;
	a += alen;
	d->zaplen = (a > lasta00 ? a - nv1 : 0)*sizeof(real);
	if (amax < a)
		amax = a;
	}

 static void
co_read(EdRead *R, cde *d, int wd)
{
	int alen;

	d->com11 = com11;
	d->n_com1 = n_com1;
	com11 += n_com1;
	n_com1 = 0;

	if (amax1 < lasta)
		amax1 = lasta;
	if (co_first) {
		co_first = 0;
		if (imap_len < lasta)
			imap_alloc();
		f_b = funnelfix(f_b);
		f_c = funnelfix(f_c);
		f_o = funnelfix(f_o);
		}
	if (!lastj) {
		lasta = lasta0;
		last_d = 0;
		}
	lastj = 0;
	last_e = 0;
	d->e = eread(R, wd);
	d->ee = last_e;
	alen = lasta - lasta0;
	if (imap_len < lasta)
		imap_alloc();
	comsubs(alen, d);
	}

 static linpart *
linpt_read(EdRead *R, int nlin)
{
	int (*Xscanf)(EdRead*, const char*, ...);
	linpart *L, *rv;

	if (nlin <= 0)
		return 0;
	Xscanf = xscanf;
	L = rv = (linpart *)mem(nlin*sizeof(linpart));
	do {
		if (Xscanf(R, "%d %lf", &L->v.i, &L->fac) != 2)
			badline(R);
		L++;
		}
		while(--nlin > 0);
	return rv;
	}

 static int
funnelkind(cexp *ce, int *ip)
{
	int i, j, k, nzc0, rv;
	int *vr, *vre;

	ce->vref = 0;
	if (!(nzc0 = nzc))
		return 0;
	for(i = k = rv = 0; i < nzc; i++)
		if ((j = zci[i]) < nv0) {
			if (k >= maxfwd)
				goto done;
			vrefx[k++] = j;
			}
		else  {
			if (!(vr = cexps[j-nv0].vref))
				goto done;
			vre = vr + *vr;
			while(++vr <= vre)
				if (!zc[*vr]++)
					zci[nzc++] = *vr;
			}
	if (k >= nvref) {
		nvref = (maxfwd + 1)*(ncom_togo < vrefGulp
					? ncom_togo : vrefGulp);
		vrefnext = (int *)M1alloc(nvref*Sizeof(int));
		}
	vr = ce->vref = vrefnext;
	*vr = k;
	vrefnext += i = k + 1;
	nvref -= i;
	for(i = 0; i < k; i++)
		*++vr = vrefx[i];
	if (nderp > 3*k && !nocopy) {
		*ip = k;
		return 2;
		}
	else {
 done:
		if (nocopy || nderp > 3*nzc0)
			rv = 1;
		}
	while(nzc > nzc0)
		zc[zci[--nzc]] = 0;
	return rv;
	}

 static void
cexp_read(EdRead *R, int k, int nlin)
{
	int a, fk, i, j, la0, na;
	funnel *f, **fp;
	linpart *L, *Le;
	expr *e;
	cplist *cl, *cl0;
	cexp *ce;

	ce = cexps + k - nv0;
	L = ce->L = linpt_read(R, ce->nlin = nlin);
	nocopy = 0;
	last_d = 0;
	last_e = 0;
	ce->z.i = la0 = lasta;
	nderps += nderp;
	nderp = 0;
	e = ce->e = eread(R, 1);
	if (la0 == lasta) {
		a = lasta++;
		if (e->op != f_OPNUM)
			new_derp(e->a, a, &edagread_one);
		}
	else
		a = e->a;
	var_e[k].a = a;
	ce->zlen = lasta - la0;
	for(Le = L + nlin; L < Le; L++) {
		new_derp(i = L->v.i, a, &L->fac);
		if (!zc[i]++)
			zci[nzc++] = i;
		}
	i = 0; /* only needed to shut up an erroneous warning */
	if ((fk = funnelkind(ce, &i))) {
		/* arrange to funnel */
		fp = k < nv0b ? &f_b : k < nv0c ? &f_c : &f_o;
		ce->funneled = f = (funnel *)mem(sizeof(funnel));
		f->next = *fp;
		*fp = f;
		f->ce = ce;
		if (imap_len < lasta)
			imap_alloc();
		if (fk == 1) {
			f->fulld = last_d;
			a = lasta00;
			for(i = nzc; --i >= 0; )
				if ((j = zci[i]) >= nv0)
					imap[var_e[j].a] = a++;
			if ((na = ce->zlen) || a > lasta00)
				na += a - nv1;
			f->fcde.zaplen = na*sizeof(real);
			i = nzc;
			derpadjust(last_d, a, 0);
			}
		else {
			f->fulld = 0;
			f->fcde.e = e;
			comsubs(ce->zlen, &f->fcde);
			memcpy(zci, vrefx, i*sizeof(int));
			}
		last_d = 0;
		cl0 = 0;
		while(i > 0)
			if ((a = var_e[zci[--i]].a) != nv1) {
				new_derp(a, lasta0, 0);
				cl = (cplist *)mem(sizeof(cplist));
				cl->next = cl0;
				cl0 = cl;
				cl->ca.i = imap[last_d->a.i];
				cl->cfa = last_d->c.rp =
					(real *)mem(sizeof(real));
				}
		f->cl = cl0;
		var_e[k].a = lasta0++;
		lasta = lasta0;
		}
	lasta0 = lasta;
	ce->d = last_d;
	ce->ee = last_e;
	ce->cref = crefs();
	--ncom_togo;
	}

 static void
cexp1_read(EdRead *R, int j, int k, int nlin)
{
	linpart *L, *Le;
	cexp1 *ce = cexps1 + (k - nv01);
	expr *e;
	expr_v *v;
	int la0;

	n_com1++;
	L = ce->L = linpt_read(R, ce->nlin = nlin);

	if (!lastj) {
		last_d = 0;
		if (amax1 < lasta)
			amax1 = lasta;
		lasta = lasta0;
		lastj = j;
		}
	last_e = 0;
	la0 = lasta;
	e = ce->e = eread(R, 1);
	ce->ee = last_e;
	v = var_e + k;
	if (la0 == lasta) {
		j = lasta++;
		if (e->op != f_OPNUM)
			new_derp(e->a, j, &edagread_one);
		}
	else
		j = e->a;
	v->a = j;
	v->dO.r = 0;
	for(Le = L + nlin; L < Le; L++)
		new_derp(L->v.i, j, &L->fac);
	}

 static expr *
hvadjust(expr *e)
{
	expr *e0;

	for(e0 = 0; e; e = e->bak) {
		e->fwd = e0;
		e0 = e;
		e->a = e->dO.i;
		}
	return e0;
	}

 static void
co_adjust(cde *d, int n)
{
	cde *de1;

	for(de1 = d + n; d < de1; d++)
		d->ef = hvadjust(d->ee);
	}

 static void
ifadjust(expr_if *e)
{
	for(; e; e = e->next) {
		e->Tv.rp = &adjoints[e->Tv.i];
		e->Fv.rp = &adjoints[e->Fv.i];
		e->Tf = hvadjust(e->Te);
		e->Ff = hvadjust(e->Fe);
		}
	}

 static void
vargadjust(expr_va *e)
{
	de *d;

	for(; e; e = e->next) {
		for(d = e->L.d; d->e; d++) {
			d->dv.rp = &adjoints[d->dv.i];
			d->ef = hvadjust(d->ee);
			}
		}
	}

 static void
funneladj1(funnel *f)
{
	real	*a	= adjoints;
	derp	*d;
	cplist	*cl;

	for(a = adjoints; f; f = f->next) {
		if ((d = f->fulld)) {
			f->fcde.d = d;
			do {
				d->a.rp = &a[d->a.i];
				d->b.rp = &a[d->b.i];
				}
				while((d = d->next));
			}
		for(cl = f->cl; cl; cl = cl->next)
			cl->ca.rp = &a[cl->ca.i];
		}
	}

 static void
funneladjust(VOID)
{
	cexp *c, *ce;
	linpart *L, *Le;
	c = cexps;
	for(ce = c + ncom0; c < ce; c++) {
		if ((L = c->L))
			for(Le = L + c->nlin; L < Le; L++)
				L->v.vp = (Char*)&var_e[L->v.i];
		c->ef = hvadjust(c->ee);
		}

	funneladj1(f_b);
	funneladj1(f_c);
	funneladj1(f_o);
	}

 static void
com1adjust(VOID)
{
	cexp1 *c, *ce;
	linpart *L, *Le;

	for(c = cexps1, ce = c + ncom1; c < ce; c++) {
		for(L = c->L, Le = L + c->nlin; L < Le; L++)
			L->v.vp = (Char*)&var_e[L->v.i];
		c->ef = hvadjust(c->ee);
		}
	}

 static void
adjust_compl_rhs(VOID)
{
	cde *C;
	expr *e;
	int *Cvar, i, j, nc, stride;
	real *L, *U, t, t1;

	L = LUrhs;
	if ((U = Urhsx))
		stride = 1;
	else {
		U = L + 1;
		stride = 2;
		}
	C = con_de;
	Cvar = cvar;
	nc = n_con;
	for(i = nlc; i < nc; i++)
		if (Cvar[i] && (e = C[i].e) && e->op == f_OPNUM
		&& (t = ((expr_n*)e)->v) != 0.) {
			t1 = t;
			if (L[j = stride*i] > negInfinity) {
				L[j] -= t;
				t1 = 0.;
				}
			if (U[j] < Infinity) {
				U[j] -= t;
				t1 = 0.;
				}
			((expr_n*)e)->v = t1;
			}
	}

 static void
adjust(ASL *a0, int flags)
{
	derp *d, **dp;
	real *a = adjoints;
	relo *r;

	for(r = relolist; r; r = r->next) {
		dp = &r->D;
		while((d = *dp)) {
			d->a.rp = &a[d->a.i];
			d->b.rp = &a[d->b.i];
			dp = &d->next;
			}
		*dp = r->Dnext;
		}
	ifadjust(iflist);
	vargadjust(varglist);
	if (ncom0)
		funneladjust();
	com1adjust();
	co_adjust(con_de, n_con);
	co_adjust(obj_de, n_obj);
	if (k_seen) {
		if (!A_vals)
			goff_comp_ASL(a0);
		else if (Fortran)
			colstart_inc_ASL(a0);
		}
	if (n_cc > nlcc && nlc < n_con
	 && !(flags & ASL_no_linear_cc_rhs_adjust))
		adjust_compl_rhs();
	}

 static void
br_read(EdRead *R, int nc, real *L, real *U, int *Cvar, int nv)
{
	ASL *asl;
	int i, inc, j, k;
	int (*Xscanf)(EdRead*, const char*, ...);

	asl = R->asl;

	if (U)
		inc = 1;
	else {
		U = L + 1;
		inc = 2;
		}
	Xscanf = xscanf;
	Xscanf(R, ""); /* purge line */
	for(i = 0; i < nc; i++, L += inc, U += inc) {
		switch(edag_peek(R) - '0') {
		  case 0:
			if (Xscanf(R,"%lf %lf",L,U)!= 2)
				badline(R);
			break;
		  case 1:
			if (Xscanf(R, "%lf", U) != 1)
				badline(R);
			*L = negInfinity;
			break;
		  case 2:
			if (Xscanf(R, "%lf", L) != 1)
				badline(R);
			*U = Infinity;
			break;
		  case 3:
			*L = negInfinity;
			*U = Infinity;
			Xscanf(R, ""); /* purge line */
			break;
		  case 4:
			if (Xscanf(R, "%lf", L) != 1)
				badline(R);
			*U = *L;
			break;
		  case 5:
			if (Cvar) {
				if (Xscanf(R, "%d %d", &k, &j) == 2
				 && j > 0 && j <= nv) {
					Cvar[i] = j;
					*L = k & 2 ? negInfinity : 0.;
					*U = k & 1 ?    Infinity : 0.;
					break;
					}
				}
		  default:
			badline(R);
		  }
		}
	}

 static expr *
aholread(EdRead *R)
{
	int i, k;
	expr_h *rvh;
	char *s1;
	FILE *nl = R->nl;

	k = getc(nl);
	if (k < '1' || k > '9')
		badline(R);
	i = k - '0';
	while((k = getc(nl)) != ':') {
		if (k < '0' || k > '9')
			badline(R);
		i = 10*i + k - '0';
		}
	rvh = (expr_h *)mem(memadj(sizeof(expr_h) + i));
	for(s1 = rvh->sym;;) {
		if ((k = getc(nl)) < 0) {
			fprintf(Stderr,
				 "Premature end of file in aholread, line %ld of %s\n",
					R->Line, R->asl->i.filename_);
				exit_ASL(R,1);
			}
		if (k == '\n') {
			R->Line++;
			if (!i)
				break;
			}
		if (--i < 0)
			badline(R);
		*s1++ = k;
		}
	*s1 = 0;
	rvh->op = (efunc*)f_OPHOL;
	rvh->a = nv1;
	return (expr *)rvh;
	}

 static expr *
bholread(EdRead *R)
{
	int i;
	expr_h *rvh;
	char *s;

	if (xscanf(R, "%d", &i) != 1)
		badline(R);
	rvh = (expr_h *)mem(memadj(sizeof(expr_h) + i));
	s = rvh->sym;
	if (fread(s, i, 1, R->nl) != 1)
		badline(R);
	s[i] = 0;
	rvh->op = (efunc*)f_OPHOL;
	rvh->a = nv1;
	for(;;) switch(*s++) {
			case 0: goto break2; /* so we return at end of fcn */
			case '\n': R->Line++;
			}
 break2:
	return (expr *)rvh;
	}

 static void
hv2init_ignore(ASL *asl, int hid_limit, int nobj, real *ow, real *y) {}

 int
fgh_read_ASL(ASL *a, FILE *nl, int flags)
{
	EdRead ER, *R;
	Jmp_buf JB;
	cgrad *cg, **cgp;
	char fname[128];
	expr_v *e;
	func_info *fi;
	int i, j, k, *ka, maxfwd1, nc, nc0, nc1, nco, ncom, nlcon, nlin;
	int no, nv, nvc, nvo, nvr, nxv, readall;
	int (*Xscanf)(EdRead*, const char*, ...);
	ograd *og, **ogp;
	real t;
	size_t *kaz, nz;
	unsigned x;

	ASL_CHECK(a, ASL_read_fgh, "fgh_read");
	flagsave_ASL(a, flags); /* includes allocation of LUv, LUrhs, A_vals or Cgrad, etc. */
	asl = (ASL_fgh*)a;
#define asl ((ASL_fgh*)a)

	ed_reset(a);
	R = EdReadInit_ASL(&ER, a, nl, 0);
	if (flags & ASL_return_read_err) {
		a->i.err_jmp_ = &JB;
		i = setjmp(JB.jb);
		if (i) {
			a->i.err_jmp_ = 0;
			return i;
			}
		}
	nlcon = a->i.n_lcon_;
	if (nlcon && !(flags & ASL_allow_CLP)) {
		if (a->i.err_jmp_)
			return ASL_readerr_CLP;
		sorry_CLP(R, "logical constraints");
		}

	readall = flags & ASL_keep_all_suffixes;
	if (nfunc)
		func_add(a);
	if (binary_nl)
		holread = bholread;
	else
		holread = aholread;

	ncom = comb + comc + como + comc1 + como1;
	nc0 = n_con;
	nc = nc1 = nc0 + a->i.nsufext[ASL_Sufkind_con];
	no = n_obj;
	nvc = c_vars;
	nvo = o_vars;
	nco = nc + no + nlcon;
	if (no < 0 || nco <= 0)
		scream(R, ASL_readerr_corrupt,
			"ed2read: nc = %d, no = %d, nlcon = %d\n",
			nc0, no, nlcon);
	nxv = a->i.nsufext[ASL_Sufkind_var];
	nvr = n_var; /* nv for reading */
	nv0 = nv1 = nvr + nxv;
	max_var = nv = nv1 + ncom;
	combc = comb + comc;
	ncom0 = ncom_togo = combc + como;
	nzclim = ncom0 >> 3;
	ncom1 = comc1 + como1;
	nv0b = nv1 + comb;
	nv0c = nv0b + comc;
	nv01 = nv1 + ncom0;
	amax = lasta = lasta0 = lasta00 = nv1 + 1;
	ka = 0;
	kaz = 0;
	if ((maxfwd1 = maxfwd + 1) > 1)
		nvref = maxfwd1*((ncom0 < vrefGulp ? ncom0 : vrefGulp) + 1);
	x = nco*sizeof(cde) + no*sizeof(ograd *)
		+ nv*(sizeof(expr_v) + 2*sizeof(int))
		+ ncom0*sizeof(cexp)
		+ ncom1*sizeof(cexp1)
		+ nfunc*sizeof(func_info *)
		+ nvref*sizeof(int)
		+ no;
	nvar0 = a->i.n_var0;
	if (!(nvinc = a->i.nvinc = a->i.n_var_ - nvar0 + nxv))
		nvar0 += ncom0 + ncom1;
	if (pi0) {
		memset(pi0, 0, nc*sizeof(real));
		if (havepi0)
			memset(havepi0, 0, nc);
		}
	if (X0)
		memset(X0, 0, nvr*sizeof(real));
	if (havex0)
		memset(havex0, 0, nvr);
	e = var_e = (expr_v *)M1zapalloc(x);
	var_ex = e + nv0;
	var_ex1 = var_ex + ncom0;
	con_de = (cde *)(e + nv);
	lcon_de = con_de + nc;
	for(k = 0; k < nv; e++) {
		e->op = f_OPVARVAL;
		e->a = k++;
		}
	obj_de = lcon_de + nlcon;
	Ograd = (ograd **)(obj_de + no);
	cexps = (cexp *)(Ograd + no);
	cexpsc = cexps + comb;
	cexpso = cexpsc + comc;
	cexps1 = (cexp1 *)(cexpse = cexps + ncom0);
	funcs = (func_info **)(cexps1 + ncom1);
	zc = (int *)(funcs + nfunc);
	zci = zc + nv;
	vrefx = zci + nv;
	objtype = (char *)(vrefx + nvref);
	if (nvref) {
		vrefnext = vrefx + maxfwd1;
		nvref -= maxfwd1;
		}
	if (n_cc && !cvar)
		cvar = (int*)M1alloc(nc*sizeof(int));
	if (cvar)
		memset(cvar, 0, nc*sizeof(int));
	last_d = 0;
	nz = 0;
	Xscanf = xscanf;
	for(;;) {
		ER.can_end = 1;
		i = edag_peek(R);
		if (i == EOF) {
			free(imap);
			adjoints = (real *)M1zapalloc(amax*Sizeof(real));
			adjoints_nv1 = &adjoints[nv1];
			adjust(a, flags);
			nzjac = nz;
			if (!Lastx)
				Lastx = (real *)M1alloc(nv0*sizeof(real));
			fclose(nl);
			nderps += nderp;
			a->p.Objval  = a->p.Objval_nomap  = obj2val_ASL;
			a->p.Objgrd  = a->p.Objgrd_nomap  = obj2grd_ASL;
			a->p.Conval  = con2val_ASL;
			a->p.Jacval  = jac2val_ASL;
			a->p.Hvcomp  = a->p.Hvcomp_nomap  = hv2comp_ASL;
			a->p.Hvcompd = hv2compd_ASL;
			a->p.Hvcomps = hv2comps_ASL;
			a->p.Hvinit  = a->p.Hvinit_nomap  = hv2init_ignore;
			a->p.Conival = con2ival_ASL;
			a->p.Conival_nomap = con2ival_nomap_ASL;
			a->p.Congrd  = con2grd_ASL;
			a->p.Congrd_nomap = con2grd_nomap_ASL;
			a->p.Lconval = lcon2val_ASL;
			a->p.Xknown  = x2known_ASL;
			return prob_adj_ASL(a);
			}
		ER.can_end = 0;
		k = -1;
		switch(i) {
			case 'C':
				Xscanf(R, "%d", &k);
				if (k < 0 || k >= nc0)
					badline(R);
				co_read(R, con_de + k, 1);
				break;
			case 'F':
				if (Xscanf(R, "%d %d %d %127s",
						&i, &j, &k, fname) != 4
				|| i < 0 || i >= nfunc)
					badline(R);
				if ((fi = func_lookup(a, fname,0))) {
					if (fi->nargs != k && fi->nargs >= 0
					 && (k >= 0 || fi->nargs < -(k+1)))
						scream(R, ASL_readerr_argerr,
				"function %s: disagreement of nargs: %d and %d\n",
					 		fname,fi->nargs, k);
					}
				else {
					fi = (func_info *)mem(sizeof(func_info));
					fi->ftype = j;
					fi->nargs = k;
					fi->funcp = 0;
					fi->name = (Const char *)strcpy((char*)
						mem(memadj(strlen(fname)+1)),
						fname);
					}
				if (!fi->funcp && !(fi->funcp = dynlink(fname)))
					scream(R, ASL_readerr_unavail,
						"function %s not available\n",
						fname);
				funcs[i] = fi;
				break;
			case 'G':
				if (Xscanf(R, "%d %d", &j, &k) != 2
				|| j < 0 || j >= no || k <= 0 || k > nvo)
					badline(R);
				ogp = Ograd + j;
				while(k--) {
					*ogp = og = (ograd *)mem(sizeof(ograd));
					ogp = &og->next;
					if (Xscanf(R, "%d %lf", &i, &og->coef) != 2)
						badline(R);
					og->varno = i;
					}
				*ogp = 0;
				break;
			case 'J':
				if (Xscanf(R, "%d %d", &j, &k) != 2
				|| j < 0 || j >= nc0 || k <= 0 || k > nvc)
					badline(R);
				nz += k;
				if (A_vals) {
					j += Fortran;
					if (ka) {
						while(k--) {
							if (Xscanf(R, "%d %lf",
								&i, &t) != 2)
								badline(R);
							i = ka[i]++;
							A_vals[i] = t;
							A_rownos[i] = j;
							}
						}
					else {
						while(k--) {
							if (Xscanf(R, "%d %lf",
								&i, &t) != 2)
								badline(R);
							i = kaz[i]++;
							A_vals[i] = t;
							A_rownos[i] = j;
							}
						}
					break;
					}
				cgp = Cgrad + j;
				j = 0;
				while(k--) {
					*cgp = cg = (cgrad *)mem(sizeof(cgrad));
					cgp = &cg->next;
					if (k_seen) {
						if (Xscanf(R, "%d %lf", &i, &cg->coef) != 2)
							badline(R);
						}
					else
						if (Xscanf(R, "%d %d %lf", &i, &j,
					    		    &cg->coef) != 3)
							badline(R);
					cg->varno = i;
					cg->goff = j;
					}
				*cgp = 0;
				break;
			case 'L':
				Xscanf(R, "%d", &k);
				if (k < 0 || k >= nlcon)
					badline(R);
				co_read(R, lcon_de + k, 0);
				break;
			case 'O':
				if (Xscanf(R, "%d %d", &k, &j) != 2
				 || k < 0 || k >= no)
					badline(R);
				objtype[k] = j;
				co_read(R, obj_de + k, 1);
				break;
			case 'V':
				if (Xscanf(R, "%d %d %d", &k, &nlin, &j) != 3)
					badline(R);
				if (k >= nvar0)
					k += nvinc;
				if (k < nv0 || k >= nv)
					badline(R);
				if (j)
					cexp1_read(R, j, k, nlin);
				else
					cexp_read(R, k, nlin);
				break;
			case 'S':
				Suf_read_ASL(R, readall);
				break;
			case 'r':
				br_read(R, asl->i.n_con0, LUrhs, Urhsx, cvar, nv0);
				break;
			case 'b':
				br_read(R, asl->i.n_var0, LUv, Uvx, 0, 0);
				break;
			case 'K':
			case 'k':
				k_seen++;
				if (ka_read_ASL(a, R, i, &ka, &kaz))
					badline(R);
				break;
			case 'x':
				if (!Xscanf(R,"%d",&k)
				|| k < 0 || k > nvr)
					badline(R);
				if (!X0 && want_xpi0 & 1) {
					x = nv1*sizeof(real);
					if (want_xpi0 & 4)
						x += nv1;
					X0 = (real *)M1zapalloc(x);
					if (want_xpi0 & 4)
						havex0 = (char*)(X0 + nv1);
					}
				while(k--) {
					if (Xscanf(R, "%d %lf", &j, &t) != 2)
						badline(R);
					if (X0) {
						X0[j] = t;
						if (havex0)
							havex0[j] = 1;
						}
					}
				break;
			case 'd':
				if (!Xscanf(R,"%d",&k)
				|| k < 0 || k > nc0)
					badline(R);
				if (!pi0 && want_xpi0 & 2) {
					x = nc*sizeof(real);
					if (want_xpi0 & 4)
						x += nc;
					pi0 = (real *)M1zapalloc(x);
					if (want_xpi0 & 4)
						havepi0 = (char*)(pi0 + nc);
					}
				while(k--) {
					if (Xscanf(R, "%d %lf", &j, &t) != 2
					 || j < 0 || j >= nc0)
						badline(R);
					if (pi0) {
						pi0[j] = t;
						if (havepi0)
							havepi0[j] = 1;
						}
					}
				break;
			default:
				badline(R);
			}
		}
	}
#ifdef __cplusplus
	}
#endif