/* nag_specfun_hyperg_gauss_real_scaled (s22bfc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.2, 2023.
 */

#include <math.h>
#include <nag.h>
#include <stdio.h>

int main(void) {
  /* Scalars */
  Integer exit_status = 0;
  Integer k, imax, scf;
  double ani, adr, bni, bdr, cni, cdr, delta, frf, x;
  /* Arrays */
  double frfv[2];
  Integer scfv[2];
  /* Nag Types */
  Nag_Boolean finite_solutions;
  NagError fail;

  imax = X02BLC;
  printf(
      "nag_specfun_hyperg_gauss_real_scaled (s22bfc) Example Program Results\n\n");

  ani = -10.0;
  bni = 2.0;
  cni = -5.0;
  delta = 1.0E-4;
  adr = delta;
  bdr = -delta;
  cdr = delta;
  x = 0.45;
  finite_solutions = Nag_TRUE;
  printf("%11s%11s%11s%11s%14s%7s%14s\n", "a", "b", "c", "x", "frf", "scf",
         "2F1(a,b;c;x)");
  for (k = 0; k < 2; k++) {
    INIT_FAIL(fail);
    /* Compute the real Gauss hypergeometric function 2F1(a,b;c;x) in scaled
     * form using nag_specfun_hyperg_gauss_real_scaled (s22bfc).
     */
    nag_specfun_hyperg_gauss_real_scaled(ani, adr, bni, bdr, cni, cdr, x, &frf,
                                         &scf, &fail);
    switch (fail.code) {
    case NE_NOERROR:
    case NW_UNDERFLOW_WARN:
    case NW_SOME_PRECISION_LOSS:
      /* A finite result has been returned. */
      if (scf < imax)
        printf(" %10.4f %10.4f %10.4f %10.4f %13.5e %6" NAG_IFMT " %13.5e\n",
               ani + adr, bni + bdr, cni + cdr, x, frf, scf,
               frf * pow(2.0, scf));
      else
        printf(" %10.4f %10.4f %10.4f %10.4f %13.5e %6" NAG_IFMT " %17s\n",
               ani + adr, bni + bdr, cni + cdr, x, frf, scf,
               "Not Representable");
      frfv[k] = frf;
      scfv[k] = scf;
      break;
    case NE_INFINITE:
      /* The result is analytically infinite. */
      finite_solutions = Nag_FALSE;
      if (frf >= 0.0)
        printf(" %10.4f %10.4f %10.4f %10.4f %13s %6" NAG_IFMT " %13s\n",
               ani + adr, bni + bdr, cni + cdr, x, "Inf", scf, "Inf");
      else
        printf(" %10.4f %10.4f %10.4f %10.4f %13s %6" NAG_IFMT " %13s\n",
               ani + adr, bni + bdr, cni + cdr, x, "-Inf", scf, "-Inf");
      break;
    case NW_OVERFLOW_WARN:
    case NE_OVERFLOW:
      /* The final result has overflowed. */
      finite_solutions = Nag_FALSE;
      if (frf >= 0.0)
        printf(" %10.4f %10.4f %10.4f %10.4f %13.5e %6s %13s\n", ani + adr,
               bni + bdr, cni + cdr, x, frf, "imax", ">pow(2,imax)");
      else
        printf(" %10.4f %10.4f %10.4f %10.4f %13.5e %6s %13s\n", ani + adr,
               bni + bdr, cni + cdr, x, frf, "imax", "<-pow(2,imax)");
      break;
    case NE_CANNOT_CALCULATE:
      /* An internal calculation resulted in an undefined result. */
      finite_solutions = Nag_FALSE;
      printf(" %10.4f %10.4f %10.4f %10.4f %13s %6" NAG_IFMT " %13s\n",
             ani + adr, bni + bdr, cni + cdr, x, "NaN", scf, "NaN");
      break;
    default:
      /* An input error has been detected. */
      printf(" %10.4f %10.4f %10.4f %10.4f %17s\n", ani + adr, bni + bdr,
             cni + cdr, x, "FAILED");
      exit_status = 1;
      goto END;
      break;
    }
    adr = -adr;
    bdr = -bdr;
    cdr = -cdr;
  }
  if (finite_solutions) {
    /* Calculate the product M1*M2. */
    frf = frfv[0] * frfv[1];
    scf = scfv[0] + scfv[1];
    printf("\n");
    if (scf < imax)
      printf("%-34s%13.5e %6" NAG_IFMT " %13.5e\n", " Solution product", frf,
             scf, frf * pow(2.0, scf));
    else
      printf("%-34s%13.5e %6" NAG_IFMT "%17s\n", " Solution product", frf, scf,
             "Not Representable");

    /* Calculate the ratio M1/M2. */
    if (frfv[1] != 0.0) {
      frf = frfv[0] / frfv[1];
      scf = scfv[0] - scfv[1];
      printf("\n");
      if (scf < imax)
        printf("%-34s%13.5e %6" NAG_IFMT " %13.5e\n", " Solution ratio", frf,
               scf, frf * pow(2.0, scf));
      else
        printf("%-34s%13.5e %6" NAG_IFMT "%17s\n", " Solution ratio", frf, scf,
               "Not Representable");
    }
  }
END:
  return exit_status;
}