Listing 3: Function XVerFun calculates verification
statistics

  1: /* ============ */
  2: /* xverfun.h    */
  3: /* ============ */
  4: // After Cody & Waite and Plauger
  5: #include <stdio.h>
  6: #include <stdlib.h>
  7: //#include <math.h>
  8: #include <time.h>
  9: #include "qfloat.h"
 10: #include "machars.h"
 11: #include "xtstdefs.h"
 12: #include "getmachx.h"
 13: /* ===================================================== */
 14: /* XVerFun - Compares function in standard arithmetic    */
 14: /* with extended                                         */
 15: /* ===================================================== */
 16: void
 17: XVerFun(MACHAR_STRU &MachData,
 18:     const int NumCpr,
 19:     const qfloat LoLim,
 20:     const qfloat HiLim,
 21:     qfloat (*LocalFun)(qfloat, qfloat &),
 22:     char *Label)
 23: {
 24:     int     k, NumHi = 0, NumLo = 0;
 25:     qfloat  XArg, XDelta, XHi, XLo, XRoot, XSqr;
 26:     qfloat  XArgMax, XErrMax, XErrSqr, XRelErr;
 27:
 28:     printf("%s\n", Label);
 29:     //xprintf("", LoLim, 40, "\nto\n");
 30:     //xprintf("", HiLim, 40, "\n");
 31:     fflush(NULL);
 32:     /* ----------------------------------------- */
 33:     /* Calculate Step for Arguments to Be Tested */
 34:     /* ----------------------------------------- */
 35:
 36:     XDelta = (HiLim - LoLim);
 37:     XDelta /= NumCpr;
 38:     //xprintf("XDelta = ", XDelta, 20, "\n");
 39:
 40:     XArgMax = 0;
 41:     XErrMax = 0;
 42:     XErrSqr = 0;
 43:
 44:     //qsrand((unsigned int)time(NULL));
 45:
 46:     for (k = 0; k < NumCpr; ++k)
 47:     {
 48:     qfloat    XAns1, XAns2, XArg;
 49:
 50:     /* Calculate Next Argument */
 51:
 52:     XArg = XDelta * xrand() +  XDelta * k + LoLim;
 53:
 54: //    printf("\nCall # %d\n", k+1);
 55: //    xprintf("Arg                 = ", XArg,  20, "\n");
 56:     XAns1 = (*LocalFun)(XArg, XAns2);
 57:
 58: //    xprintf("Extended Function   = ", XAns2, 20, "\n");
 59: //    xprintf("Function Under Test = ", XAns1, 20, "\n");
 60:
 61:     if (XAns1 == 0)
 62:     {
 63:         XRelErr = 1;
 64:     }
 65:     else
 66:     {
 67:         XRelErr = (XAns1 - XAns2) / XAns1;
 68:     }
 69:
 70:     if (XRelErr < 0)
 71:     {
 72:         XRelErr = -XRelErr;
 73:     }
 74:     if (XAns1 >= 0)
 75:     {
 76:         XAns1 = XRelErr;
 77:     }
 78:     else
 79:     {
 80:         XAns1 = -XRelErr;
 81:     }
 82:
 83:     if (XAns1 < 0)
 84:         ++NumLo;
 85:
 86:     else if (XAns1 > 0)
 87:         ++NumHi;
 88:
 89:     if (XErrMax < XRelErr)
 90:         XErrMax = XRelErr, XArgMax = (TYPE)(xtold(XArg));
 91:
 92:     XErrSqr += XRelErr * XRelErr;
 93:
 94:     if (((k+1) % 100) == 0 && isatty(fileno(stderr)))
 95:     {
 96:         fprintf(stderr, "\tArgument # %4d\r", k+1);
 97:     }
 98:     }
 99:
100:     if (isatty(fileno(stderr)))
101:     {
102:     fprintf(stderr, "\n");
103:     fflush(NULL);
104:     }
105:
106:     printf("Result was smaller %d times, equal %d times, "
107:        "and larger %d times\n",
108:         NumLo, NumCpr - NumLo - NumHi, NumHi);
109:
110:     if (XErrSqr == 0)
111:     {
112:     printf("The maximum relative error = %d ^ -(INF)\n",
113:         MachData.Radix);
114:     printf("Estimated loss of base %d significant digits "
115:         " is 0.00\n", MachData.Radix);
116:     printf("The root-mean-square relative is %d ^ -(INF)\n",
117:         MachData.Radix);
118:     printf("Estimated loss of base %d significant digits "
119:         " is 0.00\n", MachData.Radix);
120:     }
121:     else
122:     {
123:     LDBL    ErrDigs, ErrMax, LossDigs, MeanSqrErr;
124:
125:     ErrDigs = xtold(xlog(XErrMax) /
126:           xlog(MachData.Radix));
127:     ErrMax = xtold(XErrMax);
128:
129:     printf("The maximum relative error of %.10LE"
130:         " = 2 ^ %.2Lf\n\toccurred for x = %.20LG\n",
131:         ErrMax, ErrDigs, xtold(XArgMax));
132:
133:     LossDigs = ErrDigs + MachData.FracDigs;
134:
135:     if (LossDigs <= 0)
136:     {
137:         LossDigs = 0;
138:     }
139:     printf("Estimated loss of base %d significant digits "
140:         " is %.2Lf\n", MachData.Radix, LossDigs);
141:
142:     MeanSqrErr = xtold(xsqrt(XErrSqr)) / (LDBL)NumCpr;
143:     ErrDigs = xtold(xlog(MeanSqrErr) /
144:           xlog(MachData.Radix));
145:
146:     printf("The root-mean-square relative error was "
147:         "%.10LG = %d ^ %.2Lf\n",
148:         MeanSqrErr, MachData.Radix, ErrDigs);
149:
150:     LossDigs = ErrDigs + MachData.FracDigs;
151:
152:     if (LossDigs <= 0)
153:     {
154:         LossDigs = 0;
155:     }
156:     printf("Estimated loss of base %d significant digits "
157:         " is %.2Lf\n", MachData.Radix, LossDigs);
158:     }
159:     printf("\n");
160: }
/* End of File */