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