Generator take_while based on difference between elements

Question

I'm trying to approximate Hughes functional version of Newton-Raphson square roots algorithm from his 1989 paper "Why Functional Programming Matters".

I appreciate any advice for alternative approaches: more is better. My current approach uses Niebler's range-v3. You'll see in the code snippet that I've created a generator to create the successive iterations and insert them in a stream. My problem is the termination condition. I need to detect when the difference between consecutive floats in the stream falls below a threshold:

#include <iostream>
#include <range/v3/all.hpp>

using namespace ranges::view;

int main() {

  auto sqrt_stream = generate([ x = 1.f, n = 3.f ]() mutable {
    auto prevx = x;
    x = (x + n / x) / 2;
    return prevx;
  });

  std::vector<float> sequence = sqrt_stream | take_while([](int x) { ??? });

  return 0;
}

I'm not sure how to get any purchase to compare consecutive elements, such as a simple forward difference. And even if I did, once I transform the stream into a forward difference, how would I recover the appropriate element from the generator?

You could have the `take_while` be a closure with a variable storing the previous value — Passer By, Jun 01 '17 at 07:04
Yes, absolutely, I could. That's going to be my fallback solution. I'm trying to avoid any explicit handling of state (this is for a presentation), but it isn't a huge deal if I cannot. As you can see, I already have a bit of explicit local state in the generator. The point of the talk is how the use of C++ forces an interesting balance between functional and imperative design. — Timtro, Jun 01 '17 at 14:08

score 1 · Answer 1 · answered Jun 01 '17 at 16:51

The simplest solution I can think of is to generate pairs of sequence elements and filter by successive differences (live):

  auto pair_stream = generate([ x = 1., n = 7. ]() mutable {
    auto prevx = x;
    x = (x + n / x) / 2;
    return std::make_pair(prevx, x);
  });

  auto rng = pair_stream | take_while([](auto&& pair) {
    return std::abs(pair.second - pair.first) > epsilon;
  }) | transform([](auto&& x) { return x.second; });

which necessarily omits the first estimate from the sequence.

Timtro · Accepted Answer · 2017-06-02T01:41:31.663

I got a suggestion from Ivan Čukić, author of the upcoming book Functional Programming in C++. He gave me the same suggestion as +Casey, to use pairs, although I fleshed out the implementation a little differently, and made it generic. I used zip to make pairs (tuples of 2). The full solution is below, for the future reference of anyone curious.

I'm not happy with the code. It certainly isn't pretty, and I feel like there's a lot of room for clear-headed thinking to simplify the it. But I've spent my time investment. Further improvements welcomed.

#include <iostream>
#include <list>
#include <range/v3/all.hpp>

using namespace ranges::view;

// create a lazy stream given a function f and a seed value a0.
// stream will contain [a0, f(a0), f(f(a0)), f(f(f(a0))),...]
template <typename T> //
auto repeat_stream(std::function<T(T)> f, T a0) {
  return generate([ f, a = a0 ]() mutable {
    auto prev = a;
    a = f(a);
    return prev;
  });
}

// Consumes a stream until cosnecutive values differ within eps, and then the
// latest value is returned.
template <typename E, typename T>
E within(E eps, T v) {
  std::list<std::tuple<E, E>> converging_list =
      zip(v, v | drop(1))
        | take_while([eps](std::tuple<E, E> x) {
            return std::abs(std::get<0>(x) - std::get<1>(x)) > eps;
        });
  return std::get<0>(converging_list.back());
}

// Newton-Raphson Square Roots
// A la Hughes 1989 "Why Functional Programming Matters"
// http://www.cs.utexas.edu/~shmat/courses/cs345/whyfp.pdf
float nr_sqrt(float x, float x0) {
  return within(
      1E-15,
      repeat_stream<float>([n = x](float a) { return (a + n / a) / 2; }, x0));
}

int main() {

  std::cout << nr_sqrt(9, 4) << std::endl;

  return 0;
}

Edit: It bothered me that Hughes functional program, as transcribed by me into C++, was perhaps the dumbest program ever written to compute a square root. I decided to take another hack at the problem.

This time, I abstracted out the process of finding a fixed-point of the function. I used a recursive lambda to create the iteration loop, and std::bind to partially apply N in the Newton-Raphson step:

#include <cmath>
#include <functional>
#include <iostream>

using namespace std::placeholders;

template <typename T>
T fixedpoint(std::function<T(T)> f, T start, T eps) {

  std::function<T(T, T)> iter = [&iter, eps, &f](T old, T nw) -> T {
    if (std::abs(old - nw) < eps)
      return nw;
    else
      return iter(nw, f(nw));
  };

  return iter(start, f(start));
}

auto nr_step(double N, double a) { return (a + N / a) / 2; }

int main() {

  std::cout << fixedpoint<double>(std::bind(nr_step, 3, _1), 1.2, 1E-10)
            << std::endl;

  return 0;
}

If anyone knows how to resolve std::function<T(T)> with automatic type deduction so that I don't have to specify fixedpoint<double> I would be happy to hear about it.

Generator take_while based on difference between elements

2 Answers2