/// @file     power_regression_approximation.cpp
/// @author   Vo Hoang Anh, <spyofgame200@gmail.com>
/// @brief    This is a simple implementation of the power regression approximation
///           With given array (x[]) and (y[]) at same size (n)
///           It will find constant (a), (b) that provides
///           Approximation function f(x) = a * x[i] ^ b ~ y[i] 
///           
/// @license  Public domain.
 
#include <iostream>
#include <iomanip>
#include <cmath>
 
using namespace std;
 
/// Size limitation
const int LIM = 1e6 + 16;
 
/// Given array
int n;
long double x[LIM], y[LIM];
 
/// Constant need to find
long double a, b;
 
/// set it to (-1) if you only need the approximation function
/// Number of digits after decimal points of (a) and (b)
int w = 15;
 
/// First min(k, n) number to be comparing f(x[i]) vs (y[i])
int k = 0;
 
/// Funtion: a * x^b ~ y
long double f(long double x) { return a * pow(x, b); }
int main()
{
    /// Input array
    cin >> n;
    for (int i = 1; i <= n; ++i) cin >> x[i];
    for (int i = 1; i <= n; ++i) cin >> y[i];
 
    /// Detail information
    cin >> w >> k;
 
    /// Precalculating
    long double sigma_ln_x = 0;
    long double sigma_ln_y = 0;
    long double sigma_ln_xx = 0;
    long double sigma_ln_xy = 0;
    for (int i = 1; i <= n; ++i)
    {
        long double tx = log(x[i]);
        long double ty = log(y[i]);
        sigma_ln_x += tx;
        sigma_ln_y += ty;
        sigma_ln_xx += tx * tx;
        sigma_ln_xy += tx * ty;
    }
 
    /// Power approximation constant
    b = ((n * sigma_ln_xy) - (sigma_ln_x * sigma_ln_y))
      / ((n * sigma_ln_xx) - (sigma_ln_x * sigma_ln_x));
 
    /// Base approximation constant
    a = exp((sigma_ln_y - b * sigma_ln_x) / n);
 
    /// Output such function
    cout << setprecision(w) << "Approximation Function: f(x) =   " << a << " * x ^(" << b << ")" << endl;
    
    /// Incase you dont need error
    if (w == -1) return 0;

    /// Calculating error
    long double Sigma_error_value_down = 0;
    long double Sigma_error_value_up   = 0;
    long double Sigma_error_percent_down = 0;
    long double Sigma_error_percent_up   = 0;
    long double Sigma_error_value   = 0;
    long double Sigma_error_percent = 0;

    k = min(k, n);
    if (k) cout << "\nComparing f(x[i]) vs y[i]:\n";
    for (int i = 1; i <= n; ++i)
    {
        long double Error_value = y[i] - f(x[i]);
        long double Error_percent = Error_value / y[i] * 100;
        if (i <= k) /// Comparing f(x[i]) vs (y[i])
        {
            cout << "At: " << setw(6) << i << " | ";
            cout << setprecision(w) << fixed << "<" << f(x[i]) << " -> " << y[i] << ">   =   ";
            cout << setprecision(w) << fixed << Error_value << " ~ " << Error_percent << " %" << endl;
        } 
 
        if (Error_value < 0)   Sigma_error_value_down   += Error_value;
        if (Error_value > 0)   Sigma_error_value_up     += Error_value;
        Sigma_error_value   += fabs(Error_value);
 
        if (Error_percent < 0) Sigma_error_percent_down += Error_percent;
        if (Error_percent > 0) Sigma_error_percent_up   += Error_percent;
        Sigma_error_percent += fabs(Error_percent);
    }
 
    /// Average error
    long double Average_error_value_down   = Sigma_error_value_down   / n;
    long double Average_error_value_up     = Sigma_error_value_up     / n;
    long double Average_error_value        = Sigma_error_value        / n;
    long double Average_error_percent_down = Sigma_error_percent_down / n;
    long double Average_error_percent_up   = Sigma_error_percent_up   / n;
    long double Average_error_percent      = Sigma_error_percent      / n;
 
    /// Output error
    cout << "\nFunction Errors: " << setprecision(4) << Average_error_percent << "%\n";
    cout << "Sigma_error_value_down      " << " = " << setprecision(w) << fixed << Sigma_error_value_down     << endl;
    cout << "Sigma_error_value_up        " << " = " << setprecision(w) << fixed <<  Sigma_error_value_up       << endl;
    cout << "Sigma_error_value           " << " = " << setprecision(w) << fixed <<  Sigma_error_value          << endl;
    cout << "Average_error_value_down    " << " = " << setprecision(w) << fixed <<  Average_error_value_down   << endl;
    cout << "Average_error_value_up      " << " = " << setprecision(w) << fixed <<  Average_error_value_up     << endl;
    cout << "Average_error_value         " << " = " << setprecision(w) << fixed <<  Average_error_value        << endl;
    cout << "Sigma_error_percent_down    " << " = " << setprecision(w) << fixed <<  Sigma_error_percent_down   << "%" << endl;
    cout << "Sigma_error_percent_up      " << " = " << setprecision(w) << fixed <<  Sigma_error_percent_up     << "%" << endl;
    cout << "Sigma_error_percent         " << " = " << setprecision(w) << fixed <<  Sigma_error_percent        << "%" << endl;
    cout << "Average_error_percent_down  " << " = " << setprecision(w) << fixed <<  Average_error_percent_down << "%" << endl;
    cout << "Average_error_percent_up    " << " = " << setprecision(w) << fixed <<  Average_error_percent_up   << "%" << endl;
    cout << "Average_error_percent       " << " = " << setprecision(w) << fixed <<  Average_error_percent      << "%" << endl;
    return 0;
}