D Programming Language 2.0

- Comments
- D
- Search
- Home

Last update Wed Nov 4 22:39:15 2009

# std.numeric

This module is a port of a growing fragment of the numeric header in Alexander Stepanov's Standard Template Library, with a few additions.

License:
Boost License 1.0.

Authors:
Andrei Alexandrescu, Don Clugston

Copyright Andrei Alexandrescu 2008 - 2009. Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at ) http:
//www.boost.org/LICENSE_1_0.txt

enum CustomFloatFlags;
  Flags for custom floating-point formats. In the case of IEEE754 types, all of these flags are applied.

  storeNormalized
    Store values in normalized form by default, i.e. the fraction is assumed to be 1.nnnn, not 0.nnnn. Normalization is the default in all IEE754 types. If a floating-point type is defined to store only numbers inside (-1, 1), normalization may not be useful.

  allowDenorm
    Allow denormalized values (a la IEEE754 denormalized). If this flag is on, a value with an exponent containing all bits zero is considered to contain a denormalized value.

  infinity
    Support for infinity values (a la IEEE754 infinity). If this flag is on, a value with an exponent containing all bits one and fraction zero is assumed to contain a (valid) allowDenormalized value. The sign, if applicable, is decided by the sign bit.

  nan
    Support for NaN (Not a Number) values (a la IEEE754 Not a Number). If this flag is on, a value with an exponent containing all bits one and a nonzero fraction is considered to contain a allowDenormalized value.

  all
    Include all of the above options.

struct CustomFloat(bool signBit,uint fractionBits,uint exponentBits,uint bias = (1u << exponentBits - 1) - 1,CustomFloatFlags flags = CustomFloatFlags.all);
  Allows user code to define custom floating-point formats. These formats are for storage only; to do operations on custom floating-point types, one must extract the value into a float or double. After the operation is completed the result can be deposited in a custom floating-point value by using its assignment operator.

  Example:

      // Define a 16-bit floating point type
      alias CustomFloat!(1, 5, 10) HalfFloat;
      auto x = HalfFloat(0.125);
      assert(halfFloat.get!float == 0.125);

  this(float input);
    Initialize from a float.

  this(double input);
    Initialize from a double.

  template opAssign(F) if (staticIndexOf!(Unqual!(F),float,double) >= 0)
    Assigns from either a float or a double.

    void opAssign(F input);
      Assigns from either a float or a double.

  template get(F)
    Fetches the stored value either as a float or double.

    F get();
      Fetches the stored value either as a float or double.

template FPTemporary(F) if (isFloatingPoint!(F))
  Defines the fastest type to use when storing temporaries of a calculation intended to ultimately yield a result of type F (where F must be one of float, double, or real). When doing a multi-step computation, you may want to store intermediate results as FPTemporary!F.

  Example:

      // Average numbers in an array
      double avg(in double[] a)
      {
          if (a.length == 0) return 0;
          FPTemporary!double result = 0;
          foreach (e; a) result += e;
          return result / a.length;
      }

  The necessity of FPTemporary stems from the optimized floating-point operations and registers present in virtually all processors. When adding numbers in the example above, the addition may in fact be done in real precision internally. In that case, storing the intermediate result in double format is not only less precise, it is also (surprisingly) slower, because a conversion from real to double is performed every pass through the loop. This being a lose-lose situation, FPTemporary!F has been defined as the fastest type to use for calculations at precision F. There is no need to define a type for the most accurate calculations, as that is always real.

  Finally, there is no guarantee that using FPTemporary!F will always be fastest, as the speed of floating-point calculations depends on very many factors.

template secantMethod(alias fun)
  Implements the secant method for finding a root of the function fun starting from points [xn_1, x_n] (ideally close to the root). Num may be float, double, or real.

  Example:

      float f(float x) {
          return cos(x) - x*x*x;
      }
      auto x = secantMethod!(f)(0f, 1f);
      assert(approxEqual(x, 0.865474));

T findRoot(T, R)(R delegate(T) f, T a, T b);
  Find a real root of a real function f(x) via bracketing.

  Given a function f and a range [a..b] such that f(a) and f(b) have opposite signs, returns the value of x in the range which is closest to a root of f(x). If f(x) has more than one root in the range, one will be chosen arbitrarily. If f(x) returns NaN, NaN will be returned; otherwise, this algorithm is guaranteed to succeed.

  Uses an algorithm based on TOMS748, which uses inverse cubic interpolation whenever possible, otherwise reverting to parabolic or secant interpolation. Compared to TOMS748, this implementation improves worst-case performance by a factor of more than 100, and typical performance by a factor of 2. For 80-bit reals, most problems require 8 to 15 calls to f(x) to achieve full machine precision. The worst-case performance (pathological cases) is approximately twice the number of bits.

  References:
  "On Enclosing Simple Roots of Nonlinear Equations", G. Alefeld, F.A. Potra, Yixun Shi, Mathematics of Computation 61, pp733-744 (1993). Fortran code available from www.netlib.org as algorithm TOMS478.

Tuple!(T,T,R,R) findRoot(T, R)(R delegate(T) f, T ax, T bx, R fax, R fbx, bool delegate(T lo, T hi) tolerance);
  Find root of a real function f(x) by bracketing, allowing the termination condition to be specified.

  Parameters:

  f | Function to be analyzed
  ax | Left bound of initial range of f known to contain the root.
  bx | Right bound of initial range of f known to contain the root.
  fax | Value of f(ax).
  fbx | Value of f(bx). (f(ax) and f(bx) are commonly known in advance.)
  tolerance | Defines an early termination condition. Receives the current upper and lower bounds on the root. The delegate must return true when these bounds are acceptable. If this function always returns false, full machine precision will be achieved.

  Returns:
  A tuple consisting of two ranges. The first two elements are the range (in x) of the root, while the second pair of elements are the corresponding function values at those points. If an exact root was found, both of the first two elements will contain the root, and the second pair of elements will be 0.

CommonType!(ElementType!(Range1),ElementType!(Range2)) euclideanDistance(Range1, Range2)(Range1 a, Range2 b);
CommonType!(ElementType!(Range1),ElementType!(Range2)) euclideanDistance(Range1, Range2, F)(Range1 a, Range2 b, F limit);
  Computes Euclidean distance between input ranges a and b. The two ranges must have the same length. The three-parameter version stops computation as soon as the distance is greater than or equal to limit (this is useful to save computation if a small distance is sought).

CommonType!(ElementType!(Range1),ElementType!(Range2)) dotProduct(Range1, Range2)(Range1 a, Range2 b);
Unqual!(CommonType!(F1,F2)) dotProduct(F1, F2)(in F1[] avector, in F2[] bvector);
  Computes the dot product of input ranges a and b. The two ranges must have the same length. If both ranges define length, the check is done once; otherwise, it is done at each iteration.

CommonType!(ElementType!(Range1),ElementType!(Range2)) cosineSimilarity(Range1, Range2)(Range1 a, Range2 b);
  Computes the cosine similarity of input ranges a and b. The two ranges must have the same length. If both ranges define length, the check is done once; otherwise, it is done at each iteration. If either range has all-zero elements, return 0.

bool normalize(R)(R range, ElementType!(R) sum = 1);
  Normalizes values in range by multiplying each element with a number chosen such that values sum up to sum. If elements in range sum to zero, assigns sum / range.length to all. Normalization makes sense only if all elements in range are positive. normalize assumes that is the case without checking it.

  Returns:
  true if normalization completed normally, false if all elements in range were zero or if range is empty.

ElementType!(Range) entropy(Range)(Range r);
ElementType!(Range) entropy(Range, F)(Range r, F max);
  Computes entropy of input range r in bits. This function assumes (without checking) that the values in r are all in [0, 1]. For the entropy to be meaningful, often r should be normalized too (i.e., its values should sum to 1). The two-parameter version stops evaluating as soon as the intermediate result is greater than or equal to max.

CommonType!(ElementType!(Range1),ElementType!(Range2)) kullbackLeiblerDivergence(Range1, Range2)(Range1 a, Range2 b);
  Computes the Kullback-Leibler divergence between input ranges a and b, which is the sum ai * log(ai / bi). The base of logarithm is 2. The ranges are assumed to contain elements in [0, 1]. Usually the ranges are normalized probability distributions, but this is not required or checked by kullbackLeiblerDivergence. If any element bi is zero and the corresponding element ai nonzero, returns infinity. (Otherwise, if ai == 0 && bi == 0, the term ai * log(ai / bi) is considered zero.) If the inputs are normalized, the result is positive.

CommonType!(ElementType!(Range1),ElementType!(Range2)) jensenShannonDivergence(Range1, Range2)(Range1 a, Range2 b);
CommonType!(ElementType!(Range1),ElementType!(Range2)) jensenShannonDivergence(Range1, Range2, F)(Range1 a, Range2 b, F limit);
  Computes the Jensen-Shannon divergence between a and b, which is the sum (ai * log(2 * ai / (ai + bi)) + bi * log(2 * bi / (ai + bi))) / 2. The base of logarithm is 2. The ranges are assumed to contain elements in [0, 1]. Usually the ranges are normalized probability distributions, but this is not required or checked by jensenShannonDivergence. If the inputs are normalized, the result is bounded within [0, 1]. The three-parameter version stops evaluations as soon as the intermediate result is greater than or equal to limit.

F gapWeightedSimilarity(alias comp = "a == b", R1, R2, F)(R1 s, R2 t, F lambda);
  The so-called "all-lengths gap-weighted string kernel" computes a similarity measure between s and t based on all of their common subsequences of all lengths. Gapped subsequences are also included.

  To understand what gapWeightedSimilarity(s, t, lambda) computes, consider first the case lambda = 1 and the strings s = ["Hello", "brave", "new", "world"] and t = ["Hello", "new", "world"]. In that case, gapWeightedSimilarity counts the following matches:

  1. three matches of length 1, namely "Hello", "new", and "world";
  2. three matches of length 2, namely ("Hello", "new"), ("Hello", "world"), and ("new", "world");
  3. one match of length 3, namely ("Hello", "new", "world").

  The call gapWeightedSimilarity(s, t, 1) simply counts all of these matches and adds them up, returning 7.

      string[] s = ["Hello", "brave", "new", "world"];
      string[] t = ["Hello", "new", "world"];
      assert(gapWeightedSimilarity(s, t, 1) == 7);

  Note how the gaps in matching are simply ignored, for example ("Hello", "new") is deemed as good a match as ("new", "world"). This may be too permissive for some applications. To eliminate gapped matches entirely, use lambda = 0:

      string[] s = ["Hello", "brave", "new", "world"];
      string[] t = ["Hello", "new", "world"];
      assert(gapWeightedSimilarity(s, t, 0) == 4);

  The call above eliminated the gapped matches ("Hello", "new"), ("Hello", "world"), and ("Hello", "new", "world") from the tally. That leaves only 4 matches.

  The most interesting case is when gapped matches still participate in the result, but not as strongly as ungapped matches. The result will be a smooth, fine-grained similarity measure between the input strings. This is where values of lambda between 0 and 1 enter into play: gapped matches are exponentially penalized with the number of gaps with base lambda. This means that an ungapped match adds 1 to the return value; a match with one gap in either string adds lambda to the return value; ...; a match with a total of n gaps in both strings adds pow(lambda, n) to the return value. In the example above, we have 4 matches without gaps, 2 matches with one gap, and 1 match with three gaps. The latter match is ("Hello", "world"), which has two gaps in the first string and one gap in the second string, totaling to three gaps. Summing these up we get 4 + 2 * lambda + pow(lambda, 3).

      string[] s = ["Hello", "brave", "new", "world"];
      string[] t = ["Hello", "new", "world"];
      assert(gapWeightedSimilarity(s, t, 0.5) == 4 + 0.5 * 2 + 0.125);

  gapWeightedSimilarity is useful wherever a smooth similarity measure between sequences allowing for approximate matches is needed. The examples above are given with words, but any sequences with elements comparable for equality are allowed, e.g. characters or numbers. gapWeightedSimilarity uses a highly optimized dynamic programming implementation that needs 16 * min(s.length, t.length) extra bytes of memory and Ο(s.length * t.length) time to complete.

Select!(isFloatingPoint!(F),F,double) gapWeightedSimilarityNormalized(alias comp = "a == b", R1, R2, F)(R1 s, R2 t, F lambda, F sSelfSim = F.init, F tSelfSim = F.init);
  The similarity per gapWeightedSimilarity has an issue in that it grows with the lengths of the two strings, even though the strings are not actually very similar. For example, the range ["Hello", "world"] is increasingly similar with the range ["Hello", "world", "world", "world",...] as more instances of "world" are appended. To prevent that, gapWeightedSimilarityNormalized computes a normalized version of the similarity that is computed as gapWeightedSimilarity(s, t, lambda) / sqrt(gapWeightedSimilarity(s, t, lambda) * gapWeightedSimilarity(s, t, lambda)). The function gapWeightedSimilarityNormalized (a so-called normalized kernel) is bounded in [0, 1], reaches 0 only for ranges that don't match in any position, and 1 only for identical ranges.

  Example:

      string[] s = ["Hello", "brave", "new", "world"];
      string[] t = ["Hello", "new", "world"];
      assert(gapWeightedSimilarity(s, s, 1) == 15);
      assert(gapWeightedSimilarity(t, t, 1) == 7);
      assert(gapWeightedSimilarity(s, t, 1) == 7);
      assert(gapWeightedSimilarityNormalized(s, t, 1) == 7. / sqrt(15. * 7));

  The optional parameters sSelfSim and tSelfSim are meant for avoiding duplicate computation. Many applications may have already computed gapWeightedSimilarity(s, s, lambda) and/or gapWeightedSimilarity(t, t, lambda). In that case, they can be passed as sSelfSim and tSelfSim, respectively.

struct GapWeightedSimilarityIncremental(Range,F = double) if (isRandomAccessRange!(Range) && hasLength!(Range));
GapWeightedSimilarityIncremental!(R,F) gapWeightedSimilarityIncremental(R, F)(R r1, R r2, F penalty);
  Similar to gapWeightedSimilarity, just works in an incremental manner by first revealing the matches of length 1, then gapped matches of length 2, and so on. The memory requirement is Ο(s.length * t.length). The time complexity is Ο(s.length * t.length) time for computing each step. Continuing on the previous example:

      string[] s = ["Hello", "brave", "new", "world"];
      string[] t = ["Hello", "new", "world"];
      auto simIter = gapWeightedSimilarityIncremental(s, t, 1);
      assert(simIter.front == 3); // three 1-length matches
      simIter.popFront;
      assert(simIter.front == 3); // three 2-length matches
      simIter.popFront;
      assert(simIter.front == 1); // one 3-length match
      simIter.popFront;
      assert(simIter.empty);     // no more match

  The implementation is based on the pseudocode in Fig. 4 of the paper "Efﬁcient Computation of Gapped Substring Kernels on Large Alphabets" by Rousu et al., with additional algorithmic and systems-level optimizations.

  this(Range s, Range t, F lambda);
    Constructs an object given two ranges s and t and a penalty lambda. Constructor completes in Ο(s.length * t.length) time and computes all matches of length 1.

  GapWeightedSimilarityIncremental opSlice();
    Returns this.

  void popFront();
    Computes the match of the popFront length. Completes in Ο(s.length * t.length) time.

  F front();
    Returns the gapped similarity at the current match length (initially 1, grows with each call to popFront).

  bool empty();
    Returns whether there are more matches.

T gcd(T)(T a, T b);
  Computes the greatest common divisor of a and b by using Euler's algorithm.

Copyright Andrei Alexandrescu 2008 - 2009. | Page generated by Ddoc.