/*
================================================================================
  RPL/2 (R) version 4.1.32
  Copyright (C) 1989-2020 Dr. BERTRAND Joël

  This file is part of RPL/2.

  RPL/2 is free software; you can redistribute it and/or modify it
  under the terms of the CeCILL V2 License as published by the french
  CEA, CNRS and INRIA.
 
  RPL/2 is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  FITNESS FOR A PARTICULAR PURPOSE.  See the CeCILL V2 License
  for more details.
 
  You should have received a copy of the CeCILL License
  along with RPL/2. If not, write to info@cecill.info.
================================================================================
*/


#include "rpl-conv.h"


/*
================================================================================
  Fonction 'integrale_romberg'
================================================================================
  Entrées : pointeur sur une structure struct_processus
--------------------------------------------------------------------------------
  Sorties :
--------------------------------------------------------------------------------
  Effets de bord : néant
================================================================================
*/

void
integrale_romberg(struct_processus *s_etat_processus,
		struct_objet *s_expression, unsigned char *variable,
		real8 a, real8 b, real8 precision)
{
	real8						erreur;
	real8						h;
	real8						**t;
	real8						**t_tampon;
	real8						*vecteur;
	real8						*vecteur_precedent;
	real8						x;

	logical1					validite;
	logical1					erreur_memoire;

	struct_objet				*s_objet;

	struct_variable				s_variable;

	integer8					i;
	integer8					j;
	integer8					k;
	integer8					n;
	integer8					p;
	integer8					nombre_elements;
	integer8					nombre_termes;
	integer8					taille_vecteur_precedent;

	/*
	 * Création d'une variable locale représentant la variable d'intégration
	 */

	(*s_etat_processus).niveau_courant++;
	s_variable.niveau = (*s_etat_processus).niveau_courant;

	if ((s_variable.nom = (unsigned char *) malloc((strlen(variable) + 1) *
			sizeof(unsigned char))) == NULL)
	{
		(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
		return;
	}

	strcpy(s_variable.nom, variable);

	if ((s_variable.objet = allocation(s_etat_processus, REL)) == NULL)
	{
		(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
		return;
	}

	if (creation_variable(s_etat_processus, &s_variable, 'V', 'P') == d_erreur)
	{
		return;
	}

	/*
	 * Algorithme de Romberg
	 */

	if ((vecteur = (real8 *) malloc(2 * sizeof(real8))) == NULL)
	{
		(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
		return;
	}

	evaluation_romberg(s_etat_processus, s_expression, variable, &a,
			&(vecteur[nombre_termes = 0]), &validite);
	nombre_termes += (validite == d_vrai) ? 1 : 0;

	evaluation_romberg(s_etat_processus, s_expression, variable, &b,
			&(vecteur[nombre_termes]), &validite);
	nombre_termes += (validite == d_vrai) ? 1 : 0;

	if ((t = (real8 **) malloc(sizeof(real8 *))) == NULL)
	{
		(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
		return;
	}

	if ((t[0] = (real8 *) malloc(sizeof(real8))) == NULL)
	{
		(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
		return;
	}

	t[0][n = 0] = ((b - a) / 2) * sommation_vecteur_reel(vecteur,
			&nombre_termes, &erreur_memoire);

	if (erreur_memoire == d_vrai)
	{
		(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
		return;
	}

	free(vecteur);

	vecteur_precedent = NULL;
	taille_vecteur_precedent = 0;

	do
	{
		h = (b - a) / ((real8) (nombre_elements = (((integer8) 1) << (++n))));

		if (nombre_elements == 0)
		{
			// Dépassement de capacité
			n--;
			break;
		}

		x = a;

		/*
		 * Nouvelle allocation de t
		 */

		t_tampon = t;

		if ((t = (real8 **) malloc((((size_t) n) + 1) * sizeof(real8 *)))
				== NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return;
		}

		for(i = 0; i <= n; i++)
		{
			if ((t[i] = (real8 *) malloc((((size_t) n) + 1) * sizeof(real8)))
					== NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if (i != n)
			{
				for(t[i][n] = 0, j = 0; j < n; j++)
				{
					t[i][j] = t_tampon[i][j];
				}
			}
			else
			{
				for(j = 0; j <= n; t[i][j++] = 0);
			}
		}

		for(i = 0; i < n; free(t_tampon[i++]));
		free(t_tampon);

		/*
		 * Boucle principale
		 */

		if ((vecteur = (real8 *) malloc((((size_t) nombre_elements) + 1) *
				sizeof(real8))) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return;
		}


		if (vecteur_precedent == NULL)
		{
			/*
			 * Première itération
			 */

			evaluation_romberg(s_etat_processus, s_expression, variable, &a,
					&(vecteur[nombre_termes = 0]), &validite);
			nombre_termes += (validite == d_vrai) ? 1 : 0;

			evaluation_romberg(s_etat_processus, s_expression, variable, &b,
					&(vecteur[nombre_termes]), &validite);
			nombre_termes += (validite == d_vrai) ? 1 : 0;

			if (nombre_termes == 2)
			{
				vecteur[0] += vecteur[1];
				vecteur[0] /= 2;
				nombre_termes--;
			}

			for(i = 1; i < nombre_elements; i++)
			{
				x += h;

				evaluation_romberg(s_etat_processus, s_expression, variable, &x,
						&(vecteur[nombre_termes]), &validite);
				nombre_termes += (validite == d_vrai) ? 1 : 0;
			}
		}
		else
		{
			/*
			 * Pour les itérations suivantes, la moitié des points
			 * a déjà été calculée.
			 */

			for(i = 0; i < taille_vecteur_precedent; i++)
			{
				vecteur[i] = vecteur_precedent[i];
			}

			free(vecteur_precedent);

			for(nombre_termes = taille_vecteur_precedent, i = 1;
					i < nombre_elements; x += h, i += 2)
			{
				x += h;

				evaluation_romberg(s_etat_processus, s_expression, variable, &x,
						&(vecteur[nombre_termes]), &validite);
				nombre_termes += (validite == d_vrai) ? 1 : 0;
			}
		}

		t[0][n] = h * sommation_vecteur_reel(vecteur, &nombre_termes,
				&erreur_memoire);

		if (erreur_memoire == d_vrai)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return;
		}

		vecteur_precedent = vecteur;
		taille_vecteur_precedent = nombre_termes;

		for(p = 1; p <= n; p++)
		{
			k = n - p;
			t[p][k] = ((pow(4, (real8) p) * t[p - 1][k + 1]) - t[p - 1][k]) /
					((real8) (pow(4, (real8) p) - 1));
		}
	} while(((erreur = fabs(t[n][0] - t[n - 1][0])) > precision) &&
			((*s_etat_processus).var_volatile_requete_arret == 0));

	free(vecteur_precedent);

	/*
	 * Empilement du résultat et liberation de t
	 */

	if ((s_objet = allocation(s_etat_processus, REL)) == NULL)
	{
		(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
		return;
	}

	(*((real8 *) (*s_objet).objet)) = t[n][0];

	if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
			s_objet) == d_erreur)
	{
		return;
	}

	if ((s_objet = allocation(s_etat_processus, REL)) == NULL)
	{
		(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
		return;
	}

	(*((real8 *) (*s_objet).objet)) = erreur;

	if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
			s_objet) == d_erreur)
	{
		return;
	}

	for(i = 0; i <= n; free(t[i++]));
	free(t);

	/*
	 * Libération de la variable locale
	 */

	(*s_etat_processus).niveau_courant--;

	if (retrait_variable(s_etat_processus, s_variable.nom, 'L')
			== d_erreur)
	{
		if ((*s_etat_processus).erreur_systeme != d_es)
		{
			if ((*s_etat_processus).erreur_systeme ==
					d_es_variable_introuvable)
			{
				(*s_etat_processus).erreur_systeme = d_es;
			}
			else
			{
				/*
				 * Erreur système
				 */

				return;
			}
		}

		(*s_etat_processus).erreur_execution = d_ex_variable_non_definie;
		return;
	}

	return;
}


void
evaluation_romberg(struct_processus *s_etat_processus,
		struct_objet *s_expression, unsigned char *variable, real8 *point,
		real8 *valeur, logical1 *validite)
{
	integer8							hauteur_pile;

	struct_liste_pile_systeme			*l_position_normale;

	struct_objet						*s_objet;

	/*
	 * Vérification du type de la variable
	 */

	(*validite) = d_faux;
	(*valeur) = 0;

	if (recherche_variable(s_etat_processus, variable) == d_vrai)
	{
		if ((*(*s_etat_processus).pointeur_variable_courante)
				.variable_verrouillee == d_vrai)
		{
			(*s_etat_processus).erreur_execution = d_ex_variable_verrouillee;
			return;
		}

		if ((*(*(*s_etat_processus).pointeur_variable_courante)
				.objet).type != REL)
		{
			liberation(s_etat_processus, (*(*s_etat_processus)
					.pointeur_variable_courante).objet);

			if (((*(*s_etat_processus).pointeur_variable_courante).objet =
					allocation(s_etat_processus, REL)) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}
		}

		(*((real8 *) (*(*(*s_etat_processus).pointeur_variable_courante)
				.objet).objet)) = (*point);
	}
	else
	{
		/*
		 * La variable d'intégration étant locale, l'utilisateur ne peut
		 * pas l'effacer. On ne doit donc jamais passer par ici. Si
		 * c'est néanmoins le cas, une erreur système est générée
		 * provoquant l'arrêt du programme même dans une
		 * structure IFERR.
		 */

		return;
	}

	hauteur_pile = (*s_etat_processus).hauteur_pile_operationnelle;
	l_position_normale = (*s_etat_processus).l_base_pile_systeme;

	(*s_etat_processus).erreur_execution = d_ex;
	(*s_etat_processus).exception = d_ep;

	if (evaluation(s_etat_processus, s_expression, 'N') == d_erreur)
	{
		if ((*s_etat_processus).erreur_systeme != d_es)
		{
			/*
			 * Erreur système
			 */

			return;
		}
		else
		{
			/*
			 * Restauration de la pile initiale
			 */

			while((*s_etat_processus).hauteur_pile_operationnelle >
					hauteur_pile)
			{
				if (depilement(s_etat_processus, &((*s_etat_processus)
						.l_base_pile), &s_objet) == d_erreur)
				{
					(*s_etat_processus).erreur_execution =
							d_ex_manque_argument;
					break;
				}

				liberation(s_etat_processus, s_objet);
			}

			/*
			 * Restauration de la pile système
			 */

			while((*s_etat_processus).l_base_pile_systeme !=
					l_position_normale)
			{
				depilement_pile_systeme(s_etat_processus);

				if ((*s_etat_processus).erreur_systeme != d_es)
				{
					/*
					 * Une pile vide provoque une erreur système
					 */

					return;
				}
			}
		}

		(*s_etat_processus).erreur_execution = d_ex;
		(*s_etat_processus).exception = d_ep;
	}
	else
	{
		/*
		 * Retrait de la valeur de la pile
		 */

		if (depilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
				&s_objet) == d_erreur)
		{
			(*s_etat_processus).erreur_execution = d_ex_manque_argument;
			return;
		}

		if ((*s_objet).type == INT)
		{
			(*valeur) = (real8) (*((integer8 *) (*s_objet).objet));
			(*validite) = d_vrai;
		}
		else if ((*s_objet).type == REL)
		{
			(*valeur) = (*((real8 *) (*s_objet).objet));
			(*validite) = d_vrai;
		}

		liberation(s_etat_processus, s_objet);
	}

	return;
}

// vim: ts=4