/* nag_specfun_2f1_real_scaled (s22bfc) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 24, 2013.
 */

#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nags.h>
#include <nagx02.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_2f1_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_2f1_real_scaled (s22bfc).
       */
      nag_specfun_2f1_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 %6ld %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 %6ld %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 %6ld %13s\n",
                 ani+adr, bni+bdr, cni+cdr, x, "Inf", scf, "Inf");
        else
          printf(" %10.4f %10.4f %10.4f %10.4f %13s %6ld %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 %6ld %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 %6ld %13.5e\n",
               " Solution product", frf, scf, frf*pow(2.0, scf));
      else
        printf("%-34s%13.5e %6ld%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 %6ld %13.5e\n",
                   " Solution ratio", frf, scf, frf*pow(2.0, scf));
          else
            printf("%-34s%13.5e %6ld%17s\n",
                   " Solution ratio", frf, scf, "Not Representable");
        }
    }
 END:
  return exit_status;
}