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

# include "toms743.h"

int main ( );
void test01 ( int nbits );
void test02 ( int nbits, double dx, int n );
void test03 ( int nbits, double xmin, double xmax, int n );

/******************************************************************************/

int main ( )

/******************************************************************************/
/*
  Purpose:

    MAIN is the main program for TOMS743_TEST.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    16 June 2014

  Author:

    Original FORTRAN77 version by Andrew Barry, S. J. Barry, 
    Patricia Culligan-Hensley.
    This C version by John Burkardt.

  Reference:

    Andrew Barry, S. J. Barry, Patricia Culligan-Hensley,
    Algorithm 743: WAPR - A Fortran routine for calculating real 
    values of the W-function,
    ACM Transactions on Mathematical Software,
    Volume 21, Number 2, June 1995, pages 172-181.
*/
{
  double dx;
  int n;
  int nbits;
  double xmax;
  double xmin;

  timestamp ( );
  printf ( "\n" );
  printf ( "TOMS743_TEST\n" );
  printf ( "  C version\n" );
  printf ( "  Test the TOMS743 library.\n" );

  nbits = nbits_compute ( );
  printf ( "\n" );
  printf ( "  Number of bits in mantissa - 1 = %d\n", nbits );

  test01 ( nbits );

  dx = + 1.0E-09;
  n = 10;
  test02 ( nbits, dx, n );

  xmin = 0.0;
  xmax = 1.0E+20;
  n = 20;
  test03 ( nbits, xmin, xmax, n );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "TOMS743_TEST\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

  return 0;
}
/******************************************************************************/

void test01 ( int nbits )

/******************************************************************************/
/*
  Purpose:

    TEST01 compares WAPR to stored values.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    17 June 2014

  Author:

    Original FORTRAN77 version by Andrew Barry, S. J. Barry, 
    Patricia Culligan-Hensley.
    This C version by John Burkardt.

  Parameters:

    Input, int NBITS, the number of bits in the mantissa.
*/
{
  static double dx1[68] = {
    1.E-40,2.E-40,3.E-40,4.E-40,5.E-40,6.E-40,7.E-40,8.E-40, 
    9.E-40,1.E-39,1.E-30,2.E-30,3.E-30,4.E-30,5.E-30,6.E-30, 
    7.E-30,8.E-30,9.E-30,1.E-29,1.E-20,2.E-20,3.E-20,4.E-20, 
    5.E-20,6.E-20,7.E-20,8.E-20,9.E-20,1.E-19,1.E-10,2.E-10, 
    3.E-10,4.E-10,5.E-10,6.E-10,7.E-10,8.E-10,9.E-10,1.E-9,1.E-5, 
    2.E-5,3.E-5,4.E-5,5.E-5,6.E-5,7.E-5,8.E-5,9.E-5,1.E-4,2.E-4, 
    3.E-4,4.E-4,5.E-4,6.E-4,7.E-4,8.E-4,9.E-4,1.E-3,2.E-3,3.E-3, 
    4.E-3,5.E-3,6.E-3,7.E-3,8.E-3,9.E-3,1.E-2 };
  double em;
  int i;
  int nd;
  int nerror;
  double w;
  static double wm1[68] = {
           -1.000000000000000000023316439815971242034,
           -1.000000000000000000032974425414002562937,
           -1.000000000000000000040385258412884114431,
           -1.000000000000000000046632879631942484068,
           -1.000000000000000000052137144421794383842,
           -1.000000000000000000057113380167442850121,
           -1.000000000000000000061689501212464534966,
           -1.000000000000000000065948850828005125875,
           -1.000000000000000000069949319447913726103,
           -1.000000000000000000073733056744706376448,
           -1.000000000000002331643981597126015551422,
           -1.000000000000003297442541400259918073073,
           -1.000000000000004038525841288416879599233,
           -1.000000000000004663287963194255655478615,
           -1.00000000000000521371444217944744506897,
           -1.000000000000005711338016744295885146596,
           -1.000000000000006168950121246466181805002,
           -1.000000000000006594885082800527084897688,
           -1.000000000000006994931944791388919781579,
           -1.000000000000007373305674470655766501228,
           -1.000000000233164398177834299194683872234,
           -1.000000000329744254176269387087995047343,
           -1.00000000040385258418320678088280150237,
           -1.000000000466328796391912356113774800963,
           -1.000000000521371444308553232716574598644,
           -1.00000000057113380178315977436875260197,
           -1.000000000616895012251498501679602643872,
           -1.000000000659488508425026289634430354159,
           -1.000000000699493194642234170768892572882,
           -1.000000000737330567628282553087415117543,
           -1.000023316621036696460620295453277856456,
           -1.000032974787857057404928311684626421503,
           -1.000040385802079313048521250390482335902,
           -1.00004663360452259016530661221213359488,
           -1.0000521380505373899860247162705706318,
           -1.000057114467508637629346787348726350223,
           -1.000061690769779852545970707346938004879,
           -1.000065950300622136103491886720203687457,
           -1.00006995095046930174807034225078340539,
           -1.000073734868993836022551364404248563489,
           -1.007391489031309264813153180819941418531,
           -1.010463846806564696239430620915659099361,
           -1.012825626038880105597738761474228363542,
           -1.014819592594577564927398399257428492725,
           -1.016578512742400177512255407989807698099,
           -1.018170476517182636083407324035097024753,
           -1.019635932177973215702948823070039007361,
           -1.021001233780440980009663527189756124983,
           -1.022284686760270309618528459224732558767,
           -1.023499619082082348038906498637836105447,
           -1.033342436522109918536891072083243897233,
           -1.040939194524944844012076988438386306843,
           -1.047373634492196231421755878017895964061,
           -1.053065496629607111615572634090884988897,
           -1.058230030703619902820337106917657189605,
           -1.06299509412938704964298950857825124135,
           -1.067443986111355120560366087010834293872,
           -1.071634561663924136735470389832541413862,
           -1.07560894118662498941494486924522316597,
           -1.108081880631165502564629660944418191031,
           -1.133487001006868638317076487349933855181,
           -1.155245851821528613609784258821055266176,
           -1.174682608817289477552149783867901714817,
           -1.192475850408615960644596781702366026321,
           -1.209028378276581220769059281765172085749,
           -1.224602449817731587352403997390766826335,
           -1.239380103200799714836392811991400433357,
           -1.253493791367214516100457405907304877145};
  static double wm2[68] = {
           -96.67475603368003636615083422832414231073,
           -95.97433737593292677679699834708774152264,
           -95.56459382507349364043974513871359837914,
           -95.27386489130628866261760496192716897551,
           -95.0483515329550645558378163981346085731,
           -94.86408948075132603599669916724611501067,
           -94.70829516116125928735505687861553800851,
           -94.57333777268864473984718625104978190798,
           -94.45429521137271454134108055166168787618,
           -94.34780665137385269060032461028588648619,
           -73.37311031382297679706747875812087452918,
           -72.67033891766978907253811160121558649554,
           -72.25920015786413889986246462168541066725,
           -71.9674726772681844325410195162521230103,
           -71.74117979478106456261839111929684980709,
           -71.55627755942675731851469544735603279342,
           -71.39993966508440988906136771384270319452,
           -71.26450969134836299738230265916604879247,
           -71.14504894849287026061378869987876478469,
           -71.03818524971357411174259539994036186653,
           -49.96298427667447244531514297262540669957,
           -49.25557728489066973476436802404294348267,
           -48.84167348449764278180827692232244878061,
           -48.54795966722336777861228992424053274064,
           -48.32011181512544381639923088992262632623,
           -48.13392971864323755677656581326964166718,
           -47.97650308277095858785998266172487372203,
           -47.84012504158555569017852884133421231303,
           -47.71982419568730714141619502661365943996,
           -47.61220592218922310708330388890925734186,
           -26.29523881924692569411012882185491823773,
           -25.57429135222126159950976461862116397347,
           -25.15218334705420805339928870686463192335,
           -24.85251554543232259250342343156454440592,
           -24.61997095867949438248843689371454503831,
           -24.42989922074834124324823589226341161866,
           -24.26914663885402405126372567664280388966,
           -24.12985944288624210229238972590881794092,
           -24.00697058168597098928369714836882130512,
           -23.8970195845316574350263109196222825525,
           -14.16360081581018300910955630361089957762,
           -13.41624453595298662833544556875899262976,
           -12.97753279184081358418630625949303360266,
           -12.66551396826200331850774017793451747947,
           -12.42304039760186078066171146072124458063,
           -12.22461776385387453853455424320739669321,
           -12.05663003490708840623665404674007291018,
           -11.91094134143842011964821167497982287763,
           -11.78229922740701885487699061601349928173,
           -11.66711453256635441837882744697047370583,
           -10.90655739570090676132157335673785028979,
           -10.45921112040100393534625826514848865968,
           -10.14059243262036578763968437893562720385,
           -9.892699522704254067620287857665824159861,
           -9.689637966382397752838283301312347921626,
           -9.517569762038614935107630230444563521109,
           -9.368222172408836799233763466046500781388,
           -9.236251966692597369166416348621131600216,
           -9.11800647040274012125833718204681427427,
           -8.335081377982507150789361715143483020265,
           -7.872521380098708883395239767904984410792,
           -7.541940416432904084217222998374802941589,
           -7.283997135099081646930521042317118095276,
           -7.072162048994701667487346245044653243434,
           -6.892241486671583156187212318718730022068,
           -6.735741661607793269808533725369490789074,
           -6.597171733627119347342347717832724288261,
           -6.472775124394004694741057892724488037104};
  static double wp1[68] = {
           -.9999999999999999999766835601840287579665,
           -.9999999999999999999670255745859974370634,
           -.9999999999999999999596147415871158855702,
           -.9999999999999999999533671203680575159335,
           -.9999999999999999999478628555782056161596,
           -.9999999999999999999428866198325571498809,
           -.9999999999999999999383104987875354650364,
           -.9999999999999999999340511491719948741275,
           -.9999999999999999999300506805520862739007,
           -.9999999999999999999262669432552936235556,
           -.9999999999999976683560184028776088243496,
           -.9999999999999967025574585997473306784697,
           -.9999999999999959614741587115939935280812,
           -.9999999999999953367120368057588420244704,
           -.9999999999999947862855578205706768098859,
           -.9999999999999942886619832557258611080314,
           -.9999999999999938310498787535591888253966,
           -.999999999999993405114917199501910108482,
           -.9999999999999930050680552086436996003626,
           -.9999999999999926266943255293804772564852,
           -.9999999997668356018584094585181033975713,
           -.9999999996702557458962181283375794922681,
           -.9999999995961474159255244922555603070484,
           -.9999999995336712037530626747373742782638,
           -.9999999994786285558726655558473617503914,
           -.9999999994288661984343027719079710168757,
           -.9999999993831049880022078023099082447845,
           -.9999999993405114918649237720678678043098,
           -.9999999993005068056839596486461928553989,
           -.9999999992626694327341550240404575805522,
           -.9999766837414008807143234266407434345965,
           -.9999670259370180970391011806287847685011,
           -.9999596152852334187587360603177913882165,
           -.9999533678452277190993205651165793258207,
           -.9999478637616504968301143759542621752943,
           -.9999428877071168268324462680980110604599,
           -.999938311767283189655618359697592754013,
           -.999934052598878483932035240348113191731,
           -.9999300523114688962155368933174163936848,
           -.9999262687553819399632780281900349826094,
           -.9926447551971221136721993073029112268763,
           -.9896086425917686478635208903220735023288,
           -.9872831094708759013315476674998231771112,
           -.9853253899681719161468126266199947992874,
           -.9836027178149637071691226667555243369797,
           -.9820470029764667038452666345865058694192,
           -.9806177971936827573257045283891513709368,
           -.97928874641099293421931043027104578327,
           -.9780415451927629881943028498821429186059,
           -.9768628655744219140604871252425961901255,
           -.967382626983074241885253344632666448927,
           -.9601485420712594199373375259293860324633,
           -.9540768694875733222057908617314370634111,
           -.9487478690765183543410579996573536771348,
           -.9439462911219380338176477772712853402016,
           -.9395442782590063946376623684916441100361,
           -.9354585313336439336066341889767099608018,
           -.9316311953818583253420613278986500351794,
           -.9280201500545670487600430252549212247489,
           -.8991857663963733198571950343133631348347,
           -.8774287170665477623395641312875506084256,
           -.8593275036837387237312746018678000939451,
           -.8435580020488052057849697812109882706542,
           -.8294416857114015557682843481727604063558,
           -.8165758053803078481644781849709953302847,
           -.8046981564792468915744561969751509934994,
           -.7936267540949175059651534957734689407879,
           -.7832291989812967764330746819464532478021};
  static double wp2[40] = {
           9.999999990000000014999999973333333385417E-10, 
           1.9999999960000000119999999573333335E-9, 
           2.999999991000000040499999784000001265625E-9, 
           3.999999984000000095999999317333338666667E-9, 
           4.999999975000000187499998333333349609375E-9, 
           5.999999964000000323999996544000040499999E-9, 
           6.99999995100000051449999359733342086979E-9, 
           7.999999936000000767999989077333503999997E-9, 
           8.999999919000001093499982504000307546869E-9, 
           9.999999900000001499999973333333854166656E-9,
           9.901473843595011885336326816570107953628E-3, 
           1.961158933740562729168248268298370977812E-2, 
           2.913845916787001265458568152535395243296E-2, 
           3.848966594197856933287598180923987047561E-2, 
           4.767230860012937472638890051416087074706E-2, 
           5.669304377414432493107872588796066666126E-2, 
           6.555812274442272075701853672870305774479E-2, 
           7.427342455278083997072135190143718509109E-2, 
           8.284448574644162210327285639805993759142E-2, 
           9.127652716086226429989572142317956865312E-2,
           -1.000000001000000001500000002666666671875E-9, 
           -2.000000004000000012000000042666666833333E-9, 
           -3.000000009000000040500000216000001265625E-9, 
           -4.000000016000000096000000682666672E-9, 
           -5.000000025000000187500001666666682942709E-9, 
           -6.000000036000000324000003456000040500001E-9, 
           -7.000000049000000514500006402666754203126E-9, 
           -8.000000064000000768000010922666837333336E-9, 
           -9.000000081000001093500017496000307546881E-9, 
           -1.000000010000000150000002666666718750001E-8,
           -1.010152719853875327292018767138623973671E-2, 
           -2.041244405580766725973605390749548004159E-2, 
           -3.094279498284817939791038065611524917276E-2, 
           -4.170340843648447389872733812553976786256E-2, 
           -5.270598355154634795995650617915721289428E-2, 
           -6.396318935617251019529498180168867456393E-2, 
           -7.548877886579220591933915955796681153525E-2, 
           -8.729772086157992404091975866027313992649E-2, 
           -9.940635280454481474353567186786621057821E-2, 
           -1.118325591589629648335694568202658422726E-1};
  static double wp3[10] = {
           1.745528002740699383074301264875389911535,
           3.385630140290050184888244364529726867492,
           5.249602852401596227126056319697306282521,
           7.231846038093372706475618500141253883968,
           9.284571428622108983205132234759581939317,
           11.38335808614005262200015678158500428903,
           13.5143440103060912090067238511621580283,
           15.66899671545096218719628189389457073619,
           17.84172596742146918254060066535711011039,
           20.02868541330495078123430607181488729749};
  static double x2[20] = {
    1.E-9,2.E-9,3.E-9,4.E-9,5.E-9,6.E-9,7.E-9,8.E-9,9.E-9,1.E-8, 
    1.E-2,2.E-2,3.E-2,4.E-2,5.E-2,6.E-2,7.E-2,8.E-2,9.E-2,1.E-1};
  static double x3[10] = {
    1.E1,1.E2,1.E3,1.E4,1.E5,1.E6,1.E7,1.E8,1.E9,1.E10};
/*
  Compare the approximations of WAPR with the given exact values.
*/
  printf ( "\n" );
  printf ( "TEST01\n" );
  printf ( "  Compare WAPR(X) to stored values.\n" );
/*
  Wp results for x near -exp(-1).
*/
  printf ( "\n" );
  printf ( "  Wp results for x near -exp(-1)\n" );
  printf ( "\n" );
  printf ( "   Offset x    W(x) (WAPR)" );
  printf ( "   W(x) (EXACT)   Digits Correct\n" );
  printf ( "\n" );

  for ( i = 0; i < 68; i++ )
  {
    if ( dx1[i] == 0.0 )
    {
      em = - exp ( -1.0 );
      w = wapr ( em, 0, &nerror, 0 );
    }
    else
    {
      w = wapr ( dx1[i], 0, &nerror, 1 );
    }

    if ( w == wp1[i] )
    {
      nd = ( int ) ( log10 ( pow ( 2.0, nbits ) ) + 0.5 );
    }
    else
    {
      nd = ( int ) ( log10 ( fabs ( wp1[i] / ( w - wp1[i] ))) + 0.5 );
    }
    printf ( "%17.8e%17.8e%17.8e      %3d\n", dx1[i], w, wp1[i], nd );
  }
/*
  Wp results for x near 0.
*/
  printf ( "\n" );
  printf ( "  Wp results for x near 0\n" );
  printf ( "\n" );
  printf ( "      x    W(x) (WAPR)" );
  printf ( "     W(x) (EXACT)   Digits Correct\n" );
  printf ( "\n" );

  for ( i = 0; i < 20; i++ )
  {
    w = wapr ( x2[i], 0, &nerror, 0 );
    if ( w == wp2[i] )
    {
      nd = ( int ) ( log10 ( pow ( 2.0, nbits ) ) + 0.5 );
    }
    else
    {
      nd = ( int ) ( log10 ( fabs ( wp2[i] / ( w - wp2[i] ))) + 0.5 );
    }
    printf ( "%17.8e%17.8e%17.8e      %3d\n", x2[i], w, wp2[i], nd );
  }

  for ( i = 0; i < 20; i++ )
  {
    w = wapr ( -x2[i], 0, &nerror, 0 );
    if ( w == wp2[20+i] )
    {
      nd = ( int ) ( log10 ( pow ( 2.0, nbits ) ) + 0.5 );
    }
    else
    {
      nd = ( int ) ( log10 ( fabs ( wp2[20+i] / ( w - wp2[20+i] ) ) ) + 0.5 );
    }
    printf ( "%17.8e%17.8e%17.8e      %3d\n", -x2[i], w, wp2[20+i], nd );
  }
/*
  Other Wp results.
*/
  printf ( "\n" );
  printf ( "  Other Wp results\n" );
  printf ( "\n" );
  printf ( "      x     W(x) (WAPR)" );
  printf ( "     W(x) (EXACT)   Digits Correct\n" );
  printf ( "\n" );

  for ( i = 0; i < 10; i++ )
  {
    w = wapr ( x3[i], 0, &nerror, 0 );
    if ( w == wp3[i] )
    {
      nd = ( int ) ( log10 ( pow ( 2.0, nbits ) ) + 0.5 );
    }
    else
    {
      nd = ( int ) ( log10 ( fabs ( wp3[i] / ( w - wp3[i] ))) + 0.5 );
    }
    printf ( "%17.8e%17.8e%17.8e      %3d\n", x3[i], w, wp3[i], nd );
  }
/*
  Wm results for x near -exp(-1).
*/
  printf ( "\n" );
  printf ( "  Wm results for x near 0\n" );
  printf ( "\n" );
  printf ( "      x    W(x) (WAPR)" );
  printf ( "     W(x) (EXACT)   Digits Correct\n" );
  printf ( "\n" );

  for ( i = 0; i < 68; i++ )
  {
    w = wapr ( dx1[i], 1, &nerror, 1 );
    if ( w == wm1[i] )
    {
      nd = ( int ) ( log10 ( pow ( 2.0, nbits ) ) + 0.5 );
    }
    else
    {
      nd = ( int ) ( log10 ( fabs ( wm1[i] / ( w - wm1[i] ))) + 0.5 );
    }
    printf ( "%17.8e%17.8e%17.8e      %3d\n", dx1[i], w, wm1[i], nd );
  }
/*
  Wm results for x near 0.
  Check for underflow.
*/
  printf ( "\n" );
  printf ( "  Wm results for x near -exp(-1)\n" );
  printf ( "\n" );
  printf ( "   Offset x    W(x) (WAPR)" );
  printf ( "     W(x) (EXACT)   Digits Correct\n" );
  printf ( "\n" );

  for ( i = 0; i < 68; i++ )
  {
    if ( 0.0 < dx1[i] )
    {
      w = wapr ( -dx1[i], 1, &nerror, 0 );
      if ( w == wm2[i] )
      {
        nd = ( int ) ( log10 ( pow ( 2.0, nbits ) ) + 0.5 );
      }
      else
      {
        nd = ( int ) ( log10 ( fabs ( wm2[i] / ( w - wm2[i] ))) + 0.5 );
      }
      printf ( "%17.8e%17.8e%17.8e      %3d\n", -dx1[i], w, wm2[i], nd );
    }
  }
  return;
}
/******************************************************************************/

void test02 ( int nbits, double dx, int n )

/******************************************************************************/
/*
  Purpose:

    TEST02 tests WAPR(X) when X is the offset of the argument from -exp(-1).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    17 June 2014

  Author:

    Original FORTRAN77 version by Andrew Barry, S. J. Barry, 
    Patricia Culligan-Hensley.
    This C version by John Burkardt.

  Parameters:

    Input, int NBITS, the number of bits in the mantissa.

    Input, double DX, the initial offset.

    Input, int N, the number of offset arguments to generate.
*/
{
  int i;
  int ifmt;
  int iw;
  int l;
  int nd;
  int ner;
  int nerror;
  double w;
  double we;
  double x;
  double xmax;
  double xmin;

  printf ( "\n" );
  printf ( "TEST02\n" );
  printf ( "  Input X is the offset from -exp(-1).\n" );

  l = 1;
  ifmt = 0;
  xmax = n * dx - exp ( -1.0 );
  xmin = 0.0;

  if ( xmax <= 0.0 )
  {
    iw = 1;
    printf ( "  Both branches of the W function will be checked.\n" );
  }
  else
  {
    iw = 0;
    printf ( "  Wp has been selected (maximum x is > 0)\n" );
  }

  printf ( "\n" );
  printf ( "  Results for Wp(x):\n" );

  if ( ifmt == 0 )
  {
    printf ( "\n" );
    printf ( "   Offset x    W(x) (WAPR)" );
    printf ( "     W(x) (BISECT)  Digits Correct\n" );
  }
  else
  {
    printf ( "\n" );
    printf ( "     x     W(x) (WAPR)" );
    printf ( "     W(x) (BISECT)  Digits Correct\n" );
  }

  printf ( "\n" );

  for ( i = 1; i <= n + 1; i++ )
  {
    x = xmin + i * dx;
    w = wapr ( x, 0, &nerror, l );

    if ( nerror == 1 )
    {
      printf ( "  The value of X = %g is out of range.\n", x );
    }
    else
    {
      we = bisect ( x, 0, &ner, l );

      if ( ner == 1 )
      {
        printf ( "\n" );
        printf ( " BISECT did not converge for x = %g\n", x );
        printf ( "  Try reducing NBITS.\n" );
      }

      if ( w == we )
      {
        nd = ( int ) ( log10 ( pow ( 2.0, nbits ) ) + 0.5 );
      }
      else
      {
        nd = ( int ) ( log10 ( fabs ( we / ( w - we ) ) ) + 0.5 );
      }
      printf ( "%17.8e%17.8e%17.8e      %3d\n", x, w, we, nd );
    }
  }

  if ( iw == 1 )
  {
    printf ( "\n" );
    printf ( "  Results for Wm(x):\n" );

    if ( ifmt == 0 )
    {
      printf ( "\n" );
      printf ( "   Offset x    W(x) (WAPR)" );
      printf ( "     W(x) (BISECT)  Digits Correct\n" );
    }
    else
    {
      printf ( "\n" );
      printf ( "     x     W(x) (WAPR)" );
      printf ( "     W(x) (BISECT)  Digits Correct\n" );
    }

    printf ( "\n" );

    for ( i = 1; i <= n + 1; i++ )
    {
      x = xmin + i * dx;
      w = wapr ( x, 1, &nerror, l );

      if ( nerror == 1 )
      {
        printf ( "  The value of X = %g is out of range.\n", x );
      }
      else
      {
        we = bisect ( x, 1, &ner, l );

        if ( ner == 1 )
        {
          printf ( "\n" );
          printf ( " BISECT did not converge for x = %g\n", x );
          printf ( "  Try reducing NBITS.\n" );
        }
  
        if ( w == we )
        {
          nd = ( int ) ( log10 ( pow ( 2.0, nbits ) ) + 0.5 );
        }
        else
        {
          nd = ( int ) ( log10 ( fabs ( we / ( w - we ) ) ) + 0.5 );
        }
        printf ( "%17.8e%17.8e%17.8e      %3d\n", x, w, we, nd );
      }
    }
  }

  return;
}
/******************************************************************************/

void test03 ( int nbits, double xmin, double xmax, int n )

/******************************************************************************/
/*
  Purpose:

    TEST03 tests WAPR(X) when X is the argument.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    16 June 2014

  Author:

    Original FORTRAN77 version by Andrew Barry, S. J. Barry, 
    Patricia Culligan-Hensley.
    This C version by John Burkardt.

  Parameters:

    Input, int NBITS, the number of bits in the mantissa.

    Input, double XMIN, XMAX, the range.

    Input, int N, the number of equally spaced values
    in the range at which arguments are to be chosen.
*/
{
  double dx;
  int i;
  int ifmt;
  int iw;
  int l;
  int nd;
  int ner;
  int nerror;
  double temp;
  double w;
  double we;
  double x;

  printf ( "\n" );
  printf ( "TEST03\n" );
  printf ( "  Input X is the argument.\n" );

  l = 0;
  ifmt = 1;

  if ( xmax < xmin )
  {
    temp = xmin;
    xmin = xmax;
    xmax = temp;
  }

  dx = ( xmax - xmin ) / ( double ) ( n );
  xmin = xmin - dx;

  if ( xmax <= 0.0 )
  {
    iw = 1;
    printf ( "  Both branches of the W function will be checked.\n" );
  }
  else
  {
    iw = 0;
    printf ( "  Wp has been selected (maximum x is > 0)\n" );
  }

  printf ( "\n" );
  printf ( "  Results for Wp(x):\n" );

  if ( ifmt == 0 )
  {
    printf ( "\n" );
    printf ( "   Offset x    W(x) (WAPR)" );
    printf ( "     W(x) (BISECT)  Digits Correct\n" );
  }
  else
  {
    printf ( "\n" );
    printf ( "     x     W(x) (WAPR)" );
    printf ( "     W(x) (BISECT)  Digits Correct\n" );
  }
  printf ( "\n" );

  for ( i = 1; i <= n + 1; i++ )
  {
    x = xmin + i * dx;
    w = wapr ( x, 0, &nerror, l );

    if ( nerror == 1 )
    {
      printf ( "  The value of X = %g is out of range.\n", x );
    }
    else
    {
      we = bisect ( x, 0, &ner, l );

      if ( ner == 1 )
      {
        printf ( "\n" );
        printf ( " BISECT did not converge for x = %g\n", x );
        printf ( "  Try reducing NBITS.\n" );
      }

      if ( w == we )
      {
        nd = ( int ) ( log10 ( pow ( 2.0, nbits ) ) + 0.5 );
      }
      else
      {
        nd = ( int ) ( log10 ( fabs ( we / ( w - we ))) + 0.5 );
      }
      printf ( "%17.8e%17.8e%17.8e      %3d\n", x, w, we, nd );
    }
  }

  if ( iw == 1 )
  {
    printf ( "\n" );
    printf ( "  Results for Wm(x):\n" );

    if ( ifmt == 0 )
    {
      printf ( "\n" );
      printf ( "   Offset x    W(x) (WAPR)" );
      printf ( "     W(x) (BISECT)  Digits Correct\n" );
    }
    else
    {
      printf ( "\n" );
      printf ( "     x     W(x) (WAPR)" );
      printf ( "     W(x) (BISECT)  Digits Correct\n" );
    }

    for ( i = 1; i <= n + 1; i++ )
    {
      x = xmin + i * dx;
      w = wapr ( x, 1, &nerror, l );

      if ( nerror == 1 )
      {
        printf ( "  The value of X = %g is out of range.\n", x );
      }
      else
      {
        we = bisect ( x, 1, &ner, l );

        if ( ner == 1 )
        {
          printf ( "\n" );
          printf ( " BISECT did not converge for x = %g\n", x );
          printf ( "  Try reducing NBITS.\n" );
        }

        if ( w == we )
        {
          nd = ( int ) ( log10 ( pow ( 2.0, nbits ) ) + 0.5 );
        }
        else
        {
          nd = ( int ) ( log10 ( fabs ( we / ( w - we ))) + 0.5 );
        }
        printf ( "%17.8e%17.8e%17.8e      %3d\n", x, w, we, nd );
      }
    }
  }
  return;
}