/*
================================================================================
  RPL/2 (R) version 4.1.36
  Copyright (C) 1989-2024 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 réalisation le tri d'un vecteur de réels du plus petit au plus
  grand en valeur absolue (algorithme de tri dit de Shell-Metzner)
================================================================================
  Entrées : pointeur sur une structure struct_processus
--------------------------------------------------------------------------------
  Sorties :
--------------------------------------------------------------------------------
  Effets de bord : néant
================================================================================
*/

void
tri_vecteur(real8 *vecteur, integer8 taille)
{
	logical1					terminaison_boucle_1;
	logical1					terminaison_boucle_2;
	logical1					terminaison_boucle_3;

	real8						registre;

	integer8					ecartement;
	integer8					indice_i;
	integer8					indice_j;
	integer8					indice_k;
	integer8					indice_l;

	ecartement = taille;
	terminaison_boucle_1 = d_faux;

	do
	{
		ecartement = ecartement / 2;

		if (ecartement >= 1)
		{
			indice_j = 0;
			indice_k = taille - ecartement;
			terminaison_boucle_2 = d_faux;

			do
			{
				indice_i = indice_j;
				terminaison_boucle_3 = d_faux;

				do
				{
					indice_l = indice_i + ecartement;

					if ((indice_i > 0) && (indice_l > 0))
					{
						if (abs(vecteur[indice_i - 1]) >
								abs(vecteur[indice_l - 1]))
						{
							registre = vecteur[indice_i - 1];
							vecteur[indice_i - 1] = vecteur[indice_l - 1];
							vecteur[indice_l - 1] = registre;

							indice_i -= ecartement;

							if (indice_i < 1)
							{
								terminaison_boucle_3 = d_vrai;
							}
						}
						else
						{
							terminaison_boucle_3 = d_vrai;
						}
					}
					else
					{
						terminaison_boucle_3 = d_vrai;
					}
				} while(terminaison_boucle_3 == d_faux);

				indice_j++;

				if (indice_j > indice_k)
				{
					terminaison_boucle_2 = d_vrai;
				}
			} while(terminaison_boucle_2 == d_faux);
		}
		else
		{
			terminaison_boucle_1 = d_vrai;
		}
	} while(terminaison_boucle_1 == d_faux);
}


/*
================================================================================
  Fonction réalisation la sommation d'un vecteur de réel sans perte
  de précision
================================================================================
  Entrées : pointeur sur une structure struct_processus
--------------------------------------------------------------------------------
  Sorties :
--------------------------------------------------------------------------------
  Effets de bord : néant
================================================================================
*/

real8
sommation_vecteur_reel(real8 *vecteur, integer8 *taille,
		logical1 *erreur_memoire)
{
#if 0
	unsigned long				nombre_elements;
	unsigned long				pointeur;

	/*
	 * Sommation des termes en commençant par le plus petit.
	 * Algorithme optimal mais NP-complet...
	 */

	nombre_elements = (*taille);
	(*erreur_memoire) = d_faux;

	while(nombre_elements != 1)
	{
		pointeur = (*taille) - nombre_elements;
		tri_vecteur(&(vecteur[pointeur]), nombre_elements);
		vecteur[pointeur + 1] += vecteur[pointeur];
		nombre_elements--;
	}

	return(vecteur[(*taille) - 1]);
#else
	real8				erreur;
	real8				somme;
	real8				registre;
	real8				tampon;

	integer8			i;

	somme = 0;
	erreur = 0;

	(*erreur_memoire) = d_faux;

	for(i = 0; i < (*taille); i++)
	{
		tampon = somme;
		registre = vecteur[i] + erreur;
		somme = tampon + registre;
		erreur = (tampon - somme) + registre;
	}

	return(somme);
#endif
}


/*
================================================================================
  Fonction réalisation la sommation d'un vecteur de complexes sans perte
  de précision
================================================================================
  Entrées : pointeur sur une structure struct_processus
--------------------------------------------------------------------------------
  Sorties :
--------------------------------------------------------------------------------
  Effets de bord : néant
================================================================================
*/

complex16
sommation_vecteur_complexe(complex16 *vecteur, integer8 *taille,
		logical1 *erreur_memoire)
{
	complex16					cumul;

	real8						*tampon;

	integer8					i;
	integer8					nombre_elements;

 	if ((tampon = sys_malloc(((size_t) (*taille)) * sizeof(real8))) == NULL)
	{
		(*erreur_memoire) = d_vrai;

		cumul.partie_reelle = 0;
		cumul.partie_imaginaire = 0;

		return(cumul);
	}

	(*erreur_memoire) = d_faux;

	/*
	 * Sommation des termes en commençant par le plus petit
	 */

	for(i = 0, nombre_elements = (*taille); i < nombre_elements;
			tampon[i] = vecteur[i].partie_reelle, i++);
			
	cumul.partie_reelle = sommation_vecteur_reel(tampon, taille,
			erreur_memoire);

	/*
	 * Même traitement, mais sur la partie imaginaire
	 */

	for(i = 0, nombre_elements = (*taille); i < nombre_elements;
			tampon[i] = vecteur[i].partie_imaginaire, i++);
			
	cumul.partie_imaginaire = sommation_vecteur_reel(tampon, taille,
			erreur_memoire);

	sys_free(tampon);

	return(cumul);
}

// vim: ts=4