template unsigned int gamma()

in include/half.hpp [1937:2032]

• 91 lines of code
• 31 McCabe index (conditional complexity)

		template<std::float_round_style R,bool L> unsigned int gamma(unsigned int arg)
		{
/*			static const double p[] ={ 2.50662827563479526904, 225.525584619175212544, -268.295973841304927459, 80.9030806934622512966, -5.00757863970517583837, 0.0114684895434781459556 };
			double t = arg + 4.65, s = p[0];
			for(unsigned int i=0; i<5; ++i)
				s += p[i+1] / (arg+i);
			return std::log(s) + (arg-0.5)*std::log(t) - t;
*/			static const f31 pi(0xC90FDAA2, 1), lbe(0xB8AA3B29, 0);
			unsigned int abs = arg & 0x7FFF, sign = arg & 0x8000;
			bool bsign = sign != 0;
			f31 z(abs), x = sign ? (z+f31(0x80000000, 0)) : z, t = x + f31(0x94CCCCCD, 2), s =
				f31(0xA06C9901, 1) + f31(0xBBE654E2, -7)/(x+f31(0x80000000, 2)) + f31(0xA1CE6098, 6)/(x+f31(0x80000000, 1))
				+ f31(0xE1868CB7, 7)/x - f31(0x8625E279, 8)/(x+f31(0x80000000, 0)) - f31(0xA03E158F, 2)/(x+f31(0xC0000000, 1));
			int i = (s.exp>=2) + (s.exp>=4) + (s.exp>=8) + (s.exp>=16);
			s = f31((static_cast<uint32>(s.exp)<<(31-i))+(log2(s.m>>1, 28)>>i), i) / lbe;
			if(x.exp != -1 || x.m != 0x80000000)
			{
				i = (t.exp>=2) + (t.exp>=4) + (t.exp>=8);
				f31 l = f31((static_cast<uint32>(t.exp)<<(31-i))+(log2(t.m>>1, 30)>>i), i) / lbe;
				s = (x.exp<-1) ? (s-(f31(0x80000000, -1)-x)*l) : (s+(x-f31(0x80000000, -1))*l);
			}
			s = x.exp ? (s-t) : (t-s);
			if(bsign)
			{
				if(z.exp >= 0)
				{
					sign &= (L|((z.m>>(31-z.exp))&1)) - 1;
					for(z=f31((z.m<<(1+z.exp))&0xFFFFFFFF, -1); z.m<0x80000000; z.m<<=1,--z.exp) ;
				}
				if(z.exp == -1)
					z = f31(0x80000000, 0) - z;
				if(z.exp < -1)
				{
					z = z * pi;
					z.m = sincos(z.m>>(1-z.exp), 30).first;
					for(z.exp=1; z.m<0x80000000; z.m<<=1,--z.exp) ;
				}
				else
					z = f31(0x80000000, 0);
			}
			if(L)
			{
				if(bsign)
				{
					f31 l(0x92868247, 0);
					if(z.exp < 0)
					{
						uint32 m = log2((z.m+1)>>1, 27);
						z = f31(-((static_cast<uint32>(z.exp)<<26)+(m>>5)), 5);
						for(; z.m<0x80000000; z.m<<=1,--z.exp) ;
						l = l + z / lbe;
					}
					sign = static_cast<unsigned>(x.exp&&(l.exp<s.exp||(l.exp==s.exp&&l.m<s.m))) << 15;
					s = sign ? (s-l) : x.exp ? (l-s) : (l+s);
				}
				else
				{
					sign = static_cast<unsigned>(x.exp==0) << 15;
					if(s.exp < -24)
						return underflow<R>(sign);
					if(s.exp > 15)
						return overflow<R>(sign);
				}
			}
			else
			{
				s = s * lbe;
				uint32 m;
				if(s.exp < 0)
				{
					m = s.m >> -s.exp;
					s.exp = 0;
				}
				else
				{
					m = (s.m<<s.exp) & 0x7FFFFFFF;
					s.exp = (s.m>>(31-s.exp));
				}
				s.m = exp2(m, 27);
				if(!x.exp)
					s = f31(0x80000000, 0) / s;
				if(bsign)
				{
					if(z.exp < 0)
						s = s * z;
					s = pi / s;
					if(s.exp < -24)
						return underflow<R>(sign);
				}
				else if(z.exp > 0 && !(z.m&((1<<(31-z.exp))-1)))
					return ((s.exp+14)<<10) + (s.m>>21);
				if(s.exp > 15)
					return overflow<R>(sign);
			}
			return fixed2half<R,31,false,false,true>(s.m, s.exp+14, sign);
		}