/* nag_specfun_hyperg_confl_real_scaled (s22bbc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.2, 2023.
 */

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

int main(void) {
  /* Scalars */
  Integer exit_status = 0;
  Integer k, maxexponent, scm;
  double ani, adr, bni, bdr, delta, frm, x;
  /* Arrays */
  double frmv[2];
  Integer scmv[2];
  /* Nag Types */
  NagError fail;

  maxexponent = X02BLC;
  printf(
      "nag_specfun_hyperg_confl_real_scaled (s22bbc) Example Program Results\n\n");

  ani = -10.0;
  bni = 30.0;
  delta = 1.0E-4;
  adr = delta;
  bdr = -delta;
  x = 25.0;

  printf("%9s%10s%10s%12s%9s%12s\n", "a", "b", "x", "frm", "scm", "M(a,b,x)");
  for (k = 0; k < 2; k++) {
    INIT_FAIL(fail);
    /* Compute the real confluent hypergeometric function M(a,b,x) in scaled
     * form using nag_specfun_hyperg_confl_real_scaled (s22bbc).
     */
    nag_specfun_hyperg_confl_real_scaled(ani, adr, bni, bdr, x, &frm, &scm,
                                         &fail);
    switch (fail.code) {
    case NE_NOERROR:
    case NW_UNDERFLOW_WARN:
    case NW_SOME_PRECISION_LOSS: {
      if (scm < maxexponent)
        printf(" %9.4f %9.4f %9.4f %13.4e %5" NAG_IFMT " %13.4e\n", ani + adr,
               bni + bdr, x, frm, scm, frm * pow(2.0, scm));
      else
        printf(" %9.4f %9.4f %9.4f %13.4e %5" NAG_IFMT " %17s\n", ani + adr,
               bni + bdr, x, frm, scm, "Not Representable");
      frmv[k] = frm;
      scmv[k] = scm;
      break;
    }
    default: {
      /* Either the result has overflowed, no accuracy may be assumed,
       * or an input error has been detected.
       */
      printf(" %9.4f %9.4f %9.4f %17s\n", ani + adr, bni + bdr, x, "FAILED");
      exit_status = 1;
      goto END;
    }
    }
    adr = -adr;
    bdr = -bdr;
  }

  /* Calculate the product M1*M2 */
  frm = frmv[0] * frmv[1];
  scm = scmv[0] + scmv[1];
  printf("\n");
  if (scm < maxexponent)
    printf("%-30s%12.4e%6" NAG_IFMT "%12.4e\n", "Solution product", frm, scm,
           frm * pow(2.0, scm));
  else
    printf("%-30s%12.4e%6" NAG_IFMT "%17s\n", "Solution product", frm, scm,
           "Not Representable");

  /* Calculate the ratio M1/M2 */
  if (frmv[1] != 0.0) {
    frm = frmv[0] / frmv[1];
    scm = scmv[0] - scmv[1];
    printf("\n");
    if (scm < maxexponent)
      printf("%-30s%12.4e%6" NAG_IFMT "%12.4e\n", "Solution ratio", frm, scm,
             frm * pow(2.0, scm));
    else
      printf("%-30s%12.4e%6" NAG_IFMT "%17s\n", "Solution ratio", frm, scm,
             "Not Representable");
  }

END:
  return exit_status;
}