#include <stdlib.h>
#include <stdio.h>
#include <mpir.h>
#include <flint.h>
#include <fmpz.h>
#include <arith.h>

#include <mpfr.h>


int main(int argc, char* argv[])
{
	/* Variablen n und result erzeugen */
	fmpz_t n;

	fmpz_t sigma_n;
	mpz_t sigma_n_mpz;

	mpq_t harmonic_n;
	mpfr_t harmonic_n_mpfr;

	mpfr_t exp_harmonic_n_mpfr;
	mpfr_t ln_harmonic_n_mpfr;

	mpfr_t lagarias_result;


	int i;

	fmpz_init(n);
	fmpz_init(sigma_n);

	mpz_init(sigma_n_mpz);

	mpq_init(harmonic_n);
	mpfr_init(harmonic_n_mpfr);
	mpfr_init(ln_harmonic_n_mpfr);
	mpfr_init(exp_harmonic_n_mpfr);

	mpfr_init(lagarias_result);

	for(i=0;i<1000000;i++)
	{
		fmpz_set_si(n, i);
		puts("Teilersumme");
		/* Berechne die Teilersumme */
		fmpz_divisor_sigma(sigma_n, n, 1);
		/* In MPIR konvertieren */
		fmpz_get_mpz(sigma_n_mpz, sigma_n);

		puts("Harmonic");
		/* Berechne das Element i der Harmonischen Reihe */
		mpq_harmonic(harmonic_n, i);
		/* In eine MPFR-Gleitkommazahl umwandeln */
		mpfr_set_q(harmonic_n_mpfr, harmonic_n, MPFR_RNDN);

		puts("exp");
		/* e^H_n berechnen */
		mpfr_exp(exp_harmonic_n_mpfr, harmonic_n_mpfr, MPFR_RNDN);
		/* ln(H_n) berechnen */
		mpfr_log(exp_harmonic_n_mpfr, harmonic_n_mpfr, MPFR_RNDN);


		puts("mul");
		mpfr_mul(lagarias_result, exp_harmonic_n_mpfr, exp_harmonic_n_mpfr, MPFR_RNDN);
		puts("add");
		mpfr_add_q(lagarias_result, lagarias_result, harmonic_n, MPFR_RNDN);
		puts("sub");
		mpfr_sub_z(lagarias_result, lagarias_result, sigma_n_mpz, MPFR_RNDN);


		mpfr_printf("Ergebnis: %.128Rf\n", lagarias_result);
	}

	fmpz_clear(n);
	fmpz_clear(sigma_n);

	return 0;
}