phobos/std/numeric.d
2007-11-27 20:28:40 +00:00

85 lines
2.3 KiB
D

// Written in the D programming language.
/**
This module is a port of a growing fragment of the $(D_PARAM numeric)
header in Alexander Stepanov's $(LINK2 http://sgi.com/tech/stl,
Standard Template Library), with a few additions.
Macros:
WIKI = Phobos/StdNumeric
Author:
Andrei Alexandrescu
*/
/*
* Copyright (C) 2004-2006 by Digital Mars, www.digitalmars.com
* Written by Andrei Alexandrescu, www.erdani.org
*
* This software is provided 'as-is', without any express or implied
* warranty. In no event will the authors be held liable for any damages
* arising from the use of this software.
*
* Permission is granted to anyone to use this software for any purpose,
* including commercial applications, and to alter it and redistribute it
* freely, subject to the following restrictions:
*
* o The origin of this software must not be misrepresented; you must not
* claim that you wrote the original software. If you use this software
* in a product, an acknowledgment in the product documentation would be
* appreciated but is not required.
* o Altered source versions must be plainly marked as such, and must not
* be misrepresented as being the original software.
* o This notice may not be removed or altered from any source
* distribution.
*/
module std.numeric;
import std.math;
import std.stdio;
/**
Implements the $(LINK2 http://tinyurl.com/2zb9yr,secant method) for
finding a root of the function $(D_PARAM f) starting from points
[xn_1, x_n] (ideally close to the root). $(D_PARAM Num) may be
$(D_PARAM float), $(D_PARAM double), or $(D_PARAM real).
Example:
----
float f(float x) {
return cos(x) - x*x*x;
}
final x = secantMethod(&f, 0f, 1f);
assert(approxEqual(x, 0.865474));
----
*/
template secantMethod(alias F)
{
Num secantMethod(Num)(Num xn_1, Num xn) {
auto fxn = F(xn_1), d = xn_1 - xn;
typeof(fxn) fxn_1;
xn = xn_1;
while (!approxEqual(d, 0) && isfinite(d)) {
xn_1 = xn;
xn -= d;
fxn_1 = fxn;
fxn = F(xn);
d *= -fxn / (fxn - fxn_1);
}
return xn;
}
}
unittest
{
scope(failure) writeln(stderr, "Failure testing secantMethod");
float f(float x) {
return cos(x) - x*x*x;
}
final x = secantMethod!(f)(0f, 1f);
assert(approxEqual(x, 0.865474));
}