/* nag_prob_lin_non_central_chi_sq (g01jcc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */

#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg01.h>

int main(void)
{
  /* Initialized data */
  Integer maxit = 500;
  double tol = 1e-4;

  /* Scalars */
  double c, p, pdf;
  Integer exit_status, i, n;

  NagError fail;

  /* Arrays */
  double *a = 0, *rlamda = 0;
  Integer *mult = 0;

  INIT_FAIL(fail);

  exit_status = 0;
  printf("nag_prob_lin_non_central_chi_sq (g01jcc) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");

  printf("\n        A     MULT  RLAMDA\n");
  while (scanf("%" NAG_IFMT "%lf%*[^\n] ", &n, &c) != EOF)
  {
    /* Allocate memory */
    if (!(a = NAG_ALLOC(n, double)) ||
        !(rlamda = NAG_ALLOC(n, double)) || !(mult = NAG_ALLOC(n, Integer)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    printf("\n");
    for (i = 1; i <= n; ++i)
      scanf("%lf", &a[i - 1]);
    scanf("%*[^\n] ");

    for (i = 1; i <= n; ++i)
      scanf("%" NAG_IFMT "", &mult[i - 1]);
    scanf("%*[^\n] ");
    for (i = 1; i <= n; ++i)
      scanf("%lf", &rlamda[i - 1]);
    scanf("%*[^\n] ");

    /* nag_prob_lin_non_central_chi_sq (g01jcc).
     * Computes probability for a positive linear combination of
     * chi^2 variables
     */
    nag_prob_lin_non_central_chi_sq(a, mult, rlamda, n, c, &p, &pdf, tol,
                                    maxit, &fail);
    if (fail.code == NE_NOERROR || fail.code == NE_ACCURACY ||
        fail.code == NE_PROB_BOUNDARY) {
      for (i = 1; i <= n; ++i)
        printf(" %10.2f%6" NAG_IFMT "%9.2f\n", a[i - 1], mult[i - 1],
               rlamda[i - 1]);
      printf("c = %6.2f    Prob = %6.4f\n", c, p);
    }
    else {
      printf("Error from nag_normal_scores_exact (g01dac).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

    NAG_FREE(a);
    NAG_FREE(rlamda);
    NAG_FREE(mult);
  }

END:
  NAG_FREE(a);
  NAG_FREE(rlamda);
  NAG_FREE(mult);
  return exit_status;
}