Se ha calculado la recta de regresión lineal a partir de los datos de edad y fallecimientos por COVID, obteniendo una pendiente de 25.15 y una ordenada en el origen de -749.91.
A partir de esta recta, se estiman los valores de fallecimientos para cada edad sustituyendo en la ecuación 
𝑦 = 25.15
𝑥 − 749.91
y= 25.15x − 749.91
De este modo, se obtienen las predicciones del modelo para edades consideradas de 5 en 5 años.



La salida muestra los coeficientes ai, bi y ci de cada spline cuadrático calculados por el alumno y los compara con los coeficientes de referencia. Como todos coinciden, la construcción de los splines es correcta.
Después se evalúa el spline en distintos valores de la variable edad, comparando de nuevo el resultado obtenido con el valor esperado. Al coincidir todos los valores, se verifica también que la evaluación del spline se ha implementado correctamente.

Los splines de grado 1 son lineas rectas entre los puntos,
Los splines de grado 2 aseguran solo continuidad de la primera derivada, produciendo discontinuidades en la curvatura y por tanto, cambios geometricos abruptos en la forma de la curva
Los splines de grado 3 garantizan continuidad de la segunda derivada, lo que implica continuidad en la curvatura

/// INCLUSION DE LIBRERIAS NECESARIAS
#include <stdio.h>
#include "EjerciciosBasicos.h"

/// P11.1 FUNCIÓN QUE CALCULA LA MEDIA DE UN VECTOR
real mn_media(Array1D< real > &u){
    // En resumen, inicializamos una variable a 0.0 para ir sumando el contenido del Vector u y el resultado de esta suma lo dividimos entre el número de elementos del Vector, lo cual da la media.
    real answer = 0.0;
    for(int i = 0; i 
	< u.dim(); i++){
        answer += u[i];
    }
    answer = answer / u.dim();
    return answer;
}

/// P11.2 FUNCIÓN QUE CALCULA EL MAXIMO DE UN VECTOR
real mn_max(Array1D< real > &u){
    // Tomamos el primer valor dentro del Vector u como el más grande hasta que recorriendo dicho Vector desde la segunda posición encontremos uno mayor y lo guardamos en la variable answer para devolverlo una vez recorrido el Vector completo.
    real answer = u[0];
    for(int i = 1; i < u.dim(); i++){
        if(u[i] > answer){
            answer = u[i];
        }
    }
    return answer;
}

/// P11.3 FUNCIÓN QUE CALCULA EL MINIMO DE UN VECTOR
real mn_min(Array1D< real > &u){
    // Tomamos el primer valor dentro del Vector u como el más pequeño hasta que recorriendo dicho Vector desde la segunda posición encontremos uno más pequeño que este y lo devolvemos una vez hayamos recorrido el Vector completo.
    real answer = u[0];
    for(int i = 1; i < u.dim(); i++){
        if(u[i] < answer){
            answer = u[i];
        }
    }
    return answer;
}

/// P11.4 FUNCIÓN QUE ORDENA UN VECTOR DE MENOR A MAYOR
void mn_ordenar(Array1D< real > &u){
    // Creamos una variable auxiliar para hacer el intercambio de valores en el Vector de forma correcta
    // Comparamos si el i-esimo valor del Vector es mayor al i-esimo + 1, en caso afirmativo, intercambiamos estos dos numeros con la ayuda de la variable auxiliar.
    // En caso negativo, seguimos recorriendo hasta llegar al penultimo valor del Vector que se va a comparar con el ultimo, una vez realizada esa comparacion terminamos la comparación
    // ya que la siguiente iteración no es posible para la variable "j", ya que se sale fuera del tamaño del Vector.
    real auxiliar = 0.;
    for(int i = 0; i < u.dim(); i++){
        for(int j = i + 1; j < u.dim(); j++){
            if(u[i] > u[j]){
                auxiliar = u[i];
                u[i] = u[j];
                u[j] = auxiliar;
            }
        }
    }
}

/// P11.5 FUNCIÓN PARA MULTIPLICAR UNA MATRIZ POR UN VECTOR
Array1D< real > mn_multiplicacion_matriz_vector(Array2D< real > &A,Array1D< real > &u){
    // Inicializamos un Vector del mismo numero de columnas que la Matriz A, ya que multiplicar una Matriz por un Vector, nos da un Vector del mismo numero de columnas que la Matriz original.
    // Recorremos la matriz A y guardamos la suma de las multiplicaciones de cada fila por la única columna del Vector u

    Array1D<real> answer(A.dim1());
    for(int i = 0; i < A.dim1(); i++){
        for(int j = 0; j < A.dim2(); j++){
            answer[i] = answer[i] + A[i][j] * u[j];
        }
    }
    return answer;
}

/// P11.6 FUNCIÓN QUE DETERMINA SI UN NÚMERO ENTERO ES PRIMO
bool mn_es_primo(int i){
    if(i <= 0){
        return false;
    }
    if(i == 1){
        return true;
    }
    for(int k = 2; k < i; k++){
        if(i % k == 0){
            return false;
        }
    }
    return true;
}

/// P11.7 FUNCIÓN QUE CALCULA EL FACTORIAL DE UN NÚMERO NATURAL
real mn_factorial(int n){
    real answer = 1.0;
    if(n == 0 || n == 1) return 1.;
    for(int k = n; k > 0; k--){
        answer *= k;
    }
    return answer;
}


/// P11.8 FUNCIÓN QUE CALCULA UNA POTENCIA CON UN NÚMERO NATURAL
/// NO SE PUEDE USAR LA FUNCIÓN pow()
real mn_potencia(real x,int n){
    real answer = 1.0;
    if(n == 0) return 1.;
    for(int i = 0; i < n; i++){
        answer *= x;
    }
    return answer;
}

/// P11.9 FUNCIÓN QUE CALCULA EL DESARROLLO DE TAYLOR DE e^x
/// e^x = 1 + x + x^2/2! + ...... +x^n/n!
real mn_exp(real x,int n){
    real term = 1.0;
    real sum = 1.0;
    for(int i = 1; i <= n; i++){
        term = term * x / i;
        sum += term;
    }
    return sum;
}

/// P11.10 FUNCIÓN QUE CALCULA EL DESARROLLO DE TAYLOR DE cos(x)
///  cos(x) = 1 - x^2/2! + x^4/4! - x^6/6!+...... +- x^(2n)/(2n)!
real mn_cos(real x,int n){
    real term = 1.0;
    real sum = 1.0;
    for(int i = 1; i <= n; i++){
        term = term * (-x * x) / ((2*i - 1)*(2*i));
        sum += term;
    }
    return sum;
}

/// P11.11 FUNCIÓN QUE CALCULA EL DESARROLLO DE TAYLOR DE sin(x)
///  sin(x) = x - x^3/3! + x^5/5! - x^7/7!+...... +- x^(2n+1)/(2n+1)!
real mn_sin(real x,int n){
    real term = x;
    real sum = x;
    for(int i = 1; i <= n; i++){
        term = term * (-x * x) / ((2*i)*(2*i+1));
        sum += term;
    }
    return sum;
}

/// P11.12 FUNCIÓN QUE CALCULA EL DESARROLLO DE TAYLOR DE ln(x)
/// ln(x) = (x-1) - ((x-1)^2)/2 + ((x-1)^3)/3 - ((x-1)^4)/4+...... +- ((x-1)^n)/n
real mn_ln(real x,int n){
    real y = (x - 1);
    real term = y;
    real sum = y;
    for(int i = 1; i < n; i++){
        term = term * (-y) * i / (i+1);
        sum += term;
    }
    return sum;
}

/// P11.13 FUNCIÓN QUE CALCULA y^x DONDE y,x SON NÚMERO REALES
/// USAR LAS FUNCIONES IMPLEMENTADAS mn_exp() y mn_ln() TENIENDO EN CUENTA y^x=e^(x*ln(y))
real mn_pow(real y,real x,int n){
    if(y <= 0){
        return 0.;
    }
    real a = mn_ln(y, n);
    real b = x * a;
    real c = mn_exp(b, n);
    return c;
}

/// P11.13 FUNCIÓN QUE CALCULA EL LIMITE DE LA SECUENCIA  yn=(1.+1./n).^n CUANDO n TIENDE A
/// INFINITO EL ALGORITMO PARA CUANDO LA DIFERENCIA EN VALOR ABSOLUTO DE LA DIFERENCIA
/// ENTRE UN VALOR DE LA SECUENCIA Y EL ANTERIOR ES INFERIOR AL PARAMETRO tolerancia
/// EL LIMITE DE LA SECUENCIA ES EL NUMERO e=2.71828182846
/// IMPORTANTE : PARA QUE LAS CONSTANTES LAS TRATE COMO NÚMEROS REALES HAY QUE AÑADIR UN .,
/// ES DECIR, POR EJEMPLO  1. (EN LUGAR DE 1). SI HACEMOS 1/2 EL RESULTADO ES CERO PORQUE HACE
/// LA DIVISIÓN EN PRECISIÓN ENTERA. SIN EMBARGO  1./2.=1./2=1/2.=0.5
real mn_limite1(real tolerancia){
    int n = 1;
    real y = 0.0;
    real y_previo = 0.0;

    do {
        y_previo = y;
        y = mn_potencia(1.0 + 1.0 / n, n);
        n++;
    } while (mn_abs(y - y_previo) > tolerancia);
    return y;
}

/// P11.14 FUNCIÓN QUE CALCULA EL LIMITE DE LA FUNCIÓN f(x)=sin(x)/x CUANDO x TIENDE HACIA 0.
/// EL PARAMETRO tolerancia SE UTILIZA PARA PARAR EL ALGORITMO CUANDO ESTAMOS CERCA DEL LÍMITE
/// EL VALOR DEL LÍMITE ES 1.
real mn_limite2(real tolerancia){
    real x = 1.0;
    real y = 0.0;

    do {
        y = mn_sin(x, 10) / x;
        x = x/2;
    } while (mn_abs(y - 1.0) > tolerancia);
    return y;
}

/// P11.15 FUNCIÓN QUE CALCULA EL LIMITE DE LA SECUENCIA  yn=X(n+1)/X(n) DONDE X(n) ES LA
/// SUCESIÓN DE FIBONACCI DEFINIDA COMO X(n+1)=X(n)+X(n-1) EMPEZANDO POR X(1)=X(2)=1
/// EL ALGORITMO PARA CUANDO LA DIFERENCIA EN VALOR ABSOLUTO
/// ENTRE UN VALOR DE LA SECUENCIA yn Y EL ANTERIOR ES INFERIOR AL PARAMETRO tolerancia
/// EL LIMITE DE LA SECUENCIA yn ES EL NÚMERO AÚREO IGUAL A (1+SQRT(5))/2 = 1.618033988....
real mn_limite3(real tolerancia){
    real x1 = 1.0;
    real x2 = 1.0;
    real x3;

    real y_previo = 0.0;
    real y = x2/x1;
    do {
        y_previo = y;
        x3 = x2 + x1;
        y = x3/x2;
        x1 = x2;
        x2 = x3;
    } while (mn_abs(y - y_previo) > tolerancia);
    return y;
}

/// P11.16 CÁLCULO DEL NÚMERO PI POR EL MÉTODO DE MONTECARLO. EL ÁREA DEL CÍRCULO DE RADIO
/// 1 ES PI. Y EL AREA DEL CUADRADO DE LADO 2 DONDE SE INSCRIBE EL CÍRCULO ES 4. POR TANTO
/// SI SE ELIGE UN PUNTO AL AZAR EN EL CUADRADO, LA PROBABILIDAD DE QUE CAIGA EN
/// EL CÍRCULO ES PI/4. EL MÉTODO DE MONTECARLO APROXIMA PI COGIENDO PUNTOS AL AZAR EN
/// EL CUADRADO [-1,1]x[-1,1] Y VIENDO QUE PROPORCIÓN CAE EN EL CÍRCULO.
/// NOTA : LA FUNCIÓN rand() DEVUELVE UN VALOR ENTERO ALEATORIO ENTRE 0 Y RAND_MAX
real calculo_pi_montecarlo(int Nintentos){
    int dentro = 0;
    for(int i = 0; i < Nintentos; i++){
        real x = 2.0 * rand() / RAND_MAX - 1.0;
        real y = 2.0 * rand() / RAND_MAX - 1.0;
        if(x*x + y*y <= 1.0){
            dentro++;
        }
    }
    return 4.0 * dentro/Nintentos;
}


##################################################################################################


#include "mn_calculo_de_ceros.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
using namespace std;

/// PARAMETROS DE LA DISTRIBUCION GAMMA
real alfa,beta_param,d;

/// DISTRIBUCIÓN GAMMA
real Gamma(real x){
  if(x<0) return 0;
  return d*pow(x,alfa-1.)*exp(-beta_param*x);
}

/// DERIVADA DE LA DISTRIBUCIÓN GAMMA
real Gammap(real x){
  if(x<0) return 0;
  return d*(alfa-1)*pow(x,alfa-2.)*exp(-beta_param*x)-d*beta_param*pow(x,alfa-1.)*exp(-beta_param*x);
}

/// CALCULO DE LOS PARÁMETROS DE LA DISTRIBUCIÓN GAMMA A PARTIR DE LA
/// MEDIA Y VARIANZA MUESTRAL
void calculo_parametros_Gamma(real media, real varianza){
  /// calculo de alfa y beta_param
  beta_param=media/varianza;
  alfa=media*beta_param;

  /// calculo de d a través de la integral (se verá en el tema 5)
  real h=0.001;
  real suma=0;
  for(real x=0;x<100;x+=h) suma+=h*pow(x,alfa-1.)*exp(-beta_param*x);
  d=1./suma;

  cout << "\nalfa = " << alfa << " beta_param = " << beta_param << " d = " << d << "\n";

}


/// FUNCION QUE IMPLEMENTA EL METODO DE LA BISECCION
real mn_biseccion (
real (*f)(const real), /// función a la cual se calcula un cero
real &a, real &b, /// intervalo inicial para buscar la raíz
const real TOL,  /// tolerancia para parar las iteraciones del algoritmo
int &Niter) /// número de iteraciones realizadas por el método
            ///       Si Niter=-1 la función ha terminado mal
{
    // Evaluar la funcion en los extremos del intervalo
    real fa = f(a);
    real fb = f(b);

    // Comprobar que existe cambio de signo en el intervalo
    if(fa * fb > 0){
        Niter = -1;
        return -1.;
    }

    // Inicializamos contador de iteraciones
    Niter = 0;

    // Repetir mientras el tamaño del intervalo sea mayor a la tolerancia pasada por parámetro.
    while(mn_distancia(a, b) >= TOL){
        // Calculamos el punto medio del intervalo
        real c = (a + b) * 0.5;
        // Evaluamos la funcion en el punto medio
        real fc = f(c);
        // Si encontramos exactamente la raiz, terminamos.
        if(fc == 0.) return(c);

        // Seleccionamos el subintervalo donde hay cambio de signo
        if(fa * fc < 0.){
            b = c;
            fb = fc;
        } else {
            a = c;
            fa = fc;
        }
        // Incrementamos el numero de iteraciones
        Niter++;
    }
    // Devolvemos el punto medio como aproximacion de la raiz
    return (a + b) * 0.5;
}


##################################################################################################

#include "mn_calculo_de_ceros.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
using namespace std;

/// PARAMETROS DE LA DISTRIBUCION GAMMA
real alfa,beta_param,d;

/// DISTRIBUCIÓN GAMMA
real Gamma(real x){
  if(x<0) return 0;
  return d*pow(x,alfa-1.)*exp(-beta_param*x);
}

/// DERIVADA DE LA DISTRIBUCIÓN GAMMA
real Gammap(real x){
  if(x<0) return 0;
  return d*(alfa-1)*pow(x,alfa-2.)*exp(-beta_param*x)-d*beta_param*pow(x,alfa-1.)*exp(-beta_param*x);
}

/// CALCULO DE LOS PARÁMETROS DE LA DISTRIBUCIÓN GAMMA A PARTIR DE LA
/// MEDIA Y VARIANZA MUESTRAL
void calculo_parametros_Gamma(real media, real varianza){
  /// calculo de alfa y beta_param
  beta_param=media/varianza;
  alfa=media*beta_param;

  /// calculo de d a través de la integral (se verá en el tema 5)
  real h=0.001;
  real suma=0;
  for(real x=0;x<100;x+=h) suma+=h*pow(x,alfa-1.)*exp(-beta_param*x);
  d=1./suma;

  cout << "\nalfa = " << alfa << " beta_param = " << beta_param << " d = " << d << "\n";

}


/// FUNCION QUE IMPLEMENTA EL METODO DE LA REGULA FALSI  Y DEVUELVE EL NÚMERO DE ITERACIONES
/// SI ALGO VA MAL DEVUELVE -1.  LA RAÍZ SE DEVUELVE COMO PARÁMETRO
int mn_regula_falsi (
real (*f)( real), 	/// función a la cual se calcula un cero
real &a, real &b, 	/// intervalo inicial para buscar la raíz
real &x, 			/// valor de salida de la raíz
real TOL,  			/// tolerancia para parar las iteraciones del algoritmo
int NiterMax) 		/// número máximo de iteraciones permitidas
{
	/// Evaluamos la función en los extremos del intervalo inicial
	real fa=f(a);
	real fb=f(b);
	/// Si f(a) y f(b) tienen el mismo signo, no se garantiza
	/// que exista una raíz en el intervalo -> error
	if(fa * fb>0) return -1;
	
	/// Calculamos la primera aproximación usando la fórmula
	/// de la regula falsi (intersección de la recta secante con el eje x)
	x=a-(b-a)*fa/(fb-fa);
	
	/// Evaluamos la función en la aproximación
	real fx=f(x);
	
	/// Bucle principal del método
	for(int k=0;k<NiterMax;k++){
		/// Si encontramos una raíz exacta terminamos
		if(fx==0.) return k;
		
		/// Determinamos en qué subintervalo está la raíz
		/// comprobando el cambio de signo
		if(fa*fx<0){
			b =x;
			fb=fx;
		}
		else{
		  a=x;
		  fa=fx;
		}
		/// Calculamos la siguiente aproximación mediante regula falsi
		real xnew=a-(b-a)*fa/(fb-fa);
		
		/// Si la diferencia entre dos aproximaciones consecutivas
		/// es menor que la tolerancia, consideramos que ha convergido
		if(mn_distancia(x,xnew)<=TOL){
			x=xnew;
			return k;
		}
		/// Actualizamos la aproximación actual
		x=xnew;
		/// Evaluamos la función en el nuevo punto
		fx=f(x);
	}
	/// Si se alcanza el número máximo de iteraciones
	/// sin cumplir el criterio de convergencia
	return -1;
}



####################################################################################

int mn_newton_raphson (
real (*f)( real), /// función sobre la que se calcula el cero
real &x0,         /// aproximación inicial de la raíz (se irá actualizando)
int NiterMax,     /// número máximo de iteraciones permitidas
real TOL)         /// tolerancia para detener el algoritmo
{
	/// Evaluamos la función en el punto inicial
	real f_x0 = f(x0);

	/// Bucle principal del método
	for(int i = 0; i < NiterMax; i++){

		/// Si encontramos una raíz exacta terminamos
		if(f_x0 == 0.){
			return(i);
		}

		/// Calculamos la derivada de la función en x0
		/// (usando una función auxiliar que calcula la derivada numérica)
		real derivada = mn_derivada1(f, x0);

		/// Si la derivada es cero, el método no puede continuar
		/// porque habría división por cero
		if(derivada == 0.) return(-1);

		/// Fórmula del método de Newton-Raphson
		/// x_{n+1} = x_n - f(x_n)/f'(x_n)
		real x1 = x0 - f_x0 / derivada;

		/// Si la diferencia entre dos aproximaciones consecutivas
		/// es menor que la tolerancia, consideramos que ha convergido
		if(mn_distancia(x1, x0) < TOL){
			x0 = x1;      /// actualizamos la raíz final
			return(i + 1); /// devolvemos el número de iteraciones realizadas
		}

		/// Actualizamos la aproximación para la siguiente iteración
		x0 = x1;
        fx0 = f(x0);
	}

	/// Si se alcanza el número máximo de iteraciones sin converger
	return(-1);
}

####################################################################################

/// FUNCION QUE IMPLEMENTA EL METODO DE LA SECANTE
/// LA FUNCIÓN DEVUELVE EL NÚMERO DE ITERACIONES REALIZADAS SI TERMINA BIEN
/// Y DEVUELVE -1 EN CASO CONTRARIO
int mn_secante (
real (*f)( real),   /// función sobre la que se calcula el cero
real &x0,           /// primera aproximación de la raíz (se irá actualizando)
real &x1,           /// segunda aproximación de la raíz (también se actualiza)
int NiterMax,       /// número máximo de iteraciones permitidas
real TOL)           /// tolerancia para detener el algoritmo
{
    /// Evaluamos la función en las dos aproximaciones iniciales
    real f0 = f(x0);
    real f1 = f(x1);

    /// Bucle principal del método
    for(int iter = 0; iter < NiterMax; ++iter)
    {
        /// Si encontramos una raíz exacta terminamos
        if(f1 == 0.0) return iter;

        /// Calculamos la distancia entre las dos aproximaciones
        /// esto representa el denominador de la pendiente
        real deltaX = x1 - x0;

        /// Si los dos puntos coinciden no podemos calcular la secante
        if(deltaX == 0.0) return -1;

        /// Calculamos la pendiente de la recta secante
        /// que pasa por (x0,f(x0)) y (x1,f(x1))
        real pendiente = (f1 - f0) / deltaX;

        /// Si la pendiente es cero no podemos continuar
        /// porque habría división por cero en la fórmula
        if(pendiente == 0.0) return -1;

        /// Calculamos la nueva aproximación usando la fórmula
        /// del método de la secante
        real x2 = x1 - f1 / pendiente;

        /// Si la diferencia entre aproximaciones consecutivas
        /// es menor que la tolerancia, el método ha convergido
        if(mn_distancia(x2, x1) < TOL)
        {
            x1 = x2;  /// guardamos la aproximación final
            return iter + 1; /// devolvemos el número de iteraciones realizadas
        }

        /// Actualizamos los valores para la siguiente iteración
        /// desplazamos los puntos de la secante
        x0 = x1;
        x1 = x2;

        /// Actualizamos también los valores de la función
        f0 = f1;
        f1 = f(x1);
    }

    /// Si se alcanza el número máximo de iteraciones
    /// sin cumplir el criterio de convergencia
    return -1;
}


####################################################################################


/// FUNCION QUE IMPLEMENTA EL METODO DE MULLER
/// LA FUNCIÓN DEVUELVE EL NÚMERO DE ITERACIONES SI TERMINA BIEN
/// Y DEVUELVE -1 EN CASO CONTRARIO
int mn_muller (
real (*f)( real), /// función sobre la que se calcula el cero
real &x0,         /// aproximación inicial de la raíz (se va actualizando)
int NiterMax,     /// número máximo de iteraciones permitidas
real TOL)         /// tolerancia para detener el algoritmo
{
  /// Evaluamos la función en la aproximación inicial
  real f_x0 = f(x0);

  /// Bucle principal del método
  for(int iter = 0; iter < NiterMax; ++iter)
  {
      /// Si encontramos una raíz exacta terminamos
      if(f_x0 == 0.0) return iter;

      /// Variables donde se guardarán las soluciones de la ecuación cuadrática
      real z1, z2;

      /// Calculamos derivadas en el punto actual x0
      real fp  = mn_derivada1(f, x0); // primera derivada f'(x0)
      real fpp = mn_derivada2(f, x0); // segunda derivada f''(x0)

      /// Aproximamos la función mediante un polinomio cuadrático:
      /// f(x0 + z) ≈ f(x0) + f'(x0)z + (f''(x0)/2)z²
      ///
      /// Para encontrar la raíz aproximada imponemos:
      /// f(x0) + f'(x0)z + (f''(x0)/2)z² = 0
      ///
      /// Esto es una ecuación cuadrática en z

      int Nroots = mn_ceros_pol_grado_2(
                        fpp / 2.0,   // coeficiente de z²
                        fp,          // coeficiente de z
                        f_x0,        // término independiente
                        z1, z2);     // soluciones de la cuadrática

      /// Si no hay raíces reales el método no puede continuar
      if(Nroots == 0) return -1;

      /// Elegimos una de las soluciones (z1) como desplazamiento
      /// desde el punto actual
      real x1 = x0 + z1;

      /// Comprobamos si la nueva aproximación está suficientemente cerca
      /// de la anterior (criterio de convergencia)
      if(mn_distancia(x1, x0) < TOL)
      {
          x0 = x1;       /// guardamos la raíz aproximada
          return iter + 1; /// devolvemos el número de iteraciones
      }

      /// Actualizamos el punto actual para la siguiente iteración
      x0 = x1;

      /// Evaluamos la función en el nuevo punto
      f_x0 = f(x0);
  }

  /// Si no se alcanza convergencia en NiterMax iteraciones
  return -1;
}


####################################################################################

void mn_evaluar_polinomio_horner(
Array1D< real > &a 	/** coeficientes del polinomio */,
real x 				/** valor donde se evalúa el polinomio */,
real &Px 			/** valor del polinomio P(x) */,
real &PPx 			/** valor de la derivada P'(x) */)
{
	/// Inicializamos Px con el coeficiente de mayor grado del polinomio
	/// (a_n si el polinomio es a0 + a1 x + ... + an x^n)
	Px = a[a.dim()-1];

	/// Inicializamos PPx para comenzar el cálculo de la derivada
	PPx = a[a.dim()-1];

	/// Recorremos los coeficientes desde grado n-1 hasta grado 1
	/// aplicando el esquema de Horner
	for(int k = a.dim()-2; k > 0; k--) {
		/// Evaluación del polinomio mediante Horner:
		/// Px = (...((an)x + an-1)x + ... )
		Px = Px * x + a[k];

		/// Evaluación simultánea de la derivada usando
		/// la extensión del método de Horner
		PPx = PPx * x + Px;
	}

	/// Añadimos el término independiente a0 para completar P(x)
	Px = Px * x + a[0];
}

Array1D< real > mn_calcular_derivada_polinomio(
Array1D< real > &a /** coeficientes del polinomio */)
{
	/// Si el polinomio es constante (solo tiene un coeficiente),
	/// su derivada es 0. En este caso devolvemos un array vacío.
	if(a.dim() == 1) 
	return Array1D< real >();

	/// Creamos un nuevo array b para guardar los coeficientes
	/// del polinomio derivado. Su tamaño es uno menos que el original.
	Array1D< real > b(a.dim() - 1);

	/// Calculamos los coeficientes de la derivada.
	/// Si el polinomio es:
	/// P(x) = a0 + a1 x + a2 x^2 + ... + an x^n
	///
	/// entonces su derivada es:
	/// P'(x) = a1 + 2 a2 x + 3 a3 x^2 + ... + n an x^(n-1)
	///
	/// Por lo tanto:
	/// b[k] = (k+1) * a[k+1]
	for(int k = 0; k < b.dim(); k++) {
		b[k] = (k + 1.) * a[k + 1];
	}

	/// Devolvemos el polinomio derivado
	return b;
}


/// ****************************************************
///  CALCULO POLINOMIO DERIVADA
///  DEVUELVE UN VECTOR CON EL POLINOMIO DERIVADA
/// ****************************************************
Array1D< real > calcular_derivada_polinomio(
Array1D< real > &a /** coeficientes del polinomio */){

	/// Si el polinomio es constante (solo un coeficiente),
	/// su derivada es 0 → devolvemos un vector vacío
	if(a.dim()==1) return Array1D< real >();

	/// Creamos el vector que contendrá los coeficientes
	/// del polinomio derivado (tiene un grado menos)
	Array1D< real > b(a.dim()-1);

	/// Calculamos los coeficientes de la derivada
	/// Si P(x) = a0 + a1 x + a2 x² + ... + an x^n
	/// entonces P'(x) = a1 + 2a2 x + 3a3 x² + ... + n an x^(n-1)
	for(int k=0;k<b.dim();k++){
		b[k]=(k+1.)*a[k+1];
	}

	/// Devolvemos el polinomio derivado
	return b;
}



/// *************************************************************************************
///  CALCULO CEROS DE UN POLINOMIO A PARTIR DE LOS CEROS DE SU DERIVADA
///  DEVUELVE UN VECTOR CON LOS CEROS DEL POLINOMIO.
///  UTILIZA LA FUNCIÓN calculo_ceros_en_intervalo() PARA CALCULAR EL CERO EN CADA INTERVALO
/// *************************************************************************************
Array1D< real > ceros_polinomio_desde_ceros_derivada(
Array1D< real > &a /** coeficientes polinomio */,
Array1D< real > &d /** ceros polinomio derivada */,
real TOL){

	/// Si el coeficiente de mayor grado es 0, el polinomio no tiene
	/// el grado esperado → devolvemos directamente los ceros conocidos
	if(a[a.dim()-1]==0.) return d;

	/// Calculamos el coeficiente de mayor valor absoluto
	real max_c=mn_abs(a[0]);
	for(int k=1;k<a.dim()-1;k++){
		if(mn_abs(a[k])>max_c) max_c=mn_abs(a[k]);
	}

	/// Calculamos una cota A tal que todas las raíces reales
	/// del polinomio están dentro del intervalo [-A,A]
	real A=1+max_c/mn_abs(a[a.dim()-1]);

	/// Vector donde guardaremos los ceros encontrados
	Array1D< real > ceros;

	/// Construimos el vector de puntos que delimitan los intervalos
	/// empezando por -A
	Array1D< real > I(1,-A);

	/// Añadimos los ceros de la derivada
	for(int k=0;k<d.dim();k++) I.push_back(d[k]);

	/// Añadimos el límite derecho A
	I.push_back(A);

	/// Recorremos cada intervalo consecutivo
	for(int k=0;k<I.dim()-1;k++){

		/// Si hay cambio de signo del polinomio en el intervalo
		/// entonces existe al menos una raíz en él
		if(evaluar_polinomio(a,I[k])*evaluar_polinomio(a,I[k+1])<=0.){

		  /// Calculamos la raíz dentro del intervalo usando
		  /// un método numérico (por ejemplo bisección)
		  ceros.push_back(calculo_cero_en_intervalo(a,I[k],I[k+1],TOL));
		}
	}

	/// Devolvemos las raíces encontradas
	return ceros;

}



/// ****************************************************
///  CALCULO RAÍCES REALES DE UN POLINOMIO
///  DEVUELVE UN VECTOR CON LAS RAÍCES ENCONTRADAS
/// ****************************************************
Array1D< real > ceros_polinomio(
Array1D< real > &a /** coeficientes polinomio */,
real TOL /** Tolerancia para el método de cálculo de ceros en un intervalo */){

	/// Vector donde se almacenarán las raíces encontradas
	Array1D< real > ceros;

	/// Recorremos desde el grado n-1 hasta 0
	/// calculando derivadas sucesivas del polinomio
	for(int k=a.dim()-2;k>=0;k--){

		/// Copiamos el polinomio original
		Array1D< real > d=a.copy();

		/// Calculamos la derivada k veces
		/// (esto genera derivadas de orden creciente)
		for(int n=0;n<k;n++) d=calcular_derivada_polinomio(d);

		/// Calculamos los ceros del polinomio actual
		/// usando los ceros de la derivada anterior
		ceros=ceros_polinomio_desde_ceros_derivada(d,ceros,TOL);
	}

	/// Devolvemos las raíces reales encontradas
	return ceros;

}



#include "mn_aritmeticas.h"
#include "mn_interpolacion.h"

/// INTERPOLACION VECINO
real mn_interpolacion_vecino(
Array1D<real> &x, 	/// puntos de interpolación
Array1D<real> &f, 	/// valores de función en los puntos de interpolación
real x0)  			/// punto donde se evalua la función interpolada
{
    // Recorremos los puntos desde el final hasta el principio con la posibilidad de comparar pares consecutivos [n-2]
    for(int i = f.dim() - 2; i >= 0; i--){
        // Si x0 esta a la derecha de x[i], significa que x0 esta entre x[i] y x[i+1]
        if( x[i] < x0 ){
            // Cual esta mas cerca, x[i] o x[i+1] ?? En cualquier caso, devolvemos su valor de la funcion
            if( (x0 - x[i]) <= mn_abs(x0 - x[i+1]) ) return f[i];
            else return f[i+1];
        }
    }
    // Si el bucle no encuentra ningun x[i] < x0, significa que x0 <= x[0]
    return f[0];
}


#include "mn_aritmeticas.h"
#include "mn_interpolacion.h"

/// INTERPOLACION LINEAL
real mn_interpolacion_lineal(
Array1D<real> &x,   /// puntos de interpolación
Array1D<real> &f,   /// valores de función en los puntos de interpolación
real x0)            /// punto donde se evalua la función interpolada
{
    // Recorremos los puntos desde el final - 1 para poder comparar los pares consecutivos x[i] y x[i+1]
    for(int i = f.dim()-2; i>=0; i--){
        // Si x[i] < x0 significa que x0 se encuentra entre x[i] y x[i+1]
        if( x[i] < x0 ) {
            // Formula de interpolacion lineal: valor de f[i] + pendiente(x0-x[i]) * distancia
            return f[i] + (x0-x[i]) * (f[i+1] - f[i]) / (x[i+1] - x[i]);
        }
    }
    return f[0]+(x0-x[0])*(f[1]-f[0])/(x[1]-x[0]);
}



#include "mn_aritmeticas.h"
#include "mn_interpolacion.h"

/**
 LA INTERPOLACION POR SPLINES DE GRADO 2 SEGÚN SE EXPLICÓ EN CLASE.
 EN ESTE CASO SE SUPONE CONOCIDO c[0] Y SE CALCULA EL RESTO DE COEFICIENTES
 LA FUNCION DEVUELVE -1 SI ENCUENTRA ALGÚN PROBLEMA Y 0 EN CASO CONTRARIO
*/
int mn_calculo_splines_2(
    Array1D< real > &x,         /// vector con los puntos de interpolación
    Array1D< real > &f,         /// vector con los valores de la función en esos puntos
    real &c0,                   /// valor inicial de c[0] (condición inicial del spline)
    Array1D< real > &a,         /// salida: coeficientes a[i] del spline
    Array1D< real > &b,         /// salida: coeficientes b[i] del spline
    Array1D< real > &c)         /// salida: coeficientes c[i] del spline
{
  /// Comprobamos que los vectores x y f tengan la misma dimensión
  /// y que haya al menos dos puntos para interpolar
  if(f.dim()!=x.dim() || f.dim()<2 ) return -1;

  /// Verificamos que los puntos x estén estrictamente ordenados
  for(int k=1;k<x.dim();k++){
      if(x[k]<=x[k-1]) return -1;
  }

  /// Inicializamos los vectores de coeficientes del spline
  /// Habrá un polinomio cuadrático por cada intervalo
  a=Array1D< real >(x.dim()-1);
  b=Array1D< real >(x.dim()-1);
  c=Array1D< real >(x.dim()-1);

  /// Primer coeficiente del primer intervalo
  a[0] = f[0];

  /// Calculamos b[0] usando la condición de interpolación en el segundo punto
  b[0] = (f[1]-f[0]-c0*(x[1]-x[0])*(x[1]-x[0]))/(x[1]-x[0]);

  /// Asignamos el valor inicial de c[0]
  c[0] = c0;

  /// Calculamos los coeficientes de los siguientes intervalos
  for(int i=1;i<a.dim();i++){

      /// El coeficiente a[i] es simplemente el valor de la función en x[i]
      a[i]=f[i];

      /// Continuidad de la derivada entre splines consecutivos
      b[i]=b[i-1]+2*c[i-1]*(x[i]-x[i-1]);

      /// Calculamos c[i] imponiendo que el spline pase por el punto siguiente
      c[i]=(f[i+1]-a[i]-b[i]*(x[i+1]-x[i]))/((x[i+1]-x[i])*(x[i+1]-x[i]));
  }

  /// Todo ha ido bien
  return 0;
}


/**
  EVALUACIÓN SPLINE DE GRADO DOS. SE DEVUELVE EL VALOR DE LA
  EVALUACIÓN EN x0
*/
real mn_evaluar_splines_2(
Array1D< real > &x,
Array1D< real > &a,
Array1D< real > &b,
Array1D< real > &c,
real x0 ){

    /// Buscamos el intervalo al que pertenece x0
    /// recorriendo los splines desde el final
    for(int i = a.dim()-1 ;i>=0;i--){

        /// Si x0 está a la derecha de x[i], usamos ese spline
        if(x[i] < x0)

            /// Evaluamos el polinomio cuadrático del spline
            /// S_i(x) = a[i] + b[i](x-x[i]) + c[i](x-x[i])²
            return a[i]+(b[i]+c[i]*(x0-x[i]))*(x0-x[i]);
    }

    /// Si x0 está antes del primer nodo usamos el primer spline
    return a[0]+(b[0]+c[0]*(x0-x[0]))*(x0-x[0]);
}


int mn_calculo_splines_3(
    const Array1D< real > &x, /// vector con los puntos de interpolación
    const Array1D< real > &f, /// vector con los valores de la función
    const real c0,            /// primer valor del vector c[] asignado como parámetro
    const real cN,            /// último valor del vector c[] asignado como parámetro
    Array1D< real > &a,       /// vector de salida con los coeficientes a[i] del spline de grado 3
    Array1D< real > &b,       /// vector de salida con los coeficientes b[i] del spline de grado 3
    Array1D< real > &c,       /// vector de salida con los coeficientes c[i] del spline de grado 3
    Array1D< real > &d)       /// vector de salida con los coeficientes d[i] del spline de grado 3
{
    /// Vector de longitudes de cada subintervalo:
    /// h[k] = x[k+1] - x[k]
    Array1D< real > h(x.dim() - 1);

    /// Vectores del sistema tridiagonal:
    /// L = diagonal inferior
    /// M = diagonal principal
    /// U = diagonal superior
    /// B = término independiente
    Array1D< real > M(x.dim(), 1.);
    Array1D< real > B(x.dim());
    Array1D< real > U(x.dim() - 1, 0.);
    Array1D< real > L(x.dim() - 1, 0.);

    /// Calculamos las longitudes de los intervalos
    for(int k = 0; k < h.dim(); k++)
        h[k] = x[k + 1] - x[k];

    /// Condiciones en los extremos para los coeficientes c[0] y c[n]
    B[0] = c0;
    B[B.dim() - 1] = cN;

    /// Construimos el sistema tridiagonal para calcular los coeficientes c[k]
    for(int k = 1; k < M.dim() - 1; k++) {
        L[k - 1] = h[k - 1];
        M[k] = 2 * (h[k - 1] + h[k]);
        U[k] = h[k];

        /// Término independiente asociado a la continuidad de derivadas
        B[k] = 3. * (f[k + 1] - f[k]) / h[k]
             - 3. * (f[k] - f[k - 1]) / h[k - 1];
    }

    /// Resolvemos el sistema para obtener los coeficientes c
    c = solucion_sistema(L, M, U, B);

    /// Reservamos memoria para los coeficientes a, b y d
    /// Hay un spline cúbico por cada intervalo
    a = Array1D< real >(x.dim() - 1);
    b = Array1D< real >(x.dim() - 1);
    d = Array1D< real >(x.dim() - 1);

    /// Calculamos los coeficientes de cada polinomio cúbico
    for(int k = 0; k < a.dim(); k++) {
        /// Término independiente: valor de la función en el nodo izquierdo
        a[k] = f[k];

        /// Coeficiente del término cúbico
        d[k] = (c[k + 1] - c[k]) / (3 * h[k]);

        /// Coeficiente del término lineal
        b[k] = (f[k + 1] - f[k]) / h[k]
             - h[k] * (2 * c[k] + c[k + 1]) / 3.;
    }

    return 0;
}

/**
  EVALUACIÓN SPLINE DE GRADO 3. SE DEVUELVE EL VALOR DE LA
  EVALUACIÓN EN x0
*/
real mn_evaluar_splines_3(
    const Array1D< real > &x,
    Array1D< real > &a,
    Array1D< real > &b,
    Array1D< real > &c,
    Array1D< real > &d,
    real x0)
{
    /// Recorremos los intervalos desde el final para encontrar
    /// el spline al que pertenece x0
    for(int i = a.dim() - 1; i >= 0; i--) {
        if(x[i] < x0)
            /// Evaluamos el polinomio cúbico asociado al intervalo i:
            /// S_i(x) = a[i] + b[i](x-x[i]) + c[i](x-x[i])^2 + d[i](x-x[i])^3
            return ((d[i] * (x0 - x[i]) + c[i]) * (x0 - x[i]) + b[i]) * (x0 - x[i]) + a[i];
    }

    /// Si x0 está antes del primer nodo, usamos el primer spline
    return ((d[0] * (x0 - x[0]) + c[0]) * (x0 - x[0]) + b[0]) * (x0 - x[0]) + a[0];
}


#include "mn_aritmeticas.h"
#include "mn_interpolacion.h"

///* FUNCION QUE CALCULA LA RECTA DE REGRESION LINEAL y=ax+b PARA APROXIMAR UNA
/// NUBE DE PUNTOS (x[k],y[k]). LA FUNCION DEVUELVE 0 SI TERMINA BIEN Y -1 SI
/// TERMINA MAL
/// */
int mn_regresion_lineal(
    Array1D< real > &x,     /// VECTOR DE COORDENADAS x DE LA NUBE DE PUNTOS
    Array1D< real > &y,     /// VECTOR DE COORDENADAS y DE LA NUBE DE PUNTOS
    real &a,                /// COMPONENTE a DE LA RECTA DE REGRESIÓN. (PARÁMETRO DE SALIDA)
    real &b)                /// COMPONENTE b DE LA RECTA DE REGRESIÓN. (PARÁMETRO DE SALIDA)
{
    /// Número de puntos de la nube
    int N = x.dim();

    /// Comprobamos que ambos vectores tengan la misma dimensión
    /// y que haya al menos dos puntos para ajustar una recta
    if(y.dim() != N || N < 2) return(-1);

    /// Inicializamos las sumas necesarias para la fórmula
    /// de mínimos cuadrados
    real suma_x = 0, suma_y = 0, suma_xy = 0, suma_x2 = 0;

    /// Recorremos todos los puntos acumulando:
    /// suma_x  = Σ x[k]
    /// suma_y  = Σ y[k]
    /// suma_xy = Σ x[k]·y[k]
    /// suma_x2 = Σ x[k]^2
    for(int k = 0; k < N; k++) {
        suma_x += x[k];
        suma_y += y[k];
        suma_xy += x[k] * y[k];
        suma_x2 += x[k] * x[k];
    }

    /// Calculamos el denominador común de las fórmulas
    /// de la pendiente a y la ordenada en el origen b
    real denominador = N * suma_x2 - suma_x * suma_x;

    /// Si el denominador es 0, no se puede calcular la recta
    /// (por ejemplo, si todos los valores x son iguales)
    if(denominador == 0.) return(-1);

    /// Calculamos la pendiente de la recta de regresión:
    /// a = (N Σxy - Σx Σy) / (N Σx² - (Σx)²)
    a = (N * suma_xy - suma_x * suma_y) / denominador;

    /// Calculamos la ordenada en el origen:
    /// b = (Σx² Σy - Σxy Σx) / (N Σx² - (Σx)²)
    b = (suma_x2 * suma_y - suma_xy * suma_x) / denominador;

    /// Todo ha ido bien
    return(0);
}


#include "mn_aritmeticas.h"
#include "mn_interpolacion.h"

/// FUNCION QUE CALCULA EL POLINOMIO INTERPOLADOR POR LAS DIFERENCIAS DE NEWTON
/// DEVUELVE UN VECTOR CON LOS COEFICIENTES DEL POLINOMIO O UN VECTOR VACÍO SI ALGO VA MAL
Array1D< real > mn_construir_polinomio_interpolador(
    const Array1D< real > &X, /// VECTOR CON LOS VALORES DE LOS PUNTOS DE INTERPOLACION
    const Array1D< real > &F) /// VECTOR DE VALORES DE LA FUNCION EN LOS PUNTOS DE INTERPOLACION
{
    /// Si los vectores de nodos y valores no tienen la misma dimensión,
    /// no se puede construir el polinomio interpolador
    if(X.dim() != F.dim()) return(Array1D< real >());

    /// Vector donde guardaremos los coeficientes del polinomio de Newton
    Array1D< real > A(X.dim());

    /// Grado del polinomio interpolador
    int N = X.dim() - 1;

    /// Vector auxiliar para ir calculando las diferencias divididas
    Array1D< real > B(A.dim());

    /// Inicializamos B con los valores de la función en los nodos
    for(int k = 0; k <= N; k++) {
        B[k] = F[k];
    }

    /// El primer coeficiente del polinomio de Newton es F[0]
    A[0] = F[0];

    /// Calculamos sucesivamente las diferencias divididas
    for(int k = 1; k <= N; k++) {

        /// En cada paso se actualiza B con las diferencias divididas
        /// de orden k
        for(int l = 0; l <= (N - k); l++) {

            /// Si dos nodos coinciden, no se puede dividir
            /// y el polinomio no está bien definido
            if(X[k + l] == X[l]) {
                return(Array1D< real >());
            }

            /// Fórmula de las diferencias divididas:
            /// f[x_l, ..., x_{l+k}] =
            /// (f[x_{l+1}, ..., x_{l+k}] - f[x_l, ..., x_{l+k-1}]) / (x_{l+k} - x_l)
            B[l] = (B[l + 1] - B[l]) / (X[k + l] - X[l]);
        }

        /// El primer elemento de B es el coeficiente de orden k
        A[k] = B[0];
    }

    /// Devolvemos los coeficientes del polinomio interpolador de Newton
    return(A);
}

/// FUNCION QUE EVALUA EL POLINOMIO INTERPOLADOR DE NEWTON EN UN PUNTO
real mn_evaluar_polinomio_interpolador(
    const Array1D< real > &A, /// VECTOR CON LOS COEFICIENTES DEL POLINOMIO
    const Array1D< real > &X, /// VECTOR CON LOS VALORES DE LOS PUNTOS DE INTERPOLACION
    const real x0)            /// VALOR DONDE SE INTERPOLA EL POLINOMIO
{
    /// Grado del polinomio
    int N = A.dim() - 1;

    /// Inicializamos la evaluación con el coeficiente de mayor grado
    real E = A[N];

    /// Evaluamos el polinomio de Newton usando una forma anidada
    /// similar al método de Horner:
    /// P(x) = A0 + (x-X0)[A1 + (x-X1)[A2 + ... ]]
    for(int k = N - 1; k >= 0; k--) {
        E = E * (x0 - X[k]) + A[k];
    }

    /// Devolvemos el valor del polinomio en x0
    return E;
}


#include "mn_aritmeticas.h"
#include "mn_interpolacion.h"

/// FUNCIÓN PARA CALCULAR UN ZOOM DE FACTOR z A UN Array2D< real >
/// USANDO INTERPOLACIÓN POR EL VECINO MÁS CERCANO
/// SE CONSIDERA QUE LA DISTANCIA ENTRE LOS PIXELS DEL Array2D ES 1
/// DADO UN PUNTO EN LA IMAGEN DE ÍNDICES (i,j) SI AL ACCEDER AL
/// PUNTO (i+1,j) NOS SALIMOS DE LA IMAGEN, TOMAMOS EN SU LUGAR EL PUNTO (i,j)
/// LO MISMO SI NOS SALIMOS AL ACCEDER A (i,j+1) O (i+1,j+1)
Array2D< real > zoom_vecino(
    Array2D< real > &A 	/** IMAGEN ORIGINAL */,
    real z 				/** FACTOR DE ZOOM */)
{
    /// Calculamos las nuevas dimensiones de la imagen tras aplicar el zoom
    int dim1New = A.dim1() * z;
    int dim2New = A.dim2() * z;

    /// Si las nuevas dimensiones no son válidas, devolvemos una imagen vacía
    if(dim1New <= 0 || dim2New <= 0) return Array2D< real >();

    /// Creamos la imagen de salida con las nuevas dimensiones
    Array2D< real > B(dim1New, dim2New);

    /// Recorremos todos los píxeles de la imagen ampliada/reducida
    for(int ip = 0; ip < B.dim1(); ip++) {
        for(int jp = 0; jp < B.dim2(); jp++) {

            /// Calculamos la posición correspondiente en la imagen original
            real x = ip / z;
            real y = jp / z;

            /// Tomamos la parte entera como punto base
            int i = x;
            int j = y;

            /// Si la parte decimal de x supera 0.5 y existe el píxel siguiente,
            /// escogemos el vecino de abajo
            if((x - i) > 0.5 && i + 1 < A.dim1()) i = i + 1;

            /// Si la parte decimal de y supera 0.5 y existe el píxel siguiente,
            /// escogemos el vecino de la derecha
            if((y - j) > 0.5 && j + 1 < A.dim2()) j = j + 1;

            /// Asignamos a la nueva imagen el valor del píxel más cercano
            B[ip][jp] = A[i][j];
        }
    }

    /// Devolvemos la imagen transformada
    return B;
}


#include "mn_aritmeticas.h"
#include "mn_interpolacion.h"

/// FUNCIÓN PARA CALCULAR UN ZOOM DE FACTOR z A UN Array2D< real >
/// USANDO INTERPOLACIÓN BILINEAL
/// SE CONSIDERA QUE LA DISTANCIA ENTRE LOS PIXELS DEL Array2D ES 1
/// DADO UN PUNTO EN LA IMAGEN DE ÍNDICES (i,j) SI AL ACCEDER AL
/// PUNTO (i+1,j) NOS SALIMOS DE LA IMAGEN, TOMAMOS EN SU LUGAR EL PUNTO (i,j)
/// LO MISMO SI NOS SALIMOS AL ACCEDER A (i,j+1) O (i+1,j+1)
Array2D< real > zoom_bilineal(
    Array2D< real > &A /** IMAGEN ORIGINAL */,
    real z /** FACTOR DE ZOOM */)
{
    /// Calculamos las nuevas dimensiones de la imagen tras aplicar el zoom
    int dim1New = A.dim1() * z;
    int dim2New = A.dim2() * z;

    /// Si las nuevas dimensiones no son válidas, devolvemos una imagen vacía
    if(dim1New <= 0 || dim2New <= 0) return Array2D< real >();

    /// Creamos la imagen resultante con las nuevas dimensiones
    Array2D< real > B(dim1New, dim2New);

    /// Recorremos todos los píxeles de la nueva imagen
    for(int ip = 0; ip < B.dim1(); ip++) {
        for(int jp = 0; jp < B.dim2(); jp++) {

            /// Calculamos la posición real correspondiente en la imagen original
            real x = ip / z;
            real y = jp / z;

            /// Obtenemos la parte entera de la posición
            /// que corresponde al píxel superior izquierdo
            int i = x;
            int j = y;

            /// Calculamos los índices del píxel inferior y derecho
            /// sin salirnos de la imagen original
            int i1 = min(i + 1, A.dim1() - 1);
            int j1 = min(j + 1, A.dim2() - 1);

            /// Parte decimal: distancia relativa dentro del cuadrado
            real dx = x - i;
            real dy = y - j;

            /// Interpolación bilineal:
            /// primero interpolamos en dirección x entre los dos puntos superiores
            /// y entre los dos puntos inferiores, y luego interpolamos en y
            B[ip][jp] =
                (1. - dy) * ( (1. - dx) * A[i][j]  + dx * A[i1][j] ) +
                     dy  * ( (1. - dx) * A[i][j1] + dx * A[i1][j1] );
        }
    }

    /// Devolvemos la imagen redimensionada
    return B;
}