--* From postmaster%watson.vnet.ibm.com@yktvmv.watson.ibm.com  Thu Jul 14 14:50:37 1994
--* Received: from yktvmv-ob.watson.ibm.com by asharp.watson.ibm.com (AIX 3.2/UCB 5.64/930311)
--*           id AA29983; Thu, 14 Jul 1994 14:50:37 -0400
--* Received: from watson.vnet.ibm.com by yktvmv.watson.ibm.com (IBM VM SMTP V2R3)
--*    with BSMTP id 4981; Thu, 14 Jul 94 14:50:40 EDT
--* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id <A.RMC.NOTE.VAGENT2.9167.Jul.14.14:50:39.-0400>
--*           for asbugs@watson; Thu, 14 Jul 94 14:50:40 -0400
--* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id 9161; Thu, 14 Jul 1994 14:50:39 EDT
--* Received: from spadserv.watson.ibm.com by yktvmv.watson.ibm.com
--*    (IBM VM SMTP V2R3) with TCP; Thu, 14 Jul 94 14:50:39 EDT
--* Received: by spadserv.watson.ibm.com (AIX 3.2/UCB 5.64/900524)
--*           id AA17171; Thu, 14 Jul 1994 14:50:47 -0400
--* Date: Thu, 14 Jul 1994 14:50:47 -0400
--* From: rmc@spadserv.watson.ibm.com
--* X-External-Networks: yes
--* Message-Id: <9407141850.AA17171@spadserv.watson.ibm.com>
--* To: asbugs@watson.ibm.com
--* Subject: [3] Intermittent core dump (2) [compact.as][36.0]

--@ Fixed  by:  SSD   Thu Jul 28 16:52:29 EDT 1994 
--@ Tested by:  none 
--@ Summary:    Intermittent bug was not reproducible. File seems to compile and run just fine. 


-- Compact finite differences

#include "aslib"

AF ==> Array F;		SI ==> SingleInteger;

derivatives(F:Field with {sqrt:%->%}, funvals:AF, h:F):(AF) == {
  	import from SI;
  	n:SI   := # funvals;

  	(n < 5) => error("Too few function values for boundary formulas");

	-- initialize local arrays
  	b  :AF := new(n,0);  	eta :AF := new(n,0);
  	x  :AF := new(n,0);  	y   :AF := new(n,0);

	--*-*-*-* define integer<-->field operations *-*-*-*--

  	(n:SI) * (x:F):F == n::F * x;	(x:F) * (n:SI):F == n::F * x;
  	(n:SI) + (x:F):F == n::F + x;  	(x:F) + (n:SI):F == n::F + x;
  	(n:SI) - (x:F):F == n::F - x;  	(x:F) - (n:SI):F == x - n::F;
  	(n:SI) / (x:F):F == n::F / x;  	(x:F) / (n:SI):F == x/(n::F);

	--*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**-*-*-*-*-*-*-*--

	-- Constants needed for exact factorization of matrix A
	c:F := 2::F + sqrt(3::F);	cinv:= 2::F - sqrt(3::F);
	
	-- Prepare the right-hand side differences

	b(1) := (-(25*c+3)/12*funvals(1) + (24*c-5)/6*funvals(2)
		-3*(2*c-1@F)/2*funvals(3) + (8*c-3)/6*funvals(4)
		-(3*c-1@F)/12*funvals(5) )/h;

  	for i:SI in 2..(n-1) repeat
    		b(i) := 3*(funvals(i+1) - funvals(i-1))/h;

	b(n) := (103*funvals(n) - 182*funvals(n-1) + 126*funvals(n-2)
        	- 58*funvals(n-3) + 11*funvals(n-4))
          		/(12*h);

	-- Solve Ly = b

	y(1)   := b(1);
  	for i:SI in 2..n repeat
    		y(i) := b(i) - cinv*y(i-1);

	-- Solve Ux = y

  	x(n) := y(n)*cinv;
  	for i:SI in (n-1)..1 by -1 repeat
    		x(i) := cinv*(y(i) - x(i+1));

  	x
}

extend DoubleFloat: with {
	sqrt: % -> %
} == add {
	import from DoubleFloatElementaryFunctions;
	sqrt(x:%):% == (sqrt$DoubleFloatElementaryFunctions)(x);
}