```(* CAPO - Computational Analysis Platform for Oberon - by Alan Freed and Felix Friedrich. *) (* Version 1, Update 2 *) MODULE CalcD3; (** AUTHOR "adf"; PURPOSE "Computes integer-order, i.e. classical, derivatives"; *) IMPORT NbrInt, NbrRe, NbrCplx, MathRe, CalcFn; CONST (** Admissible parameters to be passed for establishing the differencing scheme used to compute a derivative. *) Forward* = 9; Central* = 10; Backward* = 11; VAR epsilon, zero: NbrRe.Real; (* Force the argument in and out of addressable memory to minimize round-off error. *) PROCEDURE DoNothing( x: NbrRe.Real ); END DoNothing; PROCEDURE DoCplxNothing( z: NbrCplx.Complex ); END DoCplxNothing; (** Computes d3f(x)/dx3 *) PROCEDURE Solve*( f: CalcFn.ReArg; atX: NbrRe.Real; differencing: NbrInt.Integer ): NbrRe.Real; VAR h, h2, h3, hOpt, hMin, power, result, temp: NbrRe.Real; BEGIN (* Select an optimum step size. See v5.7 on Numerical Derivatives in Press et al., Numerical Recipes. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 5; hOpt := ABS( atX ) * MathRe.Power( epsilon, power ); h := NbrRe.Max( hOpt, hMin ); (* Refine h so that x + h and x differ by an exactly representable number in memory. *) temp := atX + h; DoNothing( temp ); h := temp - atX; h2 := h * h; h3 := h * h2; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atX + 3 * h ); result := result - 3 * f( atX + 2 * h ); result := result + 3 * f( atX + h ); result := (result - f( atX )) / h3 ELSIF differencing = Backward THEN result := f( atX ); result := result - 3 * f( atX - h ); result := result + 3 * f( atX - 2 * h ); result := (result - f( atX - 3 * h )) / h3 ELSE (* differencing = Central *) result := f( atX + 2 * h ); result := result - 2 * f( atX + h ); result := result + 2 * f( atX - h ); result := (result - f( atX - 2 * h )) / (2 * h3) END; RETURN result END Solve; (** Computes d3f(z)/dz3 *) PROCEDURE SolveCplx*( f: CalcFn.CplxArg; atZ: NbrCplx.Complex; differencing: NbrInt.Integer ): NbrCplx.Complex; VAR h, hOpt, hMin, power: NbrRe.Real; ch, ch2, ch3, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 5; hOpt := NbrCplx.Abs( atZ ) * MathRe.Power( epsilon, power ); h := NbrRe.Max( hOpt, hMin ); NbrCplx.Set( h, h, ch ); (* Refine h so that z + ch and z differ by an exactly representable number in memory. *) temp := atZ + ch; DoCplxNothing( temp ); ch := temp - atZ; ch2 := ch * ch; ch3 := ch2 * ch; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + 3 * ch ); result := result - 3 * f( atZ + 2 * ch ); result := result + 3 * f( atZ + ch ); result := (result - f( atZ )) / ch3 ELSIF differencing = Backward THEN result := f( atZ ); result := result - 3 * f( atZ - ch ); result := result + 3 * f( atZ - 2 * ch ); result := (result - f( atZ - 3 * ch )) / ch3 ELSE (* differencing = Central *) result := f( atZ + 2 * ch ); result := result - 2 * f( atZ + ch ); result := result + 2 * f( atZ - ch ); result := (result - f( atZ - 2 * ch )) / (2 * ch3) END; RETURN result END SolveCplx; (** Computes 63f(z)/6x3, z = x + i y *) PROCEDURE SolveCplxRe*( f: CalcFn.CplxArg; atZ: NbrCplx.Complex; differencing: NbrInt.Integer ): NbrCplx.Complex; VAR h, hOpt, hMin, power: NbrRe.Real; ch, ch2, ch3, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 5; hOpt := NbrCplx.Abs( atZ ) * MathRe.Power( epsilon, power ); h := NbrRe.Max( hOpt, hMin ); NbrCplx.Set( h, zero, ch ); (* Refine h so that z + ch and z differ by an exactly representable number in memory. *) temp := atZ + ch; DoCplxNothing( temp ); ch := temp - atZ; ch2 := ch * ch; ch3 := ch2 * ch; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + 3 * ch ); result := result - 3 * f( atZ + 2 * ch ); result := result + 3 * f( atZ + ch ); result := (result - f( atZ )) / ch3 ELSIF differencing = Backward THEN result := f( atZ ); result := result - 3 * f( atZ - ch ); result := result + 3 * f( atZ - 2 * ch ); result := (result - f( atZ - 3 * ch )) / ch3 ELSE (* differencing = Central *) result := f( atZ + 2 * ch ); result := result - 2 * f( atZ + ch ); result := result + 2 * f( atZ - ch ); result := (result - f( atZ - 2 * ch )) / (2 * ch3) END; RETURN result END SolveCplxRe; (** Computes 63f(z)/6y3, z = x + i y *) PROCEDURE SolveCplxIm*( f: CalcFn.CplxArg; atZ: NbrCplx.Complex; differencing: NbrInt.Integer ): NbrCplx.Complex; VAR h, hOpt, hMin, power: NbrRe.Real; ch, ch2, ch3, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 5; hOpt := NbrCplx.Abs( atZ ) * MathRe.Power( epsilon, power ); h := NbrRe.Max( hOpt, hMin ); NbrCplx.Set( zero, h, ch ); (* Refine h so that z + ch and z differ by an exactly representable number in memory. *) temp := atZ + ch; DoCplxNothing( temp ); ch := temp - atZ; ch2 := ch * ch; ch3 := ch2 * ch; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + 3 * ch ); result := result - 3 * f( atZ + 2 * ch ); result := result + 3 * f( atZ + ch ); result := (result - f( atZ )) / ch3 ELSIF differencing = Backward THEN result := f( atZ ); result := result - 3 * f( atZ - ch ); result := result + 3 * f( atZ - 2 * ch ); result := (result - f( atZ - 3 * ch )) / ch3 ELSE (* differencing = Central *) result := f( atZ + 2 * ch ); result := result - 2 * f( atZ + ch ); result := result + 2 * f( atZ - ch ); result := (result - f( atZ - 2 * ch )) / (2 * ch3) END; RETURN result END SolveCplxIm; (** Computes 63f(z)/6r3, z = r exp( i f ) *) PROCEDURE SolveCplxAbs*( f: CalcFn.CplxArg; atZ: NbrCplx.Complex; differencing: NbrInt.Integer ): NbrCplx.Complex; VAR h, hOpt, hMin, power: NbrRe.Real; ch, ch2, ch3, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 5; hOpt := NbrCplx.Abs( atZ ) * MathRe.Power( epsilon, power ); h := NbrRe.Max( hOpt, hMin ); NbrCplx.SetPolar( h, zero, ch ); (* Refine h so that z + ch and z differ by an exactly representable number in memory. *) temp := atZ + ch; DoCplxNothing( temp ); ch := temp - atZ; ch2 := ch * ch; ch3 := ch2 * ch; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + 3 * ch ); result := result - 3 * f( atZ + 2 * ch ); result := result + 3 * f( atZ + ch ); result := (result - f( atZ )) / ch3 ELSIF differencing = Backward THEN result := f( atZ ); result := result - 3 * f( atZ - ch ); result := result + 3 * f( atZ - 2 * ch ); result := (result - f( atZ - 3 * ch )) / ch3 ELSE (* differencing = Central *) result := f( atZ + 2 * ch ); result := result - 2 * f( atZ + ch ); result := result + 2 * f( atZ - ch ); result := (result - f( atZ - 2 * ch )) / (2 * ch3) END; RETURN result END SolveCplxAbs; (** Computes 63f(z)/6f3, z = r exp( i f ) *) PROCEDURE SolveCplxArg*( f: CalcFn.CplxArg; atZ: NbrCplx.Complex; differencing: NbrInt.Integer ): NbrCplx.Complex; VAR h, hOpt, hMin, power: NbrRe.Real; ch, ch2, ch3, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 5; hOpt := NbrCplx.Arg( atZ ) * MathRe.Power( epsilon, power ); h := NbrRe.Max( hOpt, hMin ); NbrCplx.SetPolar( zero, h, ch ); (* Refine h so that z + ch and z differ by an exactly representable number in memory. *) temp := atZ + ch; DoCplxNothing( temp ); ch := temp - atZ; ch2 := ch * ch; ch3 := ch2 * ch; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + 3 * ch ); result := result - 3 * f( atZ + 2 * ch ); result := result + 3 * f( atZ + ch ); result := (result - f( atZ )) / ch3 ELSIF differencing = Backward THEN result := f( atZ ); result := result - 3 * f( atZ - ch ); result := result + 3 * f( atZ - 2 * ch ); result := (result - f( atZ - 3 * ch )) / ch3 ELSE (* differencing = Central *) result := f( atZ + 2 * ch ); result := result - 2 * f( atZ + ch ); result := result + 2 * f( atZ - ch ); result := (result - f( atZ - 2 * ch )) / (2 * ch3) END; RETURN result END SolveCplxArg; BEGIN epsilon := 100 * NbrRe.Epsilon; zero := 0 END CalcD3.```