blog dds

2004.08.16

The hypot() Mystery

I was writing a section for the Code Reading followup volume, and wanted to demonstrate the pitfalls of using homebrewn mathematical functions instead of the library ones. As an example, I chose to compare the C library hypot(x, y) function, against sqrt(x * x, y * y). I created a plot of "unit in last place" (ulp) error values between the two functions, which demonstrated how the error increased for larger values of y.

Errors in the hypot implementation

However, I was feeling uneasy. The error distribution was not what I was expecting. First I thought the culrpit was the implementation of the ulp function. I was using David M. Gay's ulp implementation supplied with the FreeBSD gdtoa library, and was not sure how it behaved outside the library's domain. I combined it with a small test stub, and found that the values it returned were indeed correct.

The next step involved educating myself on the intricacies of mathematical function implementation on William Kahan's web page, and try some simple tests on hypot. The first test involved writing a test program to test hypot on arguments that formed Pythagorean triplets (Vasilis Capouleas advised me on the existence of formulas to derive them).

Amazingly, the function that returned a non-integer result was not the sqrt(x * x, y * y) function, but hypot(). The table at the end of this entry summarizes the wrong values for triplets derived from numbers up to 100. I was able to reproduce the errors with both Microsoft's Visual C compiler and the GNU C (Win32 mingw) compiler, both with, and without the "improve floating point consistency" switches (-Op and -ffloat-store). However, I was unable to reproduce the results on any non-Windows machine: Linux, and several versions of FreeBSD. Preplexed, I searched for the imports of the Windows executable files, and saw that they were importing hypot() from MSVCRT.DLL (version 7.0.2600.1106 - xpsp1.020828-1920). Mystery solved.

The C99 standard (Section 5.2.4.2.2) states that the accuracy of floating point functions is implementation-defined, so having hypot() return a 1 ulp error is allowed by the standard. Still, I find it surprising that with Sun's implementation, which guarantees a return value < 0.97 ulp, freely available since 1993 Microsoft is still distributing a less accurate implementation.

Indicative MSVCRT.DLL hypot() Errors

For a2 + b2 = c2 (integer a, b, c) the following values generate an error of 1 ulp.

abchypot(a, b)
9920101100.99999999999999
45108117116.99999999999999
16552173173.00000000000003
153104185185.00000000000003
22160229229.00000000000003
57176185184.99999999999997
171140221220.99999999999997
209120241240.99999999999997
21220221221.00000000000003
63216225225.00000000000003
35776365365.00000000000006
39940401401.00000000000006
69260269268.99999999999994
299180349348.99999999999994
391120409409.00000000000006
48344485485.00000000000006
425168457456.99999999999994
81360369368.99999999999994
189340389389.00000000000006
459220509508.99999999999994
567144585584.99999999999989
67552677677.00000000000011
261380461461.00000000000006
319360481480.99999999999994
78356785785.00000000000011
31480481481.00000000000006
93476485484.99999999999994
89960901901.00000000000011
231520569568.99999999999989
891180909908.99999999999989
665432793793.00000000000011
11556811571156.9999999999998
481600769768.99999999999989
925372997996.99999999999989
114720411651164.9999999999998
585648873873.00000000000011
741580941940.99999999999989
113134011811181.0000000000002
123836845844.99999999999989
205828853852.99999999999989
451780901901.00000000000011
697696985985.00000000000011
110747612051205.0000000000002
15998016011600.9999999999998
43924925924.99999999999989
129920929928.99999999999989
215912937937.00000000000011
301900949949.00000000000011
387884965964.99999999999989
90370411451145.0000000000002
141938014691469.0000000000002
150531215371536.9999999999998
45101210131013.0000000000001
130559214331432.9999999999998
139553214931492.9999999999998
148546815571557.0000000000002
157540016251624.9999999999998
166532816971696.9999999999998
19358819371937.0000000000002
141110011091108.9999999999998
329108011291129.0000000000002
164549217171717.0000000000002
173942017891788.9999999999998
192726419451945.0000000000002
441116012411241.0000000000002
122588815131513.0000000000002
191144019611961.0000000000002
23039623052304.9999999999995
137793616651665.0000000000002
188761619851985.0000000000002
249910025012501.0000000000005
689132014891488.9999999999998
1325109217171716.9999999999998
55151215131512.9999999999998
825140016251625.0000000000002
1265124817771776.9999999999998
1375120018251824.9999999999998
214575222732272.9999999999995
247550025252524.9999999999995
258540826172617.0000000000005
291510829172917.0000000000005
399160016491648.9999999999998
513158416651664.9999999999998
627156416851685.0000000000002
855151217371736.9999999999998
1425131219371936.9999999999998
1539126019891988.9999999999998
279342428252825.0000000000005
302122030293028.9999999999995
177173617451745.0000000000002
767165618251824.9999999999998
2183105624252425.0000000000005
289154029412941.0000000000005
300944030413040.9999999999995
312733631453144.9999999999995
324522832533252.9999999999995
336311633653365.0000000000005
549182019011901.0000000000002
671180019211921.0000000000002
1159168020412040.9999999999998
274584828732873.0000000000005
323345632653265.0000000000005
335534833733373.0000000000005
347723634853485.0000000000005
359912036013600.9999999999995
189198019891989.0000000000002
315197219971996.9999999999998
441196020092008.9999999999998
567194420252025.0000000000002
693192420452044.9999999999998
2709106029092908.9999999999995
308778431853184.9999999999995
321368432853285.0000000000005
359136036093609.0000000000005
384312438453845.0000000000005
195210821172116.9999999999995
331581234133413.0000000000005
344570835173516.9999999999995
396525239733972.9999999999995
469222022692268.9999999999995
871216023292328.9999999999995
1139210023892389.0000000000005
2211170027892789.0000000000005
2479156029292928.9999999999995
2613148430053004.9999999999995
3149114033493349.0000000000005
355184036493648.9999999999995
395350439853985.0000000000005
435513243574357.0000000000009
207237623852385.0000000000005
621234024212420.9999999999995
2553169630653065.0000000000005
2691162031413140.9999999999995
2967145633053304.9999999999995
3519108036813681.0000000000005
365797637853785.0000000000005
379586838933893.0000000000005
393375640054004.9999999999995
497249625452544.9999999999995
781246025812581.0000000000005
1349234027012701.0000000000005
3195150835333533.0000000000005
3337141636253624.9999999999995
3621122038213820.9999999999995
3763111639253924.9999999999995
461540846334633.0000000000009
475727647654765.0000000000009
73266426652665.0000000000005
657262427052704.9999999999995
803260427252724.9999999999995
2409212032093209.0000000000005
2555205232773277.0000000000005
2993182435053505.0000000000005
3285165236773676.9999999999995
3869126040694069.0000000000005
459968046494648.9999999999991
1275266829572957.0000000000005
1425263229932992.9999999999995
3525170839173917.0000000000005
3675161240134013.0000000000005
385295229772976.9999999999995
1617274431853185.0000000000005
1771270032293228.9999999999995
2387248434453445.0000000000005
2695235235773576.9999999999995
3157212438053805.0000000000005
3311204038893889.0000000000005
4081156043694369.0000000000009
4543122447054705.0000000000009
4697110448254825.0000000000009
562130056295629.0000000000009
869306031813181.0000000000005
1027303632053205.0000000000005
1501294033013300.9999999999995
2133275634853484.9999999999995
2923243638053804.9999999999995
4977113651055104.9999999999991
560960056415641.0000000000009
567325633053304.9999999999995
729324033213320.9999999999995
1377313634253424.9999999999995
1701306035013501.0000000000005
1863301635453544.9999999999995
2025296835933593.0000000000005
4455176847934793.0000000000009
5103129652655264.9999999999991
607546860936093.0000000000009
623731662456244.9999999999991
747340434853485.0000000000005
1909318037093709.0000000000005
4897170451855185.0000000000009
5063158453055305.0000000000009
589392459655965.0000000000009
639148064096409.0000000000009
655732465656564.9999999999991
672316467256725.0000000000009
935355236733672.9999999999995
1445346837573756.9999999999995
2295324839773976.9999999999995
3145292842974297.0000000000009
3825260046254625.0000000000009
5525150057255725.0000000000009
5695136858575856.9999999999991
620594862776277.0000000000009
671549267336733.0000000000009
87378437853784.9999999999995
609376038093808.9999999999995
957372438453845.0000000000005
1131370038693869.0000000000005
2001352040494049.0000000000005
4263258449854985.0000000000009
4437248450855085.0000000000009
4785227252975296.9999999999991
5307192456455645.0000000000009
5655167258975896.9999999999991
5829154060296028.9999999999991
6351112064496449.0000000000009
669982067496749.0000000000009
687366469056905.0000000000009
704750470657065.0000000000009
722134072297228.9999999999991
739517273977397.0000000000009
89396039613961.0000000000005
267395639653964.9999999999995
445394839733973.0000000000005
979390040214021.0000000000005
3649312048014800.9999999999991
4361276051615161.0000000000009
4717255653655365.0000000000009
5429210058215820.9999999999991
5785184860736073.0000000000009
6141158063416340.9999999999991
685399669256924.9999999999991
720968072417241.0000000000009
738751674057404.9999999999991
756534875737573.0000000000009
774317677457745.0000000000009
1001408042014201.0000000000009
5369240058815880.9999999999991
6097189663856384.9999999999991
6825132869536952.9999999999991
755369675857585.0000000000009
809918081018101.0000000000009
465431243374336.9999999999991
4185331253375336.9999999999991
4371322054295429.0000000000009
4929292057295728.9999999999991
6603180468456844.9999999999991
7161136072897288.9999999999991
7347120474457444.9999999999991
771988077697768.9999999999991
790571279377937.0000000000009
809154081098109.0000000000009
1045445245734572.9999999999991
4275350055255524.9999999999991
5035310859175916.9999999999991
5225300060256024.9999999999991
5605277262536253.0000000000009
5985252864976496.9999999999991
6365226867576756.9999999999991
7125170073257325.0000000000009
7315154874777477.0000000000009
7695123277937792.9999999999991
807590081258125.0000000000009
1649456048494849.0000000000009
2037448449254924.9999999999991
4947340460056004.9999999999991
5335319262176217.0000000000009
5917284465656565.0000000000009
6111272066896689.0000000000009
6499246069496949.0000000000009
6693232470857085.0000000000009
7857142479857985.0000000000009
99490049014901.0000000000009
297489649054905.0000000000009
2079468051215120.9999999999991
3069442053815380.9999999999991
4059406057415741.0000000000009
4455388859135912.9999999999991
5049360062016201.0000000000009
6633265671457145.0000000000009
7227223675657564.9999999999991
7623193678657864.9999999999991
8019162081818180.9999999999991
960319696059604.9999999999982

Read and post comments, or share through   


Creative Commons License Last modified: Monday, August 16, 2004 7:05 pm
Unless otherwise expressly stated, all original material on this page created by Diomidis Spinellis is licensed under a Creative Commons Attribution-Share Alike 3.0 Greece License.