rns@tortuga.SanDiego.NCR.COM (Rick Schubert) (01/16/91)
In article <1991Jan14.214648.14752@SanDiego.NCR.COM> randall@thor.SanDiego.NCR.COM (Randall Rathbun) writes: >Does the sqrt function in the c programming language math library have >a problem finding the square root of squares, even though their values >are well within the precision of the mantissa? (2^56-1 for double >precision, I believe, on most machines) >When the following short program is run, the output indicates an >error: > >sqrt function - (c should be integer and e = 0 ) >a,b,c = 14111955.000000000 96473184.000000000 97499859.000000015000000 >d,e = 97499859.000000000000000 0.000000014901161 > >hypot seems to be afflicted with the same problem also. >This seems to indicate the sqrt function has a problem. Has >anyone done a study of the sqrt function and its accuracy, or done >a thorough detailed error analysis of this function? >================ bad sqrt function? c program ================ >#include<stdio.h> >#include<math.h> >main () >{ >double a,b,c,d,e; > >a=14111955.0; >b=96473184.0; >c=sqrt(a*a+b*b); >d=floor(c); >e=c-d; >(void) printf("sqrt function - \n"); >(void) printf("a,b,c = %17.9f %17.9f %24.15f\n",a,b,c); >(void) printf("d,e = %24.15f %17.15f\n",d,e); >c=hypot(a,b); >d=floor(c); >e=c-d; >(void) printf("hypot function - \n"); >(void) printf("a,b,c = %17.9f %17.9f %24.15f\n",a,b,c); >(void) printf("d,e = %24.15f %17.15f\n",d,e); >(void) fflush(stdout); >return (0); >} Normally I wouldn't answer a question like this because I'm sure that there would be a number of people who answer (and many of them usually have answered before I've even read the article, although their answers may not have arrived here yet), but since this question comes from the same building as I'm in, I may have a chance to answer and have my article arrive at other sites soon after the original post does, thus preemting a certain number of other follow-ups. At least it's worth a try. The trouble isn't with the `sqrt' function; it is with the addition operator: Try running the following program: #include<stdio.h> main () { double a,b; a=14111955.0; b=96473184.0; printf("a = %f, a*a = %f\nb = %f, b*b = %f\na*a + b*b = %f\n", a,a*a,b,b*b,a*a+b*b ); (void) fflush(stdout); return (0); } This program produces the output: a = 14111955.000000, a*a = 199147273922025.000000 b = 96473184.000000, b*b = 9307075231097856.000000 a*a + b*b = 9506222505019882.000000 Now, 199147273922025 + 9307075231097856 != 9506222505019882 a*a b*b a*a + b*b So the problem must be that addition doesn't work on this machine. :-) in case you didn't know. The problem is with the limited precision of `double'. The number a*a + b*b cannot be represented exactly in double precision on most machines (it requires 54 bits of precision; I think your assumption of 56 bits in a double mantissa doesn't hold). Try this program: #include<stdio.h> main () { double a,b,c,d,e; double x,y,z; int i; a=14111955.0; b=96473184.0; x = a*a + b*b; for (i = 1; i < 10; ++i) printf("%f + %d = %f\n",x,i,x+i); (void) fflush(stdout); return (0); } It produces the output: 9506222505019882.000000 + 1 = 9506222505019884.000000 9506222505019882.000000 + 2 = 9506222505019884.000000 9506222505019882.000000 + 3 = 9506222505019886.000000 9506222505019882.000000 + 4 = 9506222505019886.000000 9506222505019882.000000 + 5 = 9506222505019888.000000 9506222505019882.000000 + 6 = 9506222505019888.000000 9506222505019882.000000 + 7 = 9506222505019890.000000 9506222505019882.000000 + 8 = 9506222505019890.000000 9506222505019882.000000 + 9 = 9506222505019892.000000 From this we can conclude that 1==2, 3==4, 5==6, 7==9. -- Rick Schubert (rns@tortuga.SanDiego.NCR.COM)