```(* CAPO - Computational Analysis Platform for Oberon - by Alan Freed and Felix Friedrich. *) (* Version 1, Update 2 *) MODULE CalcD1; (** AUTHOR "adf"; PURPOSE "Computes a first-order derivative"; *) 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 df(x)/dx *) PROCEDURE Solve*( f: CalcFn.ReArg; atX: NbrRe.Real; differencing: NbrInt.Integer ): NbrRe.Real; VAR h, 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 / 3; hOpt := NbrRe.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; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atX + h ); result := (result - f( atX )) / h ELSIF differencing = Backward THEN result := f( atX ); result := (result - f( atX - h )) / h ELSE (* differencing = Central *) result := f( atX + h ); result := (result - f( atX - h )) / (2 * h) END; RETURN result END Solve; (** Computes df(z)/dz *) PROCEDURE SolveCplx*( f: CalcFn.CplxArg; atZ: NbrCplx.Complex; differencing: NbrInt.Integer ): NbrCplx.Complex; VAR h, hOpt, hMin, power: NbrRe.Real; ch, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 3; 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; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + ch ); result := (result - f( atZ )) / ch ELSIF differencing = Backward THEN result := f( atZ ); result := (result - f( atZ - ch )) / ch ELSE (* differencing = Central *) result := f( atZ + ch ); result := (result - f( atZ - ch )) / (2 * ch) END; RETURN result END SolveCplx; (** Computes 6f(z)/6x, 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, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 3; 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; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + ch ); result := (result - f( atZ )) / ch ELSIF differencing = Backward THEN result := f( atZ ); result := (result - f( atZ - ch )) / ch ELSE (* differencing = Central *) result := f( atZ + ch ); result := (result - f( atZ - ch )) / (2 * ch) END; RETURN result END SolveCplxRe; (** Computes 6f(z)/6y, 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, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 3; 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; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + ch ); result := (result - f( atZ )) / ch ELSIF differencing = Backward THEN result := f( atZ ); result := (result - f( atZ - ch )) / ch ELSE (* differencing = Central *) result := f( atZ + ch ); result := (result - f( atZ - ch )) / (2 * ch) END; RETURN result END SolveCplxIm; (** Computes 6f(z)/6r, 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, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 3; 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; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + ch ); result := (result - f( atZ )) / ch ELSIF differencing = Backward THEN result := f( atZ ); result := (result - f( atZ - ch )) / ch ELSE (* differencing = Central *) result := f( atZ + ch ); result := (result - f( atZ - ch )) / (2 * ch) END; RETURN result END SolveCplxAbs; (** Computes 6f(z)/6f, 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, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 3; 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; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + ch ); result := (result - f( atZ )) / ch ELSIF differencing = Backward THEN result := f( atZ ); result := (result - f( atZ - ch )) / ch ELSE (* differencing = Central *) result := f( atZ + ch ); result := (result - f( atZ - ch )) / (2 * ch) END; RETURN result END SolveCplxArg; BEGIN epsilon := 100 * NbrRe.Epsilon; zero := 0 END CalcD1.```