/*
 * Calcula la funcion zeta de Riemann para una region cuadrada, esta se particiona en n cuadrados
 * donde n es el numero de CPUs de tu maquina
 * Modifica la region cambiando en el main() la variable U que contiene dos numeros complejos de la forma a+b*I U.v y U.w
 * Este codigo solo funciona en linux
 *
 * compila:
 *
 * gcc -Wall -lm -lpth estecodigo.c -o out
 *
 * ejecutas:
 *
 * ./out
 *
 * toorandom@gmail.com
 *
 * Eduardo Ruiz Duarte (beck)
 *
 */
#define _GNU_SOURCE
#include <pthread.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <sched.h>
#include <unistd.h>
#include <math.h>
#include <complex.h>

/* Este es nuestro infinito para la serie :P */
#define INF 3000000
/* tenemos que crear una red de puntos porque no podemos calcular para todo numero complejo y esta sera separa por el DELTA */
#define DELTA_RE 0.1
#define DELTA_IM 0.1
/* CPUs en tu maquina */
int num_cpus;			/* via sysconf() con _SC_NPROCESSORS_CONF */


typedef struct cuadrado
{
  long cpuid;
  complex v, w;			/* Cuadrado en el plano complejo delimitado por v y w */
} cuadrado_t;


/* funcion zeta de Riemann acotada tradicional acotada */

double complex
zeta (complex s)
{
  double n;
  double complex r;
  for (n = 1; n < INF; n++)
    r += 1 / cpow (n, s);
  return r;
}

/* Funcion que se paralelizara calculando la zeta de una region delimitada por dos numeros complejos v y w */

void *
zeta_region (void *region)
{
  long tid;			/* thread id */
  double s_im, s_re;
  complex s, r;
  cuadrado_t *C;
  cpu_set_t cpuset;
  pthread_t thread;
  C = region;
  tid = C->cpuid;
  thread = pthread_self ();
  CPU_ZERO (&cpuset);
  CPU_SET (tid, (cpu_set_t *) & cpuset);
  if (pthread_setaffinity_np (thread, sizeof (cpu_set_t), &cpuset) != 0)
    perror ("pthread_setaffinity_np");

  printf ("CPU(%ld): funcion zeta de Riemann del cuadrado (%f,%f)x(%f,%f)\n", tid,creal(C->v),cimag(C->v),creal(C->w),cimag(C->w));

  for (s_re = creal (C->v); s_re < creal (C->w); s_re += DELTA_RE)
    for (s_im = cimag (C->v); s_im < cimag (C->w); s_im += DELTA_IM)
      {
	s = s_re + s_im * I;
	r = zeta (s);
	printf ("(cpu=%d):z(%f+%fi)=(%f,%f)\n",tid,s_re,s_im,creal (r), cimag (r));
      }
  pthread_exit (NULL);
}

int
main (int argc, char *argv[])
{
  /* obtenemos primero num_cpus */
  num_cpus = sysconf (_SC_NPROCESSORS_CONF);

  pthread_t threads[num_cpus];	/* numero de threads que haremos */
  long t;			/* id de thread a crear */
  int p;                        /* numero de particiones es el numero de cpus */
  int k=0;
  double Re1,Re2;               /* particionaremos con cuadrados verticales, necesitamos dos valores reales */
  double j;                     
  cuadrado_t U;                 /* Cuadrado que particionaremos */
  p = num_cpus;                 /* particiones coinciden con el numero de cpus */
  cuadrado_t *C;                /* p cuadrados */
  C = calloc(p,sizeof(struct cuadrado));
  memset (C, 0, sizeof (struct cuadrado));
  /* Region de cuadrado GRANDE dada por dos numeros complejos de la forma a+b*I */
  U.v = 1+1*I;
  U.w = 5+5*I;
  printf ("El numero de procesadores que hay disponibles es %d\n", num_cpus);
  for (t = 0; t < num_cpus; t++)
    {
      /* Asignamos el id de procesador t a la estructura */
      C[t].cpuid = t;
      /* Particionamos en p cuadrados del cuadradote U */
      for(j=0;j<(double)p;j++) {
              Re1 = (U.v+(fabs(creal(U.v)-creal(U.w))/p)*j);
              Re2 = (U.v+(fabs(creal(U.v)-creal(U.w))/p)*(j+1));
              C[k].v = Re1+cimag(U.v)*I;
              C[k].w = Re2+cimag(U.w)*I;
              k++;
      }
      /* creamos el thread t-esimo y le mandamos el cuadradito para que lo procese es CPU */
      if (pthread_create
	  (&threads[t], NULL, zeta_region, (void *) &C[t]) != 0)
	{
	  perror ("pthread_create");
	  exit (EXIT_FAILURE);
	}
    }
  pthread_exit (NULL);
}