paranoia.c 57 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964196519661967196819691970197119721973197419751976197719781979198019811982198319841985198619871988198919901991199219931994199519961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018201920202021202220232024202520262027202820292030203120322033203420352036203720382039204020412042204320442045204620472048204920502051205220532054205520562057205820592060206120622063206420652066206720682069207020712072207320742075207620772078207920802081208220832084208520862087208820892090209120922093209420952096209720982099210021012102210321042105210621072108210921102111211221132114211521162117211821192120212121222123212421252126212721282129213021312132213321342135213621372138213921402141214221432144214521462147214821492150215121522153215421552156215721582159216021612162216321642165216621672168216921702171217221732174217521762177217821792180218121822183218421852186218721882189219021912192219321942195219621972198219922002201220222032204
  1. #undef V9
  2. #define NOPAUSE
  3. /* A C version of Kahan's Floating Point Test "Paranoia"
  4. Thos Sumner, UCSF, Feb. 1985
  5. David Gay, BTL, Jan. 1986
  6. This is a rewrite from the Pascal version by
  7. B. A. Wichmann, 18 Jan. 1985
  8. (and does NOT exhibit good C programming style).
  9. (C) Apr 19 1983 in BASIC version by:
  10. Professor W. M. Kahan,
  11. 567 Evans Hall
  12. Electrical Engineering & Computer Science Dept.
  13. University of California
  14. Berkeley, California 94720
  15. USA
  16. converted to Pascal by:
  17. B. A. Wichmann
  18. National Physical Laboratory
  19. Teddington Middx
  20. TW11 OLW
  21. UK
  22. converted to C by:
  23. David M. Gay and Thos Sumner
  24. AT&T Bell Labs Computer Center, Rm. U-76
  25. 600 Mountain Avenue University of California
  26. Murray Hill, NJ 07974 San Francisco, CA 94143
  27. USA USA
  28. with simultaneous corrections to the Pascal source (reflected
  29. in the Pascal source available over netlib).
  30. [A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.]
  31. Reports of results on various systems from all the versions
  32. of Paranoia are being collected by Richard Karpinski at the
  33. same address as Thos Sumner. This includes sample outputs,
  34. bug reports, and criticisms.
  35. You may copy this program freely if you acknowledge its source.
  36. Comments on the Pascal version to NPL, please.
  37. The C version catches signals from floating-point exceptions.
  38. If signal(SIGFPE,...) is unavailable in your environment, you may
  39. #define NOSIGNAL to comment out the invocations of signal.
  40. This source file is too big for some C compilers, but may be split
  41. into pieces. Comments containing "SPLIT" suggest convenient places
  42. for this splitting. At the end of these comments is an "ed script"
  43. (for the UNIX(tm) editor ed) that will do this splitting.
  44. By #defining Single when you compile this source, you may obtain
  45. a single-precision C version of Paranoia.
  46. The following is from the introductory commentary from Wichmann's work:
  47. The BASIC program of Kahan is written in Microsoft BASIC using many
  48. facilities which have no exact analogy in Pascal. The Pascal
  49. version below cannot therefore be exactly the same. Rather than be
  50. a minimal transcription of the BASIC program, the Pascal coding
  51. follows the conventional style of block-structured languages. Hence
  52. the Pascal version could be useful in producing versions in other
  53. structured languages.
  54. Rather than use identifiers of minimal length (which therefore have
  55. little mnemonic significance), the Pascal version uses meaningful
  56. identifiers as follows [Note: A few changes have been made for C]:
  57. BASIC C BASIC C BASIC C
  58. A J S StickyBit
  59. A1 AInverse J0 NoErrors T
  60. B Radix [Failure] T0 Underflow
  61. B1 BInverse J1 NoErrors T2 ThirtyTwo
  62. B2 RadixD2 [SeriousDefect] T5 OneAndHalf
  63. B9 BMinusU2 J2 NoErrors T7 TwentySeven
  64. C [Defect] T8 TwoForty
  65. C1 CInverse J3 NoErrors U OneUlp
  66. D [Flaw] U0 UnderflowThreshold
  67. D4 FourD K PageNo U1
  68. E0 L Milestone U2
  69. E1 M V
  70. E2 Exp2 N V0
  71. E3 N1 V8
  72. E5 MinSqEr O Zero V9
  73. E6 SqEr O1 One W
  74. E7 MaxSqEr O2 Two X
  75. E8 O3 Three X1
  76. E9 O4 Four X8
  77. F1 MinusOne O5 Five X9 Random1
  78. F2 Half O8 Eight Y
  79. F3 Third O9 Nine Y1
  80. F6 P Precision Y2
  81. F9 Q Y9 Random2
  82. G1 GMult Q8 Z
  83. G2 GDiv Q9 Z0 PseudoZero
  84. G3 GAddSub R Z1
  85. H R1 RMult Z2
  86. H1 HInverse R2 RDiv Z9
  87. I R3 RAddSub
  88. IO NoTrials R4 RSqrt
  89. I3 IEEE R9 Random9
  90. SqRWrng
  91. All the variables in BASIC are true variables and in consequence,
  92. the program is more difficult to follow since the "constants" must
  93. be determined (the glossary is very helpful). The Pascal version
  94. uses Real constants, but checks are added to ensure that the values
  95. are correctly converted by the compiler.
  96. The major textual change to the Pascal version apart from the
  97. identifiersis that named procedures are used, inserting parameters
  98. wherehelpful. New procedures are also introduced. The
  99. correspondence is as follows:
  100. BASIC Pascal
  101. lines
  102. 90- 140 Pause
  103. 170- 250 Instructions
  104. 380- 460 Heading
  105. 480- 670 Characteristics
  106. 690- 870 History
  107. 2940-2950 Random
  108. 3710-3740 NewD
  109. 4040-4080 DoesYequalX
  110. 4090-4110 PrintIfNPositive
  111. 4640-4850 TestPartialUnderflow
  112. =*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
  113. Below is an "ed script" that splits para.c into 10 files
  114. of the form part[1-8].c, subs.c, and msgs.c, plus a header
  115. file, paranoia.h, that these files require.
  116. r paranoia.c
  117. $
  118. ?SPLIT
  119. +,$w msgs.c
  120. .,$d
  121. ?SPLIT
  122. .d
  123. +d
  124. -,$w subs.c
  125. -,$d
  126. ?part8
  127. +d
  128. ?include
  129. .,$w part8.c
  130. .,$d
  131. -d
  132. ?part7
  133. +d
  134. ?include
  135. .,$w part7.c
  136. .,$d
  137. -d
  138. ?part6
  139. +d
  140. ?include
  141. .,$w part6.c
  142. .,$d
  143. -d
  144. ?part5
  145. +d
  146. ?include
  147. .,$w part5.c
  148. .,$d
  149. -d
  150. ?part4
  151. +d
  152. ?include
  153. .,$w part4.c
  154. .,$d
  155. -d
  156. ?part3
  157. +d
  158. ?include
  159. .,$w part3.c
  160. .,$d
  161. -d
  162. ?part2
  163. +d
  164. ?include
  165. .,$w part2.c
  166. .,$d
  167. ?SPLIT
  168. .d
  169. 1,/^#include/-1d
  170. 1,$w part1.c
  171. /Computed constants/,$d
  172. 1,$s/^int/extern &/
  173. 1,$s/^FLOAT/extern &/
  174. 1,$s/^char/extern &/
  175. 1,$s! = .*!;!
  176. /^Guard/,/^Round/s/^/extern /
  177. /^jmp_buf/s/^/extern /
  178. /^Sig_type/s/^/extern /
  179. s/$/\
  180. extern void sigfpe();/
  181. w paranoia.h
  182. q
  183. */
  184. #include <stdio.h>
  185. #ifndef NOSIGNAL
  186. #include <signal.h>
  187. #endif
  188. #include <setjmp.h>
  189. extern double fabs(), floor(), log(), pow(), sqrt();
  190. #ifdef Single
  191. #define FLOAT float
  192. #define FABS(x) (float)fabs((double)(x))
  193. #define FLOOR(x) (float)floor((double)(x))
  194. #define LOG(x) (float)log((double)(x))
  195. #define POW(x,y) (float)pow((double)(x),(double)(y))
  196. #define SQRT(x) (float)sqrt((double)(x))
  197. #else
  198. #define FLOAT double
  199. #define FABS(x) fabs(x)
  200. #define FLOOR(x) floor(x)
  201. #define LOG(x) log(x)
  202. #define POW(x,y) pow(x,y)
  203. #define SQRT(x) sqrt(x)
  204. #endif
  205. jmp_buf ovfl_buf;
  206. typedef void (*Sig_type)();
  207. Sig_type sigsave;
  208. #define KEYBOARD 0
  209. FLOAT Radix, BInvrse, RadixD2, BMinusU2;
  210. FLOAT Sign(), Random();
  211. /*Small floating point constants.*/
  212. FLOAT Zero = 0.0;
  213. FLOAT Half = 0.5;
  214. FLOAT One = 1.0;
  215. FLOAT Two = 2.0;
  216. FLOAT Three = 3.0;
  217. FLOAT Four = 4.0;
  218. FLOAT Five = 5.0;
  219. FLOAT Eight = 8.0;
  220. FLOAT Nine = 9.0;
  221. FLOAT TwentySeven = 27.0;
  222. FLOAT ThirtyTwo = 32.0;
  223. FLOAT TwoForty = 240.0;
  224. FLOAT MinusOne = -1.0;
  225. FLOAT OneAndHalf = 1.5;
  226. /*Integer constants*/
  227. int NoTrials = 20; /*Number of tests for commutativity. */
  228. #define False 0
  229. #define True 1
  230. /* Definitions for declared types
  231. Guard == (Yes, No);
  232. Rounding == (Chopped, Rounded, Other);
  233. Message == packed array [1..40] of char;
  234. Class == (Flaw, Defect, Serious, Failure);
  235. */
  236. #define Yes 1
  237. #define No 0
  238. #define Chopped 2
  239. #define Rounded 1
  240. #define Other 0
  241. #define Flaw 3
  242. #define Defect 2
  243. #define Serious 1
  244. #define Failure 0
  245. typedef int Guard, Rounding, Class;
  246. typedef char Message;
  247. /* Declarations of Variables */
  248. int Indx;
  249. char ch[8];
  250. FLOAT AInvrse, A1;
  251. FLOAT C, CInvrse;
  252. FLOAT D, FourD;
  253. FLOAT E0, E1, Exp2, E3, MinSqEr;
  254. FLOAT SqEr, MaxSqEr, E9;
  255. FLOAT Third;
  256. FLOAT F6, F9;
  257. FLOAT H, HInvrse;
  258. int I;
  259. FLOAT StickyBit, J;
  260. FLOAT MyZero;
  261. FLOAT Precision;
  262. FLOAT Q, Q9;
  263. FLOAT R, Random9;
  264. FLOAT T, Underflow, S;
  265. FLOAT OneUlp, UfThold, U1, U2;
  266. FLOAT V, V0, V9;
  267. FLOAT W;
  268. FLOAT X, X1, X2, X8, Random1;
  269. FLOAT Y, Y1, Y2, Random2;
  270. FLOAT Z, PseudoZero, Z1, Z2, Z9;
  271. int ErrCnt[4];
  272. int fpecount;
  273. int Milestone;
  274. int PageNo;
  275. int M, N, N1;
  276. Guard GMult, GDiv, GAddSub;
  277. Rounding RMult, RDiv, RAddSub, RSqrt;
  278. int Break, Done, NotMonot, Monot, Anomaly, IEEE,
  279. SqRWrng, UfNGrad;
  280. /* Computed constants. */
  281. /*U1 gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
  282. /*U2 gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */
  283. /* floating point exception receiver */
  284. void
  285. sigfpe(i)
  286. {
  287. fpecount++;
  288. printf("\n* * * FLOATING-POINT ERROR * * *\n");
  289. fflush(stdout);
  290. if (sigsave) {
  291. #ifndef NOSIGNAL
  292. signal(SIGFPE, sigsave);
  293. #endif
  294. sigsave = 0;
  295. longjmp(ovfl_buf, 1);
  296. }
  297. abort();
  298. }
  299. main()
  300. {
  301. #ifdef mc
  302. char *out;
  303. ieee_flags("set", "precision", "double", &out);
  304. #endif
  305. /* First two assignments use integer right-hand sides. */
  306. Zero = 0;
  307. One = 1;
  308. Two = One + One;
  309. Three = Two + One;
  310. Four = Three + One;
  311. Five = Four + One;
  312. Eight = Four + Four;
  313. Nine = Three * Three;
  314. TwentySeven = Nine * Three;
  315. ThirtyTwo = Four * Eight;
  316. TwoForty = Four * Five * Three * Four;
  317. MinusOne = -One;
  318. Half = One / Two;
  319. OneAndHalf = One + Half;
  320. ErrCnt[Failure] = 0;
  321. ErrCnt[Serious] = 0;
  322. ErrCnt[Defect] = 0;
  323. ErrCnt[Flaw] = 0;
  324. PageNo = 1;
  325. /*=============================================*/
  326. Milestone = 0;
  327. /*=============================================*/
  328. #ifndef NOSIGNAL
  329. signal(SIGFPE, sigfpe);
  330. #endif
  331. Instructions();
  332. Pause();
  333. Heading();
  334. Pause();
  335. Characteristics();
  336. Pause();
  337. History();
  338. Pause();
  339. /*=============================================*/
  340. Milestone = 7;
  341. /*=============================================*/
  342. printf("Program is now RUNNING tests on small integers:\n");
  343. TstCond (Failure, (Zero + Zero == Zero) && (One - One == Zero)
  344. && (One > Zero) && (One + One == Two),
  345. "0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2");
  346. Z = - Zero;
  347. if (Z != 0.0) {
  348. ErrCnt[Failure] = ErrCnt[Failure] + 1;
  349. printf("Comparison alleges that -0.0 is Non-zero!\n");
  350. U1 = 0.001;
  351. Radix = 1;
  352. TstPtUf();
  353. }
  354. TstCond (Failure, (Three == Two + One) && (Four == Three + One)
  355. && (Four + Two * (- Two) == Zero)
  356. && (Four - Three - One == Zero),
  357. "3 != 2+1, 4 != 3+1, 4+2*(-2) != 0, or 4-3-1 != 0");
  358. TstCond (Failure, (MinusOne == (0 - One))
  359. && (MinusOne + One == Zero ) && (One + MinusOne == Zero)
  360. && (MinusOne + FABS(One) == Zero)
  361. && (MinusOne + MinusOne * MinusOne == Zero),
  362. "-1+1 != 0, (-1)+abs(1) != 0, or -1+(-1)*(-1) != 0");
  363. TstCond (Failure, Half + MinusOne + Half == Zero,
  364. "1/2 + (-1) + 1/2 != 0");
  365. /*=============================================*/
  366. /*SPLIT
  367. part2();
  368. part3();
  369. part4();
  370. part5();
  371. part6();
  372. part7();
  373. part8();
  374. }
  375. #include "paranoia.h"
  376. part2(){
  377. */
  378. Milestone = 10;
  379. /*=============================================*/
  380. TstCond (Failure, (Nine == Three * Three)
  381. && (TwentySeven == Nine * Three) && (Eight == Four + Four)
  382. && (ThirtyTwo == Eight * Four)
  383. && (ThirtyTwo - TwentySeven - Four - One == Zero),
  384. "9 != 3*3, 27 != 9*3, 32 != 8*4, or 32-27-4-1 != 0");
  385. TstCond (Failure, (Five == Four + One) &&
  386. (TwoForty == Four * Five * Three * Four)
  387. && (TwoForty / Three - Four * Four * Five == Zero)
  388. && ( TwoForty / Four - Five * Three * Four == Zero)
  389. && ( TwoForty / Five - Four * Three * Four == Zero),
  390. "5 != 4+1, 240/3 != 80, 240/4 != 60, or 240/5 != 48");
  391. if (ErrCnt[Failure] == 0) {
  392. printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
  393. printf("\n");
  394. }
  395. printf("Searching for Radix and Precision.\n");
  396. W = One;
  397. do {
  398. W = W + W;
  399. Y = W + One;
  400. Z = Y - W;
  401. Y = Z - One;
  402. } while (MinusOne + FABS(Y) < Zero);
  403. /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/
  404. Precision = Zero;
  405. Y = One;
  406. do {
  407. Radix = W + Y;
  408. Y = Y + Y;
  409. Radix = Radix - W;
  410. } while ( Radix == Zero);
  411. if (Radix < Two) Radix = One;
  412. printf("Radix = %f .\n", Radix);
  413. if (Radix != 1) {
  414. W = One;
  415. do {
  416. Precision = Precision + One;
  417. W = W * Radix;
  418. Y = W + One;
  419. } while ((Y - W) == One);
  420. }
  421. /*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1
  422. ...*/
  423. U1 = One / W;
  424. U2 = Radix * U1;
  425. printf("Closest relative separation found is U1 = %.7e .\n\n", U1);
  426. printf("Recalculating radix and precision\n ");
  427. /*save old values*/
  428. E0 = Radix;
  429. E1 = U1;
  430. E9 = U2;
  431. E3 = Precision;
  432. X = Four / Three;
  433. Third = X - One;
  434. F6 = Half - Third;
  435. X = F6 + F6;
  436. X = FABS(X - Third);
  437. if (X < U2) X = U2;
  438. /*... now X = (unknown no.) ulps of 1+...*/
  439. do {
  440. U2 = X;
  441. Y = Half * U2 + ThirtyTwo * U2 * U2;
  442. Y = One + Y;
  443. X = Y - One;
  444. } while ( ! ((U2 <= X) || (X <= Zero)));
  445. /*... now U2 == 1 ulp of 1 + ... */
  446. X = Two / Three;
  447. F6 = X - Half;
  448. Third = F6 + F6;
  449. X = Third - Half;
  450. X = FABS(X + F6);
  451. if (X < U1) X = U1;
  452. /*... now X == (unknown no.) ulps of 1 -... */
  453. do {
  454. U1 = X;
  455. Y = Half * U1 + ThirtyTwo * U1 * U1;
  456. Y = Half - Y;
  457. X = Half + Y;
  458. Y = Half - X;
  459. X = Half + Y;
  460. } while ( ! ((U1 <= X) || (X <= Zero)));
  461. /*... now U1 == 1 ulp of 1 - ... */
  462. if (U1 == E1) printf("confirms closest relative separation U1 .\n");
  463. else printf("gets better closest relative separation U1 = %.7e .\n", U1);
  464. W = One / U1;
  465. F9 = (Half - U1) + Half;
  466. Radix = FLOOR(0.01 + U2 / U1);
  467. if (Radix == E0) printf("Radix confirmed.\n");
  468. else printf("MYSTERY: recalculated Radix = %.7e .\n", Radix);
  469. TstCond (Defect, Radix <= Eight + Eight,
  470. "Radix is too big: roundoff problems");
  471. TstCond (Flaw, (Radix == Two) || (Radix == 10)
  472. || (Radix == One), "Radix is not as good as 2 or 10");
  473. /*=============================================*/
  474. Milestone = 20;
  475. /*=============================================*/
  476. TstCond (Failure, F9 - Half < Half,
  477. "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?");
  478. X = F9;
  479. I = 1;
  480. Y = X - Half;
  481. Z = Y - Half;
  482. TstCond (Failure, (X != One)
  483. || (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0");
  484. X = One + U2;
  485. I = 0;
  486. /*=============================================*/
  487. Milestone = 25;
  488. /*=============================================*/
  489. /*... BMinusU2 = nextafter(Radix, 0) */
  490. BMinusU2 = Radix - One;
  491. BMinusU2 = (BMinusU2 - U2) + One;
  492. /* Purify Integers */
  493. if (Radix != One) {
  494. X = - TwoForty * LOG(U1) / LOG(Radix);
  495. Y = FLOOR(Half + X);
  496. if (FABS(X - Y) * Four < One) X = Y;
  497. Precision = X / TwoForty;
  498. Y = FLOOR(Half + Precision);
  499. if (FABS(Precision - Y) * TwoForty < Half) Precision = Y;
  500. }
  501. if ((Precision != FLOOR(Precision)) || (Radix == One)) {
  502. printf("Precision cannot be characterized by an Integer number\n");
  503. printf("of significant digits but, by itself, this is a minor flaw.\n");
  504. }
  505. if (Radix == One)
  506. printf("logarithmic encoding has precision characterized solely by U1.\n");
  507. else printf("The number of significant digits of the Radix is %f .\n",
  508. Precision);
  509. TstCond (Serious, U2 * Nine * Nine * TwoForty < One,
  510. "Precision worse than 5 decimal figures ");
  511. /*=============================================*/
  512. Milestone = 30;
  513. /*=============================================*/
  514. /* Test for extra-precise subepressions */
  515. X = FABS(((Four / Three - One) - One / Four) * Three - One / Four);
  516. do {
  517. Z2 = X;
  518. X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
  519. } while ( ! ((Z2 <= X) || (X <= Zero)));
  520. X = Y = Z = FABS((Three / Four - Two / Three) * Three - One / Four);
  521. do {
  522. Z1 = Z;
  523. Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
  524. + One / Two)) + One / Two;
  525. } while ( ! ((Z1 <= Z) || (Z <= Zero)));
  526. do {
  527. do {
  528. Y1 = Y;
  529. Y = (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half
  530. )) + Half;
  531. } while ( ! ((Y1 <= Y) || (Y <= Zero)));
  532. X1 = X;
  533. X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
  534. } while ( ! ((X1 <= X) || (X <= Zero)));
  535. if ((X1 != Y1) || (X1 != Z1)) {
  536. BadCond(Serious, "Disagreements among the values X1, Y1, Z1,\n");
  537. printf("respectively %.7e, %.7e, %.7e,\n", X1, Y1, Z1);
  538. printf("are symptoms of inconsistencies introduced\n");
  539. printf("by extra-precise evaluation of arithmetic subexpressions.\n");
  540. notify("Possibly some part of this");
  541. if ((X1 == U1) || (Y1 == U1) || (Z1 == U1)) printf(
  542. "That feature is not tested further by this program.\n") ;
  543. }
  544. else {
  545. if ((Z1 != U1) || (Z2 != U2)) {
  546. if ((Z1 >= U1) || (Z2 >= U2)) {
  547. BadCond(Failure, "");
  548. notify("Precision");
  549. printf("\tU1 = %.7e, Z1 - U1 = %.7e\n",U1,Z1-U1);
  550. printf("\tU2 = %.7e, Z2 - U2 = %.7e\n",U2,Z2-U2);
  551. }
  552. else {
  553. if ((Z1 <= Zero) || (Z2 <= Zero)) {
  554. printf("Because of unusual Radix = %f", Radix);
  555. printf(", or exact rational arithmetic a result\n");
  556. printf("Z1 = %.7e, or Z2 = %.7e ", Z1, Z2);
  557. notify("of an\nextra-precision");
  558. }
  559. if (Z1 != Z2 || Z1 > Zero) {
  560. X = Z1 / U1;
  561. Y = Z2 / U2;
  562. if (Y > X) X = Y;
  563. Q = - LOG(X);
  564. printf("Some subexpressions appear to be calculated extra\n");
  565. printf("precisely with about %g extra B-digits, i.e.\n",
  566. (Q / LOG(Radix)));
  567. printf("roughly %g extra significant decimals.\n",
  568. Q / LOG(10.));
  569. }
  570. printf("That feature is not tested further by this program.\n");
  571. }
  572. }
  573. }
  574. Pause();
  575. /*=============================================*/
  576. /*SPLIT
  577. }
  578. #include "paranoia.h"
  579. part3(){
  580. */
  581. Milestone = 35;
  582. /*=============================================*/
  583. if (Radix >= Two) {
  584. X = W / (Radix * Radix);
  585. Y = X + One;
  586. Z = Y - X;
  587. T = Z + U2;
  588. X = T - Z;
  589. TstCond (Failure, X == U2,
  590. "Subtraction is not normalized X=Y,X+Z != Y+Z!");
  591. if (X == U2) printf(
  592. "Subtraction appears to be normalized, as it should be.");
  593. }
  594. printf("\nChecking for guard digit in *, /, and -.\n");
  595. Y = F9 * One;
  596. Z = One * F9;
  597. X = F9 - Half;
  598. Y = (Y - Half) - X;
  599. Z = (Z - Half) - X;
  600. X = One + U2;
  601. T = X * Radix;
  602. R = Radix * X;
  603. X = T - Radix;
  604. X = X - Radix * U2;
  605. T = R - Radix;
  606. T = T - Radix * U2;
  607. X = X * (Radix - One);
  608. T = T * (Radix - One);
  609. if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)) GMult = Yes;
  610. else {
  611. GMult = No;
  612. TstCond (Serious, False,
  613. "* lacks a Guard Digit, so 1*X != X");
  614. }
  615. Z = Radix * U2;
  616. X = One + Z;
  617. Y = FABS((X + Z) - X * X) - U2;
  618. X = One - U2;
  619. Z = FABS((X - U2) - X * X) - U1;
  620. TstCond (Failure, (Y <= Zero)
  621. && (Z <= Zero), "* gets too many final digits wrong.\n");
  622. Y = One - U2;
  623. X = One + U2;
  624. Z = One / Y;
  625. Y = Z - X;
  626. X = One / Three;
  627. Z = Three / Nine;
  628. X = X - Z;
  629. T = Nine / TwentySeven;
  630. Z = Z - T;
  631. TstCond(Defect, X == Zero && Y == Zero && Z == Zero,
  632. "Division lacks a Guard Digit, so error can exceed 1 ulp\nor 1/3 and 3/9 and 9/27 may disagree");
  633. Y = F9 / One;
  634. X = F9 - Half;
  635. Y = (Y - Half) - X;
  636. X = One + U2;
  637. T = X / One;
  638. X = T - X;
  639. if ((X == Zero) && (Y == Zero) && (Z == Zero)) GDiv = Yes;
  640. else {
  641. GDiv = No;
  642. TstCond (Serious, False,
  643. "Division lacks a Guard Digit, so X/1 != X");
  644. }
  645. X = One / (One + U2);
  646. Y = X - Half - Half;
  647. TstCond (Serious, Y < Zero,
  648. "Computed value of 1/1.000..1 >= 1");
  649. X = One - U2;
  650. Y = One + Radix * U2;
  651. Z = X * Radix;
  652. T = Y * Radix;
  653. R = Z / Radix;
  654. StickyBit = T / Radix;
  655. X = R - X;
  656. Y = StickyBit - Y;
  657. TstCond (Failure, X == Zero && Y == Zero,
  658. "* and/or / gets too many last digits wrong");
  659. Y = One - U1;
  660. X = One - F9;
  661. Y = One - Y;
  662. T = Radix - U2;
  663. Z = Radix - BMinusU2;
  664. T = Radix - T;
  665. if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2)) GAddSub = Yes;
  666. else {
  667. GAddSub = No;
  668. TstCond (Serious, False,
  669. "- lacks Guard Digit, so cancellation is obscured");
  670. }
  671. if (F9 != One && F9 - One >= Zero) {
  672. BadCond(Serious, "comparison alleges (1-U1) < 1 although\n");
  673. printf(" subtraction yields (1-U1) - 1 = 0 , thereby vitiating\n");
  674. printf(" such precautions against division by zero as\n");
  675. printf(" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
  676. }
  677. if (GMult == Yes && GDiv == Yes && GAddSub == Yes) printf(
  678. " *, /, and - appear to have guard digits, as they should.\n");
  679. /*=============================================*/
  680. Milestone = 40;
  681. /*=============================================*/
  682. Pause();
  683. printf("Checking rounding on multiply, divide and add/subtract.\n");
  684. RMult = Other;
  685. RDiv = Other;
  686. RAddSub = Other;
  687. RadixD2 = Radix / Two;
  688. A1 = Two;
  689. Done = False;
  690. do {
  691. AInvrse = Radix;
  692. do {
  693. X = AInvrse;
  694. AInvrse = AInvrse / A1;
  695. } while ( ! (FLOOR(AInvrse) != AInvrse));
  696. Done = (X == One) || (A1 > Three);
  697. if (! Done) A1 = Nine + One;
  698. } while ( ! (Done));
  699. if (X == One) A1 = Radix;
  700. AInvrse = One / A1;
  701. X = A1;
  702. Y = AInvrse;
  703. Done = False;
  704. do {
  705. Z = X * Y - Half;
  706. TstCond (Failure, Z == Half,
  707. "X * (1/X) differs from 1");
  708. Done = X == Radix;
  709. X = Radix;
  710. Y = One / X;
  711. } while ( ! (Done));
  712. Y2 = One + U2;
  713. Y1 = One - U2;
  714. X = OneAndHalf - U2;
  715. Y = OneAndHalf + U2;
  716. Z = (X - U2) * Y2;
  717. T = Y * Y1;
  718. Z = Z - X;
  719. T = T - X;
  720. X = X * Y2;
  721. Y = (Y + U2) * Y1;
  722. X = X - OneAndHalf;
  723. Y = Y - OneAndHalf;
  724. if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) {
  725. X = (OneAndHalf + U2) * Y2;
  726. Y = OneAndHalf - U2 - U2;
  727. Z = OneAndHalf + U2 + U2;
  728. T = (OneAndHalf - U2) * Y1;
  729. X = X - (Z + U2);
  730. StickyBit = Y * Y1;
  731. S = Z * Y2;
  732. T = T - Y;
  733. Y = (U2 - Y) + StickyBit;
  734. Z = S - (Z + U2 + U2);
  735. StickyBit = (Y2 + U2) * Y1;
  736. Y1 = Y2 * Y1;
  737. StickyBit = StickyBit - Y2;
  738. Y1 = Y1 - Half;
  739. if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
  740. && ( StickyBit == Zero) && (Y1 == Half)) {
  741. RMult = Rounded;
  742. printf("Multiplication appears to round correctly.\n");
  743. }
  744. else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
  745. && (T < Zero) && (StickyBit + U2 == Zero)
  746. && (Y1 < Half)) {
  747. RMult = Chopped;
  748. printf("Multiplication appears to chop.\n");
  749. }
  750. else printf("* is neither chopped nor correctly rounded.\n");
  751. if ((RMult == Rounded) && (GMult == No)) notify("Multiplication");
  752. }
  753. else printf("* is neither chopped nor correctly rounded.\n");
  754. /*=============================================*/
  755. Milestone = 45;
  756. /*=============================================*/
  757. Y2 = One + U2;
  758. Y1 = One - U2;
  759. Z = OneAndHalf + U2 + U2;
  760. X = Z / Y2;
  761. T = OneAndHalf - U2 - U2;
  762. Y = (T - U2) / Y1;
  763. Z = (Z + U2) / Y2;
  764. X = X - OneAndHalf;
  765. Y = Y - T;
  766. T = T / Y1;
  767. Z = Z - (OneAndHalf + U2);
  768. T = (U2 - OneAndHalf) + T;
  769. if (! ((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) {
  770. X = OneAndHalf / Y2;
  771. Y = OneAndHalf - U2;
  772. Z = OneAndHalf + U2;
  773. X = X - Y;
  774. T = OneAndHalf / Y1;
  775. Y = Y / Y1;
  776. T = T - (Z + U2);
  777. Y = Y - Z;
  778. Z = Z / Y2;
  779. Y1 = (Y2 + U2) / Y2;
  780. Z = Z - OneAndHalf;
  781. Y2 = Y1 - Y2;
  782. Y1 = (F9 - U1) / F9;
  783. if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
  784. && (Y2 == Zero) && (Y2 == Zero)
  785. && (Y1 - Half == F9 - Half )) {
  786. RDiv = Rounded;
  787. printf("Division appears to round correctly.\n");
  788. if (GDiv == No) notify("Division");
  789. }
  790. else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
  791. && (Y2 < Zero) && (Y1 - Half < F9 - Half)) {
  792. RDiv = Chopped;
  793. printf("Division appears to chop.\n");
  794. }
  795. }
  796. if (RDiv == Other) printf("/ is neither chopped nor correctly rounded.\n");
  797. BInvrse = One / Radix;
  798. TstCond (Failure, (BInvrse * Radix - Half == Half),
  799. "Radix * ( 1 / Radix ) differs from 1");
  800. /*=============================================*/
  801. /*SPLIT
  802. }
  803. #include "paranoia.h"
  804. part4(){
  805. */
  806. Milestone = 50;
  807. /*=============================================*/
  808. TstCond (Failure, ((F9 + U1) - Half == Half)
  809. && ((BMinusU2 + U2 ) - One == Radix - One),
  810. "Incomplete carry-propagation in Addition");
  811. X = One - U1 * U1;
  812. Y = One + U2 * (One - U2);
  813. Z = F9 - Half;
  814. X = (X - Half) - Z;
  815. Y = Y - One;
  816. if ((X == Zero) && (Y == Zero)) {
  817. RAddSub = Chopped;
  818. printf("Add/Subtract appears to be chopped.\n");
  819. }
  820. if (GAddSub == Yes) {
  821. X = (Half + U2) * U2;
  822. Y = (Half - U2) * U2;
  823. X = One + X;
  824. Y = One + Y;
  825. X = (One + U2) - X;
  826. Y = One - Y;
  827. if ((X == Zero) && (Y == Zero)) {
  828. X = (Half + U2) * U1;
  829. Y = (Half - U2) * U1;
  830. X = One - X;
  831. Y = One - Y;
  832. X = F9 - X;
  833. Y = One - Y;
  834. if ((X == Zero) && (Y == Zero)) {
  835. RAddSub = Rounded;
  836. printf("Addition/Subtraction appears to round correctly.\n");
  837. if (GAddSub == No) notify("Add/Subtract");
  838. }
  839. else printf("Addition/Subtraction neither rounds nor chops.\n");
  840. }
  841. else printf("Addition/Subtraction neither rounds nor chops.\n");
  842. }
  843. else printf("Addition/Subtraction neither rounds nor chops.\n");
  844. S = One;
  845. X = One + Half * (One + Half);
  846. Y = (One + U2) * Half;
  847. Z = X - Y;
  848. T = Y - X;
  849. StickyBit = Z + T;
  850. if (StickyBit != Zero) {
  851. S = Zero;
  852. BadCond(Flaw, "(X - Y) + (Y - X) is non zero!\n");
  853. }
  854. StickyBit = Zero;
  855. if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
  856. && (RMult == Rounded) && (RDiv == Rounded)
  857. && (RAddSub == Rounded) && (FLOOR(RadixD2) == RadixD2)) {
  858. printf("Checking for sticky bit.\n");
  859. X = (Half + U1) * U2;
  860. Y = Half * U2;
  861. Z = One + Y;
  862. T = One + X;
  863. if ((Z - One <= Zero) && (T - One >= U2)) {
  864. Z = T + Y;
  865. Y = Z - X;
  866. if ((Z - T >= U2) && (Y - T == Zero)) {
  867. X = (Half + U1) * U1;
  868. Y = Half * U1;
  869. Z = One - Y;
  870. T = One - X;
  871. if ((Z - One == Zero) && (T - F9 == Zero)) {
  872. Z = (Half - U1) * U1;
  873. T = F9 - Z;
  874. Q = F9 - Y;
  875. if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) {
  876. Z = (One + U2) * OneAndHalf;
  877. T = (OneAndHalf + U2) - Z + U2;
  878. X = One + Half / Radix;
  879. Y = One + Radix * U2;
  880. Z = X * Y;
  881. if (T == Zero && X + Radix * U2 - Z == Zero) {
  882. if (Radix != Two) {
  883. X = Two + U2;
  884. Y = X / Two;
  885. if ((Y - One == Zero)) StickyBit = S;
  886. }
  887. else StickyBit = S;
  888. }
  889. }
  890. }
  891. }
  892. }
  893. }
  894. if (StickyBit == One) printf("Sticky bit apparently used correctly.\n");
  895. else printf("Sticky bit used incorrectly or not at all.\n");
  896. TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
  897. RMult == Other || RDiv == Other || RAddSub == Other),
  898. "lack(s) of guard digits or failure(s) to correctly round or chop\n(noted above) count as one flaw in the final tally below");
  899. /*=============================================*/
  900. Milestone = 60;
  901. /*=============================================*/
  902. printf("\n");
  903. printf("Does Multiplication commute? ");
  904. printf("Testing on %d random pairs.\n", NoTrials);
  905. Random9 = SQRT(3.0);
  906. Random1 = Third;
  907. I = 1;
  908. do {
  909. X = Random();
  910. Y = Random();
  911. Z9 = Y * X;
  912. Z = X * Y;
  913. Z9 = Z - Z9;
  914. I = I + 1;
  915. } while ( ! ((I > NoTrials) || (Z9 != Zero)));
  916. if (I == NoTrials) {
  917. Random1 = One + Half / Three;
  918. Random2 = (U2 + U1) + One;
  919. Z = Random1 * Random2;
  920. Y = Random2 * Random1;
  921. Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
  922. Three) * ((U2 + U1) + One);
  923. }
  924. if (! ((I == NoTrials) || (Z9 == Zero)))
  925. BadCond(Defect, "X * Y == Y * X trial fails.\n");
  926. else printf(" No failures found in %d integer pairs.\n", NoTrials);
  927. /*=============================================*/
  928. Milestone = 70;
  929. /*=============================================*/
  930. printf("\nRunning test of square root(x).\n");
  931. TstCond (Failure, (Zero == SQRT(Zero))
  932. && (- Zero == SQRT(- Zero))
  933. && (One == SQRT(One)), "Square root of 0.0, -0.0 or 1.0 wrong");
  934. MinSqEr = Zero;
  935. MaxSqEr = Zero;
  936. J = Zero;
  937. X = Radix;
  938. OneUlp = U2;
  939. SqXMinX (Serious);
  940. X = BInvrse;
  941. OneUlp = BInvrse * U1;
  942. SqXMinX (Serious);
  943. X = U1;
  944. OneUlp = U1 * U1;
  945. SqXMinX (Serious);
  946. if (J != Zero) Pause();
  947. printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
  948. J = Zero;
  949. X = Two;
  950. Y = Radix;
  951. if ((Radix != One)) do {
  952. X = Y;
  953. Y = Radix * Y;
  954. } while ( ! ((Y - X >= NoTrials)));
  955. OneUlp = X * U2;
  956. I = 1;
  957. while (I <= NoTrials) {
  958. X = X + One;
  959. SqXMinX (Defect);
  960. if (J > Zero) break;
  961. I = I + 1;
  962. }
  963. printf("Test for sqrt monotonicity.\n");
  964. I = - 1;
  965. X = BMinusU2;
  966. Y = Radix;
  967. Z = Radix + Radix * U2;
  968. NotMonot = False;
  969. Monot = False;
  970. while ( ! (NotMonot || Monot)) {
  971. I = I + 1;
  972. X = SQRT(X);
  973. Q = SQRT(Y);
  974. Z = SQRT(Z);
  975. if ((X > Q) || (Q > Z)) NotMonot = True;
  976. else {
  977. Q = FLOOR(Q + Half);
  978. if ((I > 0) || (Radix == Q * Q)) Monot = True;
  979. else if (I > 0) {
  980. if (I > 1) Monot = True;
  981. else {
  982. Y = Y * BInvrse;
  983. X = Y - U1;
  984. Z = Y + U1;
  985. }
  986. }
  987. else {
  988. Y = Q;
  989. X = Y - U2;
  990. Z = Y + U2;
  991. }
  992. }
  993. }
  994. if (Monot) printf("sqrt has passed a test for Monotonicity.\n");
  995. else {
  996. BadCond(Defect, "");
  997. printf("sqrt(X) is non-monotonic for X near %.7e .\n", Y);
  998. }
  999. /*=============================================*/
  1000. /*SPLIT
  1001. }
  1002. #include "paranoia.h"
  1003. part5(){
  1004. */
  1005. Milestone = 80;
  1006. /*=============================================*/
  1007. MinSqEr = MinSqEr + Half;
  1008. MaxSqEr = MaxSqEr - Half;
  1009. Y = (SQRT(One + U2) - One) / U2;
  1010. SqEr = (Y - One) + U2 / Eight;
  1011. if (SqEr > MaxSqEr) MaxSqEr = SqEr;
  1012. SqEr = Y + U2 / Eight;
  1013. if (SqEr < MinSqEr) MinSqEr = SqEr;
  1014. Y = ((SQRT(F9) - U2) - (One - U2)) / U1;
  1015. SqEr = Y + U1 / Eight;
  1016. if (SqEr > MaxSqEr) MaxSqEr = SqEr;
  1017. SqEr = (Y + One) + U1 / Eight;
  1018. if (SqEr < MinSqEr) MinSqEr = SqEr;
  1019. OneUlp = U2;
  1020. X = OneUlp;
  1021. for( Indx = 1; Indx <= 3; ++Indx) {
  1022. Y = SQRT((X + U1 + X) + F9);
  1023. Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;
  1024. Z = ((U1 - X) + F9) * Half * X * X / OneUlp;
  1025. SqEr = (Y + Half) + Z;
  1026. if (SqEr < MinSqEr) MinSqEr = SqEr;
  1027. SqEr = (Y - Half) + Z;
  1028. if (SqEr > MaxSqEr) MaxSqEr = SqEr;
  1029. if (((Indx == 1) || (Indx == 3)))
  1030. X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));
  1031. else {
  1032. OneUlp = U1;
  1033. X = - OneUlp;
  1034. }
  1035. }
  1036. /*=============================================*/
  1037. Milestone = 85;
  1038. /*=============================================*/
  1039. SqRWrng = False;
  1040. Anomaly = False;
  1041. RSqrt = Other; /* ~dgh */
  1042. if (Radix != One) {
  1043. printf("Testing whether sqrt is rounded or chopped.\n");
  1044. D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));
  1045. /* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
  1046. X = D / Radix;
  1047. Y = D / A1;
  1048. if ((X != FLOOR(X)) || (Y != FLOOR(Y))) {
  1049. Anomaly = True;
  1050. }
  1051. else {
  1052. X = Zero;
  1053. Z2 = X;
  1054. Y = One;
  1055. Y2 = Y;
  1056. Z1 = Radix - One;
  1057. FourD = Four * D;
  1058. do {
  1059. if (Y2 > Z2) {
  1060. Q = Radix;
  1061. Y1 = Y;
  1062. do {
  1063. X1 = FABS(Q + FLOOR(Half - Q / Y1) * Y1);
  1064. Q = Y1;
  1065. Y1 = X1;
  1066. } while ( ! (X1 <= Zero));
  1067. if (Q <= One) {
  1068. Z2 = Y2;
  1069. Z = Y;
  1070. }
  1071. }
  1072. Y = Y + Two;
  1073. X = X + Eight;
  1074. Y2 = Y2 + X;
  1075. if (Y2 >= FourD) Y2 = Y2 - FourD;
  1076. } while ( ! (Y >= D));
  1077. X8 = FourD - Z2;
  1078. Q = (X8 + Z * Z) / FourD;
  1079. X8 = X8 / Eight;
  1080. if (Q != FLOOR(Q)) Anomaly = True;
  1081. else {
  1082. Break = False;
  1083. do {
  1084. X = Z1 * Z;
  1085. X = X - FLOOR(X / Radix) * Radix;
  1086. if (X == One)
  1087. Break = True;
  1088. else
  1089. Z1 = Z1 - One;
  1090. } while ( ! (Break || (Z1 <= Zero)));
  1091. if ((Z1 <= Zero) && (! Break)) Anomaly = True;
  1092. else {
  1093. if (Z1 > RadixD2) Z1 = Z1 - Radix;
  1094. do {
  1095. NewD();
  1096. } while ( ! (U2 * D >= F9));
  1097. if (D * Radix - D != W - D) Anomaly = True;
  1098. else {
  1099. Z2 = D;
  1100. I = 0;
  1101. Y = D + (One + Z) * Half;
  1102. X = D + Z + Q;
  1103. SR3750();
  1104. Y = D + (One - Z) * Half + D;
  1105. X = D - Z + D;
  1106. X = X + Q + X;
  1107. SR3750();
  1108. NewD();
  1109. if (D - Z2 != W - Z2) Anomaly = True;
  1110. else {
  1111. Y = (D - Z2) + (Z2 + (One - Z) * Half);
  1112. X = (D - Z2) + (Z2 - Z + Q);
  1113. SR3750();
  1114. Y = (One + Z) * Half;
  1115. X = Q;
  1116. SR3750();
  1117. if (I == 0) Anomaly = True;
  1118. }
  1119. }
  1120. }
  1121. }
  1122. }
  1123. if ((I == 0) || Anomaly) {
  1124. BadCond(Failure, "Anomalous arithmetic with Integer < ");
  1125. printf("Radix^Precision = %.7e\n", W);
  1126. printf(" fails test whether sqrt rounds or chops.\n");
  1127. SqRWrng = True;
  1128. }
  1129. }
  1130. if (! Anomaly) {
  1131. if (! ((MinSqEr < Zero) || (MaxSqEr > Zero))) {
  1132. RSqrt = Rounded;
  1133. printf("Square root appears to be correctly rounded.\n");
  1134. }
  1135. else {
  1136. if ((MaxSqEr + U2 > U2 - Half) || (MinSqEr > Half)
  1137. || (MinSqEr + Radix < Half)) SqRWrng = True;
  1138. else {
  1139. RSqrt = Chopped;
  1140. printf("Square root appears to be chopped.\n");
  1141. }
  1142. }
  1143. }
  1144. if (SqRWrng) {
  1145. printf("Square root is neither chopped nor correctly rounded.\n");
  1146. printf("Observed errors run from %.7e ", MinSqEr - Half);
  1147. printf("to %.7e ulps.\n", Half + MaxSqEr);
  1148. TstCond (Serious, MaxSqEr - MinSqEr < Radix * Radix,
  1149. "sqrt gets too many last digits wrong");
  1150. }
  1151. /*=============================================*/
  1152. Milestone = 90;
  1153. /*=============================================*/
  1154. Pause();
  1155. printf("Testing powers Z^i for small Integers Z and i.\n");
  1156. N = 0;
  1157. /* ... test powers of zero. */
  1158. I = 0;
  1159. Z = -Zero;
  1160. M = 3.0;
  1161. Break = False;
  1162. do {
  1163. X = One;
  1164. SR3980();
  1165. if (I <= 10) {
  1166. I = 1023;
  1167. SR3980();
  1168. }
  1169. if (Z == MinusOne) Break = True;
  1170. else {
  1171. Z = MinusOne;
  1172. PrintIfNPositive();
  1173. N = 0;
  1174. /* .. if(-1)^N is invalid, replace MinusOne by One. */
  1175. I = - 4;
  1176. }
  1177. } while ( ! Break);
  1178. PrintIfNPositive();
  1179. N1 = N;
  1180. N = 0;
  1181. Z = A1;
  1182. M = FLOOR(Two * LOG(W) / LOG(A1));
  1183. Break = False;
  1184. do {
  1185. X = Z;
  1186. I = 1;
  1187. SR3980();
  1188. if (Z == AInvrse) Break = True;
  1189. else Z = AInvrse;
  1190. } while ( ! (Break));
  1191. /*=============================================*/
  1192. Milestone = 100;
  1193. /*=============================================*/
  1194. /* Powers of Radix have been tested, */
  1195. /* next try a few primes */
  1196. M = NoTrials;
  1197. Z = Three;
  1198. do {
  1199. X = Z;
  1200. I = 1;
  1201. SR3980();
  1202. do {
  1203. Z = Z + Two;
  1204. } while ( Three * FLOOR(Z / Three) == Z );
  1205. } while ( Z < Eight * Three );
  1206. if (N > 0) {
  1207. printf("Errors like this may invalidate financial calculations\n");
  1208. printf("\tinvolving interest rates.\n");
  1209. }
  1210. PrintIfNPositive();
  1211. N += N1;
  1212. if (N == 0) printf("... no discrepancis found.\n");
  1213. if (N > 0) Pause();
  1214. else printf("\n");
  1215. /*=============================================*/
  1216. /*SPLIT
  1217. }
  1218. #include "paranoia.h"
  1219. part6(){
  1220. */
  1221. Milestone = 110;
  1222. /*=============================================*/
  1223. printf("Seeking Underflow thresholds UfThold and E0.\n");
  1224. D = U1;
  1225. if (Precision != FLOOR(Precision)) {
  1226. D = BInvrse;
  1227. X = Precision;
  1228. do {
  1229. D = D * BInvrse;
  1230. X = X - One;
  1231. } while ( X > Zero);
  1232. }
  1233. Y = One;
  1234. Z = D;
  1235. /* ... D is power of 1/Radix < 1. */
  1236. do {
  1237. C = Y;
  1238. Y = Z;
  1239. Z = Y * Y;
  1240. } while ((Y > Z) && (Z + Z > Z));
  1241. Y = C;
  1242. Z = Y * D;
  1243. do {
  1244. C = Y;
  1245. Y = Z;
  1246. Z = Y * D;
  1247. } while ((Y > Z) && (Z + Z > Z));
  1248. if (Radix < Two) HInvrse = Two;
  1249. else HInvrse = Radix;
  1250. H = One / HInvrse;
  1251. /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
  1252. CInvrse = One / C;
  1253. E0 = C;
  1254. Z = E0 * H;
  1255. /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
  1256. do {
  1257. Y = E0;
  1258. E0 = Z;
  1259. Z = E0 * H;
  1260. } while ((E0 > Z) && (Z + Z > Z));
  1261. UfThold = E0;
  1262. E1 = Zero;
  1263. Q = Zero;
  1264. E9 = U2;
  1265. S = One + E9;
  1266. D = C * S;
  1267. if (D <= C) {
  1268. E9 = Radix * U2;
  1269. S = One + E9;
  1270. D = C * S;
  1271. if (D <= C) {
  1272. BadCond(Failure, "multiplication gets too many last digits wrong.\n");
  1273. Underflow = E0;
  1274. Y1 = Zero;
  1275. PseudoZero = Z;
  1276. Pause();
  1277. }
  1278. }
  1279. else {
  1280. Underflow = D;
  1281. PseudoZero = Underflow * H;
  1282. UfThold = Zero;
  1283. do {
  1284. Y1 = Underflow;
  1285. Underflow = PseudoZero;
  1286. if (E1 + E1 <= E1) {
  1287. Y2 = Underflow * HInvrse;
  1288. E1 = FABS(Y1 - Y2);
  1289. Q = Y1;
  1290. if ((UfThold == Zero) && (Y1 != Y2)) UfThold = Y1;
  1291. }
  1292. PseudoZero = PseudoZero * H;
  1293. } while ((Underflow > PseudoZero)
  1294. && (PseudoZero + PseudoZero > PseudoZero));
  1295. }
  1296. /* Comment line 4530 .. 4560 */
  1297. if (PseudoZero != Zero) {
  1298. printf("\n");
  1299. Z = PseudoZero;
  1300. /* ... Test PseudoZero for "phoney- zero" violates */
  1301. /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
  1302. ... */
  1303. if (PseudoZero <= Zero) {
  1304. BadCond(Failure, "Positive expressions can underflow to an\n");
  1305. printf("allegedly negative value\n");
  1306. printf("PseudoZero that prints out as: %g .\n", PseudoZero);
  1307. X = - PseudoZero;
  1308. if (X <= Zero) {
  1309. printf("But -PseudoZero, which should be\n");
  1310. printf("positive, isn't; it prints out as %g .\n", X);
  1311. }
  1312. }
  1313. else {
  1314. BadCond(Flaw, "Underflow can stick at an allegedly positive\n");
  1315. printf("value PseudoZero that prints out as %g .\n", PseudoZero);
  1316. }
  1317. TstPtUf();
  1318. }
  1319. /*=============================================*/
  1320. Milestone = 120;
  1321. /*=============================================*/
  1322. if (CInvrse * Y > CInvrse * Y1) {
  1323. S = H * S;
  1324. E0 = Underflow;
  1325. }
  1326. if (! ((E1 == Zero) || (E1 == E0))) {
  1327. BadCond(Defect, "");
  1328. if (E1 < E0) {
  1329. printf("Products underflow at a higher");
  1330. printf(" threshold than differences.\n");
  1331. if (PseudoZero == Zero)
  1332. E0 = E1;
  1333. }
  1334. else {
  1335. printf("Difference underflows at a higher");
  1336. printf(" threshold than products.\n");
  1337. }
  1338. }
  1339. printf("Smallest strictly positive number found is E0 = %g .\n", E0);
  1340. Z = E0;
  1341. TstPtUf();
  1342. Underflow = E0;
  1343. if (N == 1) Underflow = Y;
  1344. I = 4;
  1345. if (E1 == Zero) I = 3;
  1346. if (UfThold == Zero) I = I - 2;
  1347. UfNGrad = True;
  1348. switch (I) {
  1349. case 1:
  1350. UfThold = Underflow;
  1351. if ((CInvrse * Q) != ((CInvrse * Y) * S)) {
  1352. UfThold = Y;
  1353. BadCond(Failure, "Either accuracy deteriorates as numbers\n");
  1354. printf("approach a threshold = %.17e\n", UfThold);;
  1355. printf(" coming down from %.17e\n", C);
  1356. printf(" or else multiplication gets too many last digits wrong.\n");
  1357. }
  1358. Pause();
  1359. break;
  1360. case 2:
  1361. BadCond(Failure, "Underflow confuses Comparison, which alleges that\n");
  1362. printf("Q == Y while denying that |Q - Y| == 0; these values\n");
  1363. printf("print out as Q = %.17e, Y = %.17e .\n", Q, Y2);
  1364. printf ("|Q - Y| = %.17e .\n" , FABS(Q - Y2));
  1365. UfThold = Q;
  1366. break;
  1367. case 3:
  1368. X = X;
  1369. break;
  1370. case 4:
  1371. if ((Q == UfThold) && (E1 == E0)
  1372. && (FABS( UfThold - E1 / E9) <= E1)) {
  1373. UfNGrad = False;
  1374. printf("Underflow is gradual; it incurs Absolute Error =\n");
  1375. printf("(roundoff in UfThold) < E0.\n");
  1376. Y = E0 * CInvrse;
  1377. Y = Y * (OneAndHalf + U2);
  1378. X = CInvrse * (One + U2);
  1379. Y = Y / X;
  1380. IEEE = (Y == E0);
  1381. }
  1382. }
  1383. if (UfNGrad) {
  1384. printf("\n");
  1385. sigsave = sigfpe;
  1386. if (setjmp(ovfl_buf)) {
  1387. printf("Underflow / UfThold failed!\n");
  1388. R = H + H;
  1389. }
  1390. else R = SQRT(Underflow / UfThold);
  1391. sigsave = 0;
  1392. if (R <= H) {
  1393. Z = R * UfThold;
  1394. X = Z * (One + R * H * (One + H));
  1395. }
  1396. else {
  1397. Z = UfThold;
  1398. X = Z * (One + H * H * (One + H));
  1399. }
  1400. if (! ((X == Z) || (X - Z != Zero))) {
  1401. BadCond(Flaw, "");
  1402. printf("X = %.17e\n\tis not equal to Z = %.17e .\n", X, Z);
  1403. Z9 = X - Z;
  1404. printf("yet X - Z yields %.17e .\n", Z9);
  1405. printf(" Should this NOT signal Underflow, ");
  1406. printf("this is a SERIOUS DEFECT\nthat causes ");
  1407. printf("confusion when innocent statements like\n");;
  1408. printf(" if (X == Z) ... else");
  1409. printf(" ... (f(X) - f(Z)) / (X - Z) ...\n");
  1410. printf("encounter Division by Zero although actually\n");
  1411. sigsave = sigfpe;
  1412. if (setjmp(ovfl_buf)) printf("X / Z fails!\n");
  1413. else printf("X / Z = 1 + %g .\n", (X / Z - Half) - Half);
  1414. sigsave = 0;
  1415. }
  1416. }
  1417. printf("The Underflow threshold is %.17e, %s\n", UfThold,
  1418. " below which");
  1419. printf("calculation may suffer larger Relative error than ");
  1420. printf("merely roundoff.\n");
  1421. Y2 = U1 * U1;
  1422. Y = Y2 * Y2;
  1423. Y2 = Y * U1;
  1424. if (Y2 <= UfThold) {
  1425. if (Y > E0) {
  1426. BadCond(Defect, "");
  1427. I = 5;
  1428. }
  1429. else {
  1430. BadCond(Serious, "");
  1431. I = 4;
  1432. }
  1433. printf("Range is too narrow; U1^%d Underflows.\n", I);
  1434. }
  1435. /*=============================================*/
  1436. /*SPLIT
  1437. }
  1438. #include "paranoia.h"
  1439. part7(){
  1440. */
  1441. Milestone = 130;
  1442. /*=============================================*/
  1443. Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;
  1444. Y2 = Y + Y;
  1445. printf("Since underflow occurs below the threshold\n");
  1446. printf("UfThold = (%.17e) ^ (%.17e)\nonly underflow ", HInvrse, Y);
  1447. printf("should afflict the expression\n\t(%.17e) ^ (%.17e);\n", HInvrse, Y);
  1448. V9 = POW(HInvrse, Y2);
  1449. printf("actually calculating yields: %.17e .\n", V9);
  1450. if (! ((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) {
  1451. BadCond(Serious, "this is not between 0 and underflow\n");
  1452. printf(" threshold = %.17e .\n", UfThold);
  1453. }
  1454. else if (! (V9 > UfThold * (One + E9)))
  1455. printf("This computed value is O.K.\n");
  1456. else {
  1457. BadCond(Defect, "this is not between 0 and underflow\n");
  1458. printf(" threshold = %.17e .\n", UfThold);
  1459. }
  1460. /*=============================================*/
  1461. Milestone = 140;
  1462. /*=============================================*/
  1463. printf("\n");
  1464. /* ...calculate Exp2 == exp(2) == 7.389056099... */
  1465. X = Zero;
  1466. I = 2;
  1467. Y = Two * Three;
  1468. Q = Zero;
  1469. N = 0;
  1470. do {
  1471. Z = X;
  1472. I = I + 1;
  1473. Y = Y / (I + I);
  1474. R = Y + Q;
  1475. X = Z + R;
  1476. Q = (Z - X) + R;
  1477. } while(X > Z);
  1478. Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);
  1479. X = Z * Z;
  1480. Exp2 = X * X;
  1481. X = F9;
  1482. Y = X - U1;
  1483. printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17e as X -> 1.\n",
  1484. Exp2);
  1485. for(I = 1;;) {
  1486. Z = X - BInvrse;
  1487. Z = (X + One) / (Z - (One - BInvrse));
  1488. Q = POW(X, Z) - Exp2;
  1489. if (FABS(Q) > TwoForty * U2) {
  1490. N = 1;
  1491. V9 = (X - BInvrse) - (One - BInvrse);
  1492. BadCond(Defect, "Calculated");
  1493. printf(" %.17e for\n", POW(X,Z));
  1494. printf("\t(1 + (%.17e) ^ (%.17e);\n", V9, Z);
  1495. printf("\tdiffers from correct value by %.17e .\n", Q);
  1496. printf("\tThis much error may spoil financial\n");
  1497. printf("\tcalculations involving tiny interest rates.\n");
  1498. break;
  1499. }
  1500. else {
  1501. Z = (Y - X) * Two + Y;
  1502. X = Y;
  1503. Y = Z;
  1504. Z = One + (X - F9)*(X - F9);
  1505. if (Z > One && I < NoTrials) I++;
  1506. else {
  1507. if (X > One) {
  1508. if (N == 0)
  1509. printf("Accuracy seems adequate.\n");
  1510. break;
  1511. }
  1512. else {
  1513. X = One + U2;
  1514. Y = U2 + U2;
  1515. Y += X;
  1516. I = 1;
  1517. }
  1518. }
  1519. }
  1520. }
  1521. /*=============================================*/
  1522. Milestone = 150;
  1523. /*=============================================*/
  1524. printf("Testing powers Z^Q at four nearly extreme values.\n");
  1525. N = 0;
  1526. Z = A1;
  1527. Q = FLOOR(Half - LOG(C) / LOG(A1));
  1528. Break = False;
  1529. do {
  1530. X = CInvrse;
  1531. Y = POW(Z, Q);
  1532. IsYeqX();
  1533. Q = - Q;
  1534. X = C;
  1535. Y = POW(Z, Q);
  1536. IsYeqX();
  1537. if (Z < One) Break = True;
  1538. else Z = AInvrse;
  1539. } while ( ! (Break));
  1540. PrintIfNPositive();
  1541. if (N == 0) printf(" ... no discrepancies found.\n");
  1542. printf("\n");
  1543. /*=============================================*/
  1544. Milestone = 160;
  1545. /*=============================================*/
  1546. Pause();
  1547. printf("Searching for Overflow threshold:\n");
  1548. printf("This may generate an error.\n");
  1549. Y = - CInvrse;
  1550. V9 = HInvrse * Y;
  1551. sigsave = sigfpe;
  1552. if (setjmp(ovfl_buf)) { I = 0; V9 = Y; goto overflow; }
  1553. do {
  1554. V = Y;
  1555. Y = V9;
  1556. V9 = HInvrse * Y;
  1557. } while(V9 < Y);
  1558. I = 1;
  1559. overflow:
  1560. sigsave = 0;
  1561. Z = V9;
  1562. printf("Can `Z = -Y' overflow?\n");
  1563. printf("Trying it on Y = %.17e .\n", Y);
  1564. V9 = - Y;
  1565. V0 = V9;
  1566. if (V - Y == V + V0) printf("Seems O.K.\n");
  1567. else {
  1568. printf("finds a ");
  1569. BadCond(Flaw, "-(-Y) differs from Y.\n");
  1570. }
  1571. if (Z != Y) {
  1572. BadCond(Serious, "");
  1573. printf("overflow past %.17e\n\tshrinks to %.17e .\n", Y, Z);
  1574. }
  1575. if (I) {
  1576. Y = V * (HInvrse * U2 - HInvrse);
  1577. Z = Y + ((One - HInvrse) * U2) * V;
  1578. if (Z < V0) Y = Z;
  1579. if (Y < V0) V = Y;
  1580. if (V0 - V < V0) V = V0;
  1581. }
  1582. else {
  1583. V = Y * (HInvrse * U2 - HInvrse);
  1584. V = V + ((One - HInvrse) * U2) * Y;
  1585. }
  1586. printf("Overflow threshold is V = %.17e .\n", V);
  1587. if (I) printf("Overflow saturates at V0 = %.17e .\n", V0);
  1588. else printf("There is no saturation value because the system traps on overflow.\n");
  1589. V9 = V * One;
  1590. printf("No Overflow should be signaled for V * 1 = %.17e\n", V9);
  1591. V9 = V / One;
  1592. printf(" nor for V / 1 = %.17e .\n", V9);
  1593. printf("Any overflow signal separating this * from the one\n");
  1594. printf("above is a DEFECT.\n");
  1595. /*=============================================*/
  1596. Milestone = 170;
  1597. /*=============================================*/
  1598. if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V)) {
  1599. BadCond(Failure, "Comparisons involving ");
  1600. printf("+-%g, +-%g\nand +-%g are confused by Overflow.",
  1601. V, V0, UfThold);
  1602. }
  1603. /*=============================================*/
  1604. Milestone = 175;
  1605. /*=============================================*/
  1606. printf("\n");
  1607. for(Indx = 1; Indx <= 3; ++Indx) {
  1608. switch (Indx) {
  1609. case 1: Z = UfThold; break;
  1610. case 2: Z = E0; break;
  1611. case 3: Z = PseudoZero; break;
  1612. }
  1613. if (Z != Zero) {
  1614. V9 = SQRT(Z);
  1615. Y = V9 * V9;
  1616. if (Y / (One - Radix * E9) < Z
  1617. || Y > (One + Radix * E9) * Z) { /* dgh: + E9 --> * E9 */
  1618. if (V9 > U1) BadCond(Serious, "");
  1619. else BadCond(Defect, "");
  1620. printf("Comparison alleges that what prints as Z = %.17e\n", Z);
  1621. printf(" is too far from sqrt(Z) ^ 2 = %.17e .\n", Y);
  1622. }
  1623. }
  1624. }
  1625. /*=============================================*/
  1626. Milestone = 180;
  1627. /*=============================================*/
  1628. for(Indx = 1; Indx <= 2; ++Indx) {
  1629. if (Indx == 1) Z = V;
  1630. else Z = V0;
  1631. V9 = SQRT(Z);
  1632. X = (One - Radix * E9) * V9;
  1633. V9 = V9 * X;
  1634. if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) {
  1635. Y = V9;
  1636. if (X < W) BadCond(Serious, "");
  1637. else BadCond(Defect, "");
  1638. printf("Comparison alleges that Z = %17e\n", Z);
  1639. printf(" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y);
  1640. }
  1641. }
  1642. /*=============================================*/
  1643. /*SPLIT
  1644. }
  1645. #include "paranoia.h"
  1646. part8(){
  1647. */
  1648. Milestone = 190;
  1649. /*=============================================*/
  1650. Pause();
  1651. X = UfThold * V;
  1652. Y = Radix * Radix;
  1653. if (X*Y < One || X > Y) {
  1654. if (X * Y < U1 || X > Y/U1) BadCond(Defect, "Badly");
  1655. else BadCond(Flaw, "");
  1656. printf(" unbalanced range; UfThold * V = %.17e\n\t%s\n",
  1657. X, "is too far from 1.\n");
  1658. }
  1659. /*=============================================*/
  1660. Milestone = 200;
  1661. /*=============================================*/
  1662. for (Indx = 1; Indx <= 5; ++Indx) {
  1663. X = F9;
  1664. switch (Indx) {
  1665. case 2: X = One + U2; break;
  1666. case 3: X = V; break;
  1667. case 4: X = UfThold; break;
  1668. case 5: X = Radix;
  1669. }
  1670. Y = X;
  1671. sigsave = sigfpe;
  1672. if (setjmp(ovfl_buf))
  1673. printf(" X / X traps when X = %g\n", X);
  1674. else {
  1675. V9 = (Y / X - Half) - Half;
  1676. if (V9 == Zero) continue;
  1677. if (V9 == - U1 && Indx < 5) BadCond(Flaw, "");
  1678. else BadCond(Serious, "");
  1679. printf(" X / X differs from 1 when X = %.17e\n", X);
  1680. printf(" instead, X / X - 1/2 - 1/2 = %.17e .\n", V9);
  1681. }
  1682. sigsave = 0;
  1683. }
  1684. /*=============================================*/
  1685. Milestone = 210;
  1686. /*=============================================*/
  1687. MyZero = Zero;
  1688. printf("\n");
  1689. printf("What message and/or values does Division by Zero produce?\n") ;
  1690. #ifndef NOPAUSE
  1691. printf("This can interupt your program. You can ");
  1692. printf("skip this part if you wish.\n");
  1693. printf("Do you wish to compute 1 / 0? ");
  1694. fflush(stdout);
  1695. read (KEYBOARD, ch, 8);
  1696. if ((ch[0] == 'Y') || (ch[0] == 'y')) {
  1697. #endif
  1698. sigsave = sigfpe;
  1699. printf(" Trying to compute 1 / 0 produces ...");
  1700. if (!setjmp(ovfl_buf)) printf(" %.7e .\n", One / MyZero);
  1701. sigsave = 0;
  1702. #ifndef NOPAUSE
  1703. }
  1704. else printf("O.K.\n");
  1705. printf("\nDo you wish to compute 0 / 0? ");
  1706. fflush(stdout);
  1707. read (KEYBOARD, ch, 80);
  1708. if ((ch[0] == 'Y') || (ch[0] == 'y')) {
  1709. #endif
  1710. sigsave = sigfpe;
  1711. printf("\n Trying to compute 0 / 0 produces ...");
  1712. if (!setjmp(ovfl_buf)) printf(" %.7e .\n", Zero / MyZero);
  1713. sigsave = 0;
  1714. #ifndef NOPAUSE
  1715. }
  1716. else printf("O.K.\n");
  1717. #endif
  1718. /*=============================================*/
  1719. Milestone = 220;
  1720. /*=============================================*/
  1721. Pause();
  1722. printf("\n");
  1723. {
  1724. static char *msg[] = {
  1725. "FAILUREs encountered =",
  1726. "SERIOUS DEFECTs discovered =",
  1727. "DEFECTs discovered =",
  1728. "FLAWs discovered =" };
  1729. int i;
  1730. for(i = 0; i < 4; i++) if (ErrCnt[i])
  1731. printf("The number of %-29s %d.\n",
  1732. msg[i], ErrCnt[i]);
  1733. }
  1734. printf("\n");
  1735. if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]
  1736. + ErrCnt[Flaw]) > 0) {
  1737. if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[
  1738. Defect] == 0) && (ErrCnt[Flaw] > 0)) {
  1739. printf("The arithmetic diagnosed seems ");
  1740. printf("Satisfactory though flawed.\n");
  1741. }
  1742. if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)
  1743. && ( ErrCnt[Defect] > 0)) {
  1744. printf("The arithmetic diagnosed may be Acceptable\n");
  1745. printf("despite inconvenient Defects.\n");
  1746. }
  1747. if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {
  1748. printf("The arithmetic diagnosed has ");
  1749. printf("unacceptable Serious Defects.\n");
  1750. }
  1751. if (ErrCnt[Failure] > 0) {
  1752. printf("Potentially fatal FAILURE may have spoiled this");
  1753. printf(" program's subsequent diagnoses.\n");
  1754. }
  1755. }
  1756. else {
  1757. printf("No failures, defects nor flaws have been discovered.\n");
  1758. if (! ((RMult == Rounded) && (RDiv == Rounded)
  1759. && (RAddSub == Rounded) && (RSqrt == Rounded)))
  1760. printf("The arithmetic diagnosed seems Satisfactory.\n");
  1761. else {
  1762. if (StickyBit >= One &&
  1763. (Radix - Two) * (Radix - Nine - One) == Zero) {
  1764. printf("Rounding appears to conform to ");
  1765. printf("the proposed IEEE standard P");
  1766. if ((Radix == Two) &&
  1767. ((Precision - Four * Three * Two) *
  1768. ( Precision - TwentySeven -
  1769. TwentySeven + One) == Zero))
  1770. printf("754");
  1771. else printf("854");
  1772. if (IEEE) printf(".\n");
  1773. else {
  1774. printf(",\nexcept for possibly Double Rounding");
  1775. printf(" during Gradual Underflow.\n");
  1776. }
  1777. }
  1778. printf("The arithmetic diagnosed appears to be Excellent!\n");
  1779. }
  1780. }
  1781. if (fpecount)
  1782. printf("\nA total of %d floating point exceptions were registered.\n",
  1783. fpecount);
  1784. printf("END OF TEST.\n");
  1785. return 0;
  1786. }
  1787. /*SPLIT subs.c
  1788. #include "paranoia.h"
  1789. */
  1790. /* Sign */
  1791. FLOAT Sign (X)
  1792. FLOAT X;
  1793. { return X >= 0. ? 1.0 : -1.0; }
  1794. /* Pause */
  1795. Pause()
  1796. {
  1797. #ifndef NOPAUSE
  1798. char ch[8];
  1799. printf("\nTo continue, press RETURN");
  1800. fflush(stdout);
  1801. read(KEYBOARD, ch, 8);
  1802. #endif
  1803. printf("\nDiagnosis resumes after milestone Number %d", Milestone);
  1804. printf(" Page: %d\n\n", PageNo);
  1805. ++Milestone;
  1806. ++PageNo;
  1807. }
  1808. /* TstCond */
  1809. TstCond (K, Valid, T)
  1810. int K, Valid;
  1811. char *T;
  1812. { if (! Valid) { BadCond(K,T); printf(".\n"); } }
  1813. BadCond(K, T)
  1814. int K;
  1815. char *T;
  1816. {
  1817. static char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
  1818. ErrCnt [K] = ErrCnt [K] + 1;
  1819. printf("%s: %s", msg[K], T);
  1820. }
  1821. /* Random */
  1822. /* Random computes
  1823. X = (Random1 + Random9)^5
  1824. Random1 = X - FLOOR(X) + 0.000005 * X;
  1825. and returns the new value of Random1
  1826. */
  1827. FLOAT Random()
  1828. {
  1829. FLOAT X, Y;
  1830. X = Random1 + Random9;
  1831. Y = X * X;
  1832. Y = Y * Y;
  1833. X = X * Y;
  1834. Y = X - FLOOR(X);
  1835. Random1 = Y + X * 0.000005;
  1836. return(Random1);
  1837. }
  1838. /* SqXMinX */
  1839. SqXMinX (ErrKind)
  1840. int ErrKind;
  1841. {
  1842. FLOAT XA, XB;
  1843. XB = X * BInvrse;
  1844. XA = X - XB;
  1845. SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp;
  1846. if (SqEr != Zero) {
  1847. if (SqEr < MinSqEr) MinSqEr = SqEr;
  1848. if (SqEr > MaxSqEr) MaxSqEr = SqEr;
  1849. J = J + 1.0;
  1850. BadCond(ErrKind, "\n");
  1851. printf("sqrt( %.17e) - %.17e = %.17e\n", X * X, X, OneUlp * SqEr);
  1852. printf("\tinstead of correct value 0 .\n");
  1853. }
  1854. }
  1855. /* NewD */
  1856. NewD()
  1857. {
  1858. X = Z1 * Q;
  1859. X = FLOOR(Half - X / Radix) * Radix + X;
  1860. Q = (Q - X * Z) / Radix + X * X * (D / Radix);
  1861. Z = Z - Two * X * D;
  1862. if (Z <= Zero) {
  1863. Z = - Z;
  1864. Z1 = - Z1;
  1865. }
  1866. D = Radix * D;
  1867. }
  1868. /* SR3750 */
  1869. SR3750()
  1870. {
  1871. if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {
  1872. I = I + 1;
  1873. X2 = SQRT(X * D);
  1874. Y2 = (X2 - Z2) - (Y - Z2);
  1875. X2 = X8 / (Y - Half);
  1876. X2 = X2 - Half * X2 * X2;
  1877. SqEr = (Y2 + Half) + (Half - X2);
  1878. if (SqEr < MinSqEr) MinSqEr = SqEr;
  1879. SqEr = Y2 - X2;
  1880. if (SqEr > MaxSqEr) MaxSqEr = SqEr;
  1881. }
  1882. }
  1883. /* IsYeqX */
  1884. IsYeqX()
  1885. {
  1886. if (Y != X) {
  1887. if (N <= 0) {
  1888. if (Z == Zero && Q <= Zero)
  1889. printf("WARNING: computing\n");
  1890. else BadCond(Defect, "computing\n");
  1891. printf("\t(%.17e) ^ (%.17e)\n", Z, Q);
  1892. printf("\tyielded %.17e;\n", Y);
  1893. printf("\twhich compared unequal to correct %.17e ;\n",
  1894. X);
  1895. printf("\t\tthey differ by %.17e .\n", Y - X);
  1896. }
  1897. N = N + 1; /* ... count discrepancies. */
  1898. }
  1899. }
  1900. /* SR3980 */
  1901. SR3980()
  1902. {
  1903. do {
  1904. Q = (FLOAT) I;
  1905. Y = POW(Z, Q);
  1906. IsYeqX();
  1907. if (++I > M) break;
  1908. X = Z * X;
  1909. } while ( X < W );
  1910. }
  1911. /* PrintIfNPositive */
  1912. PrintIfNPositive()
  1913. {
  1914. if (N > 0) printf("Similar discrepancies have occurred %d times.\n", N);
  1915. }
  1916. /* TstPtUf */
  1917. TstPtUf()
  1918. {
  1919. N = 0;
  1920. if (Z != Zero) {
  1921. printf("Since comparison denies Z = 0, evaluating ");
  1922. printf("(Z + Z) / Z should be safe.\n");
  1923. sigsave = sigfpe;
  1924. if (setjmp(ovfl_buf)) goto very_serious;
  1925. Q9 = (Z + Z) / Z;
  1926. printf("What the machine gets for (Z + Z) / Z is %.17e .\n",
  1927. Q9);
  1928. if (FABS(Q9 - Two) < Radix * U2) {
  1929. printf("This is O.K., provided Over/Underflow");
  1930. printf(" has NOT just been signaled.\n");
  1931. }
  1932. else {
  1933. if ((Q9 < One) || (Q9 > Two)) {
  1934. very_serious:
  1935. N = 1;
  1936. ErrCnt [Serious] = ErrCnt [Serious] + 1;
  1937. printf("This is a VERY SERIOUS DEFECT!\n");
  1938. }
  1939. else {
  1940. N = 1;
  1941. ErrCnt [Defect] = ErrCnt [Defect] + 1;
  1942. printf("This is a DEFECT!\n");
  1943. }
  1944. }
  1945. sigsave = 0;
  1946. V9 = Z * One;
  1947. Random1 = V9;
  1948. V9 = One * Z;
  1949. Random2 = V9;
  1950. V9 = Z / One;
  1951. if ((Z == Random1) && (Z == Random2) && (Z == V9)) {
  1952. if (N > 0) Pause();
  1953. }
  1954. else {
  1955. N = 1;
  1956. BadCond(Defect, "What prints as Z = ");
  1957. printf("%.17e\n\tcompares different from ", Z);
  1958. if (Z != Random1) printf("Z * 1 = %.17e ", Random1);
  1959. if (! ((Z == Random2)
  1960. || (Random2 == Random1)))
  1961. printf("1 * Z == %g\n", Random2);
  1962. if (! (Z == V9)) printf("Z / 1 = %.17e\n", V9);
  1963. if (Random2 != Random1) {
  1964. ErrCnt [Defect] = ErrCnt [Defect] + 1;
  1965. BadCond(Defect, "Multiplication does not commute!\n");
  1966. printf("\tComparison alleges that 1 * Z = %.17e\n",
  1967. Random2);
  1968. printf("\tdiffers from Z * 1 = %.17e\n", Random1);
  1969. }
  1970. Pause();
  1971. }
  1972. }
  1973. }
  1974. notify(s)
  1975. char *s;
  1976. {
  1977. printf("%s test appears to be inconsistent...\n", s);
  1978. printf(" PLEASE NOTIFY KARPINKSI!\n");
  1979. }
  1980. /*SPLIT msgs.c */
  1981. /* Instructions */
  1982. msglist(s)
  1983. char **s;
  1984. { while(*s) printf("%s\n", *s++); }
  1985. Instructions()
  1986. {
  1987. static char *instr[] = {
  1988. "Lest this program stop prematurely, i.e. before displaying\n",
  1989. " `END OF TEST',\n",
  1990. "try to persuade the computer NOT to terminate execution when an",
  1991. "error like Over/Underflow or Division by Zero occurs, but rather",
  1992. "to persevere with a surrogate value after, perhaps, displaying some",
  1993. "warning. If persuasion avails naught, don't despair but run this",
  1994. "program anyway to see how many milestones it passes, and then",
  1995. "amend it to make further progress.\n",
  1996. "Answer questions with Y, y, N or n (unless otherwise indicated).\n",
  1997. 0};
  1998. msglist(instr);
  1999. }
  2000. /* Heading */
  2001. Heading()
  2002. {
  2003. static char *head[] = {
  2004. "Users are invited to help debug and augment this program so it will",
  2005. "cope with unanticipated and newly uncovered arithmetic pathologies.\n",
  2006. "Please send suggestions and interesting results to",
  2007. "\tRichard Karpinski",
  2008. "\tComputer Center U-76",
  2009. "\tUniversity of California",
  2010. "\tSan Francisco, CA 94143-0704, USA\n",
  2011. "In doing so, please include the following information:",
  2012. #ifdef Single
  2013. "\tPrecision:\tsingle;",
  2014. #else
  2015. "\tPrecision:\tdouble;",
  2016. #endif
  2017. "\tVersion:\t10 February 1989;",
  2018. "\tComputer:\n",
  2019. "\tCompiler:\n",
  2020. "\tOptimization level:\n",
  2021. "\tOther relevant compiler options:",
  2022. 0};
  2023. msglist(head);
  2024. }
  2025. /* Characteristics */
  2026. Characteristics()
  2027. {
  2028. static char *chars[] = {
  2029. "Running this program should reveal these characteristics:",
  2030. " Radix = 1, 2, 4, 8, 10, 16, 100, 256 ...",
  2031. " Precision = number of significant digits carried.",
  2032. " U2 = Radix/Radix^Precision = One Ulp",
  2033. "\t(OneUlpnit in the Last Place) of 1.000xxx .",
  2034. " U1 = 1/Radix^Precision = One Ulp of numbers a little less than 1.0 .",
  2035. " Adequacy of guard digits for Mult., Div. and Subt.",
  2036. " Whether arithmetic is chopped, correctly rounded, or something else",
  2037. "\tfor Mult., Div., Add/Subt. and Sqrt.",
  2038. " Whether a Sticky Bit used correctly for rounding.",
  2039. " UnderflowThreshold = an underflow threshold.",
  2040. " E0 and PseudoZero tell whether underflow is abrupt, gradual, or fuzzy.",
  2041. " V = an overflow threshold, roughly.",
  2042. " V0 tells, roughly, whether Infinity is represented.",
  2043. " Comparisions are checked for consistency with subtraction",
  2044. "\tand for contamination with pseudo-zeros.",
  2045. " Sqrt is tested. Y^X is not tested.",
  2046. " Extra-precise subexpressions are revealed but NOT YET tested.",
  2047. " Decimal-Binary conversion is NOT YET tested for accuracy.",
  2048. 0};
  2049. msglist(chars);
  2050. }
  2051. History()
  2052. { /* History */
  2053. /* Converted from Brian Wichmann's Pascal version to C by Thos Sumner,
  2054. with further massaging by David M. Gay. */
  2055. static char *hist[] = {
  2056. "The program attempts to discriminate among",
  2057. " FLAWs, like lack of a sticky bit,",
  2058. " Serious DEFECTs, like lack of a guard digit, and",
  2059. " FAILUREs, like 2+2 == 5 .",
  2060. "Failures may confound subsequent diagnoses.\n",
  2061. "The diagnostic capabilities of this program go beyond an earlier",
  2062. "program called `MACHAR', which can be found at the end of the",
  2063. "book `Software Manual for the Elementary Functions' (1980) by",
  2064. "W. J. Cody and W. Waite. Although both programs try to discover",
  2065. "the Radix, Precision and range (over/underflow thresholds)",
  2066. "of the arithmetic, this program tries to cope with a wider variety",
  2067. "of pathologies, and to say how well the arithmetic is implemented.",
  2068. "\nThe program is based upon a conventional radix representation for",
  2069. "floating-point numbers, but also allows logarithmic encoding",
  2070. "as used by certain early WANG machines.\n",
  2071. "BASIC version of this program (C) 1983 by Prof. W. M. Kahan;",
  2072. "see source comments for more history.",
  2073. 0};
  2074. msglist(hist);
  2075. }
  2076. double
  2077. pow(x, y) /* return x ^ y (exponentiation) */
  2078. double x, y;
  2079. {
  2080. extern double exp(), frexp(), ldexp(), log(), modf();
  2081. double xy, ye;
  2082. long i;
  2083. int ex, ey = 0, flip = 0;
  2084. if (!y) return 1.0;
  2085. if ((y < -1100. || y > 1100.) && x != -1.) return exp(y * log(x));
  2086. if (y < 0.) { y = -y; flip = 1; }
  2087. y = modf(y, &ye);
  2088. if (y) xy = exp(y * log(x));
  2089. else xy = 1.0;
  2090. /* next several lines assume >= 32 bit integers */
  2091. x = frexp(x, &ex);
  2092. if (i = ye) for(;;) {
  2093. if (i & 1) { xy *= x; ey += ex; }
  2094. if (!(i >>= 1)) break;
  2095. x *= x;
  2096. ex *= 2;
  2097. if (x < .5) { x *= 2.; ex -= 1; }
  2098. }
  2099. if (flip) { xy = 1. / xy; ey = -ey; }
  2100. return ldexp(xy, ey);
  2101. }