| 1935 | /// \exception FE_UNDERFLOW on underflows |
| 1936 | /// \exception FE_INEXACT if \a arg is not a positive integer |
| 1937 | template<std::float_round_style R,bool L> unsigned int gamma(unsigned int arg) |
| 1938 | { |
| 1939 | /* static const double p[] ={ 2.50662827563479526904, 225.525584619175212544, -268.295973841304927459, 80.9030806934622512966, -5.00757863970517583837, 0.0114684895434781459556 }; |
| 1940 | double t = arg + 4.65, s = p[0]; |
| 1941 | for(unsigned int i=0; i<5; ++i) |
| 1942 | s += p[i+1] / (arg+i); |
| 1943 | return std::log(s) + (arg-0.5)*std::log(t) - t; |
| 1944 | */ static const f31 pi(0xC90FDAA2, 1), lbe(0xB8AA3B29, 0); |
| 1945 | unsigned int abs = arg & 0x7FFF, sign = arg & 0x8000; |
| 1946 | bool bsign = sign != 0; |
| 1947 | f31 z(abs), x = sign ? (z+f31(0x80000000, 0)) : z, t = x + f31(0x94CCCCCD, 2), s = |
| 1948 | f31(0xA06C9901, 1) + f31(0xBBE654E2, -7)/(x+f31(0x80000000, 2)) + f31(0xA1CE6098, 6)/(x+f31(0x80000000, 1)) |
| 1949 | + f31(0xE1868CB7, 7)/x - f31(0x8625E279, 8)/(x+f31(0x80000000, 0)) - f31(0xA03E158F, 2)/(x+f31(0xC0000000, 1)); |
| 1950 | int i = (s.exp>=2) + (s.exp>=4) + (s.exp>=8) + (s.exp>=16); |
| 1951 | s = f31((static_cast<uint32>(s.exp)<<(31-i))+(log2(s.m>>1, 28)>>i), i) / lbe; |
| 1952 | if(x.exp != -1 || x.m != 0x80000000) |
| 1953 | { |
| 1954 | i = (t.exp>=2) + (t.exp>=4) + (t.exp>=8); |
| 1955 | f31 l = f31((static_cast<uint32>(t.exp)<<(31-i))+(log2(t.m>>1, 30)>>i), i) / lbe; |
| 1956 | s = (x.exp<-1) ? (s-(f31(0x80000000, -1)-x)*l) : (s+(x-f31(0x80000000, -1))*l); |
| 1957 | } |
| 1958 | s = x.exp ? (s-t) : (t-s); |
| 1959 | if(bsign) |
| 1960 | { |
| 1961 | if(z.exp >= 0) |
| 1962 | { |
| 1963 | sign &= (L|((z.m>>(31-z.exp))&1)) - 1; |
| 1964 | for(z=f31((z.m<<(1+z.exp))&0xFFFFFFFF, -1); z.m<0x80000000; z.m<<=1,--z.exp) ; |
| 1965 | } |
| 1966 | if(z.exp == -1) |
| 1967 | z = f31(0x80000000, 0) - z; |
| 1968 | if(z.exp < -1) |
| 1969 | { |
| 1970 | z = z * pi; |
| 1971 | z.m = sincos(z.m>>(1-z.exp), 30).first; |
| 1972 | for(z.exp=1; z.m<0x80000000; z.m<<=1,--z.exp) ; |
| 1973 | } |
| 1974 | else |
| 1975 | z = f31(0x80000000, 0); |
| 1976 | } |
| 1977 | if(L) |
| 1978 | { |
| 1979 | if(bsign) |
| 1980 | { |
| 1981 | f31 l(0x92868247, 0); |
| 1982 | if(z.exp < 0) |
| 1983 | { |
| 1984 | uint32 m = log2((z.m+1)>>1, 27); |
| 1985 | z = f31(-((static_cast<uint32>(z.exp)<<26)+(m>>5)), 5); |
| 1986 | for(; z.m<0x80000000; z.m<<=1,--z.exp) ; |
| 1987 | l = l + z / lbe; |
| 1988 | } |
| 1989 | sign = static_cast<unsigned>(x.exp&&(l.exp<s.exp||(l.exp==s.exp&&l.m<s.m))) << 15; |
| 1990 | s = sign ? (s-l) : x.exp ? (l-s) : (l+s); |
| 1991 | } |
| 1992 | else |
| 1993 | { |
| 1994 | sign = static_cast<unsigned>(x.exp==0) << 15; |