Skip to content

[new concept under development] Python's fsum and others #2991

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
9il wants to merge 1 commit into from
Closed

[new concept under development] Python's fsum and others #2991

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion posix.mak
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -159,7 +159,7 @@ STD_PACKAGES = std $(addprefix std/,\

PACKAGE_std = array ascii base64 bigint bitmanip compiler complex concurrency \
concurrencybase conv cstream csv datetime demangle encoding exception file format \
functional getopt json math mathspecial meta metastrings mmfile numeric \
functional getopt json math mathspecial meta metastrings mmfile \
outbuffer parallelism path process random signals socket socketstream stdint \
stdio stdiobase stream string syserror system traits typecons typetuple uni \
uri utf uuid variant xml zip zlib
Expand All @@ -177,6 +177,7 @@ PACKAGE_std_experimental_allocator_building_blocks = \
kernighan_ritchie null_allocator package quantizer \
region scoped_allocator segregator stats_collector
PACKAGE_std_net = curl isemail
PACKAGE_std_numeric = package summation
PACKAGE_std_range = interfaces package primitives
PACKAGE_std_regex = package $(addprefix internal/,generator ir parser \
backtracking kickstart tests thompson)
Expand Down
203 changes: 2 additions & 201 deletions std/algorithm/iteration.d
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,8 @@ import std.functional; // : unaryFun, binaryFun;
import std.range.primitives;
import std.traits;

public import std.numeric.summation: Summation, sum;

template aggregate(fun...) if (fun.length >= 1)
{
/* --Intentionally not ddoc--
Expand Down Expand Up @@ -3892,207 +3894,6 @@ if (isSomeChar!C)
}
}

// sum
/**
Sums elements of $(D r), which must be a finite
$(XREF_PACK_NAMED range,primitives,isInputRange,input range). Although
conceptually $(D sum(r)) is equivalent to $(LREF reduce)!((a, b) => a +
b)(0, r), $(D sum) uses specialized algorithms to maximize accuracy,
as follows.

$(UL
$(LI If $(D $(XREF_PACK range,primitives,ElementType)!R) is a floating-point
type and $(D R) is a
$(XREF_PACK_NAMED range,primitives,isRandomAccessRange,random-access range) with
length and slicing, then $(D sum) uses the
$(WEB en.wikipedia.org/wiki/Pairwise_summation, pairwise summation)
algorithm.)
$(LI If $(D ElementType!R) is a floating-point type and $(D R) is a
finite input range (but not a random-access range with slicing), then
$(D sum) uses the $(WEB en.wikipedia.org/wiki/Kahan_summation,
Kahan summation) algorithm.)
$(LI In all other cases, a simple element by element addition is done.)
)

For floating point inputs, calculations are made in
$(DDLINK spec/type, Types, $(D real))
precision for $(D real) inputs and in $(D double) precision otherwise
(Note this is a special case that deviates from $(D reduce)'s behavior,
which would have kept $(D float) precision for a $(D float) range).
For all other types, the calculations are done in the same type obtained
from from adding two elements of the range, which may be a different
type from the elements themselves (for example, in case of
$(DDSUBLINK spec/type,integer-promotions, integral promotion)).

A seed may be passed to $(D sum). Not only will this seed be used as an initial
value, but its type will override all the above, and determine the algorithm
and precision used for summation.

Note that these specialized summing algorithms execute more primitive operations
than vanilla summation. Therefore, if in certain cases maximum speed is required
at expense of precision, one can use $(D reduce!((a, b) => a + b)(0, r)), which
is not specialized for summation.

Params:
seed = the initial value of the summation
r = a finite input range

Returns:
The sum of all the elements in the range r.
*/
auto sum(R)(R r)
if (isInputRange!R && !isInfinite!R && is(typeof(r.front + r.front)))
{
alias E = Unqual!(ElementType!R);
static if (isFloatingPoint!E)
alias Seed = typeof(E.init + 0.0); //biggest of double/real
else
alias Seed = typeof(r.front + r.front);
return sum(r, Unqual!Seed(0));
}
/// ditto
auto sum(R, E)(R r, E seed)
if (isInputRange!R && !isInfinite!R && is(typeof(seed = seed + r.front)))
{
static if (isFloatingPoint!E)
{
static if (hasLength!R && hasSlicing!R)
return seed + sumPairwise!E(r);
else
return sumKahan!E(seed, r);
}
else
{
return reduce!"a + b"(seed, r);
}
}

// Pairwise summation http://en.wikipedia.org/wiki/Pairwise_summation
private auto sumPairwise(Result, R)(R r)
{
static assert (isFloatingPoint!Result);
switch (r.length)
{
case 0: return cast(Result) 0;
case 1: return cast(Result) r.front;
case 2: return cast(Result) r[0] + cast(Result) r[1];
default: return sumPairwise!Result(r[0 .. $ / 2]) + sumPairwise!Result(r[$ / 2 .. $]);
}
}

// Kahan algo http://en.wikipedia.org/wiki/Kahan_summation_algorithm
private auto sumKahan(Result, R)(Result result, R r)
{
static assert (isFloatingPoint!Result && isMutable!Result);
Result c = 0;
for (; !r.empty; r.popFront())
{
auto y = r.front - c;
auto t = result + y;
c = (t - result) - y;
result = t;
}
return result;
}

/// Ditto
@safe pure nothrow unittest
{
import std.range;

//simple integral sumation
assert(sum([ 1, 2, 3, 4]) == 10);

//with integral promotion
assert(sum([false, true, true, false, true]) == 3);
assert(sum(ubyte.max.repeat(100)) == 25500);

//The result may overflow
assert(uint.max.repeat(3).sum() == 4294967293U );
//But a seed can be used to change the sumation primitive
assert(uint.max.repeat(3).sum(ulong.init) == 12884901885UL);

//Floating point sumation
assert(sum([1.0, 2.0, 3.0, 4.0]) == 10);

//Floating point operations have double precision minimum
static assert(is(typeof(sum([1F, 2F, 3F, 4F])) == double));
assert(sum([1F, 2, 3, 4]) == 10);

//Force pair-wise floating point sumation on large integers
import std.math : approxEqual;
assert(iota(ulong.max / 2, ulong.max / 2 + 4096).sum(0.0)
.approxEqual((ulong.max / 2) * 4096.0 + 4096^^2 / 2));
}

@safe pure nothrow unittest
{
static assert(is(typeof(sum([cast( byte)1])) == int));
static assert(is(typeof(sum([cast(ubyte)1])) == int));
static assert(is(typeof(sum([ 1, 2, 3, 4])) == int));
static assert(is(typeof(sum([ 1U, 2U, 3U, 4U])) == uint));
static assert(is(typeof(sum([ 1L, 2L, 3L, 4L])) == long));
static assert(is(typeof(sum([1UL, 2UL, 3UL, 4UL])) == ulong));

int[] empty;
assert(sum(empty) == 0);
assert(sum([42]) == 42);
assert(sum([42, 43]) == 42 + 43);
assert(sum([42, 43, 44]) == 42 + 43 + 44);
assert(sum([42, 43, 44, 45]) == 42 + 43 + 44 + 45);
}

@safe pure nothrow unittest
{
static assert(is(typeof(sum([1.0, 2.0, 3.0, 4.0])) == double));
static assert(is(typeof(sum([ 1F, 2F, 3F, 4F])) == double));
const(float[]) a = [1F, 2F, 3F, 4F];
static assert(is(typeof(sum(a)) == double));
const(float)[] b = [1F, 2F, 3F, 4F];
static assert(is(typeof(sum(a)) == double));

double[] empty;
assert(sum(empty) == 0);
assert(sum([42.]) == 42);
assert(sum([42., 43.]) == 42 + 43);
assert(sum([42., 43., 44.]) == 42 + 43 + 44);
assert(sum([42., 43., 44., 45.5]) == 42 + 43 + 44 + 45.5);
}

@safe pure nothrow unittest
{
import std.container;
static assert(is(typeof(sum(SList!float()[])) == double));
static assert(is(typeof(sum(SList!double()[])) == double));
static assert(is(typeof(sum(SList!real()[])) == real));

assert(sum(SList!double()[]) == 0);
assert(sum(SList!double(1)[]) == 1);
assert(sum(SList!double(1, 2)[]) == 1 + 2);
assert(sum(SList!double(1, 2, 3)[]) == 1 + 2 + 3);
assert(sum(SList!double(1, 2, 3, 4)[]) == 10);
}

@safe pure nothrow unittest // 12434
{
immutable a = [10, 20];
auto s1 = sum(a); // Error
auto s2 = a.map!(x => x).sum; // Error
}

unittest
{
import std.bigint;
import std.range;

immutable BigInt[] a = BigInt("1_000_000_000_000_000_000").repeat(10).array();
immutable ulong[] b = (ulong.max/2).repeat(10).array();
auto sa = a.sum();
auto sb = b.sum(BigInt(0)); //reduce ulongs into bigint
assert(sa == BigInt("10_000_000_000_000_000_000"));
assert(sb == (BigInt(ulong.max/2) * 10));
}

// uniq
/**
Lazily iterates unique consecutive elements of the given range (functionality
Expand Down
43 changes: 3 additions & 40 deletions std/numeric.d → std/numeric/package.d
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,8 @@ Distributed under the Boost Software License, Version 1.0.
*/
module std.numeric;

public import std.numeric.summation;

import std.complex;
import std.exception;
import std.math;
Expand Down Expand Up @@ -1593,6 +1595,7 @@ unittest
}
}


/**
Normalizes values in $(D range) by multiplying each element with a
number chosen such that values sum up to $(D sum). If elements in $(D
Expand Down Expand Up @@ -1656,46 +1659,6 @@ unittest
assert(a == [ 0.5, 0.5 ]);
}

/**
Compute the sum of binary logarithms of the input range $(D r).
The error of this method is much smaller than with a naive sum of log2.
*/
ElementType!Range sumOfLog2s(Range)(Range r)
if (isInputRange!Range && isFloatingPoint!(ElementType!Range))
{
long exp = 0;
Unqual!(typeof(return)) x = 1;
foreach (e; r)
{
if (e < 0)
return typeof(return).nan;
int lexp = void;
x *= frexp(e, lexp);
exp += lexp;
if (x < 0.5)
{
x *= 2;
exp--;
}
}
return exp + log2(x);
}

///
unittest
{
assert(sumOfLog2s(new double[0]) == 0);
assert(sumOfLog2s([0.0L]) == -real.infinity);
assert(sumOfLog2s([-0.0L]) == -real.infinity);
assert(sumOfLog2s([2.0L]) == 1);
assert(sumOfLog2s([-2.0L]).isNaN());
assert(sumOfLog2s([real.nan]).isNaN());
assert(sumOfLog2s([-real.nan]).isNaN());
assert(sumOfLog2s([real.infinity]) == real.infinity);
assert(sumOfLog2s([-real.infinity]).isNaN());
assert(sumOfLog2s([ 0.25, 0.25, 0.25, 0.125 ]) == -9);
}

/**
Computes $(LUCKY _entropy) of input range $(D r) in bits. This
Copy link
Contributor

@kuettler kuettler Feb 20, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Has this been released before? Is an alias needed to maintain compatibility?

Copy link
Member Author

@9il 9il Feb 20, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It will be released in 2.067. An alias is not needed since std.numeric.summation is publicly imported in std.numeric.

function assumes (without checking) that the values in $(D r) are all
Expand Down
Loading