#include <cassert>
#include <cstddef>
#include <iostream>

/***************************************************
 * Simple linear regression
 *
 * Suppose you have data
 *
 *     (x1, y1), ..., (xn, yn)
 *
 * which, when plotted, looks like this:
 *
 * |                   x             x
 * |                             x
 * |
 * |  x         x           x
 * |
 * |       x
 * |
 * +------------------------------------
 *
 * You are asked to find the best-fitting line
 *
 *     y = ax + b
 *
 * through these points:
 *
 * |                   x             x
 * |                           **x******
 * |                  *********
 * |  x      ***x*****      x
 * |*********
 * |       x
 * |
 * +------------------------------------
 *
 * The least-squares values of a and b satisfy
 *
 *     a SS_x + b S_x = SP_xy
 *
 *     a S_x  + b n   = S_y
 *
 * where
 *
 *     S_x   is the sum of the x values,
 *
 *     SS_x  is the sum of the squares of the
 *           x values,
 *
 *     S_y   is the sum of the y values, and
 *
 *     SP_xy is the sum of the products of
 *           corresponding x and y values.
 *
 * To solve a system of linear equations
 *
 *     a A + b B = C
 *     a D + b E = F,
 *
 * calculate
 *
 *         CE - BF             AF - CD
 *     a = -------,        b = -------.
 *         AE - BD             AE - BD
 *
 * Implement a function that takes two arrays with
 * the same capacity: one containing the x values
 * and the other containing the corresponding
 * y values.
 *
 * Two parameters are passed by reference and are
 * assigned the slope and constant of the
 * best-fitting line.
 *
 * If all the x values are equal, there is no
 * unique best-fitting slope. In this case, set the
 * slope to 0 and the constant to the average of
 * the y values.
 ***************************************************/

// Function declarations
int main();

void least_squares(
    double const x[],
    double const y[],
    std::size_t const capacity,
    double &slope,
    double &constant
);

// Function definitions
int main() {
    std::size_t const N{ 7 };

    double x[N]{
        0.9, 2.1, 3.2, 4.8, 6.0, 7.1, 8.1
    };

    // To test the case when all x values are equal:
    //
    // double x[N]{
    //     0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9
    // };
    //
    // The result will be
    //
    //     y = 0x + 5.3

    double y[N]{
        4.2, 2.8, 4.0, 7.4, 4.9, 6.3, 7.5
    };

    double a{};
    double b{};

    least_squares( x, y, N, a, b );

    std::cout << "The best-fitting line is y = "
              << a << "x ";

    if ( b >= 0.0 ) {
        std::cout << "+ " << b;
    } else {
        std::cout << "- " << (-b);
    }

    std::cout << std::endl;

    return 0;
}

/***************************************************
 * least_squares
 *
 * Calculate the slope and constant of the
 * best-fitting line
 *
 *     y = slope*x + constant.
 *
 * There must be at least one data point.
 ***************************************************/
void least_squares(
    double const x[],
    double const y[],
    std::size_t const capacity,
    double &slope,
    double &constant
) {
    assert( capacity > 0 );

    slope    = 0.0;
    constant = 0.0;
}




























// Hint:
//
// First calculate
//
//     S_x, S_y, SS_x, and SP_xy.
//
// You could write four separate functions to
// calculate these quantities, although that is not
// necessary here.
//
// Then substitute these values into the formulas
// given above.
//
// Remember that mathematical formulas often number
// data from 1 through n, whereas in C++ the entries
// of an array of capacity n have indices
//
//     0, 1, ..., n - 1.
//
// There are still n entries.
//
// You must also consider what happens when all the
// x values are equal. In that case, the denominator
// in the formulas is zero.




































// Here is one solution:

/*
void least_squares(
    double const x[],
    double const y[],
    std::size_t const capacity,
    double &slope,
    double &constant
) {
    assert( capacity > 0 );

    double sum_x{ 0.0 };
    double sum_y{ 0.0 };
    double sum_squares_x{ 0.0 };
    double sum_products_xy{ 0.0 };

    bool all_equal_x{ true };

    for ( std::size_t k{ 0 }; k < capacity; ++k ) {
        if ( x[k] != x[0] ) {
            all_equal_x = false;
        }

        sum_x           += x[k];
        sum_y           += y[k];
        sum_squares_x   += x[k]*x[k];
        sum_products_xy += x[k]*y[k];
    }

    if ( all_equal_x ) {
        slope    = 0.0;
        constant = sum_y/capacity;

        return;
    }

    // There is no point in calculating this twice.
    double const denominator{
        sum_squares_x*capacity - sum_x*sum_x
    };

    slope = (
        sum_products_xy*capacity - sum_x*sum_y
    )/denominator;

    constant = (
        sum_squares_x*sum_y
            - sum_products_xy*sum_x
    )/denominator;
}
*/