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

#include <stdio.h>
#include <string.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, 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_1f1_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_1f1_real_scaled (s22bbc).
       */
      nag_specfun_1f1_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 %5ld %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 %5ld %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%6ld%12.4e\n",
           "Solution product", frm, scm, frm*pow(2.0, scm));
  else
    printf("%-30s%12.4e%6ld%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%6ld%12.4e\n",
               "Solution ratio", frm, scm, frm*pow(2.0, scm));
      else
        printf("%-30s%12.4e%6ld%17s\n",
               "Solution ratio", frm, scm, "Not Representable");
    }

 END:
  return exit_status;
}