```(* CAPO - Computational Analysis Platform for Oberon - by Alan Freed and Felix Friedrich. *) (* Version 1, Update 2 *) MODULE CalcD4; (** 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 d4f(x)/dx4 *) PROCEDURE Solve*( f: CalcFn.ReArg; atX: NbrRe.Real; differencing: NbrInt.Integer ): NbrRe.Real; VAR h, h2, h4, 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 / 6; 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; h4 := h2 * h2; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atX + 4 * h ); result := result - 4 * f( atX + 3 * h ); result := result + 6 * f( atX + 2 * h ); result := result - 4 * f( atX + h ); result := (result + f( atX )) / h4 ELSIF differencing = Backward THEN result := f( atX ); result := result - 4 * f( atX - h ); result := result + 6 * f( atX - 2 * h ); result := result - 4 * f( atX - 3 * h ); result := (result + f( atX - 4 * h )) / h4 ELSE (* differencing = Central *) result := f( atX + 2 * h ); result := result - 4 * f( atX + h ); result := result + 6 * f( atX ); result := result - 4 * f( atX - h ); result := (result + f( atX - 2 * h )) / h4 END; RETURN result END Solve; (** Computes d4f(z)/dz4 *) PROCEDURE SolveCplx*( f: CalcFn.CplxArg; atZ: NbrCplx.Complex; differencing: NbrInt.Integer ): NbrCplx.Complex; VAR h, hOpt, hMin, power: NbrRe.Real; ch, ch2, ch4, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 6; 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; ch4 := ch2 * ch2; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + 4 * ch ); result := result - 4 * f( atZ + 3 * ch ); result := result + 6 * f( atZ + 2 * ch ); result := result - 4 * f( atZ + ch ); result := (result + f( atZ )) / ch4 ELSIF differencing = Backward THEN result := f( atZ ); result := result - 4 * f( atZ - ch ); result := result + 6 * f( atZ - 2 * ch ); result := result - 4 * f( atZ - 3 * ch ); result := (result + f( atZ - 4 * ch )) / ch4 ELSE (* differencing = Central *) result := f( atZ + 2 * ch ); result := result - 4 * f( atZ + ch ); result := result + 6 * f( atZ ); result := result - 4 * f( atZ - ch ); result := (result + f( atZ - 2 * ch )) / ch4 END; RETURN result END SolveCplx; (** Computes 64f(z)/6x4, 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, ch4, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 6; 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; ch4 := ch2 * ch2; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + 4 * ch ); result := result - 4 * f( atZ + 3 * ch ); result := result + 6 * f( atZ + 2 * ch ); result := result - 4 * f( atZ + ch ); result := (result + f( atZ )) / ch4 ELSIF differencing = Backward THEN result := f( atZ ); result := result - 4 * f( atZ - ch ); result := result + 6 * f( atZ - 2 * ch ); result := result - 4 * f( atZ - 3 * ch ); result := (result + f( atZ - 4 * ch )) / ch4 ELSE (* differencing = Central *) result := f( atZ + 2 * ch ); result := result - 4 * f( atZ + ch ); result := result + 6 * f( atZ ); result := result - 4 * f( atZ - ch ); result := (result + f( atZ - 2 * ch )) / ch4 END; RETURN result END SolveCplxRe; (** Computes 64f(z)/6y4, 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, ch4, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 6; 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; ch4 := ch2 * ch2; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + 4 * ch ); result := result - 4 * f( atZ + 3 * ch ); result := result + 6 * f( atZ + 2 * ch ); result := result - 4 * f( atZ + ch ); result := (result + f( atZ )) / ch4 ELSIF differencing = Backward THEN result := f( atZ ); result := result - 4 * f( atZ - ch ); result := result + 6 * f( atZ - 2 * ch ); result := result - 4 * f( atZ - 3 * ch ); result := (result + f( atZ - 4 * ch )) / ch4 ELSE (* differencing = Central *) result := f( atZ + 2 * ch ); result := result - 4 * f( atZ + ch ); result := result + 6 * f( atZ ); result := result - 4 * f( atZ - ch ); result := (result + f( atZ - 2 * ch )) / ch4 END; RETURN result END SolveCplxIm; (** Computes 64f(z)/6r4, 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, ch4, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 6; 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; ch4 := ch2 * ch2; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + 4 * ch ); result := result - 4 * f( atZ + 3 * ch ); result := result + 6 * f( atZ + 2 * ch ); result := result - 4 * f( atZ + ch ); result := (result + f( atZ )) / ch4 ELSIF differencing = Backward THEN result := f( atZ ); result := result - 4 * f( atZ - ch ); result := result + 6 * f( atZ - 2 * ch ); result := result - 4 * f( atZ - 3 * ch ); result := (result + f( atZ - 4 * ch )) / ch4 ELSE (* differencing = Central *) result := f( atZ + 2 * ch ); result := result - 4 * f( atZ + ch ); result := result + 6 * f( atZ ); result := result - 4 * f( atZ - ch ); result := (result + f( atZ - 2 * ch )) / ch4 END; RETURN result END SolveCplxAbs; (** Computes 64f(z)/6f4, 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, ch4, result, temp: NbrCplx.Complex; BEGIN (* Select an optimum step size. *) power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 6; 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; ch4 := ch2 * ch2; (* Compute an approximate value for the derivative. *) IF differencing = Forward THEN result := f( atZ + 4 * ch ); result := result - 4 * f( atZ + 3 * ch ); result := result + 6 * f( atZ + 2 * ch ); result := result - 4 * f( atZ + ch ); result := (result + f( atZ )) / ch4 ELSIF differencing = Backward THEN result := f( atZ ); result := result - 4 * f( atZ - ch ); result := result + 6 * f( atZ - 2 * ch ); result := result - 4 * f( atZ - 3 * ch ); result := (result + f( atZ - 4 * ch )) / ch4 ELSE (* differencing = Central *) result := f( atZ + 2 * ch ); result := result - 4 * f( atZ + ch ); result := result + 6 * f( atZ ); result := result - 4 * f( atZ - ch ); result := (result + f( atZ - 2 * ch )) / ch4 END; RETURN result END SolveCplxArg; BEGIN epsilon := 100 * NbrRe.Epsilon; zero := 0 END CalcD4.```