Following the recent discussion of "reduce", I decided to try my hand at generic programming in C++.

Here's how I would write the standard deviation function in C

double stddev(double *x, int N)
{ int i;
double mean = 0;
double sumvar2 = 0;

for (i =0; i < N; i++)
mean += x[i];
mean /= N;

for (i =0; i < N; i++)
sumvar2 += (x[i] - mean)*(x[i]-mean);

return sqrt(sumvar2/(N-1));
}

Here's my attempt in C++.

#include <iostream>
#include <numeric>
#include <algorithm>
#include <vector>
#include <cmath>

template <class Iterator, typename func_t>
double stddev(Iterator first, Iterator last, func_t todouble)
{ typedef typename std::iterator_traits<Iterator>::value_type T;
double sum = std::transform_reduce(first, last, 0.0,
[](double a, double b){return a + b;}, todouble);
typename std::iterator_traits<Iterator>::difference_type N = std::distance(first, last);
double sume2 = std::transform_reduce(first, last, 0.0,
[](double a, double b){return a + b;},
[sum, N, todouble](T &value){double x = todouble(value);
return (x - sum/N)*(x-sum/N);}
);

return sqrt(sume2/(N-1));
}

A few comments. Whilst C has a reputation for efficiency, in fact
the C implementation will usually be less efficient. That's because
typically data comes in structs. You might want the standard deviation
of the employees' salaries. In C, that means allocating a temporary array,
traversing the employee array and picking out the salaries, calling stddev,
and freeing the temporary. In C++, it's easy to provide an access function.
(You can do this in C, but it makes the signature of the function too
complicated).

C++ is meant to be generic. But the standard deviation of a set of
integers is a real. How are you meant to express that? I'm sure there's
a way, but there didn't seem to be any tidy way of doing it. I resorted
to using a double as the return value. A lot of the complexity expresses
the simple idea that the function operates on the type that the iterators
point to. This is a legacy of the way C++ has grown by accretion.

transform_reduce can be parallelised. So in theory, once we've written
this function, we never need to write it again. When we get whizzy hardware
which is really good at fine-grained parallelisation, we rewrite transform_reduce() to take advantage of it.

In reality, it's hard to think of a case where calculating the standard
deviation would need to be faster than a commodity PC can do it
on a single core. In a real program, the x values have to mean something,
which usually means they are collected from the outside world. Collection
is so expensive and slow that rocessing speed is almost irrelevant.

However, of course, I chose this function because it is easy to implement
with transform_reduce.

I tried to develop this on replit. It failed to compile because the C++
version didn't have transform_reduce. So I tried my commandline g++.
Same story. g++ is mean to be a front end to clang. My IDE did compile
it. However setting up a project is heavyweight, so I cut and pasted
into another project. Finally it compiled. But this sort of messing about
with implementations greatly reduces the utility of code in the real world,
and is why I write a lot of code in a conservative subset of C89 and
the later versions. Then it almost always compiles in a new environment
without problems.

On 06/03/2022 11:17, Malcolm McLean wrote:
> Following the recent discussion of "reduce", I decided to try my hand at generic programming in C++.
>
> Here's how I would write the standard deviation function in C
>
> double stddev(double *x, int N)
> {
> int i;
> double mean = 0;
> double sumvar2 = 0;
>
> for (i =0; i < N; i++)
> mean += x[i];
> mean /= N;
>
> for (i =0; i < N; i++)
> sumvar2 += (x[i] - mean)*(x[i]-mean);
>
> return sqrt(sumvar2/(N-1));
> }
>
> Here's my attempt in C++.
....
> template <class Iterator, typename func_t>
> double stddev(Iterator first, Iterator last, func_t todouble)
> {
> typedef typename std::iterator_traits<Iterator>::value_type T;
> double sum = std::transform_reduce(first, last, 0.0,
> [](double a, double b){return a + b;}, todouble);
> typename std::iterator_traits<Iterator>::difference_type N = std::distance(first, last);
> double sume2 = std::transform_reduce(first, last, 0.0,
> [](double a, double b){return a + b;},
> [sum, N, todouble](T &value){double x = todouble(value);
> return (x - sum/N)*(x-sum/N);}
> );
>
> return sqrt(sume2/(N-1));
> }
>
> A few comments.

I don't need to read any further. Just look at the difference between
those two examples. I can easily follow the algorithm expressed as C,
but I can'make out anything in the C++. Is there even an algorithm in
that lot? All I can see is the last line that matches C's.

That this is supposed to be generic is no excuse. Generic code should
not need to look much different from non-generic code. With dynamic
typing, it would be exactly the same source code**.

According to the examples posted by David Brown, a generic version of
your the C function, written in C++, ought to look like this:

double stddev(auto *x, int N) {...}

x is a pointer to int, long long int, float or double. But I couldn't
get this to compile. If code needs to look like your C++, then it is a
failure IMO.

--

** I'm in the process of adding dynamic types to my otherwise static
systems language. The following example should be partly working by
tomorrow (when A is a simple List of variants), and fully working in a
week or two (allow any static arrays to be passed).

function stddev(var A)real =
real mean:=0, sumvar2:=0

for x in A do
mean +:= x
od

mean /:= A.len

for x in A do
sumvar2 +:= sqr(x-mean)
od

sqrt(sumvar2/(A.len-1))
end

The only downside is that expressions involving A and x need to use
dynamic type dispatch. But this does provide a generic routine if one is
needed.

In article <b2017949-b97f-4d1e-87c5-5b6ed79ecd10n@googlegroups.com>,
Malcolm McLean <malcolm.arthur.mclean@gmail.com> wrote:
>Following the recent discussion of "reduce", I decided to try my hand at generic
>programming in C++.

Off topic.

Try comp.lang.c++

--
BigBusiness types (aka, Republicans/Conservatives/Independents/Liberatarians/whatevers)
don't hate big government. They *love* big government as a means for them to get
rich, sucking off the public teat. What they don't like is *democracy* - you know,
like people actually having the right to vote and stuff like that.

Malcolm McLean <malcolm.arthur.mclean@gmail.com> writes:

> Following the recent discussion of "reduce", I decided to try my hand
> at generic programming in C++.
>
> Here's how I would write the standard deviation function in C
>
> double stddev(double *x, int N)
> {
> int i;
> double mean = 0;
> double sumvar2 = 0;
>
> for (i =0; i < N; i++)
> mean += x[i];
> mean /= N;
>
> for (i =0; i < N; i++)
> sumvar2 += (x[i] - mean)*(x[i]-mean);
>
> return sqrt(sumvar2/(N-1));
> }

If I was not worried about extreme data values, I'd write

double stddev(double *x, int N)
{
double sum = 0, sum_sq = 0;
for (int i = 0; i < N; i++) {
sum += x[i];
sum_sq += x[i] * x[i];
}
return sqrt((sum_sq - sum*sum/N) / (N-1));
}

> Here's my attempt in C++.
>
> #include <iostream>
> #include <numeric>
> #include <algorithm>
> #include <vector>
> #include <cmath>
>
> template <class Iterator, typename func_t>
> double stddev(Iterator first, Iterator last, func_t todouble)
> {
> typedef typename std::iterator_traits<Iterator>::value_type T;
> double sum = std::transform_reduce(first, last, 0.0,
> [](double a, double b){return a + b;}, todouble);
> typename std::iterator_traits<Iterator>::difference_type N = std::distance(first, last);
> double sume2 = std::transform_reduce(first, last, 0.0,
> [](double a, double b){return a + b;},
> [sum, N, todouble](T &value){double x = todouble(value);
> return (x - sum/N)*(x-sum/N);}
> );
>
> return sqrt(sume2/(N-1));
> }
>
> A few comments. Whilst C has a reputation for efficiency, in fact
> the C implementation will usually be less efficient. That's because
> typically data comes in structs. You might want the standard deviation
> of the employees' salaries. In C, that means allocating a temporary array,
> traversing the employee array and picking out the salaries, calling stddev,
> and freeing the temporary. In C++, it's easy to provide an access function.
> (You can do this in C, but it makes the signature of the function too
> complicated).
>
> C++ is meant to be generic. But the standard deviation of a set of
> integers is a real. How are you meant to express that? I'm sure there's
> a way, but there didn't seem to be any tidy way of doing it. I resorted
> to using a double as the return value. A lot of the complexity expresses
> the simple idea that the function operates on the type that the iterators
> point to. This is a legacy of the way C++ has grown by accretion.

You don't need to do that. Just tidying up a bit and giving the code
some room to breathe, I get a more readable template:

template <class Iterator, typename func_t>
double stddev(Iterator first, Iterator last, func_t todouble)
{
auto add_func = [](double a, double b){ return a + b; };

double sum = std::transform_reduce(first, last, 0.0, add_func, todouble);
auto N = std::distance(first, last);

double sume2 = std::transform_reduce(
first, last, 0.0, add_func,
[&](auto &value){
double x = todouble(value);
return (x - sum/N) * (x - sum/N);
});

return sqrt(sume2/(N-1));
}

But I don't know why you are worried about the conversion from int to
double:

template <class Iterator>
double stddevY(Iterator first, Iterator last)
{
auto add_func = [](double a, double b){ return a + b; };

double sum = std::accumulate(first, last, 0.0, add_func);
auto N = std::distance(first, last);

double sume2 = std::transform_reduce(
first, last, 0.0, add_func,
[&](double x){ return (x - sum/N) * (x - sum/N); }
);

return sqrt(sume2/(N-1));
}

But surely you could just use accumulate:

template <class Iterator>
double stddev(Iterator first, Iterator last)
{
double sum = std::accumulate(first, last, 0);
double sum_sq = std::accumulate(
first, last, 0, [](double sqs, double x){ return sqs + x*x; });
double N = std::distance(first, last);
return sqrt((sum_sq - sum*sum/N)/(N-1));
}

> transform_reduce can be parallelised.

I assume that accumulate can also be parallelised, but if not you can
just use something that can be like transform_reduce with the identity
function.

> So in theory, once we've written this function, we never need to write
> it again. When we get whizzy hardware which is really good at
> fine-grained parallelisation, we rewrite transform_reduce() to take
> advantage of it.
>
> In reality, it's hard to think of a case where calculating the standard
> deviation would need to be faster than a commodity PC can do it
> on a single core. In a real program, the x values have to mean something,
> which usually means they are collected from the outside world. Collection
> is so expensive and slow that rocessing speed is almost irrelevant.

A long time ago I was asked to write a function that (in effect)
averaged the standard deviations of every subset of a data set (with
some minimum size which I can't remember). It was, as you would expect,
very slow! But your main point stands: parallelisation of stddev would
not have helped much in this case, especially as I just kept adjusting a
running sum and sum of squares rather than calling stddev a lot.

However, my point is that can't say speed of F(X) is almost irrelevant
just because collecting X is expensive.

Also, collecting X is not always expensive. A network monitor that
collects real-time trip times and displays the standard deviation comes
to mind. Also a computer simulation of something might run at many,
many times the speed of the events in the real world and therefore end
up accumulating lots of data very fast.

> However, of course, I chose this function because it is easy to implement
> with transform_reduce.

It's an odd thing to post in comp.lang.c. You are really asking for
expert C++ advice on generics. Take what I say above with a pinch of
salt. It could be very bad C++ advice.

--
Ben.

Bart <bc@freeuk.com> writes:

> That this is supposed to be generic is no excuse. Generic code should
> not need to look much different from non-generic code. With dynamic
> typing, it would be exactly the same source code**.
>
> According to the examples posted by David Brown, a generic version of
> your the C function, written in C++, ought to look like this:
>
> double stddev(auto *x, int N) {...}

double stddev(const auto *x, int N)
{
double sum = 0, sum_sq = 0;
for (int i = 0; i < N; i++) {
sum += x[i];
sum_sq += x[i]*x[i];
}
return sqrt((sum_sq - sum*sum/N)/(N-1));
}

But then you then probably loose the automatic parallelisation that the
standard library might give you.

> x is a pointer to int, long long int, float or double. But I couldn't
> get this to compile. If code needs to look like your C++, then it is a
> failure IMO.

The trouble is that the C++ code shown was intended to do more than the
C, and not just in terms of handling generic types. Also, you can't
judge C++ by a single example from one person. For example, I don't
think

double stddev(const auto &first, const auto &last)
{
double sum = std::accumulate(first, last);
double sum_sq = std::accumulate(
first, last, [](double sqs, double x){ return sqs + x*x; });
double N = std::distance(first, last);
return sqrt((sum_sq - sum*sum/N)/(N-1));
}

is particularly unclear.

--
Ben.

On Sunday, 6 March 2022 at 15:05:04 UTC, Ben Bacarisse wrote:
> Bart <b...@freeuk.com> writes:
>
> > That this is supposed to be generic is no excuse. Generic code should
> > not need to look much different from non-generic code. With dynamic
> > typing, it would be exactly the same source code**.
> >
> > According to the examples posted by David Brown, a generic version of
> > your the C function, written in C++, ought to look like this:
> >
> > double stddev(auto *x, int N) {...}
> double stddev(const auto *x, int N)
> {
> double sum = 0, sum_sq = 0;
> for (int i = 0; i < N; i++) {
> sum += x[i];
> sum_sq += x[i]*x[i];
> }
> return sqrt((sum_sq - sum*sum/N)/(N-1));
> }
> But then you then probably loose the automatic parallelisation that the
> standard library might give you.
> > x is a pointer to int, long long int, float or double. But I couldn't
> > get this to compile. If code needs to look like your C++, then it is a
> > failure IMO.
> The trouble is that the C++ code shown was intended to do more than the
> C, and not just in terms of handling generic types. Also, you can't
> judge C++ by a single example from one person. For example, I don't
> think
>
> double stddev(const auto &first, const auto &last)
> {
> double sum = std::accumulate(first, last);
> double sum_sq = std::accumulate(
> first, last, [](double sqs, double x){ return sqs + x*x; });
> double N = std::distance(first, last);
> return sqrt((sum_sq - sum*sum/N)/(N-1));
> }
>
> is particularly unclear.
>
Yes, but that's not generic. If first and last iterator over employees, the "+'
operator has to be defined to add two employees together. That's unlikely
to make any sort of sense.
In fact we'll want that standard deviation of the salaries. So we can pass
in a trivial little lambda [](Employee &employee){return employee.salary;}
If we've not another variable associated with the employees, like hours
worked, we can get the standard deviation of that just by changing the lambda.

reduce is better than accumulate because accumulate requires left-to
right evaluation of the operands. That's not necessary for a standard deviation,
unless annoyingly, we need every bit of precision (floating point arithmetic
isn't exactly associative, though it's near enough for almost all purposes).

You've used a slightly different formula for the standard deviation., which
is more efficient and creates tighter code.

Ben Bacarisse <ben.usenet@bsb.me.uk> writes:
>Bart <bc@freeuk.com> writes:
>
>> That this is supposed to be generic is no excuse. Generic code should
>> not need to look much different from non-generic code. With dynamic
>> typing, it would be exactly the same source code**.
>>
>> According to the examples posted by David Brown, a generic version of
>> your the C function, written in C++, ought to look like this:
>>
>> double stddev(auto *x, int N) {...}
>
> double stddev(const auto *x, int N)
> {
> double sum = 0, sum_sq = 0;
> for (int i = 0; i < N; i++) {
> sum += x[i];
> sum_sq += x[i]*x[i];
> }
> return sqrt((sum_sq - sum*sum/N)/(N-1));
> }
>
>But then you then probably loose the automatic parallelisation that the
>standard library might give you.

Unless you're using one of the ubiquitous automatically vectorizing
compilers.

scott@slp53.sl.home (Scott Lurndal) writes:

> Ben Bacarisse <ben.usenet@bsb.me.uk> writes:
>>Bart <bc@freeuk.com> writes:
>>
>>> That this is supposed to be generic is no excuse. Generic code should
>>> not need to look much different from non-generic code. With dynamic
>>> typing, it would be exactly the same source code**.
>>>
>>> According to the examples posted by David Brown, a generic version of
>>> your the C function, written in C++, ought to look like this:
>>>
>>> double stddev(auto *x, int N) {...}
>>
>> double stddev(const auto *x, int N)
>> {
>> double sum = 0, sum_sq = 0;
>> for (int i = 0; i < N; i++) {
>> sum += x[i];
>> sum_sq += x[i]*x[i];
>> }
>> return sqrt((sum_sq - sum*sum/N)/(N-1));
>> }
>>
>>But then you then probably loose the automatic parallelisation that the
>>standard library might give you.
>
> Unless you're using one of the ubiquitous automatically vectorizing
> compilers.

Do they work well with type-generic code, given that they won't always
know how + and * are defined?

--
Ben.

On 06/03/2022 16:36, Ben Bacarisse wrote:
> scott@slp53.sl.home (Scott Lurndal) writes:
>
>> Ben Bacarisse <ben.usenet@bsb.me.uk> writes:
>>> Bart <bc@freeuk.com> writes:
>>>
>>>> That this is supposed to be generic is no excuse. Generic code should
>>>> not need to look much different from non-generic code. With dynamic
>>>> typing, it would be exactly the same source code**.
>>>>
>>>> According to the examples posted by David Brown, a generic version of
>>>> your the C function, written in C++, ought to look like this:
>>>>
>>>> double stddev(auto *x, int N) {...}
>>>
>>> double stddev(const auto *x, int N)
>>> {
>>> double sum = 0, sum_sq = 0;
>>> for (int i = 0; i < N; i++) {
>>> sum += x[i];
>>> sum_sq += x[i]*x[i];
>>> }
>>> return sqrt((sum_sq - sum*sum/N)/(N-1));
>>> }
>>>
>>> But then you then probably loose the automatic parallelisation that the
>>> standard library might give you.
>>
>> Unless you're using one of the ubiquitous automatically vectorizing
>> compilers.
>
> Do they work well with type-generic code, given that they won't always
> know how + and * are defined?
>

They know how + and * are defined when compiling the code - it's
compile-time polymorphism. Compilers can do exactly as much with such
templated (or implicitly templated with "auto") code as they can with
explicit types.

But I think there's a bit of mixup here between "parallelisation" and
"vectorisation" - they don't mean the same thing. Compilers can do
automatic vectorisation on code like this, but automatic parallelisation
is rarer. In theory, the C++ algorithm library functions support
parallel execution strategy - in practice they usually don't (i.e., you
can ask for parallel strategies, but you still get sequential ones).

Malcolm McLean <malcolm.arthur.mclean@gmail.com> writes:

> On Sunday, 6 March 2022 at 15:05:04 UTC, Ben Bacarisse wrote:
>> Bart <b...@freeuk.com> writes:
>>
>> > That this is supposed to be generic is no excuse. Generic code should
>> > not need to look much different from non-generic code. With dynamic
>> > typing, it would be exactly the same source code**.
>> >
>> > According to the examples posted by David Brown, a generic version of
>> > your the C function, written in C++, ought to look like this:
>> >
>> > double stddev(auto *x, int N) {...}
>> double stddev(const auto *x, int N)
>> {
>> double sum = 0, sum_sq = 0;
>> for (int i = 0; i < N; i++) {
>> sum += x[i];
>> sum_sq += x[i]*x[i];
>> }
>> return sqrt((sum_sq - sum*sum/N)/(N-1));
>> }
>> But then you then probably loose the automatic parallelisation that the
>> standard library might give you.
>> > x is a pointer to int, long long int, float or double. But I couldn't
>> > get this to compile. If code needs to look like your C++, then it is a
>> > failure IMO.
>> The trouble is that the C++ code shown was intended to do more than the
>> C, and not just in terms of handling generic types. Also, you can't
>> judge C++ by a single example from one person. For example, I don't
>> think
>>
>> double stddev(const auto &first, const auto &last)
>> {
>> double sum = std::accumulate(first, last);
>> double sum_sq = std::accumulate(
>> first, last, [](double sqs, double x){ return sqs + x*x; });
>> double N = std::distance(first, last);
>> return sqrt((sum_sq - sum*sum/N)/(N-1));
>> }
>>
>> is particularly unclear.
>>
> Yes, but that's not generic. If first and last iterator over employees, the "+'
> operator has to be defined to add two employees together. That's unlikely
> to make any sort of sense.

I think that's a good thing. Standard deviation is an arithmetic
property.

> In fact we'll want that standard deviation of the salaries. So we can pass
> in a trivial little lambda [](Employee &employee){return
> employee.salary;}

I think you should use the ranges::view types to pass iterators to that
data. Otherwise, every arithmetic operation has to acquire a little
function to do the conversion, and the most common cases need a trivial
function that does no more than the language does by default.

> If we've not another variable associated with the employees, like hours
> worked, we can get the standard deviation of that just by changing the
> lambda.

I think views can handle this just fine without complicating every
arithmetic operation. But this is a classic example of why you should
be posting where experts can reply -- I may be wrong.

> reduce is better than accumulate because accumulate requires left-to
> right evaluation of the operands. That's not necessary for a standard
> deviation, unless annoyingly, we need every bit of precision (floating
> point arithmetic isn't exactly associative, though it's near enough
> for almost all purposes).

Sure. You can just plug in reduce if you don't mid the risk of losing
control over the order:

double stddev(const auto &first, const auto &last)
{
double sum = std::reduce(first, last);
double sum_sq = std::reduce(
first, last, [](double sqs, double x){ return sqs + x*x; });
double N = std::distance(first, last);
return sqrt((sum_sq - sum*sum/N)/(N-1));
}

> You've used a slightly different formula for the standard deviation., which
> is more efficient and creates tighter code.

It comes with a risk, though, because the sum of squares grows much
faster than the sum of variances.

--
Ben.

On Sunday, 6 March 2022 at 15:54:29 UTC, David Brown wrote:
> On 06/03/2022 16:36, Ben Bacarisse wrote:
> > sc...@slp53.sl.home (Scott Lurndal) writes:
> >
> >> Ben Bacarisse <ben.u...@bsb.me.uk> writes:
> >>> Bart <b...@freeuk.com> writes:
> >>>
> >>>> That this is supposed to be generic is no excuse. Generic code should
> >>>> not need to look much different from non-generic code. With dynamic
> >>>> typing, it would be exactly the same source code**.
> >>>>
> >>>> According to the examples posted by David Brown, a generic version of
> >>>> your the C function, written in C++, ought to look like this:
> >>>>
> >>>> double stddev(auto *x, int N) {...}
> >>>
> >>> double stddev(const auto *x, int N)
> >>> {
> >>> double sum = 0, sum_sq = 0;
> >>> for (int i = 0; i < N; i++) {
> >>> sum += x[i];
> >>> sum_sq += x[i]*x[i];
> >>> }
> >>> return sqrt((sum_sq - sum*sum/N)/(N-1));
> >>> }
> >>>
> >>> But then you then probably loose the automatic parallelisation that the
> >>> standard library might give you.
> >>
> >> Unless you're using one of the ubiquitous automatically vectorizing
> >> compilers.
> >
> > Do they work well with type-generic code, given that they won't always
> > know how + and * are defined?
> >
> They know how + and * are defined when compiling the code - it's
> compile-time polymorphism. Compilers can do exactly as much with such
> templated (or implicitly templated with "auto") code as they can with
> explicit types.
>
> But I think there's a bit of mixup here between "parallelisation" and
> "vectorisation" - they don't mean the same thing. Compilers can do
> automatic vectorisation on code like this, but automatic parallelisation
> is rarer. In theory, the C++ algorithm library functions support
> parallel execution strategy - in practice they usually don't (i.e., you
> can ask for parallel strategies, but you still get sequential ones).
>
The idea is to future-proof the code. Currently most readily-available
platforms can't do the fine-grained parallelisation you'd need to make
a parallel sddev worthwhile. But that might change.

However you could argue that a simple C for loop might be easier for
a future compiler to paralleise, and thus more future-proof than the C++
transform_reduce.

David Brown <david.brown@hesbynett.no> writes:

> On 06/03/2022 16:36, Ben Bacarisse wrote:
>> scott@slp53.sl.home (Scott Lurndal) writes:
>>
>>> Ben Bacarisse <ben.usenet@bsb.me.uk> writes:
>>>> Bart <bc@freeuk.com> writes:
>>>>
>>>>> That this is supposed to be generic is no excuse. Generic code should
>>>>> not need to look much different from non-generic code. With dynamic
>>>>> typing, it would be exactly the same source code**.
>>>>>
>>>>> According to the examples posted by David Brown, a generic version of
>>>>> your the C function, written in C++, ought to look like this:
>>>>>
>>>>> double stddev(auto *x, int N) {...}
>>>>
>>>> double stddev(const auto *x, int N)
>>>> {
>>>> double sum = 0, sum_sq = 0;
>>>> for (int i = 0; i < N; i++) {
>>>> sum += x[i];
>>>> sum_sq += x[i]*x[i];
>>>> }
>>>> return sqrt((sum_sq - sum*sum/N)/(N-1));
>>>> }
>>>>
>>>> But then you then probably loose the automatic parallelisation that the
>>>> standard library might give you.
>>>
>>> Unless you're using one of the ubiquitous automatically vectorizing
>>> compilers.
>>
>> Do they work well with type-generic code, given that they won't always
>> know how + and * are defined?
>
> They know how + and * are defined when compiling the code - it's
> compile-time polymorphism. Compilers can do exactly as much with such
> templated (or implicitly templated with "auto") code as they can with
> explicit types.

But do they? If you include link-time with compile-time, sure the
compiler knows what * is, but it's usually too late to change what the
loop compiled to. Does the code generator write two forms and the
linker picks one or the other when certain conditions are met?

Or maybe I've missed your point altogether.

> But I think there's a bit of mixup here between "parallelisation" and
> "vectorisation" - they don't mean the same thing. Compilers can do
> automatic vectorisation on code like this, but automatic parallelisation
> is rarer.

I would say that vectorisation is a form of parallelisation -- the
latter being the generic term for whole range of specialist
optimisations. I may well have used one or other of them incorrectly,
but I know they are different!

> In theory, the C++ algorithm library functions support
> parallel execution strategy - in practice they usually don't (i.e., you
> can ask for parallel strategies, but you still get sequential ones).

That's a good point. Explicit loops (that Malcolm was, I think, trying
to avoid) may well be better targets for optimisation than the library
functions he was leaning towards.

--
Ben.

I use this code:

template<typename Iterator>
requires std::is_convertible_v<typename
std::iterator_traits<Iterator>::value_type, double>
double stddev( Iterator begin, Iterator end )
{ using namespace std;
using T = typename iterator_traits<Iterator>::value_type;
ptrdiff_t n;
double mean = 0.0;
if constexpr( random_access_iterator<Iterator> )
n = std::distance( begin, end ),
for_each( begin, end, [&]( T &value ) { mean += (double)value; } );
else
n = 0,
for_each( begin, end, [&]( T &value ) { ++n; mean += (double)value; } );
mean /= n;
double sqSum = 0.0;
for_each( begin, end, [&]( T &value ) { double v = value; sqSum += (v -
mean) * (v - mean); } );
return sqrt( sqSum / n );
}

I don't use reduce because the type inside the array might
be asymmetrical to double.

Am 06.03.2022 um 13:42 schrieb Kenny McCormack:
> In article <b2017949-b97f-4d1e-87c5-5b6ed79ecd10n@googlegroups.com>,
> Malcolm McLean <malcolm.arthur.mclean@gmail.com> wrote:
>> Following the recent discussion of "reduce", I decided to try my hand at generic
>> programming in C++.
>
> Off topic.
>
> Try comp.lang.c++

Comparison of C against other languages is on-topic here.

This might be better:

template<typename Iterator>
requires std::is_convertible_v<typename
std::iterator_traits<Iterator>::value_type, double> &&
std::input_iterator<Iterator>
double stddev( Iterator begin, Iterator end )
{ using namespace std;
using T = typename iterator_traits<Iterator>::value_type;
ptrdiff_t n;
double mean = 0.0;
if constexpr( random_access_iterator<Iterator> )
n = std::distance( begin, end ),
for_each( begin, end, [&]( T &value ) { mean += (double)value; } );
else
n = 0,
for_each( begin, end, [&]( T &value ) { ++n; mean += (double)value; } );
mean /= n;
double sqSum = 0.0;
for_each( begin, end, [&]( T &value ) { double diff = (double)value -
mean; sqSum += diff * diff; } );
return sqrt( sqSum / n );
}

On Sunday, 6 March 2022 at 15:59:14 UTC, Ben Bacarisse wrote:
> Malcolm McLean <malcolm.ar...@gmail.com> writes:
>
> > On Sunday, 6 March 2022 at 15:05:04 UTC, Ben Bacarisse wrote:
> >> Bart <b...@freeuk.com> writes:
> >>
> >> > That this is supposed to be generic is no excuse. Generic code should
> >> > not need to look much different from non-generic code. With dynamic
> >> > typing, it would be exactly the same source code**.
> >> >
> >> > According to the examples posted by David Brown, a generic version of
> >> > your the C function, written in C++, ought to look like this:
> >> >
> >> > double stddev(auto *x, int N) {...}
> >> double stddev(const auto *x, int N)
> >> {
> >> double sum = 0, sum_sq = 0;
> >> for (int i = 0; i < N; i++) {
> >> sum += x[i];
> >> sum_sq += x[i]*x[i];
> >> }
> >> return sqrt((sum_sq - sum*sum/N)/(N-1));
> >> }
> >> But then you then probably loose the automatic parallelisation that the
> >> standard library might give you.
> >> > x is a pointer to int, long long int, float or double. But I couldn't
> >> > get this to compile. If code needs to look like your C++, then it is a
> >> > failure IMO.
> >> The trouble is that the C++ code shown was intended to do more than the
> >> C, and not just in terms of handling generic types. Also, you can't
> >> judge C++ by a single example from one person. For example, I don't
> >> think
> >>
> >> double stddev(const auto &first, const auto &last)
> >> {
> >> double sum = std::accumulate(first, last);
> >> double sum_sq = std::accumulate(
> >> first, last, [](double sqs, double x){ return sqs + x*x; });
> >> double N = std::distance(first, last);
> >> return sqrt((sum_sq - sum*sum/N)/(N-1));
> >> }
> >>
> >> is particularly unclear.
> >>
> > Yes, but that's not generic. If first and last iterator over employees, the "+'
> > operator has to be defined to add two employees together. That's unlikely
> > to make any sort of sense.
> I think that's a good thing. Standard deviation is an arithmetic
> property.
> > In fact we'll want that standard deviation of the salaries. So we can pass
> > in a trivial little lambda [](Employee &employee){return
> > employee.salary;}
> I think you should use the ranges::view types to pass iterators to that
> data. Otherwise, every arithmetic operation has to acquire a little
> function to do the conversion, and the most common cases need a trivial
> function that does no more than the language does by default.
> > If we've not another variable associated with the employees, like hours
> > worked, we can get the standard deviation of that just by changing the
> > lambda.
> I think views can handle this just fine without complicating every
> arithmetic operation. But this is a classic example of why you should
> be posting where experts can reply -- I may be wrong.
> > reduce is better than accumulate because accumulate requires left-to
> > right evaluation of the operands. That's not necessary for a standard
> > deviation, unless annoyingly, we need every bit of precision (floating
> > point arithmetic isn't exactly associative, though it's near enough
> > for almost all purposes).
> Sure. You can just plug in reduce if you don't mid the risk of losing
> control over the order:
> double stddev(const auto &first, const auto &last)
> {
> double sum = std::reduce(first, last);
> double sum_sq = std::reduce(
> first, last, [](double sqs, double x){ return sqs + x*x; });
> double N = std::distance(first, last);
> return sqrt((sum_sq - sum*sum/N)/(N-1));
> }
>
This won't work. You want the sum of squares. But if you evaluate left to right,
the first element is not squared. That's easily fixed by passing in an init.
However if you use "reduce" instead of "accumulate" then the problem is
much more serious. Rndom pairs might not have the first operand squared.

You must use reduce_transform.

On 06/03/2022 18:03, Ben Bacarisse wrote:
> David Brown <david.brown@hesbynett.no> writes:
>
>> On 06/03/2022 16:36, Ben Bacarisse wrote:
>>> scott@slp53.sl.home (Scott Lurndal) writes:
>>>
>>>> Ben Bacarisse <ben.usenet@bsb.me.uk> writes:
>>>>> Bart <bc@freeuk.com> writes:
>>>>>
>>>>>> That this is supposed to be generic is no excuse. Generic code should
>>>>>> not need to look much different from non-generic code. With dynamic
>>>>>> typing, it would be exactly the same source code**.
>>>>>>
>>>>>> According to the examples posted by David Brown, a generic version of
>>>>>> your the C function, written in C++, ought to look like this:
>>>>>>
>>>>>> double stddev(auto *x, int N) {...}
>>>>>
>>>>> double stddev(const auto *x, int N)
>>>>> {
>>>>> double sum = 0, sum_sq = 0;
>>>>> for (int i = 0; i < N; i++) {
>>>>> sum += x[i];
>>>>> sum_sq += x[i]*x[i];
>>>>> }
>>>>> return sqrt((sum_sq - sum*sum/N)/(N-1));
>>>>> }
>>>>>
>>>>> But then you then probably loose the automatic parallelisation that the
>>>>> standard library might give you.
>>>>
>>>> Unless you're using one of the ubiquitous automatically vectorizing
>>>> compilers.
>>>
>>> Do they work well with type-generic code, given that they won't always
>>> know how + and * are defined?
>>
>> They know how + and * are defined when compiling the code - it's
>> compile-time polymorphism. Compilers can do exactly as much with such
>> templated (or implicitly templated with "auto") code as they can with
>> explicit types.
>
> But do they? If you include link-time with compile-time, sure the
> compiler knows what * is, but it's usually too late to change what the
> loop compiled to. Does the code generator write two forms and the
> linker picks one or the other when certain conditions are met?
>
> Or maybe I've missed your point altogether.
>

When you use "auto" in a function declaration, it is just like a
template (indeed it /is/ a template, written with a shorter syntax).
Code is not generated for it at all until the compiler knows the type
being used, and the details of it. Of course, if the type is a class
with operators + and *, which are declared at the point of template
instantiation but defined externally, then the compiler won't be able to
do anything special with them. But if the definitions of the operators
are known when the template is instantiated (and that's likely to be the
case for arithmetic types), then the compiler has all the information it
needs at compile time.

>> But I think there's a bit of mixup here between "parallelisation" and
>> "vectorisation" - they don't mean the same thing. Compilers can do
>> automatic vectorisation on code like this, but automatic parallelisation
>> is rarer.
>
> I would say that vectorisation is a form of parallelisation -- the
> latter being the generic term for whole range of specialist
> optimisations. I may well have used one or other of them incorrectly,
> but I know they are different!
>

I would use "vectorisation" to refer to SIMD operations in the same
thread, and "parallelisation" to refer to splitting the calculation into
parts and running them on separate threads on separate cores. There is
an overlap, but I think the distinction is clear.

>> In theory, the C++ algorithm library functions support
>> parallel execution strategy - in practice they usually don't (i.e., you
>> can ask for parallel strategies, but you still get sequential ones).
>
> That's a good point. Explicit loops (that Malcolm was, I think, trying
> to avoid) may well be better targets for optimisation than the library
> functions he was leaning towards.
>

Explicit loops are often easier for optimisation when the loop count is
small - function calls can be better when the loop count is high. (You
get the same thing for memcpy() - a memcpy with a small known count can
be turned into a few move instructions, but a library function can do
things like alignment of cache lines to optimise large movements better.)

In particular, an explicit loop /might/ be vectorised, depending of
flags and the details of the code. A function call might lead to code
that has different vectorised versions for different processors - larger
overhead, but fast operation regardless of flags.

Malcolm McLean <malcolm.arthur.mclean@gmail.com> writes:

> On Sunday, 6 March 2022 at 15:59:14 UTC, Ben Bacarisse wrote:
>> Malcolm McLean <malcolm.ar...@gmail.com> writes:
>>
>> > On Sunday, 6 March 2022 at 15:05:04 UTC, Ben Bacarisse wrote:
>> >> Bart <b...@freeuk.com> writes:
>> >>
>> >> > That this is supposed to be generic is no excuse. Generic code should
>> >> > not need to look much different from non-generic code. With dynamic
>> >> > typing, it would be exactly the same source code**.
>> >> >
>> >> > According to the examples posted by David Brown, a generic version of
>> >> > your the C function, written in C++, ought to look like this:
>> >> >
>> >> > double stddev(auto *x, int N) {...}
>> >> double stddev(const auto *x, int N)
>> >> {
>> >> double sum = 0, sum_sq = 0;
>> >> for (int i = 0; i < N; i++) {
>> >> sum += x[i];
>> >> sum_sq += x[i]*x[i];
>> >> }
>> >> return sqrt((sum_sq - sum*sum/N)/(N-1));
>> >> }
>> >> But then you then probably loose the automatic parallelisation that the
>> >> standard library might give you.
>> >> > x is a pointer to int, long long int, float or double. But I couldn't
>> >> > get this to compile. If code needs to look like your C++, then it is a
>> >> > failure IMO.
>> >> The trouble is that the C++ code shown was intended to do more than the
>> >> C, and not just in terms of handling generic types. Also, you can't
>> >> judge C++ by a single example from one person. For example, I don't
>> >> think
>> >>
>> >> double stddev(const auto &first, const auto &last)
>> >> {
>> >> double sum = std::accumulate(first, last);
>> >> double sum_sq = std::accumulate(
>> >> first, last, [](double sqs, double x){ return sqs + x*x; });
>> >> double N = std::distance(first, last);
>> >> return sqrt((sum_sq - sum*sum/N)/(N-1));
>> >> }
>> >>
>> >> is particularly unclear.
>> >>
>> > Yes, but that's not generic. If first and last iterator over employees, the "+'
>> > operator has to be defined to add two employees together. That's unlikely
>> > to make any sort of sense.
>> I think that's a good thing. Standard deviation is an arithmetic
>> property.
>> > In fact we'll want that standard deviation of the salaries. So we can pass
>> > in a trivial little lambda [](Employee &employee){return
>> > employee.salary;}
>> I think you should use the ranges::view types to pass iterators to that
>> data. Otherwise, every arithmetic operation has to acquire a little
>> function to do the conversion, and the most common cases need a trivial
>> function that does no more than the language does by default.
>> > If we've not another variable associated with the employees, like hours
>> > worked, we can get the standard deviation of that just by changing the
>> > lambda.
>> I think views can handle this just fine without complicating every
>> arithmetic operation. But this is a classic example of why you should
>> be posting where experts can reply -- I may be wrong.
>> > reduce is better than accumulate because accumulate requires left-to
>> > right evaluation of the operands. That's not necessary for a standard
>> > deviation, unless annoyingly, we need every bit of precision (floating
>> > point arithmetic isn't exactly associative, though it's near enough
>> > for almost all purposes).
>> Sure. You can just plug in reduce if you don't mid the risk of losing
>> control over the order:
>> double stddev(const auto &first, const auto &last)
>> {
>> double sum = std::reduce(first, last);
>> double sum_sq = std::reduce(
>> first, last, [](double sqs, double x){ return sqs + x*x; });
>> double N = std::distance(first, last);
>> return sqrt((sum_sq - sum*sum/N)/(N-1));
>> }
>>
> This won't work. You want the sum of squares. But if you evaluate left
> to right, the first element is not squared. That's easily fixed by
> passing in an init.

Thanks. I removed the 0 init from the accumulate version without
thinking. But...

> However if you use "reduce" instead of "accumulate" then the problem is
> much more serious. Rndom pairs might not have the first operand
> squared.

Ah yes. The binary operator must be associative and commutative.

> You must use reduce_transform.

(transform_reduce)

So I read the documentation (as I clearly should have done before!) and,
unless I've messed up again, an even neater solution presents itself:

double stddev(const auto &first, const auto &last)
{
double sum = std::reduce(first, last);
double sum_sq = std::transform_reduce(first, last, first);
double N = std::distance(first, last);
return sqrt((sum_sq - sum*sum/N)/(N-1));
}

because std::transform_reduce(first, last, first) is essentially a
parallel inner product of the sequence with itself. It's possible,
given the prevalence of inner products, that this might be special-cased
and particularly fast on some hardware.

BTW, can you do anything about your Usenet client messing up the
formatting? It's particularly unhelpful with code, but "paragraph
breaks" get lost too.

--
Ben.

On Sunday, 6 March 2022 at 19:42:41 UTC, Ben Bacarisse wrote:
> Malcolm McLean <malcolm.ar...@gmail.com> writes:
>
> > On Sunday, 6 March 2022 at 15:59:14 UTC, Ben Bacarisse wrote:
> >> Malcolm McLean <malcolm.ar...@gmail.com> writes:
> >>
> >> > On Sunday, 6 March 2022 at 15:05:04 UTC, Ben Bacarisse wrote:
> >> >> Bart <b...@freeuk.com> writes:
> >> >>
> >> >> > That this is supposed to be generic is no excuse. Generic code should
> >> >> > not need to look much different from non-generic code. With dynamic
> >> >> > typing, it would be exactly the same source code**.
> >> >> >
> >> >> > According to the examples posted by David Brown, a generic version of
> >> >> > your the C function, written in C++, ought to look like this:
> >> >> >
> >> >> > double stddev(auto *x, int N) {...}
> >> >> double stddev(const auto *x, int N)
> >> >> {
> >> >> double sum = 0, sum_sq = 0;
> >> >> for (int i = 0; i < N; i++) {
> >> >> sum += x[i];
> >> >> sum_sq += x[i]*x[i];
> >> >> }
> >> >> return sqrt((sum_sq - sum*sum/N)/(N-1));
> >> >> }
> >> >> But then you then probably loose the automatic parallelisation that the
> >> >> standard library might give you.
> >> >> > x is a pointer to int, long long int, float or double. But I couldn't
> >> >> > get this to compile. If code needs to look like your C++, then it is a
> >> >> > failure IMO.
> >> >> The trouble is that the C++ code shown was intended to do more than the
> >> >> C, and not just in terms of handling generic types. Also, you can't
> >> >> judge C++ by a single example from one person. For example, I don't
> >> >> think
> >> >>
> >> >> double stddev(const auto &first, const auto &last)
> >> >> {
> >> >> double sum = std::accumulate(first, last);
> >> >> double sum_sq = std::accumulate(
> >> >> first, last, [](double sqs, double x){ return sqs + x*x; });
> >> >> double N = std::distance(first, last);
> >> >> return sqrt((sum_sq - sum*sum/N)/(N-1));
> >> >> }
> >> >>
> >> >> is particularly unclear.
> >> >>
> >> > Yes, but that's not generic. If first and last iterator over employees, the "+'
> >> > operator has to be defined to add two employees together. That's unlikely
> >> > to make any sort of sense.
> >> I think that's a good thing. Standard deviation is an arithmetic
> >> property.
> >> > In fact we'll want that standard deviation of the salaries. So we can pass
> >> > in a trivial little lambda [](Employee &employee){return
> >> > employee.salary;}
> >> I think you should use the ranges::view types to pass iterators to that
> >> data. Otherwise, every arithmetic operation has to acquire a little
> >> function to do the conversion, and the most common cases need a trivial
> >> function that does no more than the language does by default.
> >> > If we've not another variable associated with the employees, like hours
> >> > worked, we can get the standard deviation of that just by changing the
> >> > lambda.
> >> I think views can handle this just fine without complicating every
> >> arithmetic operation. But this is a classic example of why you should
> >> be posting where experts can reply -- I may be wrong.
> >> > reduce is better than accumulate because accumulate requires left-to
> >> > right evaluation of the operands. That's not necessary for a standard
> >> > deviation, unless annoyingly, we need every bit of precision (floating
> >> > point arithmetic isn't exactly associative, though it's near enough
> >> > for almost all purposes).
> >> Sure. You can just plug in reduce if you don't mid the risk of losing
> >> control over the order:
> >> double stddev(const auto &first, const auto &last)
> >> {
> >> double sum = std::reduce(first, last);
> >> double sum_sq = std::reduce(
> >> first, last, [](double sqs, double x){ return sqs + x*x; });
> >> double N = std::distance(first, last);
> >> return sqrt((sum_sq - sum*sum/N)/(N-1));
> >> }
> >>
> > This won't work. You want the sum of squares. But if you evaluate left
> > to right, the first element is not squared. That's easily fixed by
> > passing in an init.
> Thanks. I removed the 0 init from the accumulate version without
> thinking. But...
> > However if you use "reduce" instead of "accumulate" then the problem is
> > much more serious. Rndom pairs might not have the first operand
> > squared.
> Ah yes. The binary operator must be associative and commutative.
>
> > You must use reduce_transform.
>
> (transform_reduce)
>
> So I read the documentation (as I clearly should have done before!) and,
> unless I've messed up again, an even neater solution presents itself:
> double stddev(const auto &first, const auto &last)
> {
> double sum = std::reduce(first, last);
> double sum_sq = std::transform_reduce(first, last, first);
> double N = std::distance(first, last);
> return sqrt((sum_sq - sum*sum/N)/(N-1));
> }
That's a huge readability improvement.
>
> because std::transform_reduce(first, last, first) is essentially a
> parallel inner product of the sequence with itself. It's possible,
> given the prevalence of inner products, that this might be special-cased
> and particularly fast on some hardware.
>
I'm in two minds about this. Whilst it's good to know that inner products
are particularly fast, ina sense special-casing is a cheat. The idea is that
the transorm reduce paradigm can be widely used to construct a wide
varienty of functions. I chose stddev because it lends itself to this
approach.
>
> BTW, can you do anything about your Usenet client messing up the
> formatting? It's particularly unhelpful with code, but "paragraph
> breaks" get lost too.
>
Unfortunately the machine I use is a work machine, and not my property.

On 06/03/2022 12:19, Bart wrote:
>
>
> I don't need to read any further. Just look at the difference between
> those two examples. I can easily follow the algorithm expressed as C,
> but I can'make out anything in the C++. Is there even an algorithm in
> that lot? All I can see is the last line that matches C's.

A simple C program without any bells & whistles!!!!

<====================================================================>
#include <math.h>
#include <stdio.h>
#include <time.h>
#include <stdlib.h>

float calculateSD(float data[]);
int main()
{ int i;
srand(time(NULL));
float data[10];

for (i = 0; i < 10; ++i)
{
data[i] = rand() % 99;
}
printf("The 10 numbers are: \n");
for (int i = 0; i < 10; i++)
{
printf("%4.0f", data[i]);
}
printf("\nStandard Deviation = %.6f", calculateSD(data));
return 0;
}

float calculateSD(float data[])
{ float sum = 0.0, mean, SD = 0.0;
int i;
for (i = 0; i < 10; ++i)
{
sum += data[i];
}
mean = sum / 10;
for (i = 0; i < 10; ++i)
{
SD += pow(data[i] - mean, 2);
}
return sqrt(SD / 10);
}

<====================================================================>

Am 06.03.2022 um 20:30 schrieb Manu Raju:
> On 06/03/2022 12:19, Bart wrote:
>>
>>
>> I don't need to read any further. Just look at the difference between
>> those two examples. I can easily follow the algorithm expressed as C,
>> but I can'make out anything in the C++. Is there even an algorithm in
>> that lot? All I can see is the last line that matches C's.
>
> A simple C program without any bells & whistles!!!!
>
> <====================================================================>
> #include <math.h>
> #include <stdio.h>
> #include <time.h>
> #include <stdlib.h>
>
> float calculateSD(float data[]);
> int main()
> {
>    int i;
>    srand(time(NULL));    float data[10];
>    for (i = 0; i < 10; ++i)
>    {
>        data[i] = rand() % 99;
>    }
>    printf("The 10 numbers are: \n");
>    for (int i = 0; i < 10; i++)
>    {
>        printf("%4.0f", data[i]);
>    }
>    printf("\nStandard Deviation = %.6f", calculateSD(data));
>    return 0;
> }
>
> float calculateSD(float data[])
> {
>    float sum = 0.0, mean, SD = 0.0;
>    int i;
>    for (i = 0; i < 10; ++i)
>    {
>        sum += data[i];
>    }
>    mean = sum / 10;
>    for (i = 0; i < 10; ++i)
>    {
>        SD += pow(data[i] - mean, 2);
>    }
>    return sqrt(SD / 10);
> }

Works only with float-arrays.

This works with any input-iterable data that is
convertible to double:

template<typename Iterator>
requires std::is_convertible_v<typename
std::iterator_traits<Iterator>::value_type, double> &&
std::input_iterator<Iterator>
double stddev( Iterator begin, Iterator end )
{ using namespace std;
using T = typename iterator_traits<Iterator>::value_type;
ptrdiff_t n;
double mean = 0.0;
if constexpr( random_access_iterator<Iterator> )
n = std::distance( begin, end ),
for_each( begin, end, [&]( T &value ) { mean += (double)value; } );
else
n = 0,
for_each( begin, end, [&]( T &value ) { ++n; mean +=
(double)value; } );
mean /= n;
double sqSum = 0.0;
for_each( begin, end, [&]( T &value ) { double diff = (double)value
- mean; sqSum += diff * diff; } );
return sqrt( sqSum / n );
}

Magnitudes more flexible.

>        SD += pow(data[i] - mean, 2);

#include <cmath>

double f( double d )
{ return pow( d, 2 );
}

MSVC:

movsd xmm1, QWORD PTR __real@4000000000000000
jmp pow

David Brown <david.brown@hesbynett.no> writes:

> On 06/03/2022 18:03, Ben Bacarisse wrote:
>> David Brown <david.brown@hesbynett.no> writes:
>>
>>> On 06/03/2022 16:36, Ben Bacarisse wrote:
>>>> scott@slp53.sl.home (Scott Lurndal) writes:
>>>>
>>>>> Ben Bacarisse <ben.usenet@bsb.me.uk> writes:
>>>>>> Bart <bc@freeuk.com> writes:
>>>>>>
>>>>>>> That this is supposed to be generic is no excuse. Generic code should
>>>>>>> not need to look much different from non-generic code. With dynamic
>>>>>>> typing, it would be exactly the same source code**.
>>>>>>>
>>>>>>> According to the examples posted by David Brown, a generic version of
>>>>>>> your the C function, written in C++, ought to look like this:
>>>>>>>
>>>>>>> double stddev(auto *x, int N) {...}
>>>>>>
>>>>>> double stddev(const auto *x, int N)
>>>>>> {
>>>>>> double sum = 0, sum_sq = 0;
>>>>>> for (int i = 0; i < N; i++) {
>>>>>> sum += x[i];
>>>>>> sum_sq += x[i]*x[i];
>>>>>> }
>>>>>> return sqrt((sum_sq - sum*sum/N)/(N-1));
>>>>>> }
>>>>>>
>>>>>> But then you then probably loose the automatic parallelisation that the
>>>>>> standard library might give you.
>>>>>
>>>>> Unless you're using one of the ubiquitous automatically vectorizing
>>>>> compilers.
>>>>
>>>> Do they work well with type-generic code, given that they won't always
>>>> know how + and * are defined?
>>>
>>> They know how + and * are defined when compiling the code - it's
>>> compile-time polymorphism. Compilers can do exactly as much with such
>>> templated (or implicitly templated with "auto") code as they can with
>>> explicit types.
>>
>> But do they? If you include link-time with compile-time, sure the
>> compiler knows what * is, but it's usually too late to change what the
>> loop compiled to. Does the code generator write two forms and the
>> linker picks one or the other when certain conditions are met?
>>
>> Or maybe I've missed your point altogether.
>
> When you use "auto" in a function declaration, it is just like a
> template (indeed it /is/ a template, written with a shorter syntax).
> Code is not generated for it at all until the compiler knows the type
> being used, and the details of it. Of course, if the type is a class
> with operators + and *, which are declared at the point of template
> instantiation but defined externally, then the compiler won't be able to
> do anything special with them.

Your "of course" is exactly my point -- apparently what I was saying is
obvious! Your "[t]hey know how + and * are defined when compiling the
code" is either not always true, or, in the case of whole-program
optimisation, not current enough to be much use.

>>> But I think there's a bit of mixup here between "parallelisation" and
>>> "vectorisation" - they don't mean the same thing. Compilers can do
>>> automatic vectorisation on code like this, but automatic parallelisation
>>> is rarer.
>>
>> I would say that vectorisation is a form of parallelisation -- the
>> latter being the generic term for whole range of specialist
>> optimisations. I may well have used one or other of them incorrectly,
>> but I know they are different!
>
> I would use "vectorisation" to refer to SIMD operations in the same
> thread, and "parallelisation" to refer to splitting the calculation into
> parts and running them on separate threads on separate cores. There is
> an overlap, but I think the distinction is clear.

I would call that concurrency. It's another form of parallelism. My
usage was common (at least in the UK) some time ago, but maybe the kids
use the words differently nowadays and the meaning has drifted.

--
Ben.

Malcolm McLean <malcolm.arthur.mclean@gmail.com> writes:

> On Sunday, 6 March 2022 at 19:42:41 UTC, Ben Bacarisse wrote:
>> Malcolm McLean <malcolm.ar...@gmail.com> writes:
>>
>> > On Sunday, 6 March 2022 at 15:59:14 UTC, Ben Bacarisse wrote:
>> >> Malcolm McLean <malcolm.ar...@gmail.com> writes:
>> >>
>> >> > On Sunday, 6 March 2022 at 15:05:04 UTC, Ben Bacarisse wrote:
>> >> >> Bart <b...@freeuk.com> writes:
>> >> >>
>> >> >> > That this is supposed to be generic is no excuse. Generic code should
>> >> >> > not need to look much different from non-generic code. With dynamic
>> >> >> > typing, it would be exactly the same source code**.
>> >> >> >
>> >> >> > According to the examples posted by David Brown, a generic version of
>> >> >> > your the C function, written in C++, ought to look like this:
>> >> >> >
>> >> >> > double stddev(auto *x, int N) {...}
>> >> >> double stddev(const auto *x, int N)
>> >> >> {
>> >> >> double sum = 0, sum_sq = 0;
>> >> >> for (int i = 0; i < N; i++) {
>> >> >> sum += x[i];
>> >> >> sum_sq += x[i]*x[i];
>> >> >> }
>> >> >> return sqrt((sum_sq - sum*sum/N)/(N-1));
>> >> >> }
>> >> >> But then you then probably loose the automatic parallelisation that the
>> >> >> standard library might give you.
>> >> >> > x is a pointer to int, long long int, float or double. But I couldn't
>> >> >> > get this to compile. If code needs to look like your C++, then it is a
>> >> >> > failure IMO.
>> >> >> The trouble is that the C++ code shown was intended to do more than the
>> >> >> C, and not just in terms of handling generic types. Also, you can't
>> >> >> judge C++ by a single example from one person. For example, I don't
>> >> >> think
>> >> >>
>> >> >> double stddev(const auto &first, const auto &last)
>> >> >> {
>> >> >> double sum = std::accumulate(first, last);
>> >> >> double sum_sq = std::accumulate(
>> >> >> first, last, [](double sqs, double x){ return sqs + x*x; });
>> >> >> double N = std::distance(first, last);
>> >> >> return sqrt((sum_sq - sum*sum/N)/(N-1));
>> >> >> }
>> >> >>
>> >> >> is particularly unclear.
>> >> >>
>> >> > Yes, but that's not generic. If first and last iterator over employees, the "+'
>> >> > operator has to be defined to add two employees together. That's unlikely
>> >> > to make any sort of sense.
>> >> I think that's a good thing. Standard deviation is an arithmetic
>> >> property.
>> >> > In fact we'll want that standard deviation of the salaries. So we can pass
>> >> > in a trivial little lambda [](Employee &employee){return
>> >> > employee.salary;}
>> >> I think you should use the ranges::view types to pass iterators to that
>> >> data. Otherwise, every arithmetic operation has to acquire a little
>> >> function to do the conversion, and the most common cases need a trivial
>> >> function that does no more than the language does by default.
>> >> > If we've not another variable associated with the employees, like hours
>> >> > worked, we can get the standard deviation of that just by changing the
>> >> > lambda.
>> >> I think views can handle this just fine without complicating every
>> >> arithmetic operation. But this is a classic example of why you should
>> >> be posting where experts can reply -- I may be wrong.
>> >> > reduce is better than accumulate because accumulate requires left-to
>> >> > right evaluation of the operands. That's not necessary for a standard
>> >> > deviation, unless annoyingly, we need every bit of precision (floating
>> >> > point arithmetic isn't exactly associative, though it's near enough
>> >> > for almost all purposes).
>> >> Sure. You can just plug in reduce if you don't mid the risk of losing
>> >> control over the order:
>> >> double stddev(const auto &first, const auto &last)
>> >> {
>> >> double sum = std::reduce(first, last);
>> >> double sum_sq = std::reduce(
>> >> first, last, [](double sqs, double x){ return sqs + x*x; });
>> >> double N = std::distance(first, last);
>> >> return sqrt((sum_sq - sum*sum/N)/(N-1));
>> >> }
>> >>
>> > This won't work. You want the sum of squares. But if you evaluate left
>> > to right, the first element is not squared. That's easily fixed by
>> > passing in an init.
>> Thanks. I removed the 0 init from the accumulate version without
>> thinking. But...
>> > However if you use "reduce" instead of "accumulate" then the problem is
>> > much more serious. Rndom pairs might not have the first operand
>> > squared.
>> Ah yes. The binary operator must be associative and commutative.
>>
>> > You must use reduce_transform.
>>
>> (transform_reduce)
>>
>> So I read the documentation (as I clearly should have done before!) and,
>> unless I've messed up again, an even neater solution presents itself:
>> double stddev(const auto &first, const auto &last)
>> {
>> double sum = std::reduce(first, last);
>> double sum_sq = std::transform_reduce(first, last, first);
>> double N = std::distance(first, last);
>> return sqrt((sum_sq - sum*sum/N)/(N-1));
>> }
> That's a huge readability improvement.
>>
>> because std::transform_reduce(first, last, first) is essentially a
>> parallel inner product of the sequence with itself. It's possible,
>> given the prevalence of inner products, that this might be special-cased
>> and particularly fast on some hardware.
>>
> I'm in two minds about this. Whilst it's good to know that inner products
> are particularly fast, ina sense special-casing is a cheat. The idea is that
> the transorm reduce paradigm can be widely used to construct a wide
> varienty of functions. I chose stddev because it lends itself to this
> approach.

I would not write the code that way /because/ it might be special-cased.
I would write it that way because it's, to me, clearer than using the
long form with two functions.

Maybe your objection is that the fact that

std::transform_reduce(first, last, first)

means

std::transform_reduce(first, last, first, std::plus<>(), std::multiplies<>())

is rather less than obvious? Maybe. But once you know what the three
iterator version of reduce is doing, the defaults for the functions
become obvious.

--
Ben.

On 06/03/2022 19:30, Manu Raju wrote:
> On 06/03/2022 12:19, Bart wrote:
>>
>>
>> I don't need to read any further. Just look at the difference between
>> those two examples. I can easily follow the algorithm expressed as C,
>> but I can'make out anything in the C++. Is there even an algorithm in
>> that lot? All I can see is the last line that matches C's.
>
> A simple C program without any bells & whistles!!!!
>
> <====================================================================>
> #include <math.h>
> #include <stdio.h>
> #include <time.h>
> #include <stdlib.h>
>
> float calculateSD(float data[]);
> int main()
> {
>    int i;
>    srand(time(NULL));    float data[10];
>    for (i = 0; i < 10; ++i)
>    {
>        data[i] = rand() % 99;
>    }
>    printf("The 10 numbers are: \n");
>    for (int i = 0; i < 10; i++)
>    {
>        printf("%4.0f", data[i]);
>    }
>    printf("\nStandard Deviation = %.6f", calculateSD(data));
>    return 0;
> }
>
> float calculateSD(float data[])
> {
>    float sum = 0.0, mean, SD = 0.0;
>    int i;
>    for (i = 0; i < 10; ++i)
>    {
>        sum += data[i];
>    }
>    mean = sum / 10;
>    for (i = 0; i < 10; ++i)
>    {
>        SD += pow(data[i] - mean, 2);
>    }
>    return sqrt(SD / 10);
> }
>
> <====================================================================>

This is a little /too/ simple. The size of the data should be a
parameter to calculateSD(), otherwise it will only ever work for arrays
of length 10.

> return sqrt(SD / 10);

This should probably divide by 9 (or N-1 when done properly).

Also, it there is only one such routine, making it work on doubles is
better.