home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 2 / 2334 < prev    next >
Internet Message Format  |  1990-12-28  |  56KB
  1. From: glenn@qed.physics.su.OZ.AU (Glenn Geers)
  2. Newsgroups: alt.sources
  3. Subject: Alternative 80386 math lib v2.0 part05/06
  4. Message-ID: <1990Dec16.210740.28604@metro.ucc.su.OZ.AU>
  5. Date: 16 Dec 90 21:07:40 GMT
  6.  
  7. Submitted-by: root@trantor
  8. Archive-name: mathlib2.0/part05
  9.  
  10. ---- Cut Here and feed the following to sh ----
  11. #!/bin/sh
  12. # this is mathlib.05 (part 5 of mathlib2.0)
  13. # do not concatenate these parts, unpack them in order with /bin/sh
  14. # file paranoia.c continued
  15. #
  16. if test ! -r _shar_seq_.tmp; then
  17.     echo 'Please unpack part 1 first!'
  18.     exit 1
  19. fi
  20. (read Scheck
  21.  if test "$Scheck" != 5; then
  22.     echo Please unpack part "$Scheck" next!
  23.     exit 1
  24.  else
  25.     exit 0
  26.  fi
  27. ) < _shar_seq_.tmp || exit 1
  28. if test ! -f _shar_wnt_.tmp; then
  29.     echo 'x - still skipping paranoia.c'
  30. else
  31. echo 'x - continuing file paranoia.c'
  32. sed 's/^X//' << 'SHAR_EOF' >> 'paranoia.c' &&
  33. X    X = OneAndHalf - U2;
  34. X    Y = OneAndHalf + U2;
  35. X    Z = (X - U2) * Y2;
  36. X    T = Y * Y1;
  37. X    Z = Z - X;
  38. X    T = T - X;
  39. X    X = X * Y2;
  40. X    Y = (Y + U2) * Y1;
  41. X    X = X - OneAndHalf;
  42. X    Y = Y - OneAndHalf;
  43. X    if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) {
  44. X        X = (OneAndHalf + U2) * Y2;
  45. X        Y = OneAndHalf - U2 - U2;
  46. X        Z = OneAndHalf + U2 + U2;
  47. X        T = (OneAndHalf - U2) * Y1;
  48. X        X = X - (Z + U2);
  49. X        StickyBit = Y * Y1;
  50. X        S = Z * Y2;
  51. X        T = T - Y;
  52. X        Y = (U2 - Y) + StickyBit;
  53. X        Z = S - (Z + U2 + U2);
  54. X        StickyBit = (Y2 + U2) * Y1;
  55. X        Y1 = Y2 * Y1;
  56. X        StickyBit = StickyBit - Y2;
  57. X        Y1 = Y1 - Half;
  58. X        if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
  59. X            && ( StickyBit == Zero) && (Y1 == Half)) {
  60. X            RMult = Rounded;
  61. X            printf("Multiplication appears to round correctly.\n");
  62. X            }
  63. X        else    if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
  64. X                && (T < Zero) && (StickyBit + U2 == Zero)
  65. X                && (Y1 < Half)) {
  66. X                RMult = Chopped;
  67. X                printf("Multiplication appears to chop.\n");
  68. X                }
  69. X            else printf("* is neither chopped nor correctly rounded.\n");
  70. X        if ((RMult == Rounded) && (GMult == No)) notify("Multiplication");
  71. X        }
  72. X    else printf("* is neither chopped nor correctly rounded.\n");
  73. X    /*=============================================*/
  74. X    Milestone = 45;
  75. X    /*=============================================*/
  76. X    Y2 = One + U2;
  77. X    Y1 = One - U2;
  78. X    Z = OneAndHalf + U2 + U2;
  79. X    X = Z / Y2;
  80. X    T = OneAndHalf - U2 - U2;
  81. X    Y = (T - U2) / Y1;
  82. X    Z = (Z + U2) / Y2;
  83. X    X = X - OneAndHalf;
  84. X    Y = Y - T;
  85. X    T = T / Y1;
  86. X    Z = Z - (OneAndHalf + U2);
  87. X    T = (U2 - OneAndHalf) + T;
  88. X    if (! ((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) {
  89. X        X = OneAndHalf / Y2;
  90. X        Y = OneAndHalf - U2;
  91. X        Z = OneAndHalf + U2;
  92. X        X = X - Y;
  93. X        T = OneAndHalf / Y1;
  94. X        Y = Y / Y1;
  95. X        T = T - (Z + U2);
  96. X        Y = Y - Z;
  97. X        Z = Z / Y2;
  98. X        Y1 = (Y2 + U2) / Y2;
  99. X        Z = Z - OneAndHalf;
  100. X        Y2 = Y1 - Y2;
  101. X        Y1 = (F9 - U1) / F9;
  102. X        if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
  103. X            && (Y2 == Zero) && (Y2 == Zero)
  104. X            && (Y1 - Half == F9 - Half )) {
  105. X            RDiv = Rounded;
  106. X            printf("Division appears to round correctly.\n");
  107. X            if (GDiv == No) notify("Division");
  108. X            }
  109. X        else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
  110. X            && (Y2 < Zero) && (Y1 - Half < F9 - Half)) {
  111. X            RDiv = Chopped;
  112. X            printf("Division appears to chop.\n");
  113. X            }
  114. X        }
  115. X    if (RDiv == Other) printf("/ is neither chopped nor correctly rounded.\n");
  116. X    BInvrse = One / Radix;
  117. X    TstCond (Failure, (BInvrse * Radix - Half == Half),
  118. X           "Radix * ( 1 / Radix ) differs from 1");
  119. X    /*=============================================*/
  120. X    /*SPLIT
  121. X    }
  122. #include "paranoia.h"
  123. part4(){
  124. */
  125. X    Milestone = 50;
  126. X    /*=============================================*/
  127. X    TstCond (Failure, ((F9 + U1) - Half == Half)
  128. X           && ((BMinusU2 + U2 ) - One == Radix - One),
  129. X           "Incomplete carry-propagation in Addition");
  130. X    X = One - U1 * U1;
  131. X    Y = One + U2 * (One - U2);
  132. X    Z = F9 - Half;
  133. X    X = (X - Half) - Z;
  134. X    Y = Y - One;
  135. X    if ((X == Zero) && (Y == Zero)) {
  136. X        RAddSub = Chopped;
  137. X        printf("Add/Subtract appears to be chopped.\n");
  138. X        }
  139. X    if (GAddSub == Yes) {
  140. X        X = (Half + U2) * U2;
  141. X        Y = (Half - U2) * U2;
  142. X        X = One + X;
  143. X        Y = One + Y;
  144. X        X = (One + U2) - X;
  145. X        Y = One - Y;
  146. X        if ((X == Zero) && (Y == Zero)) {
  147. X            X = (Half + U2) * U1;
  148. X            Y = (Half - U2) * U1;
  149. X            X = One - X;
  150. X            Y = One - Y;
  151. X            X = F9 - X;
  152. X            Y = One - Y;
  153. X            if ((X == Zero) && (Y == Zero)) {
  154. X                RAddSub = Rounded;
  155. X                printf("Addition/Subtraction appears to round correctly.\n");
  156. X                if (GAddSub == No) notify("Add/Subtract");
  157. X                }
  158. X            else printf("Addition/Subtraction neither rounds nor chops.\n");
  159. X            }
  160. X        else printf("Addition/Subtraction neither rounds nor chops.\n");
  161. X        }
  162. X    else printf("Addition/Subtraction neither rounds nor chops.\n");
  163. X    S = One;
  164. X    X = One + Half * (One + Half);
  165. X    Y = (One + U2) * Half;
  166. X    Z = X - Y;
  167. X    T = Y - X;
  168. X    StickyBit = Z + T;
  169. X    if (StickyBit != Zero) {
  170. X        S = Zero;
  171. X        BadCond(Flaw, "(X - Y) + (Y - X) is non zero!\n");
  172. X        }
  173. X    StickyBit = Zero;
  174. X    if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
  175. X        && (RMult == Rounded) && (RDiv == Rounded)
  176. X        && (RAddSub == Rounded) && (FLOOR(RadixD2) == RadixD2)) {
  177. X        printf("Checking for sticky bit.\n");
  178. X        X = (Half + U1) * U2;
  179. X        Y = Half * U2;
  180. X        Z = One + Y;
  181. X        T = One + X;
  182. X        if ((Z - One <= Zero) && (T - One >= U2)) {
  183. X            Z = T + Y;
  184. X            Y = Z - X;
  185. X            if ((Z - T >= U2) && (Y - T == Zero)) {
  186. X                X = (Half + U1) * U1;
  187. X                Y = Half * U1;
  188. X                Z = One - Y;
  189. X                T = One - X;
  190. X                if ((Z - One == Zero) && (T - F9 == Zero)) {
  191. X                    Z = (Half - U1) * U1;
  192. X                    T = F9 - Z;
  193. X                    Q = F9 - Y;
  194. X                    if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) {
  195. X                        Z = (One + U2) * OneAndHalf;
  196. X                        T = (OneAndHalf + U2) - Z + U2;
  197. X                        X = One + Half / Radix;
  198. X                        Y = One + Radix * U2;
  199. X                        Z = X * Y;
  200. X                        if (T == Zero && X + Radix * U2 - Z == Zero) {
  201. X                            if (Radix != Two) {
  202. X                                X = Two + U2;
  203. X                                Y = X / Two;
  204. X                                if ((Y - One == Zero)) StickyBit = S;
  205. X                                }
  206. X                            else StickyBit = S;
  207. X                            }
  208. X                        }
  209. X                    }
  210. X                }
  211. X            }
  212. X        }
  213. X    if (StickyBit == One) printf("Sticky bit apparently used correctly.\n");
  214. X    else printf("Sticky bit used incorrectly or not at all.\n");
  215. X    TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
  216. X            RMult == Other || RDiv == Other || RAddSub == Other),
  217. X        "lack(s) of guard digits or failure(s) to correctly round or chop\n\
  218. (noted above) count as one flaw in the final tally below");
  219. X    /*=============================================*/
  220. X    Milestone = 60;
  221. X    /*=============================================*/
  222. X    printf("\n");
  223. X    printf("Does Multiplication commute?  ");
  224. X    printf("Testing on %d random pairs.\n", NoTrials);
  225. X    Random9 = SQRT(3.0);
  226. X    Random1 = Third;
  227. X    I = 1;
  228. X    do  {
  229. X        X = Random();
  230. X        Y = Random();
  231. X        Z9 = Y * X;
  232. X        Z = X * Y;
  233. X        Z9 = Z - Z9;
  234. X        I = I + 1;
  235. X        } while ( ! ((I > NoTrials) || (Z9 != Zero)));
  236. X    if (I == NoTrials) {
  237. X        Random1 = One + Half / Three;
  238. X        Random2 = (U2 + U1) + One;
  239. X        Z = Random1 * Random2;
  240. X        Y = Random2 * Random1;
  241. X        Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
  242. X            Three) * ((U2 + U1) + One);
  243. X        }
  244. X    if (! ((I == NoTrials) || (Z9 == Zero)))
  245. X        BadCond(Defect, "X * Y == Y * X trial fails.\n");
  246. X    else printf("     No failures found in %d integer pairs.\n", NoTrials);
  247. X    /*=============================================*/
  248. X    Milestone = 70;
  249. X    /*=============================================*/
  250. X    printf("\nRunning test of square root(x).\n");
  251. X    TstCond (Failure, (Zero == SQRT(Zero))
  252. X           && (- Zero == SQRT(- Zero))
  253. X           && (One == SQRT(One)), "Square root of 0.0, -0.0 or 1.0 wrong");
  254. X    MinSqEr = Zero;
  255. X    MaxSqEr = Zero;
  256. X    J = Zero;
  257. X    X = Radix;
  258. X    OneUlp = U2;
  259. X    SqXMinX (Serious);
  260. X    X = BInvrse;
  261. X    OneUlp = BInvrse * U1;
  262. X    SqXMinX (Serious);
  263. X    X = U1;
  264. X    OneUlp = U1 * U1;
  265. X    SqXMinX (Serious);
  266. X    if (J != Zero) Pause();
  267. X    printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
  268. X    J = Zero;
  269. X    X = Two;
  270. X    Y = Radix;
  271. X    if ((Radix != One)) do  {
  272. X        X = Y;
  273. X        Y = Radix * Y;
  274. X        } while ( ! ((Y - X >= NoTrials)));
  275. X    OneUlp = X * U2;
  276. X    I = 1;
  277. X    while (I <= NoTrials) {
  278. X        X = X + One;
  279. X        SqXMinX (Defect);
  280. X        if (J > Zero) break;
  281. X        I = I + 1;
  282. X        }
  283. X    printf("Test for sqrt monotonicity.\n");
  284. X    I = - 1;
  285. X    X = BMinusU2;
  286. X    Y = Radix;
  287. X    Z = Radix + Radix * U2;
  288. X    NotMonot = False;
  289. X    Monot = False;
  290. X    while ( ! (NotMonot || Monot)) {
  291. X        I = I + 1;
  292. X        X = SQRT(X);
  293. X        Q = SQRT(Y);
  294. X        Z = SQRT(Z);
  295. X        if ((X > Q) || (Q > Z)) NotMonot = True;
  296. X        else {
  297. X            Q = FLOOR(Q + Half);
  298. X            if ((I > 0) || (Radix == Q * Q)) Monot = True;
  299. X            else if (I > 0) {
  300. X            if (I > 1) Monot = True;
  301. X            else {
  302. X                Y = Y * BInvrse;
  303. X                X = Y - U1;
  304. X                Z = Y + U1;
  305. X                }
  306. X            }
  307. X            else {
  308. X                Y = Q;
  309. X                X = Y - U2;
  310. X                Z = Y + U2;
  311. X                }
  312. X            }
  313. X        }
  314. X    if (Monot) printf("sqrt has passed a test for Monotonicity.\n");
  315. X    else {
  316. X        BadCond(Defect, "");
  317. X        printf("sqrt(X) is non-monotonic for X near %.7e .\n", Y);
  318. X        }
  319. X    /*=============================================*/
  320. X    /*SPLIT
  321. X    }
  322. #include "paranoia.h"
  323. part5(){
  324. */
  325. X    Milestone = 80;
  326. X    /*=============================================*/
  327. X    MinSqEr = MinSqEr + Half;
  328. X    MaxSqEr = MaxSqEr - Half;
  329. X    Y = (SQRT(One + U2) - One) / U2;
  330. X    SqEr = (Y - One) + U2 / Eight;
  331. X    if (SqEr > MaxSqEr) MaxSqEr = SqEr;
  332. X    SqEr = Y + U2 / Eight;
  333. X    if (SqEr < MinSqEr) MinSqEr = SqEr;
  334. X    Y = ((SQRT(F9) - U2) - (One - U2)) / U1;
  335. X    SqEr = Y + U1 / Eight;
  336. X    if (SqEr > MaxSqEr) MaxSqEr = SqEr;
  337. X    SqEr = (Y + One) + U1 / Eight;
  338. X    if (SqEr < MinSqEr) MinSqEr = SqEr;
  339. X    OneUlp = U2;
  340. X    X = OneUlp;
  341. X    for( Indx = 1; Indx <= 3; ++Indx) {
  342. X        Y = SQRT((X + U1 + X) + F9);
  343. X        Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;
  344. X        Z = ((U1 - X) + F9) * Half * X * X / OneUlp;
  345. X        SqEr = (Y + Half) + Z;
  346. X        if (SqEr < MinSqEr) MinSqEr = SqEr;
  347. X        SqEr = (Y - Half) + Z;
  348. X        if (SqEr > MaxSqEr) MaxSqEr = SqEr;
  349. X        if (((Indx == 1) || (Indx == 3))) 
  350. X            X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));
  351. X        else {
  352. X            OneUlp = U1;
  353. X            X = - OneUlp;
  354. X            }
  355. X        }
  356. X    /*=============================================*/
  357. X    Milestone = 85;
  358. X    /*=============================================*/
  359. X    SqRWrng = False;
  360. X    Anomaly = False;
  361. X    RSqrt = Other; /* ~dgh */
  362. X    if (Radix != One) {
  363. X        printf("Testing whether sqrt is rounded or chopped.\n");
  364. X        D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));
  365. X    /* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
  366. X        X = D / Radix;
  367. X        Y = D / A1;
  368. X        if ((X != FLOOR(X)) || (Y != FLOOR(Y))) {
  369. X            Anomaly = True;
  370. X            }
  371. X        else {
  372. X            X = Zero;
  373. X            Z2 = X;
  374. X            Y = One;
  375. X            Y2 = Y;
  376. X            Z1 = Radix - One;
  377. X            FourD = Four * D;
  378. X            do  {
  379. X                if (Y2 > Z2) {
  380. X                    Q = Radix;
  381. X                    Y1 = Y;
  382. X                    do  {
  383. X                        X1 = FABS(Q + FLOOR(Half - Q / Y1) * Y1);
  384. X                        Q = Y1;
  385. X                        Y1 = X1;
  386. X                        } while ( ! (X1 <= Zero));
  387. X                    if (Q <= One) {
  388. X                        Z2 = Y2;
  389. X                        Z = Y;
  390. X                        }
  391. X                    }
  392. X                Y = Y + Two;
  393. X                X = X + Eight;
  394. X                Y2 = Y2 + X;
  395. X                if (Y2 >= FourD) Y2 = Y2 - FourD;
  396. X                } while ( ! (Y >= D));
  397. X            X8 = FourD - Z2;
  398. X            Q = (X8 + Z * Z) / FourD;
  399. X            X8 = X8 / Eight;
  400. X            if (Q != FLOOR(Q)) Anomaly = True;
  401. X            else {
  402. X                Break = False;
  403. X                do  {
  404. X                    X = Z1 * Z;
  405. X                    X = X - FLOOR(X / Radix) * Radix;
  406. X                    if (X == One) 
  407. X                        Break = True;
  408. X                    else
  409. X                        Z1 = Z1 - One;
  410. X                    } while ( ! (Break || (Z1 <= Zero)));
  411. X                if ((Z1 <= Zero) && (! Break)) Anomaly = True;
  412. X                else {
  413. X                    if (Z1 > RadixD2) Z1 = Z1 - Radix;
  414. X                    do  {
  415. X                        NewD();
  416. X                        } while ( ! (U2 * D >= F9));
  417. X                    if (D * Radix - D != W - D) Anomaly = True;
  418. X                    else {
  419. X                        Z2 = D;
  420. X                        I = 0;
  421. X                        Y = D + (One + Z) * Half;
  422. X                        X = D + Z + Q;
  423. X                        SR3750();
  424. X                        Y = D + (One - Z) * Half + D;
  425. X                        X = D - Z + D;
  426. X                        X = X + Q + X;
  427. X                        SR3750();
  428. X                        NewD();
  429. X                        if (D - Z2 != W - Z2) Anomaly = True;
  430. X                        else {
  431. X                            Y = (D - Z2) + (Z2 + (One - Z) * Half);
  432. X                            X = (D - Z2) + (Z2 - Z + Q);
  433. X                            SR3750();
  434. X                            Y = (One + Z) * Half;
  435. X                            X = Q;
  436. X                            SR3750();
  437. X                            if (I == 0) Anomaly = True;
  438. X                            }
  439. X                        }
  440. X                    }
  441. X                }
  442. X            }
  443. X        if ((I == 0) || Anomaly) {
  444. X            BadCond(Failure, "Anomalous arithmetic with Integer < ");
  445. X            printf("Radix^Precision = %.7e\n", W);
  446. X            printf(" fails test whether sqrt rounds or chops.\n");
  447. X            SqRWrng = True;
  448. X            }
  449. X        }
  450. X    if (! Anomaly) {
  451. X        if (! ((MinSqEr < Zero) || (MaxSqEr > Zero))) {
  452. X            RSqrt = Rounded;
  453. X            printf("Square root appears to be correctly rounded.\n");
  454. X            }
  455. X        else  {
  456. X            if ((MaxSqEr + U2 > U2 - Half) || (MinSqEr > Half)
  457. X                || (MinSqEr + Radix < Half)) SqRWrng = True;
  458. X            else {
  459. X                RSqrt = Chopped;
  460. X                printf("Square root appears to be chopped.\n");
  461. X                }
  462. X            }
  463. X        }
  464. X    if (SqRWrng) {
  465. X        printf("Square root is neither chopped nor correctly rounded.\n");
  466. X        printf("Observed errors run from %.7e ", MinSqEr - Half);
  467. X        printf("to %.7e ulps.\n", Half + MaxSqEr);
  468. X        TstCond (Serious, MaxSqEr - MinSqEr < Radix * Radix,
  469. X            "sqrt gets too many last digits wrong");
  470. X        }
  471. X    /*=============================================*/
  472. X    Milestone = 90;
  473. X    /*=============================================*/
  474. X    Pause();
  475. X    printf("Testing powers Z^i for small Integers Z and i.\n");
  476. X    N = 0;
  477. X    /* ... test powers of zero. */
  478. X    I = 0;
  479. X    Z = -Zero;
  480. X    M = 3.0;
  481. X    Break = False;
  482. X    do  {
  483. X        X = One;
  484. X        SR3980();
  485. X        if (I <= 10) {
  486. X            I = 1023;
  487. X            SR3980();
  488. X            }
  489. X        if (Z == MinusOne) Break = True;
  490. X        else {
  491. X            Z = MinusOne;
  492. X            PrintIfNPositive();
  493. X            N = 0;
  494. X            /* .. if(-1)^N is invalid, replace MinusOne by One. */
  495. X            I = - 4;
  496. X            }
  497. X        } while ( ! Break);
  498. X    PrintIfNPositive();
  499. X    N1 = N;
  500. X    N = 0;
  501. X    Z = A1;
  502. X    M = FLOOR(Two * LOG(W) / LOG(A1));
  503. X    Break = False;
  504. X    do  {
  505. X        X = Z;
  506. X        I = 1;
  507. X        SR3980();
  508. X        if (Z == AInvrse) Break = True;
  509. X        else Z = AInvrse;
  510. X        } while ( ! (Break));
  511. X    /*=============================================*/
  512. X        Milestone = 100;
  513. X    /*=============================================*/
  514. X    /*  Powers of Radix have been tested, */
  515. X    /*         next try a few primes     */
  516. X    M = NoTrials;
  517. X    Z = Three;
  518. X    do  {
  519. X        X = Z;
  520. X        I = 1;
  521. X        SR3980();
  522. X        do  {
  523. X            Z = Z + Two;
  524. X            } while ( Three * FLOOR(Z / Three) == Z );
  525. X        } while ( Z < Eight * Three );
  526. X    if (N > 0) {
  527. X        printf("Errors like this may invalidate financial calculations\n");
  528. X        printf("\tinvolving interest rates.\n");
  529. X        }
  530. X    PrintIfNPositive();
  531. X    N += N1;
  532. X    if (N == 0) printf("... no discrepancis found.\n");
  533. X    if (N > 0) Pause();
  534. X    else printf("\n");
  535. X    /*=============================================*/
  536. X    /*SPLIT
  537. X    }
  538. #include "paranoia.h"
  539. part6(){
  540. */
  541. X    Milestone = 110;
  542. X    /*=============================================*/
  543. X    printf("Seeking Underflow thresholds UfThold and E0.\n");
  544. X    D = U1;
  545. X    if (Precision != FLOOR(Precision)) {
  546. X        D = BInvrse;
  547. X        X = Precision;
  548. X        do  {
  549. X            D = D * BInvrse;
  550. X            X = X - One;
  551. X            } while ( X > Zero);
  552. X        }
  553. X    Y = One;
  554. X    Z = D;
  555. X    /* ... D is power of 1/Radix < 1. */
  556. X    do  {
  557. X        C = Y;
  558. X        Y = Z;
  559. X        Z = Y * Y;
  560. X        } while ((Y > Z) && (Z + Z > Z));
  561. X    Y = C;
  562. X    Z = Y * D;
  563. X    do  {
  564. X        C = Y;
  565. X        Y = Z;
  566. X        Z = Y * D;
  567. X        } while ((Y > Z) && (Z + Z > Z));
  568. X    if (Radix < Two) HInvrse = Two;
  569. X    else HInvrse = Radix;
  570. X    H = One / HInvrse;
  571. X    /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
  572. X    CInvrse = One / C;
  573. X    E0 = C;
  574. X    Z = E0 * H;
  575. X    /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
  576. X    do  {
  577. X        Y = E0;
  578. X        E0 = Z;
  579. X        Z = E0 * H;
  580. X        } while ((E0 > Z) && (Z + Z > Z));
  581. X    UfThold = E0;
  582. X    E1 = Zero;
  583. X    Q = Zero;
  584. X    E9 = U2;
  585. X    S = One + E9;
  586. X    D = C * S;
  587. X    if (D <= C) {
  588. X        E9 = Radix * U2;
  589. X        S = One + E9;
  590. X        D = C * S;
  591. X        if (D <= C) {
  592. X            BadCond(Failure, "multiplication gets too many last digits wrong.\n");
  593. X            Underflow = E0;
  594. X            Y1 = Zero;
  595. X            PseudoZero = Z;
  596. X            Pause();
  597. X            }
  598. X        }
  599. X    else {
  600. X        Underflow = D;
  601. X        PseudoZero = Underflow * H;
  602. X        UfThold = Zero;
  603. X        do  {
  604. X            Y1 = Underflow;
  605. X            Underflow = PseudoZero;
  606. X            if (E1 + E1 <= E1) {
  607. X                Y2 = Underflow * HInvrse;
  608. X                E1 = FABS(Y1 - Y2);
  609. X                Q = Y1;
  610. X                if ((UfThold == Zero) && (Y1 != Y2)) UfThold = Y1;
  611. X                }
  612. X            PseudoZero = PseudoZero * H;
  613. X            } while ((Underflow > PseudoZero)
  614. X                && (PseudoZero + PseudoZero > PseudoZero));
  615. X        }
  616. X    /* Comment line 4530 .. 4560 */
  617. X    if (PseudoZero != Zero) {
  618. X        printf("\n");
  619. X        Z = PseudoZero;
  620. X    /* ... Test PseudoZero for "phoney- zero" violates */
  621. X    /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
  622. X           ... */
  623. X        if (PseudoZero <= Zero) {
  624. X            BadCond(Failure, "Positive expressions can underflow to an\n");
  625. X            printf("allegedly negative value\n");
  626. X            printf("PseudoZero that prints out as: %g .\n", PseudoZero);
  627. X            X = - PseudoZero;
  628. X            if (X <= Zero) {
  629. X                printf("But -PseudoZero, which should be\n");
  630. X                printf("positive, isn't; it prints out as  %g .\n", X);
  631. X                }
  632. X            }
  633. X        else {
  634. X            BadCond(Flaw, "Underflow can stick at an allegedly positive\n");
  635. X            printf("value PseudoZero that prints out as %g .\n", PseudoZero);
  636. X            }
  637. X        TstPtUf();
  638. X        }
  639. X    /*=============================================*/
  640. X    Milestone = 120;
  641. X    /*=============================================*/
  642. X    if (CInvrse * Y > CInvrse * Y1) {
  643. X        S = H * S;
  644. X        E0 = Underflow;
  645. X        }
  646. X    if (! ((E1 == Zero) || (E1 == E0))) {
  647. X        BadCond(Defect, "");
  648. X        if (E1 < E0) {
  649. X            printf("Products underflow at a higher");
  650. X            printf(" threshold than differences.\n");
  651. X            if (PseudoZero == Zero) 
  652. X            E0 = E1;
  653. X            }
  654. X        else {
  655. X            printf("Difference underflows at a higher");
  656. X            printf(" threshold than products.\n");
  657. X            }
  658. X        }
  659. X    printf("Smallest strictly positive number found is E0 = %g .\n", E0);
  660. X    Z = E0;
  661. X    TstPtUf();
  662. X    Underflow = E0;
  663. X    if (N == 1) Underflow = Y;
  664. X    I = 4;
  665. X    if (E1 == Zero) I = 3;
  666. X    if (UfThold == Zero) I = I - 2;
  667. X    UfNGrad = True;
  668. X    switch (I)  {
  669. X        case    1:
  670. X        UfThold = Underflow;
  671. X        if ((CInvrse * Q) != ((CInvrse * Y) * S)) {
  672. X            UfThold = Y;
  673. X            BadCond(Failure, "Either accuracy deteriorates as numbers\n");
  674. X            printf("approach a threshold = %.17e\n", UfThold);;
  675. X            printf(" coming down from %.17e\n", C);
  676. X            printf(" or else multiplication gets too many last digits wrong.\n");
  677. X            }
  678. X        Pause();
  679. X        break;
  680. X    
  681. X        case    2:
  682. X        BadCond(Failure, "Underflow confuses Comparison, which alleges that\n");
  683. X        printf("Q == Y while denying that |Q - Y| == 0; these values\n");
  684. X        printf("print out as Q = %.17e, Y = %.17e .\n", Q, Y2);
  685. X        printf ("|Q - Y| = %.17e .\n" , FABS(Q - Y2));
  686. X        UfThold = Q;
  687. X        break;
  688. X    
  689. X        case    3:
  690. X        X = X;
  691. X        break;
  692. X    
  693. X        case    4:
  694. X        if ((Q == UfThold) && (E1 == E0)
  695. X            && (FABS( UfThold - E1 / E9) <= E1)) {
  696. X            UfNGrad = False;
  697. X            printf("Underflow is gradual; it incurs Absolute Error =\n");
  698. X            printf("(roundoff in UfThold) < E0.\n");
  699. X            Y = E0 * CInvrse;
  700. X            Y = Y * (OneAndHalf + U2);
  701. X            X = CInvrse * (One + U2);
  702. X            Y = Y / X;
  703. X            IEEE = (Y == E0);
  704. X            }
  705. X        }
  706. X    if (UfNGrad) {
  707. X        printf("\n");
  708. X        sigsave = sigfpe;
  709. X        if (setjmp(ovfl_buf)) {
  710. X            printf("Underflow / UfThold failed!\n");
  711. X            R = H + H;
  712. X            }
  713. X        else R = SQRT(Underflow / UfThold);
  714. X        sigsave = 0;
  715. X        if (R <= H) {
  716. X            Z = R * UfThold;
  717. X            X = Z * (One + R * H * (One + H));
  718. X            }
  719. X        else {
  720. X            Z = UfThold;
  721. X            X = Z * (One + H * H * (One + H));
  722. X            }
  723. X        if (! ((X == Z) || (X - Z != Zero))) {
  724. X            BadCond(Flaw, "");
  725. X            printf("X = %.17e\n\tis not equal to Z = %.17e .\n", X, Z);
  726. X            Z9 = X - Z;
  727. X            printf("yet X - Z yields %.17e .\n", Z9);
  728. X            printf("    Should this NOT signal Underflow, ");
  729. X            printf("this is a SERIOUS DEFECT\nthat causes ");
  730. X            printf("confusion when innocent statements like\n");;
  731. X            printf("    if (X == Z)  ...  else");
  732. X            printf("  ... (f(X) - f(Z)) / (X - Z) ...\n");
  733. X            printf("encounter Division by Zero although actually\n");
  734. X            sigsave = sigfpe;
  735. X            if (setjmp(ovfl_buf)) printf("X / Z fails!\n");
  736. X            else printf("X / Z = 1 + %g .\n", (X / Z - Half) - Half);
  737. X            sigsave = 0;
  738. X            }
  739. X        }
  740. X    printf("The Underflow threshold is %.17e, %s\n", UfThold,
  741. X           " below which");
  742. X    printf("calculation may suffer larger Relative error than ");
  743. X    printf("merely roundoff.\n");
  744. X    Y2 = U1 * U1;
  745. X    Y = Y2 * Y2;
  746. X    Y2 = Y * U1;
  747. X    if (Y2 <= UfThold) {
  748. X        if (Y > E0) {
  749. X            BadCond(Defect, "");
  750. X            I = 5;
  751. X            }
  752. X        else {
  753. X            BadCond(Serious, "");
  754. X            I = 4;
  755. X            }
  756. X        printf("Range is too narrow; U1^%d Underflows.\n", I);
  757. X        }
  758. X    /*=============================================*/
  759. X    /*SPLIT
  760. X    }
  761. #include "paranoia.h"
  762. part7(){
  763. */
  764. X    Milestone = 130;
  765. X    /*=============================================*/
  766. X    Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;
  767. X    Y2 = Y + Y;
  768. X    printf("Since underflow occurs below the threshold\n");
  769. X    printf("UfThold = (%.17e) ^ (%.17e)\nonly underflow ", HInvrse, Y);
  770. X    printf("should afflict the expression\n\t(%.17e) ^ (%.17e);\n", HInvrse, Y);
  771. X    V9 = POW(HInvrse, Y2);
  772. X    printf("actually calculating yields: %.17e .\n", V9);
  773. X    if (! ((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) {
  774. X        BadCond(Serious, "this is not between 0 and underflow\n");
  775. X        printf("   threshold = %.17e .\n", UfThold);
  776. X        }
  777. X    else if (! (V9 > UfThold * (One + E9)))
  778. X        printf("This computed value is O.K.\n");
  779. X    else {
  780. X        BadCond(Defect, "this is not between 0 and underflow\n");
  781. X        printf("   threshold = %.17e .\n", UfThold);
  782. X        }
  783. X    /*=============================================*/
  784. X    Milestone = 140;
  785. X    /*=============================================*/
  786. X    printf("\n");
  787. X    /* ...calculate Exp2 == exp(2) == 7.389056099... */
  788. X    X = Zero;
  789. X    I = 2;
  790. X    Y = Two * Three;
  791. X    Q = Zero;
  792. X    N = 0;
  793. X    do  {
  794. X        Z = X;
  795. X        I = I + 1;
  796. X        Y = Y / (I + I);
  797. X        R = Y + Q;
  798. X        X = Z + R;
  799. X        Q = (Z - X) + R;
  800. X        } while(X > Z);
  801. X    Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);
  802. X    X = Z * Z;
  803. X    Exp2 = X * X;
  804. X    X = F9;
  805. X    Y = X - U1;
  806. X    printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17e as X -> 1.\n",
  807. X        Exp2);
  808. X    for(I = 1;;) {
  809. X        Z = X - BInvrse;
  810. X        Z = (X + One) / (Z - (One - BInvrse));
  811. X        Q = POW(X, Z) - Exp2;
  812. X        if (FABS(Q) > TwoForty * U2) {
  813. X            N = 1;
  814. X             V9 = (X - BInvrse) - (One - BInvrse);
  815. X            BadCond(Defect, "Calculated");
  816. X            printf(" %.17e for\n", POW(X,Z));
  817. X            printf("\t(1 + (%.17e) ^ (%.17e);\n", V9, Z);
  818. X            printf("\tdiffers from correct value by %.17e .\n", Q);
  819. X            printf("\tThis much error may spoil financial\n");
  820. X            printf("\tcalculations involving tiny interest rates.\n");
  821. X            break;
  822. X            }
  823. X        else {
  824. X            Z = (Y - X) * Two + Y;
  825. X            X = Y;
  826. X            Y = Z;
  827. X            Z = One + (X - F9)*(X - F9);
  828. X            if (Z > One && I < NoTrials) I++;
  829. X            else  {
  830. X                if (X > One) {
  831. X                    if (N == 0)
  832. X                       printf("Accuracy seems adequate.\n");
  833. X                    break;
  834. X                    }
  835. X                else {
  836. X                    X = One + U2;
  837. X                    Y = U2 + U2;
  838. X                    Y += X;
  839. X                    I = 1;
  840. X                    }
  841. X                }
  842. X            }
  843. X        }
  844. X    /*=============================================*/
  845. X    Milestone = 150;
  846. X    /*=============================================*/
  847. X    printf("Testing powers Z^Q at four nearly extreme values.\n");
  848. X    N = 0;
  849. X    Z = A1;
  850. X    Q = FLOOR(Half - LOG(C) / LOG(A1));
  851. X    Break = False;
  852. X    do  {
  853. X        X = CInvrse;
  854. X        Y = POW(Z, Q);
  855. X        IsYeqX();
  856. X        Q = - Q;
  857. X        X = C;
  858. X        Y = POW(Z, Q);
  859. X        IsYeqX();
  860. X        if (Z < One) Break = True;
  861. X        else Z = AInvrse;
  862. X        } while ( ! (Break));
  863. X    PrintIfNPositive();
  864. X    if (N == 0) printf(" ... no discrepancies found.\n");
  865. X    printf("\n");
  866. X    
  867. X    /*=============================================*/
  868. X    Milestone = 160;
  869. X    /*=============================================*/
  870. X    Pause();
  871. X    printf("Searching for Overflow threshold:\n");
  872. X    printf("This may generate an error.\n");
  873. X    Y = - CInvrse;
  874. X    V9 = HInvrse * Y;
  875. X    sigsave = sigfpe;
  876. X    if (setjmp(ovfl_buf)) { I = 0; V9 = Y; goto overflow; }
  877. X    do {
  878. X        V = Y;
  879. X        Y = V9;
  880. X        V9 = HInvrse * Y;
  881. X        } while(V9 < Y);
  882. X    I = 1;
  883. overflow:
  884. X    sigsave = 0;
  885. X    Z = V9;
  886. X    printf("Can `Z = -Y' overflow?\n");
  887. X    printf("Trying it on Y = %.17e .\n", Y);
  888. X    V9 = - Y;
  889. X    V0 = V9;
  890. X    if (V - Y == V + V0) printf("Seems O.K.\n");
  891. X    else {
  892. X        printf("finds a ");
  893. X        BadCond(Flaw, "-(-Y) differs from Y.\n");
  894. X        }
  895. X    if (Z != Y) {
  896. X        BadCond(Serious, "");
  897. X        printf("overflow past %.17e\n\tshrinks to %.17e .\n", Y, Z);
  898. X        }
  899. X    if (I) {
  900. X        Y = V * (HInvrse * U2 - HInvrse);
  901. X        Z = Y + ((One - HInvrse) * U2) * V;
  902. X        if (Z < V0) Y = Z;
  903. X        if (Y < V0) V = Y;
  904. X        if (V0 - V < V0) V = V0;
  905. X        }
  906. X    else {
  907. X        V = Y * (HInvrse * U2 - HInvrse);
  908. X        V = V + ((One - HInvrse) * U2) * Y;
  909. X        }
  910. X    printf("Overflow threshold is V  = %.17e .\n", V);
  911. X    if (I) printf("Overflow saturates at V0 = %.17e .\n", V0);
  912. X    else printf("There is no saturation value because \
  913. the system traps on overflow.\n");
  914. X    V9 = V * One;
  915. X    printf("No Overflow should be signaled for V * 1 = %.17e\n", V9);
  916. X    V9 = V / One;
  917. X    printf("                           nor for V / 1 = %.17e .\n", V9);
  918. X    printf("Any overflow signal separating this * from the one\n");
  919. X    printf("above is a DEFECT.\n");
  920. X    /*=============================================*/
  921. X    Milestone = 170;
  922. X    /*=============================================*/
  923. X    if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V)) {
  924. X        BadCond(Failure, "Comparisons involving ");
  925. X        printf("+-%g, +-%g\nand +-%g are confused by Overflow.",
  926. X            V, V0, UfThold);
  927. X        }
  928. X    /*=============================================*/
  929. X    Milestone = 175;
  930. X    /*=============================================*/
  931. X    printf("\n");
  932. X    for(Indx = 1; Indx <= 3; ++Indx) {
  933. X        switch (Indx)  {
  934. X            case 1: Z = UfThold; break;
  935. X            case 2: Z = E0; break;
  936. X            case 3: Z = PseudoZero; break;
  937. X            }
  938. X        if (Z != Zero) {
  939. X            V9 = SQRT(Z);
  940. X            Y = V9 * V9;
  941. X            if (Y / (One - Radix * E9) < Z
  942. X               || Y > (One + Radix * E9) * Z) { /* dgh: + E9 --> * E9 */
  943. X                if (V9 > U1) BadCond(Serious, "");
  944. X                else BadCond(Defect, "");
  945. X                printf("Comparison alleges that what prints as Z = %.17e\n", Z);
  946. X                printf(" is too far from sqrt(Z) ^ 2 = %.17e .\n", Y);
  947. X                }
  948. X            }
  949. X        }
  950. X    /*=============================================*/
  951. X    Milestone = 180;
  952. X    /*=============================================*/
  953. X    for(Indx = 1; Indx <= 2; ++Indx) {
  954. X        if (Indx == 1) Z = V;
  955. X        else Z = V0;
  956. X        V9 = SQRT(Z);
  957. X        X = (One - Radix * E9) * V9;
  958. X        V9 = V9 * X;
  959. X        if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) {
  960. X            Y = V9;
  961. X            if (X < W) BadCond(Serious, "");
  962. X            else BadCond(Defect, "");
  963. X            printf("Comparison alleges that Z = %17e\n", Z);
  964. X            printf(" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y);
  965. X            }
  966. X        }
  967. X    /*=============================================*/
  968. X    /*SPLIT
  969. X    }
  970. #include "paranoia.h"
  971. part8(){
  972. */
  973. X    Milestone = 190;
  974. X    /*=============================================*/
  975. X    Pause();
  976. X    X = UfThold * V;
  977. X    Y = Radix * Radix;
  978. X    if (X*Y < One || X > Y) {
  979. X        if (X * Y < U1 || X > Y/U1) BadCond(Defect, "Badly");
  980. X        else BadCond(Flaw, "");
  981. X            
  982. X        printf(" unbalanced range; UfThold * V = %.17e\n\t%s\n",
  983. X            X, "is too far from 1.\n");
  984. X        }
  985. X    /*=============================================*/
  986. X    Milestone = 200;
  987. X    /*=============================================*/
  988. X    for (Indx = 1; Indx <= 5; ++Indx)  {
  989. X        X = F9;
  990. X        switch (Indx)  {
  991. X            case 2: X = One + U2; break;
  992. X            case 3: X = V; break;
  993. X            case 4: X = UfThold; break;
  994. X            case 5: X = Radix;
  995. X            }
  996. X        Y = X;
  997. X        sigsave = sigfpe;
  998. X        if (setjmp(ovfl_buf))
  999. X            printf("  X / X  traps when X = %g\n", X);
  1000. X        else {
  1001. X            V9 = (Y / X - Half) - Half;
  1002. X            if (V9 == Zero) continue;
  1003. X            if (V9 == - U1 && Indx < 5) BadCond(Flaw, "");
  1004. X            else BadCond(Serious, "");
  1005. X            printf("  X / X differs from 1 when X = %.17e\n", X);
  1006. X            printf("  instead, X / X - 1/2 - 1/2 = %.17e .\n", V9);
  1007. X            }
  1008. X        sigsave = 0;
  1009. X        }
  1010. X    /*=============================================*/
  1011. X    Milestone = 210;
  1012. X    /*=============================================*/
  1013. X    MyZero = Zero;
  1014. X    printf("\n");
  1015. X    printf("What message and/or values does Division by Zero produce?\n") ;
  1016. #ifndef NOPAUSE
  1017. X    printf("This can interupt your program.  You can ");
  1018. X    printf("skip this part if you wish.\n");
  1019. X    printf("Do you wish to compute 1 / 0? ");
  1020. X    fflush(stdout);
  1021. X    read (KEYBOARD, ch, 8);
  1022. X    if ((ch[0] == 'Y') || (ch[0] == 'y')) {
  1023. #endif
  1024. X        sigsave = sigfpe;
  1025. X        printf("    Trying to compute 1 / 0 produces ...");
  1026. X        if (!setjmp(ovfl_buf)) printf("  %.7e .\n", One / MyZero);
  1027. X        sigsave = 0;
  1028. #ifndef NOPAUSE
  1029. X        }
  1030. X    else printf("O.K.\n");
  1031. X    printf("\nDo you wish to compute 0 / 0? ");
  1032. X    fflush(stdout);
  1033. X    read (KEYBOARD, ch, 80);
  1034. X    if ((ch[0] == 'Y') || (ch[0] == 'y')) {
  1035. #endif
  1036. X        sigsave = sigfpe;
  1037. X        printf("\n    Trying to compute 0 / 0 produces ...");
  1038. X        if (!setjmp(ovfl_buf)) printf("  %.7e .\n", Zero / MyZero);
  1039. X        sigsave = 0;
  1040. #ifndef NOPAUSE
  1041. X        }
  1042. X    else printf("O.K.\n");
  1043. #endif
  1044. X    /*=============================================*/
  1045. X    Milestone = 220;
  1046. X    /*=============================================*/
  1047. X    Pause();
  1048. X    printf("\n");
  1049. X    {
  1050. X        static char *msg[] = {
  1051. X            "FAILUREs  encountered =",
  1052. X            "SERIOUS DEFECTs  discovered =",
  1053. X            "DEFECTs  discovered =",
  1054. X            "FLAWs  discovered =" };
  1055. X        int i;
  1056. X        for(i = 0; i < 4; i++) if (ErrCnt[i])
  1057. X            printf("The number of  %-29s %d.\n",
  1058. X                msg[i], ErrCnt[i]);
  1059. X        }
  1060. X    printf("\n");
  1061. X    if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]
  1062. X            + ErrCnt[Flaw]) > 0) {
  1063. X        if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[
  1064. X            Defect] == 0) && (ErrCnt[Flaw] > 0)) {
  1065. X            printf("The arithmetic diagnosed seems ");
  1066. X            printf("Satisfactory though flawed.\n");
  1067. X            }
  1068. X        if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)
  1069. X            && ( ErrCnt[Defect] > 0)) {
  1070. X            printf("The arithmetic diagnosed may be Acceptable\n");
  1071. X            printf("despite inconvenient Defects.\n");
  1072. X            }
  1073. X        if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {
  1074. X            printf("The arithmetic diagnosed has ");
  1075. X            printf("unacceptable Serious Defects.\n");
  1076. X            }
  1077. X        if (ErrCnt[Failure] > 0) {
  1078. X            printf("Potentially fatal FAILURE may have spoiled this");
  1079. X            printf(" program's subsequent diagnoses.\n");
  1080. X            }
  1081. X        }
  1082. X    else {
  1083. X        printf("No failures, defects nor flaws have been discovered.\n");
  1084. X        if (! ((RMult == Rounded) && (RDiv == Rounded)
  1085. X            && (RAddSub == Rounded) && (RSqrt == Rounded))) 
  1086. X            printf("The arithmetic diagnosed seems Satisfactory.\n");
  1087. X        else {
  1088. X            if (StickyBit >= One &&
  1089. X                (Radix - Two) * (Radix - Nine - One) == Zero) {
  1090. X                printf("Rounding appears to conform to ");
  1091. X                printf("the proposed IEEE standard P");
  1092. X                if ((Radix == Two) &&
  1093. X                     ((Precision - Four * Three * Two) *
  1094. X                      ( Precision - TwentySeven -
  1095. X                       TwentySeven + One) == Zero)) 
  1096. X                    printf("754");
  1097. X                else printf("854");
  1098. X                if (IEEE) printf(".\n");
  1099. X                else {
  1100. X                    printf(",\nexcept for possibly Double Rounding");
  1101. X                    printf(" during Gradual Underflow.\n");
  1102. X                    }
  1103. X                }
  1104. X            printf("The arithmetic diagnosed appears to be Excellent!\n");
  1105. X            }
  1106. X        }
  1107. X    if (fpecount)
  1108. X        printf("\nA total of %d floating point exceptions were registered.\n",
  1109. X            fpecount);
  1110. X    printf("END OF TEST.\n");
  1111. X
  1112. #ifdef TEST
  1113. X    ieee_retrospective((FILE *)NULL);
  1114. #endif
  1115. X    asm("fnclex");
  1116. X        return(0);
  1117. X    }
  1118. X
  1119. /*SPLIT subs.c
  1120. #include "paranoia.h"
  1121. */
  1122. X
  1123. /* Sign */
  1124. X
  1125. FLOAT Sign (X)
  1126. FLOAT X;
  1127. { return X >= 0. ? 1.0 : -1.0; }
  1128. X
  1129. /* Pause */
  1130. X
  1131. Pause()
  1132. {
  1133. #ifndef NOPAUSE
  1134. X    char ch[8];
  1135. X
  1136. X    printf("\nTo continue, press RETURN");
  1137. X    fflush(stdout);
  1138. X    read(KEYBOARD, ch, 8);
  1139. #endif
  1140. X    printf("\nDiagnosis resumes after milestone Number %d", Milestone);
  1141. X    printf("          Page: %d\n\n", PageNo);
  1142. X    ++Milestone;
  1143. X    ++PageNo;
  1144. X    }
  1145. X
  1146. X /* TstCond */
  1147. X
  1148. TstCond (K, Valid, T)
  1149. int K, Valid;
  1150. char *T;
  1151. { if (! Valid) { BadCond(K,T); printf(".\n"); } }
  1152. X
  1153. BadCond(K, T)
  1154. int K;
  1155. char *T;
  1156. {
  1157. X    static char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
  1158. X
  1159. X    ErrCnt [K] = ErrCnt [K] + 1;
  1160. X    printf("%s:  %s", msg[K], T);
  1161. X    }
  1162. X
  1163. /* Random */
  1164. /*  Random computes
  1165. X     X = (Random1 + Random9)^5
  1166. X     Random1 = X - FLOOR(X) + 0.000005 * X;
  1167. X   and returns the new value of Random1
  1168. */
  1169. X
  1170. FLOAT Random()
  1171. {
  1172. X    FLOAT X, Y;
  1173. X    
  1174. X    X = Random1 + Random9;
  1175. X    Y = X * X;
  1176. X    Y = Y * Y;
  1177. X    X = X * Y;
  1178. X    Y = X - FLOOR(X);
  1179. X    Random1 = Y + X * 0.000005;
  1180. X    return(Random1);
  1181. X    }
  1182. X
  1183. /* SqXMinX */
  1184. X
  1185. SqXMinX (ErrKind)
  1186. int ErrKind;
  1187. {
  1188. X    FLOAT XA, XB;
  1189. X    
  1190. X    XB = X * BInvrse;
  1191. X    XA = X - XB;
  1192. X    SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp;
  1193. X    if (SqEr != Zero) {
  1194. X        if (SqEr < MinSqEr) MinSqEr = SqEr;
  1195. X        if (SqEr > MaxSqEr) MaxSqEr = SqEr;
  1196. X        J = J + 1.0;
  1197. X        BadCond(ErrKind, "\n");
  1198. X        printf("sqrt( %.17e) - %.17e  = %.17e\n", X * X, X, OneUlp * SqEr);
  1199. X        printf("\tinstead of correct value 0 .\n");
  1200. X        }
  1201. X    }
  1202. X
  1203. /* NewD */
  1204. X
  1205. NewD()
  1206. {
  1207. X    X = Z1 * Q;
  1208. X    X = FLOOR(Half - X / Radix) * Radix + X;
  1209. X    Q = (Q - X * Z) / Radix + X * X * (D / Radix);
  1210. X    Z = Z - Two * X * D;
  1211. X    if (Z <= Zero) {
  1212. X        Z = - Z;
  1213. X        Z1 = - Z1;
  1214. X        }
  1215. X    D = Radix * D;
  1216. X    }
  1217. X
  1218. /* SR3750 */
  1219. X
  1220. SR3750()
  1221. {
  1222. X    if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {
  1223. X        I = I + 1;
  1224. X        X2 = SQRT(X * D);
  1225. X        Y2 = (X2 - Z2) - (Y - Z2);
  1226. X        X2 = X8 / (Y - Half);
  1227. X        X2 = X2 - Half * X2 * X2;
  1228. X        SqEr = (Y2 + Half) + (Half - X2);
  1229. X        if (SqEr < MinSqEr) MinSqEr = SqEr;
  1230. X        SqEr = Y2 - X2;
  1231. X        if (SqEr > MaxSqEr) MaxSqEr = SqEr;
  1232. X        }
  1233. X    }
  1234. X
  1235. /* IsYeqX */
  1236. X
  1237. IsYeqX()
  1238. {
  1239. X    if (Y != X) {
  1240. X        if (N <= 0) {
  1241. X            if (Z == Zero && Q <= Zero)
  1242. X                printf("WARNING:  computing\n");
  1243. X            else BadCond(Defect, "computing\n");
  1244. X            printf("\t(%.17e) ^ (%.17e)\n", Z, Q);
  1245. X            printf("\tyielded %.17e;\n", Y);
  1246. X            printf("\twhich compared unequal to correct %.17e ;\n",
  1247. X                X);
  1248. X            printf("\t\tthey differ by %.17e .\n", Y - X);
  1249. X            }
  1250. X        N = N + 1; /* ... count discrepancies. */
  1251. X        }
  1252. X    }
  1253. X
  1254. /* SR3980 */
  1255. X
  1256. SR3980()
  1257. {
  1258. X    do {
  1259. X        Q = (FLOAT) I;
  1260. X        Y = POW(Z, Q);
  1261. X        IsYeqX();
  1262. X        if (++I > M) break;
  1263. X        X = Z * X;
  1264. X        } while ( X < W );
  1265. X    }
  1266. X
  1267. /* PrintIfNPositive */
  1268. X
  1269. PrintIfNPositive()
  1270. {
  1271. X    if (N > 0) printf("Similar discrepancies have occurred %d times.\n", N);
  1272. X    }
  1273. X
  1274. /* TstPtUf */
  1275. X
  1276. TstPtUf()
  1277. {
  1278. X    N = 0;
  1279. X    if (Z != Zero) {
  1280. X        printf("Since comparison denies Z = 0, evaluating ");
  1281. X        printf("(Z + Z) / Z should be safe.\n");
  1282. X        sigsave = sigfpe;
  1283. X        if (setjmp(ovfl_buf)) goto very_serious;
  1284. X        Q9 = (Z + Z) / Z;
  1285. X        printf("What the machine gets for (Z + Z) / Z is  %.17e .\n",
  1286. X            Q9);
  1287. X        if (FABS(Q9 - Two) < Radix * U2) {
  1288. X            printf("This is O.K., provided Over/Underflow");
  1289. X            printf(" has NOT just been signaled.\n");
  1290. X            }
  1291. X        else {
  1292. X            if ((Q9 < One) || (Q9 > Two)) {
  1293. very_serious:
  1294. X                N = 1;
  1295. X                ErrCnt [Serious] = ErrCnt [Serious] + 1;
  1296. X                printf("This is a VERY SERIOUS DEFECT!\n");
  1297. X                }
  1298. X            else {
  1299. X                N = 1;
  1300. X                ErrCnt [Defect] = ErrCnt [Defect] + 1;
  1301. X                printf("This is a DEFECT!\n");
  1302. X                }
  1303. X            }
  1304. X        sigsave = 0;
  1305. X        V9 = Z * One;
  1306. X        Random1 = V9;
  1307. X        V9 = One * Z;
  1308. X        Random2 = V9;
  1309. X        V9 = Z / One;
  1310. X        if ((Z == Random1) && (Z == Random2) && (Z == V9)) {
  1311. X            if (N > 0) Pause();
  1312. X            }
  1313. X        else {
  1314. X            N = 1;
  1315. X            BadCond(Defect, "What prints as Z = ");
  1316. X            printf("%.17e\n\tcompares different from  ", Z);
  1317. X            if (Z != Random1) printf("Z * 1 = %.17e ", Random1);
  1318. X            if (! ((Z == Random2)
  1319. X                || (Random2 == Random1)))
  1320. X                printf("1 * Z == %g\n", Random2);
  1321. X            if (! (Z == V9)) printf("Z / 1 = %.17e\n", V9);
  1322. X            if (Random2 != Random1) {
  1323. X                ErrCnt [Defect] = ErrCnt [Defect] + 1;
  1324. X                BadCond(Defect, "Multiplication does not commute!\n");
  1325. X                printf("\tComparison alleges that 1 * Z = %.17e\n",
  1326. X                    Random2);
  1327. X                printf("\tdiffers from Z * 1 = %.17e\n", Random1);
  1328. X                }
  1329. X            Pause();
  1330. X            }
  1331. X        }
  1332. X    }
  1333. X
  1334. notify(s)
  1335. char *s;
  1336. {
  1337. X    printf("%s test appears to be inconsistent...\n", s);
  1338. X    printf("   PLEASE NOTIFY KARPINKSI!\n");
  1339. X    }
  1340. X
  1341. /*SPLIT msgs.c */
  1342. X
  1343. /* Instructions */
  1344. X
  1345. msglist(s)
  1346. char **s;
  1347. { while(*s) printf("%s\n", *s++); }
  1348. X
  1349. Instructions()
  1350. {
  1351. X  static char *instr[] = {
  1352. X    "Lest this program stop prematurely, i.e. before displaying\n",
  1353. X    "    `END OF TEST',\n",
  1354. X    "try to persuade the computer NOT to terminate execution when an",
  1355. X    "error like Over/Underflow or Division by Zero occurs, but rather",
  1356. X    "to persevere with a surrogate value after, perhaps, displaying some",
  1357. X    "warning.  If persuasion avails naught, don't despair but run this",
  1358. X    "program anyway to see how many milestones it passes, and then",
  1359. X    "amend it to make further progress.\n",
  1360. X    "Answer questions with Y, y, N or n (unless otherwise indicated).\n",
  1361. X    0};
  1362. X
  1363. X    msglist(instr);
  1364. X    }
  1365. X
  1366. /* Heading */
  1367. X
  1368. Heading()
  1369. {
  1370. X  static char *head[] = {
  1371. X    "Users are invited to help debug and augment this program so it will",
  1372. X    "cope with unanticipated and newly uncovered arithmetic pathologies.\n",
  1373. X    "Please send suggestions and interesting results to",
  1374. X    "\tRichard Karpinski",
  1375. X    "\tComputer Center U-76",
  1376. X    "\tUniversity of California",
  1377. X    "\tSan Francisco, CA 94143-0704, USA\n",
  1378. X    "In doing so, please include the following information:",
  1379. #ifdef Single
  1380. X    "\tPrecision:\tsingle;",
  1381. #else
  1382. X    "\tPrecision:\tdouble;",
  1383. #endif
  1384. X    "\tVersion:\t10 February 1989;",
  1385. X    "\tComputer:\n",
  1386. X    "\tCompiler:\n",
  1387. X    "\tOptimization level:\n",
  1388. X    "\tOther relevant compiler options:",
  1389. X    0};
  1390. X
  1391. X    msglist(head);
  1392. X    }
  1393. X
  1394. /* Characteristics */
  1395. X
  1396. Characteristics()
  1397. {
  1398. X    static char *chars[] = {
  1399. X     "Running this program should reveal these characteristics:",
  1400. X    "     Radix = 1, 2, 4, 8, 10, 16, 100, 256 ...",
  1401. X    "     Precision = number of significant digits carried.",
  1402. X    "     U2 = Radix/Radix^Precision = One Ulp",
  1403. X    "\t(OneUlpnit in the Last Place) of 1.000xxx .",
  1404. X    "     U1 = 1/Radix^Precision = One Ulp of numbers a little less than 1.0 .",
  1405. X    "     Adequacy of guard digits for Mult., Div. and Subt.",
  1406. X    "     Whether arithmetic is chopped, correctly rounded, or something else",
  1407. X    "\tfor Mult., Div., Add/Subt. and Sqrt.",
  1408. X    "     Whether a Sticky Bit used correctly for rounding.",
  1409. X    "     UnderflowThreshold = an underflow threshold.",
  1410. X    "     E0 and PseudoZero tell whether underflow is abrupt, gradual, or fuzzy.",
  1411. X    "     V = an overflow threshold, roughly.",
  1412. X    "     V0  tells, roughly, whether  Infinity  is represented.",
  1413. X    "     Comparisions are checked for consistency with subtraction",
  1414. X    "\tand for contamination with pseudo-zeros.",
  1415. X    "     Sqrt is tested.  Y^X is not tested.",
  1416. X    "     Extra-precise subexpressions are revealed but NOT YET tested.",
  1417. X    "     Decimal-Binary conversion is NOT YET tested for accuracy.",
  1418. X    0};
  1419. X
  1420. X    msglist(chars);
  1421. X    }
  1422. X
  1423. History()
  1424. X
  1425. { /* History */
  1426. X /* Converted from Brian Wichmann's Pascal version to C by Thos Sumner,
  1427. X    with further massaging by David M. Gay. */
  1428. X
  1429. X  static char *hist[] = {
  1430. X    "The program attempts to discriminate among",
  1431. X    "   FLAWs, like lack of a sticky bit,",
  1432. X    "   Serious DEFECTs, like lack of a guard digit, and",
  1433. X    "   FAILUREs, like 2+2 == 5 .",
  1434. X    "Failures may confound subsequent diagnoses.\n",
  1435. X    "The diagnostic capabilities of this program go beyond an earlier",
  1436. X    "program called `MACHAR', which can be found at the end of the",
  1437. X    "book  `Software Manual for the Elementary Functions' (1980) by",
  1438. X    "W. J. Cody and W. Waite. Although both programs try to discover",
  1439. X    "the Radix, Precision and range (over/underflow thresholds)",
  1440. X    "of the arithmetic, this program tries to cope with a wider variety",
  1441. X    "of pathologies, and to say how well the arithmetic is implemented.",
  1442. X    "\nThe program is based upon a conventional radix representation for",
  1443. X    "floating-point numbers, but also allows logarithmic encoding",
  1444. X    "as used by certain early WANG machines.\n",
  1445. X    "BASIC version of this program (C) 1983 by Prof. W. M. Kahan;",
  1446. X    "see source comments for more history.",
  1447. X    0};
  1448. X
  1449. X    msglist(hist);
  1450. X    }
  1451. SHAR_EOF
  1452. echo 'File paranoia.c is complete' &&
  1453. chmod 0644 paranoia.c ||
  1454. echo 'restore of paranoia.c failed'
  1455. Wc_c="`wc -c < 'paranoia.c'`"
  1456. test 57395 -eq "$Wc_c" ||
  1457.     echo 'paranoia.c: original size 57395, current size' "$Wc_c"
  1458. rm -f _shar_wnt_.tmp
  1459. fi
  1460. # ============= CHANGELOG ==============
  1461. if test -f 'CHANGELOG' -a X"$1" != X"-c"; then
  1462.     echo 'x - skipping CHANGELOG (File already exists)'
  1463.     rm -f _shar_wnt_.tmp
  1464. else
  1465. > _shar_wnt_.tmp
  1466. echo 'x - extracting CHANGELOG (Text)'
  1467. sed 's/^X//' << 'SHAR_EOF' > 'CHANGELOG' &&
  1468. Thu Nov 15 07:05:39 EDT 1990
  1469. Shar'd beta version for release to testers.
  1470. X
  1471. Sat Nov 17 17:06:17 EDT 1990
  1472. Fixed floor and ceil so that they address their local variables at -ve
  1473. offsets relative to the frame pointer.
  1474. X
  1475. Fixed hypot.s so that unnecessary register loads are avoided. I wasn't aware
  1476. of the correct assembler mnemonic to multiply %st(0) by the contents of a memory
  1477. location.
  1478. X
  1479. Sat Nov 17 18:49:31 EDT 1990
  1480. Copysign rewritten in assembler.
  1481. X
  1482. Sun Nov 18 11:21:24 EDT 1990
  1483. Sinh and cosh now in assembler.
  1484. Tanh now in assembler.
  1485. Fixed nasty bug in acos.s.
  1486. exp__E no longer needed - source deleted.
  1487. X
  1488. Tue Nov 20 21:24:26 EDT 1990
  1489. Sinh and cosh now multiply by 0.5 instead of dividing by 2. It's faster.
  1490. Added setinternal for convenience in some cases. Using this makes the
  1491. output of paranoia closer to that obtained on a Sun 4.
  1492. Added setcont.
  1493. Modified the Makefile so that test and paranoia are synonymous.
  1494. Modified atan2.c to avoid a floating point register load.
  1495. Fixed bug in atan2 - was giving incorrect result when y/x = +/-inf.
  1496. X
  1497. Thu Nov 22 06:52:47 EDT 1990
  1498. Fixed local variable handling in assembler files. Stack needs to be 
  1499. allocated for local variables.
  1500. All assembler routines are now .align 4.
  1501. X
  1502. Thu Nov 22 19:24:09 EDT 1990
  1503. Added MASK_ALL to fpumath.h
  1504. X
  1505. Fri Nov 23 21:11:33 EDT 1990
  1506. Inverse hyperbolics missing from fpumath.h - fixed.
  1507. Some definitions missing from atanh.c - didn't work on some systems.
  1508. Coded inverse hyperbolics in assembler.
  1509. Coded atan2 in assembler.
  1510. X
  1511. Sat Nov 24 17:12:34 EDT 1990
  1512. Modified hypot.s to avoid use of ffree.
  1513. atan2 was not working correctly - fixed
  1514. Ready for netwide release.
  1515. X
  1516. Mon Nov 26 19:25:56 EDT 1990
  1517. Put in explicit call to cc for compiling lgamma. This means that the
  1518. library can link to code not compiled with gcc (not quite true pow.s has to be
  1519. changed).
  1520. X
  1521. Wed Nov 28 22:57:43 EDT 1990
  1522. Coded the is??? functions required for IEEE conformance
  1523. Coded ieee_retrospective() which prints a summary of IEEE exceptions that are on
  1524. when the function is called.
  1525. X
  1526. Fri Nov 30 17:35:29 EDT 1990
  1527. Removed the last vestige of GNU specifics. You can now use the library with 'cc'
  1528. although you need gcc/gas to create it.
  1529. X
  1530. Sat Dec  1 15:48:19 EDT 1990
  1531. Found a bug in pow (pow(-x,x) shoudn't work for x non-integral) fixed.
  1532. X
  1533. Sat Dec  1 21:00:48 EST 1990
  1534. Rewrote copysign in a more efficient way.
  1535. X
  1536. Sun Dec  2 06:51:59 EDT 1990
  1537. Added infinity().
  1538. Amended README to show the current state.
  1539. Changed fpumath.h so as to define HUGE = infinity() if you compile (user
  1540. code) with IEEE defined.
  1541. X
  1542. Sun Dec  2 08:26:20 EDT 1990
  1543. Modified log functions and inverse hyperbolics for greater accuracy.
  1544. X
  1545. Sun Dec  2 10:08:47 EDT 1990
  1546. Added MASK_ALL1 to fpumath.h. This enables 64 bit precision significands.
  1547. Coded sqrtp to ensure that sqrt is rounded to 53 bits - used in
  1548. paranoia.c.
  1549. Put sqrtp in fpumath.h and in README.
  1550. Made very minor changes to pow.s, floor.s and ceil.s.
  1551. ieee_retrospective was printing its leading info for a LOS exception -
  1552. fixed.
  1553. Rewrote finite so that it does not use the floating point unit and so does
  1554. not raise any exceptions.
  1555. Produced d2dcomb.summ.
  1556. X
  1557. Sun Dec  2 18:40:47 EDT 1990
  1558. Fixed exp, exp10, expm1, exp2, sinh, cosh, tanh to gain more speed.
  1559. X
  1560. Wed Dec  5 17:25:50 EDT 1990
  1561. Coded max_..., min_..., etc.
  1562. X
  1563. Wed Dec  5 19:38:56 EDT 1990
  1564. Fixed some local parameter handling errors in a few of the assembler
  1565. files.
  1566. X
  1567. Wed Dec  5 20:54:37 EDT 1990
  1568. Define MAX(MIN)DOUBLE in terms of max(min)_(sub)normal in fpumath.h.
  1569. X
  1570. Sat Dec  8 06:09:01 EDT 1990
  1571. Fixed bug in isnan 8(%ebp) not $8(%ebp) - typo!
  1572. Fixed bug in isinf - not doing correct test.
  1573. X
  1574. Sat Dec  8 06:36:07 EDT 1990
  1575. Coded nextafter.
  1576. X
  1577. Sat Dec  8 06:46:44 EDT 1990
  1578. MAXDOUBLE is wrong in math.h it is correct in fpumath.h
  1579. Released to net. Posted to alt.sources.
  1580. X
  1581. Sun Dec  9 08:26:23 EDT 1990
  1582. Added a couple more special cases to nextafter.
  1583. X
  1584. Sun Dec  9 12:16:04 EDT 1990
  1585. Slight modification to Makefile
  1586. X
  1587. Sun Dec  9 14:33:13 EDT 1990
  1588. Coded float version - still testing. You must have an ANSI compiler (gcc is
  1589. ok) to use this version.
  1590. X
  1591. Mon Dec 10 17:21:16 EDT 1990
  1592. Float version done.
  1593. Fully prototyped all functions in fpumath.h.
  1594. Version 2.0 done.
  1595. X
  1596. Tue Dec 11 06:22:51 EDT 1990
  1597. Added a couple of extra cases to nextafter[f].
  1598. Made a couple of minor changes to fpumath.h.
  1599. Removed AT&T code (j0.c, j1.c, jn.c, erf.c, lgamma.c). These were
  1600. avaialble with BSD 4.3 so I had no way of knowing...
  1601. X
  1602. Wed Dec 12 05:50:31 EDT 1990
  1603. Finally got nextafter[f] working correctly.
  1604. X
  1605. Thu Dec 13 09:44:37 EDT 1990
  1606. sinh[f]() now gives good results for all argument ranges. Compute using a 
  1607. combination of exp() and a ratio of polynomials.
  1608. X
  1609. Thu Dec 13 18:48:55 EDT 1990
  1610. Compute expm1() as expm1() = exp(x)*(1-exp(-x)). This is much better than
  1611. before but still not perfect. I need Hart & Cheney...
  1612. X
  1613. Thu Dec 13 19:38:31 EDT 1990
  1614. cosh() was not producing good results for small x. Fixed.
  1615. X
  1616. Fri Dec 14 08:19:56 EDT 1990
  1617. Berkeley provided j0.c, j1.c, erf.c, gamma.c and lgamma.c (no jn.c yet) now
  1618. being used. This is *not* BSD code and is not covered by the BSD licence.
  1619. ***No Berkeley code used.***
  1620. mathimpl.h deleted.
  1621. README updated.
  1622. Defined _getsw() (in ieee_retrospective.c) as a short.
  1623. X
  1624. Fri Dec 14 11:43:30 EDT 1990
  1625. Added shar to Makefile.
  1626. Version 2.0 ready for distribution.
  1627. SHAR_EOF
  1628. chmod 0644 CHANGELOG ||
  1629. echo 'restore of CHANGELOG failed'
  1630. Wc_c="`wc -c < 'CHANGELOG'`"
  1631. test 5231 -eq "$Wc_c" ||
  1632.     echo 'CHANGELOG: original size 5231, current size' "$Wc_c"
  1633. rm -f _shar_wnt_.tmp
  1634. fi
  1635. # ============= COPYING ==============
  1636. if test -f 'COPYING' -a X"$1" != X"-c"; then
  1637.     echo 'x - skipping COPYING (File already exists)'
  1638.     rm -f _shar_wnt_.tmp
  1639. else
  1640. > _shar_wnt_.tmp
  1641. echo 'x - extracting COPYING (Text)'
  1642. sed 's/^X//' << 'SHAR_EOF' > 'COPYING' &&
  1643. X
  1644. X            GNU GENERAL PUBLIC LICENSE
  1645. X             Version 1, February 1989
  1646. X
  1647. X Copyright (C) 1989 Free Software Foundation, Inc.
  1648. X                    675 Mass Ave, Cambridge, MA 02139, USA
  1649. X Everyone is permitted to copy and distribute verbatim copies
  1650. X of this license document, but changing it is not allowed.
  1651. X
  1652. X                Preamble
  1653. X
  1654. X  The license agreements of most software companies try to keep users
  1655. at the mercy of those companies.  By contrast, our General Public
  1656. License is intended to guarantee your freedom to share and change free
  1657. software--to make sure the software is free for all its users.  The
  1658. General Public License applies to the Free Software Foundation's
  1659. software and to any other program whose authors commit to using it.
  1660. You can use it for your programs, too.
  1661. X
  1662. X  When we speak of free software, we are referring to freedom, not
  1663. price.  Specifically, the General Public License is designed to make
  1664. sure that you have the freedom to give away or sell copies of free
  1665. software, that you receive source code or can get it if you want it,
  1666. that you can change the software or use pieces of it in new free
  1667. programs; and that you know you can do these things.
  1668. X
  1669. X  To protect your rights, we need to make restrictions that forbid
  1670. anyone to deny you these rights or to ask you to surrender the rights.
  1671. These restrictions translate to certain responsibilities for you if you
  1672. distribute copies of the software, or if you modify it.
  1673. X
  1674. X  For example, if you distribute copies of a such a program, whether
  1675. gratis or for a fee, you must give the recipients all the rights that
  1676. you have.  You must make sure that they, too, receive or can get the
  1677. source code.  And you must tell them their rights.
  1678. X
  1679. X  We protect your rights with two steps: (1) copyright the software, and
  1680. (2) offer you this license which gives you legal permission to copy,
  1681. distribute and/or modify the software.
  1682. X
  1683. X  Also, for each author's protection and ours, we want to make certain
  1684. that everyone understands that there is no warranty for this free
  1685. software.  If the software is modified by someone else and passed on, we
  1686. want its recipients to know that what they have is not the original, so
  1687. that any problems introduced by others will not reflect on the original
  1688. authors' reputations.
  1689. X
  1690. X  The precise terms and conditions for copying, distribution and
  1691. modification follow.
  1692. X
  1693. X            GNU GENERAL PUBLIC LICENSE
  1694. X   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
  1695. X
  1696. X  0. This License Agreement applies to any program or other work which
  1697. contains a notice placed by the copyright holder saying it may be
  1698. distributed under the terms of this General Public License.  The
  1699. "Program", below, refers to any such program or work, and a "work based
  1700. on the Program" means either the Program or any work containing the
  1701. Program or a portion of it, either verbatim or with modifications.  Each
  1702. licensee is addressed as "you".
  1703. X
  1704. X  1. You may copy and distribute verbatim copies of the Program's source
  1705. code as you receive it, in any medium, provided that you conspicuously and
  1706. appropriately publish on each copy an appropriate copyright notice and
  1707. disclaimer of warranty; keep intact all the notices that refer to this
  1708. General Public License and to the absence of any warranty; and give any
  1709. other recipients of the Program a copy of this General Public License
  1710. along with the Program.  You may charge a fee for the physical act of
  1711. transferring a copy.
  1712. X
  1713. X  2. You may modify your copy or copies of the Program or any portion of
  1714. it, and copy and distribute such modifications under the terms of Paragraph
  1715. 1 above, provided that you also do the following:
  1716. X
  1717. X    a) cause the modified files to carry prominent notices stating that
  1718. X    you changed the files and the date of any change; and
  1719. X
  1720. X    b) cause the whole of any work that you distribute or publish, that
  1721. X    in whole or in part contains the Program or any part thereof, either
  1722. X    with or without modifications, to be licensed at no charge to all
  1723. X    third parties under the terms of this General Public License (except
  1724. X    that you may choose to grant warranty protection to some or all
  1725. X    third parties, at your option).
  1726. X
  1727. X    c) If the modified program normally reads commands interactively when
  1728. X    run, you must cause it, when started running for such interactive use
  1729. X    in the simplest and most usual way, to print or display an
  1730. X    announcement including an appropriate copyright notice and a notice
  1731. X    that there is no warranty (or else, saying that you provide a
  1732. X    warranty) and that users may redistribute the program under these
  1733. X    conditions, and telling the user how to view a copy of this General
  1734. X    Public License.
  1735. X
  1736. X    d) You may charge a fee for the physical act of transferring a
  1737. X    copy, and you may at your option offer warranty protection in
  1738. X    exchange for a fee.
  1739. X
  1740. Mere aggregation of another independent work with the Program (or its
  1741. derivative) on a volume of a storage or distribution medium does not bring
  1742. the other work under the scope of these terms.
  1743. X
  1744. X  3. You may copy and distribute the Program (or a portion or derivative of
  1745. it, under Paragraph 2) in object code or executable form under the terms of
  1746. Paragraphs 1 and 2 above provided that you also do one of the following:
  1747. X
  1748. X    a) accompany it with the complete corresponding machine-readable
  1749. X    source code, which must be distributed under the terms of
  1750. X    Paragraphs 1 and 2 above; or,
  1751. X
  1752. X    b) accompany it with a written offer, valid for at least three
  1753. X    years, to give any third party free (except for a nominal charge
  1754. X    for the cost of distribution) a complete machine-readable copy of the
  1755. X    corresponding source code, to be distributed under the terms of
  1756. X    Paragraphs 1 and 2 above; or,
  1757. X
  1758. X    c) accompany it with the information you received as to where the
  1759. X    corresponding source code may be obtained.  (This alternative is
  1760. X    allowed only for noncommercial distribution and only if you
  1761. X    received the program in object code or executable form alone.)
  1762. X
  1763. Source code for a work means the preferred form of the work for making
  1764. modifications to it.  For an executable file, complete source code means
  1765. all the source code for all modules it contains; but, as a special
  1766. exception, it need not include source code for modules which are standard
  1767. libraries that accompany the operating system on which the executable
  1768. file runs, or for standard header files or definitions files that
  1769. accompany that operating system.
  1770. X
  1771. X  4. You may not copy, modify, sublicense, distribute or transfer the
  1772. Program except as expressly provided under this General Public License.
  1773. Any attempt otherwise to copy, modify, sublicense, distribute or transfer
  1774. the Program is void, and will automatically terminate your rights to use
  1775. the Program under this License.  However, parties who have received
  1776. copies, or rights to use copies, from you under this General Public
  1777. License will not have their licenses terminated so long as such parties
  1778. remain in full compliance.
  1779. X
  1780. X  5. By copying, distributing or modifying the Program (or any work based
  1781. on the Program) you indicate your acceptance of this license to do so,
  1782. and all its terms and conditions.
  1783. X
  1784. X  6. Each time you redistribute the Program (or any work based on the
  1785. Program), the recipient automatically receives a license from the original
  1786. licensor to copy, distribute or modify the Program subject to these
  1787. terms and conditions.  You may not impose any further restrictions on the
  1788. recipients' exercise of the rights granted herein.
  1789. X
  1790. X  7. The Free Software Foundation may publish revised and/or new versions
  1791. of the General Public License from time to time.  Such new versions will
  1792. be similar in spirit to the present version, but may differ in detail to
  1793. address new problems or concerns.
  1794. X
  1795. Each version is given a distinguishing version number.  If the Program
  1796. specifies a version number of the license which applies to it and "any
  1797. later version", you have the option of following the terms and conditions
  1798. either of that version or of any later version published by the Free
  1799. Software Foundation.  If the Program does not specify a version number of
  1800. the license, you may choose any version ever published by the Free Software
  1801. Foundation.
  1802. X
  1803. X  8. If you wish to incorporate parts of the Program into other free
  1804. programs whose distribution conditions are different, write to the author
  1805. to ask for permission.  For software which is copyrighted by the Free
  1806. Software Foundation, write to the Free Software Foundation; we sometimes
  1807. make exceptions for this.  Our decision will be guided by the two goals
  1808. of preserving the free status of all derivatives of our free software and
  1809. of promoting the sharing and reuse of software generally.
  1810. X
  1811. X                NO WARRANTY
  1812. X
  1813. X  9. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
  1814. FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
  1815. OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
  1816. PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
  1817. OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
  1818. MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
  1819. TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
  1820. PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
  1821. REPAIR OR CORRECTION.
  1822. X
  1823. X  10. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
  1824. WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
  1825. REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
  1826. INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
  1827. OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
  1828. TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
  1829. YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
  1830. PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
  1831. POSSIBILITY OF SUCH DAMAGES.
  1832. X
  1833. X             END OF TERMS AND CONDITIONS
  1834. X
  1835. X    Appendix: How to Apply These Terms to Your New Programs
  1836. X
  1837. X  If you develop a new program, and you want it to be of the greatest
  1838. possible use to humanity, the best way to achieve this is to make it
  1839. free software which everyone can redistribute and change under these
  1840. terms.
  1841. X
  1842. X  To do so, attach the following notices to the program.  It is safest to
  1843. attach them to the start of each source file to most effectively convey
  1844. the exclusion of warranty; and each file should have at least the
  1845. "copyright" line and a pointer to where the full notice is found.
  1846. X
  1847. X    <one line to give the program's name and a brief idea of what it does.>
  1848. X    Copyright (C) 19yy  <name of author>
  1849. SHAR_EOF
  1850. true || echo 'restore of COPYING failed'
  1851. fi
  1852. echo 'End of mathlib2.0 part 5'
  1853. echo 'File COPYING is continued in part 6'
  1854. echo 6 > _shar_seq_.tmp
  1855. exit 0
  1856. --
  1857. Glenn Geers                       | "So when it's over, we're back to people.
  1858. Department of Theoretical Physics |  Just to prove that human touch can have
  1859. The University of Sydney          |  no equal."
  1860. Sydney NSW 2006 Australia         |  - Basia Trzetrzelewska, 'Prime Time TV'
  1861.