+++ /dev/null
-/** @file\r
- Compute the logrithm of x.\r
-\r
- Copyright (c) 2010 - 2011, Intel Corporation. All rights reserved.<BR>\r
- This program and the accompanying materials are licensed and made available under\r
- the terms and conditions of the BSD License that accompanies this distribution.\r
- The full text of the license may be found at\r
- http://opensource.org/licenses/bsd-license.\r
-\r
- THE PROGRAM IS DISTRIBUTED UNDER THE BSD LICENSE ON AN "AS IS" BASIS,\r
- WITHOUT WARRANTIES OR REPRESENTATIONS OF ANY KIND, EITHER EXPRESS OR IMPLIED.\r
-\r
- * ====================================================\r
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\r
- *\r
- * Developed at SunPro, a Sun Microsystems, Inc. business.\r
- * Permission to use, copy, modify, and distribute this\r
- * software is freely granted, provided that this notice\r
- * is preserved.\r
- * ====================================================\r
-\r
- e_sqrt.c 5.1 93/09/24\r
- NetBSD: e_sqrt.c,v 1.12 2002/05/26 22:01:52 wiz Exp\r
-**/\r
-#include <LibConfig.h>\r
-#include <sys/EfiCdefs.h>\r
-\r
-#include <errno.h>\r
-#include "math.h"\r
-#include "math_private.h"\r
-\r
-#if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */\r
-// potential divide by 0 -- near line 129, (x-x)/(x-x) is on purpose\r
-#pragma warning ( disable : 4723 )\r
-#endif\r
-\r
-/* __ieee754_sqrt(x)\r
- * Return correctly rounded sqrt.\r
- * ------------------------------------------\r
- * | Use the hardware sqrt if you have one |\r
- * ------------------------------------------\r
- * Method:\r
- * Bit by bit method using integer arithmetic. (Slow, but portable)\r
- * 1. Normalization\r
- * Scale x to y in [1,4) with even powers of 2:\r
- * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then\r
- * sqrt(x) = 2^k * sqrt(y)\r
- * 2. Bit by bit computation\r
- * Let q = sqrt(y) truncated to i bit after binary point (q = 1),\r
- * i 0\r
- * i+1 2\r
- * s = 2*q , and y = 2 * ( y - q ). (1)\r
- * i i i i\r
- *\r
- * To compute q from q , one checks whether\r
- * i+1 i\r
- *\r
- * -(i+1) 2\r
- * (q + 2 ) <= y. (2)\r
- * i\r
- * -(i+1)\r
- * If (2) is false, then q = q ; otherwise q = q + 2 .\r
- * i+1 i i+1 i\r
- *\r
- * With some algebric manipulation, it is not difficult to see\r
- * that (2) is equivalent to\r
- * -(i+1)\r
- * s + 2 <= y (3)\r
- * i i\r
- *\r
- * The advantage of (3) is that s and y can be computed by\r
- * i i\r
- * the following recurrence formula:\r
- * if (3) is false\r
- *\r
- * s = s , y = y ; (4)\r
- * i+1 i i+1 i\r
- *\r
- * otherwise,\r
- * -i -(i+1)\r
- * s = s + 2 , y = y - s - 2 (5)\r
- * i+1 i i+1 i i\r
- *\r
- * One may easily use induction to prove (4) and (5).\r
- * Note. Since the left hand side of (3) contain only i+2 bits,\r
- * it does not necessary to do a full (53-bit) comparison\r
- * in (3).\r
- * 3. Final rounding\r
- * After generating the 53 bits result, we compute one more bit.\r
- * Together with the remainder, we can decide whether the\r
- * result is exact, bigger than 1/2ulp, or less than 1/2ulp\r
- * (it will never equal to 1/2ulp).\r
- * The rounding mode can be detected by checking whether\r
- * huge + tiny is equal to huge, and whether huge - tiny is\r
- * equal to huge for some floating point number "huge" and "tiny".\r
- *\r
- * Special cases:\r
- * sqrt(+-0) = +-0 ... exact\r
- * sqrt(inf) = inf\r
- * sqrt(-ve) = NaN ... with invalid signal\r
- * sqrt(NaN) = NaN ... with invalid signal for signaling NaN\r
- *\r
- * Other methods : see the appended file at the end of the program below.\r
- *---------------\r
- */\r
-\r
-static const double one = 1.0, tiny=1.0e-300;\r
-\r
-double\r
-__ieee754_sqrt(double x)\r
-{\r
- double z;\r
- int32_t sign = (int)0x80000000;\r
- int32_t ix0,s0,q,m,t,i;\r
- u_int32_t r,t1,s1,ix1,q1;\r
-\r
- EXTRACT_WORDS(ix0,ix1,x);\r
-\r
- /* take care of Inf and NaN */\r
- if((ix0&0x7ff00000)==0x7ff00000) {\r
- return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf\r
- sqrt(-inf)=sNaN */\r
- }\r
- /* take care of zero */\r
- if(ix0<=0) {\r
- if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */\r
- else if(ix0<0) {\r
- errno = EDOM;\r
- return (x-x)/(x-x); /* sqrt(-ve) = sNaN */\r
- }\r
- }\r
- /* normalize x */\r
- m = (ix0>>20);\r
- if(m==0) { /* subnormal x */\r
- while(ix0==0) {\r
- m -= 21;\r
- ix0 |= (ix1>>11); ix1 <<= 21;\r
- }\r
- for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;\r
- m -= i-1;\r
- ix0 |= (ix1>>(32-i));\r
- ix1 <<= i;\r
- }\r
- m -= 1023; /* unbias exponent */\r
- ix0 = (ix0&0x000fffff)|0x00100000;\r
- if(m&1){ /* odd m, double x to make it even */\r
- ix0 += ix0 + ((ix1&sign)>>31);\r
- ix1 += ix1;\r
- }\r
- m >>= 1; /* m = [m/2] */\r
-\r
- /* generate sqrt(x) bit by bit */\r
- ix0 += ix0 + ((ix1&sign)>>31);\r
- ix1 += ix1;\r
- q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */\r
- r = 0x00200000; /* r = moving bit from right to left */\r
-\r
- while(r!=0) {\r
- t = s0+r;\r
- if(t<=ix0) {\r
- s0 = t+r;\r
- ix0 -= t;\r
- q += r;\r
- }\r
- ix0 += ix0 + ((ix1&sign)>>31);\r
- ix1 += ix1;\r
- r>>=1;\r
- }\r
-\r
- r = sign;\r
- while(r!=0) {\r
- t1 = s1+r;\r
- t = s0;\r
- if((t<ix0)||((t==ix0)&&(t1<=ix1))) {\r
- s1 = t1+r;\r
- if(((t1&sign)==(u_int32_t)sign)&&(s1&sign)==0) s0 += 1;\r
- ix0 -= t;\r
- if (ix1 < t1) ix0 -= 1;\r
- ix1 -= t1;\r
- q1 += r;\r
- }\r
- ix0 += ix0 + ((ix1&sign)>>31);\r
- ix1 += ix1;\r
- r>>=1;\r
- }\r
-\r
- /* use floating add to find out rounding direction */\r
- if((ix0|ix1)!=0) {\r
- z = one-tiny; /* trigger inexact flag */\r
- if (z>=one) {\r
- z = one+tiny;\r
- if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}\r
- else if (z>one) {\r
- if (q1==(u_int32_t)0xfffffffe) q+=1;\r
- q1+=2;\r
- } else\r
- q1 += (q1&1);\r
- }\r
- }\r
- ix0 = (q>>1)+0x3fe00000;\r
- ix1 = q1>>1;\r
- if ((q&1)==1) ix1 |= sign;\r
- ix0 += (m <<20);\r
- INSERT_WORDS(z,ix0,ix1);\r
- return z;\r
-}\r
-\r
-/*\r
-Other methods (use floating-point arithmetic)\r
--------------\r
-(This is a copy of a drafted paper by Prof W. Kahan\r
-and K.C. Ng, written in May, 1986)\r
-\r
- Two algorithms are given here to implement sqrt(x)\r
- (IEEE double precision arithmetic) in software.\r
- Both supply sqrt(x) correctly rounded. The first algorithm (in\r
- Section A) uses newton iterations and involves four divisions.\r
- The second one uses reciproot iterations to avoid division, but\r
- requires more multiplications. Both algorithms need the ability\r
- to chop results of arithmetic operations instead of round them,\r
- and the INEXACT flag to indicate when an arithmetic operation\r
- is executed exactly with no roundoff error, all part of the\r
- standard (IEEE 754-1985). The ability to perform shift, add,\r
- subtract and logical AND operations upon 32-bit words is needed\r
- too, though not part of the standard.\r
-\r
-A. sqrt(x) by Newton Iteration\r
-\r
- (1) Initial approximation\r
-\r
- Let x0 and x1 be the leading and the trailing 32-bit words of\r
- a floating point number x (in IEEE double format) respectively\r
-\r
- 1 11 52 ...widths\r
- ------------------------------------------------------\r
- x: |s| e | f |\r
- ------------------------------------------------------\r
- msb lsb msb lsb ...order\r
-\r
-\r
- ------------------------ ------------------------\r
- x0: |s| e | f1 | x1: | f2 |\r
- ------------------------ ------------------------\r
-\r
- By performing shifts and subtracts on x0 and x1 (both regarded\r
- as integers), we obtain an 8-bit approximation of sqrt(x) as\r
- follows.\r
-\r
- k := (x0>>1) + 0x1ff80000;\r
- y0 := k - T1[31&(k>>15)]. ... y ~ sqrt(x) to 8 bits\r
- Here k is a 32-bit integer and T1[] is an integer array containing\r
- correction terms. Now magically the floating value of y (y's\r
- leading 32-bit word is y0, the value of its trailing word is 0)\r
- approximates sqrt(x) to almost 8-bit.\r
-\r
- Value of T1:\r
- static int T1[32]= {\r
- 0, 1024, 3062, 5746, 9193, 13348, 18162, 23592,\r
- 29598, 36145, 43202, 50740, 58733, 67158, 75992, 85215,\r
- 83599, 71378, 60428, 50647, 41945, 34246, 27478, 21581,\r
- 16499, 12183, 8588, 5674, 3403, 1742, 661, 130,};\r
-\r
- (2) Iterative refinement\r
-\r
- Apply Heron's rule three times to y, we have y approximates\r
- sqrt(x) to within 1 ulp (Unit in the Last Place):\r
-\r
- y := (y+x/y)/2 ... almost 17 sig. bits\r
- y := (y+x/y)/2 ... almost 35 sig. bits\r
- y := y-(y-x/y)/2 ... within 1 ulp\r
-\r
-\r
- Remark 1.\r
- Another way to improve y to within 1 ulp is:\r
-\r
- y := (y+x/y) ... almost 17 sig. bits to 2*sqrt(x)\r
- y := y - 0x00100006 ... almost 18 sig. bits to sqrt(x)\r
-\r
- 2\r
- (x-y )*y\r
- y := y + 2* ---------- ...within 1 ulp\r
- 2\r
- 3y + x\r
-\r
-\r
- This formula has one division fewer than the one above; however,\r
- it requires more multiplications and additions. Also x must be\r
- scaled in advance to avoid spurious overflow in evaluating the\r
- expression 3y*y+x. Hence it is not recommended uless division\r
- is slow. If division is very slow, then one should use the\r
- reciproot algorithm given in section B.\r
-\r
- (3) Final adjustment\r
-\r
- By twiddling y's last bit it is possible to force y to be\r
- correctly rounded according to the prevailing rounding mode\r
- as follows. Let r and i be copies of the rounding mode and\r
- inexact flag before entering the square root program. Also we\r
- use the expression y+-ulp for the next representable floating\r
- numbers (up and down) of y. Note that y+-ulp = either fixed\r
- point y+-1, or multiply y by nextafter(1,+-inf) in chopped\r
- mode.\r
-\r
- I := FALSE; ... reset INEXACT flag I\r
- R := RZ; ... set rounding mode to round-toward-zero\r
- z := x/y; ... chopped quotient, possibly inexact\r
- If(not I) then { ... if the quotient is exact\r
- if(z=y) {\r
- I := i; ... restore inexact flag\r
- R := r; ... restore rounded mode\r
- return sqrt(x):=y.\r
- } else {\r
- z := z - ulp; ... special rounding\r
- }\r
- }\r
- i := TRUE; ... sqrt(x) is inexact\r
- If (r=RN) then z=z+ulp ... rounded-to-nearest\r
- If (r=RP) then { ... round-toward-+inf\r
- y = y+ulp; z=z+ulp;\r
- }\r
- y := y+z; ... chopped sum\r
- y0:=y0-0x00100000; ... y := y/2 is correctly rounded.\r
- I := i; ... restore inexact flag\r
- R := r; ... restore rounded mode\r
- return sqrt(x):=y.\r
-\r
- (4) Special cases\r
-\r
- Square root of +inf, +-0, or NaN is itself;\r
- Square root of a negative number is NaN with invalid signal.\r
-\r
-\r
-B. sqrt(x) by Reciproot Iteration\r
-\r
- (1) Initial approximation\r
-\r
- Let x0 and x1 be the leading and the trailing 32-bit words of\r
- a floating point number x (in IEEE double format) respectively\r
- (see section A). By performing shifs and subtracts on x0 and y0,\r
- we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.\r
-\r
- k := 0x5fe80000 - (x0>>1);\r
- y0:= k - T2[63&(k>>14)]. ... y ~ 1/sqrt(x) to 7.8 bits\r
-\r
- Here k is a 32-bit integer and T2[] is an integer array\r
- containing correction terms. Now magically the floating\r
- value of y (y's leading 32-bit word is y0, the value of\r
- its trailing word y1 is set to zero) approximates 1/sqrt(x)\r
- to almost 7.8-bit.\r
-\r
- Value of T2:\r
- static int T2[64]= {\r
- 0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,\r
- 0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,\r
- 0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,\r
- 0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,\r
- 0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,\r
- 0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,\r
- 0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,\r
- 0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,};\r
-\r
- (2) Iterative refinement\r
-\r
- Apply Reciproot iteration three times to y and multiply the\r
- result by x to get an approximation z that matches sqrt(x)\r
- to about 1 ulp. To be exact, we will have\r
- -1ulp < sqrt(x)-z<1.0625ulp.\r
-\r
- ... set rounding mode to Round-to-nearest\r
- y := y*(1.5-0.5*x*y*y) ... almost 15 sig. bits to 1/sqrt(x)\r
- y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)\r
- ... special arrangement for better accuracy\r
- z := x*y ... 29 bits to sqrt(x), with z*y<1\r
- z := z + 0.5*z*(1-z*y) ... about 1 ulp to sqrt(x)\r
-\r
- Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that\r
- (a) the term z*y in the final iteration is always less than 1;\r
- (b) the error in the final result is biased upward so that\r
- -1 ulp < sqrt(x) - z < 1.0625 ulp\r
- instead of |sqrt(x)-z|<1.03125ulp.\r
-\r
- (3) Final adjustment\r
-\r
- By twiddling y's last bit it is possible to force y to be\r
- correctly rounded according to the prevailing rounding mode\r
- as follows. Let r and i be copies of the rounding mode and\r
- inexact flag before entering the square root program. Also we\r
- use the expression y+-ulp for the next representable floating\r
- numbers (up and down) of y. Note that y+-ulp = either fixed\r
- point y+-1, or multiply y by nextafter(1,+-inf) in chopped\r
- mode.\r
-\r
- R := RZ; ... set rounding mode to round-toward-zero\r
- switch(r) {\r
- case RN: ... round-to-nearest\r
- if(x<= z*(z-ulp)...chopped) z = z - ulp; else\r
- if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;\r
- break;\r
- case RZ:case RM: ... round-to-zero or round-to--inf\r
- R:=RP; ... reset rounding mod to round-to-+inf\r
- if(x<z*z ... rounded up) z = z - ulp; else\r
- if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;\r
- break;\r
- case RP: ... round-to-+inf\r
- if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else\r
- if(x>z*z ...chopped) z = z+ulp;\r
- break;\r
- }\r
-\r
- Remark 3. The above comparisons can be done in fixed point. For\r
- example, to compare x and w=z*z chopped, it suffices to compare\r
- x1 and w1 (the trailing parts of x and w), regarding them as\r
- two's complement integers.\r
-\r
- ...Is z an exact square root?\r
- To determine whether z is an exact square root of x, let z1 be the\r
- trailing part of z, and also let x0 and x1 be the leading and\r
- trailing parts of x.\r
-\r
- If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0\r
- I := 1; ... Raise Inexact flag: z is not exact\r
- else {\r
- j := 1 - [(x0>>20)&1] ... j = logb(x) mod 2\r
- k := z1 >> 26; ... get z's 25-th and 26-th\r
- fraction bits\r
- I := i or (k&j) or ((k&(j+j+1))!=(x1&3));\r
- }\r
- R:= r ... restore rounded mode\r
- return sqrt(x):=z.\r
-\r
- If multiplication is cheaper than the foregoing red tape, the\r
- Inexact flag can be evaluated by\r
-\r
- I := i;\r
- I := (z*z!=x) or I.\r
-\r
- Note that z*z can overwrite I; this value must be sensed if it is\r
- True.\r
-\r
- Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be\r
- zero.\r
-\r
- --------------------\r
- z1: | f2 |\r
- --------------------\r
- bit 31 bit 0\r
-\r
- Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd\r
- or even of logb(x) have the following relations:\r
-\r
- -------------------------------------------------\r
- bit 27,26 of z1 bit 1,0 of x1 logb(x)\r
- -------------------------------------------------\r
- 00 00 odd and even\r
- 01 01 even\r
- 10 10 odd\r
- 10 00 even\r
- 11 01 even\r
- -------------------------------------------------\r
-\r
- (4) Special cases (see (4) of Section A).\r
-\r
- */\r
-\r
projects / mirror_edk2.git / blobdiff