🔒 Locked

Template Recursion

Started by Rattrap Aug 31, 2005 at 11:17 PM 14 replies 2.2k views

Original Post

Rattrap

125

Rattrap

Aug 31, 2005 at 11:17 PM

I've been playing with template recursion to compute constants at compile time. I have successfully implemented a factorial based one. I was trying to write one that computed e. e is defined as 1 + 1/1! + 1/2! + 1/3! + 1/n! So I tried


namespace KS
{
	namespace Math
	{
		template <unsigned int Limit, typename T>
		struct E
		{
			//static const T	result;
		};

		template <unsigned int Limit>
		struct E <Limit, float>
		{
			static const float	result;
		};
		template <>
		const float E <0, float>::result = 1.0f;
		//  Assume that KS::Math::Factoral does return a succesful factorial, which tests show that it does
		template <unsigned int Limit>
		const float E <Limit, float>::result = ((KS::Math::E<(Limit - 1), float>::result) + (1.0f / KS::Math::Factorial<(Limit)>::result));
	}
}


	//  Here is my test code
	std::cout << KS::Math::E<0, float>::result << std::endl;
	std::cout << KS::Math::E<1, float>::result << std::endl;
	std::cout << KS::Math::E<2, float>::result << std::endl;
	std::cout << KS::Math::E<6, float>::result << std::endl;

Here is the result of the code. 1 (This would be doing 1 which would be a straight constant, so corect) 2 (This should be 1/1! + 1, so again correct. Note that if I didn't actually put the std::cout << KS::Math::E<0, float>::result << std::endl; before it, it resolves to just 1, which would just be 1/1!) 2.5 (This should be 1/2! + 1/1! + 1, so again correct. Note that if I didn't actually put the std::cout << KS::Math::E<0, float>::result << std::endl; & std::cout << KS::Math::E<1, float>::result << std::endl; before it, it resolves to just 0.5, which would just be 1/2!) 0.00138889 (This is just resolving to 1/6!, instead of doing any of the recursion, as the above 3 show, I would guess that if I manually do E<3>, E<4>, & E<5> before it, they would all resolve correctly. So what it looks like it is doing is if a E is on its own, it just does the 1 / some number! math. If all of the neccessary calculations are predefined before it is needed, it is successfully solving for it. Is there just a problem with my logic somewhere? Note: You can't go higher than E<12, float>, because the factorial would overflow an int. [edit] This is all purely academic. So please, do not recommend using someone else's libraries ala boost or just using a constant. [/edit]

jflanglois

Sep 1, 2005 at 12:15 AM

I get the same behavior if I don't use the previous precisions.

This is my code (compiled with VC++.NET 2003):

// factorial.h -----------------------------------------------------#pragma oncenamespace Test {	template< unsigned i > struct factorial {		static const unsigned int result = i * Test::factorial< i - 1 >::result;	};	template<> struct factorial< 1 > {		static const unsigned int result = 1;	};	template<> struct factorial< 0 > {		static const unsigned int result = 1;	};}// e.h ------------------------------------------------------------#pragma once#include "factorial.h"namespace Test {	template< unsigned int precision > struct e {		static const float result;	};	template< unsigned int precision >	const float e< precision >::result =		Test::e< precision - 1 >::result + 			( 1.0f / Test::factorial< precision >::result );	template<> struct e< 0 > {		static const float result;	};	const float e< 0 >::result = 1.0f;}// main.cpp --------------------------------------------------------#include <iostream>#include <iomanip>#include "factorial.h"#include "e.h"using std::cout;using std::endl;int main() {	cout.precision( 20 );	cout << Test::factorial< 4 >::result << endl;	cout << Test::e< 1 >::result << endl		<< Test::e< 2 >::result << endl		<< Test::e< 3 >::result << endl		<< Test::e< 4 >::result << endl;	return 0;}

The result is:

2422.52.66666674613952642.7083334922790527

Which seems fine.

If I remove the first three Test::e templates, I get this:

240.041666667908430099

So, the same problem, it is only doing 1 / precision!.

I too would be interested in an academic solution.

[edit] slight changes to code
[edit2] Tested on:
G++ 3.2.3 - same problem
G++ 3.4.4 - works

jfl.

[Edited by - jflanglois on September 1, 2005 12:15:59 AM]

dragongame

Sep 1, 2005 at 12:46 AM

#include <iostream>namespace Test {        template< unsigned i > struct factorial {                static const float result = i*Test::factorial< i - 1 >::result;        };        template<> struct factorial< 1 > {                static const float result = 1;        };        template<> struct factorial< 0 > {                static const float result = 1;        };        template <unsigned int Limit>        struct E        {                static const float      result;        };        template <> const float E <0>::result = 1.0f;        template <unsigned int Limit> const float E <Limit>::result = ((E<(Limit - 1)>::result) + (1.0f / factorial<(Limit)>::result));}int main(){        std::cout << Test::E<1>::result << std::endl;        std::cout << Test::E<3>::result << std::endl;        std::cout << Test::E<6>::result << std::endl;        std::cout << Test::E<20>::result << std::endl;        std::cout << Test::E<25>::result << std::endl;        return 0;}

Gives me the following output:
2
2.66667
2.71806
2.71828
2.71828

This is with
powerpc-apple-darwin8-g++-4.0.0 (GCC) 4.0.0 20041026 (Apple Computer, Inc. build 4061) compiled.

Rattrap

Author

125

Rattrap

Sep 1, 2005 at 7:20 AM

I tested it with gcc in cygwin with 3.4.4. And it looks like it resolves correctly. So this looks like a problem with Visual C++ (I am currently compiling it in VS2005 beta 2).

Thanks for the help guys.

Rattrap

Author

125

Rattrap

Sep 2, 2005 at 12:07 PM

Just for an update, here is the official answer from Microsoft

Quote:

Resolved as Postponed by Microsoft on 2005-09-02 at 09:18:22

The order of initialization of static data members of class templates is not defined by the standard.

Please refer to section 3.6.2 of the standard. You can also examine the DR link below:

http://www.open-std.org/jtc1/sc22/wg21/docs/cwg_defects.html#270

The implementation does not require us to deal with this, and our compiler does not. We are not in the position that we can examine such a change for the current ship cycle, as we are only taking critical fixes. We will reexamine this after VC2005 ships.

For managed classes, the order is well defined and will produce the desired results.

Thank You,

Dean Wills
Visual C++ Team.

Sneftel

296

Sneftel

Sep 2, 2005 at 12:14 PM

template<bool TemplateTricksAreSilly>struct E{    static const float result;};template<>const float E<true>::result = 2.718281828f;

[grin]

Rattrap

Author

125

Rattrap

Sep 2, 2005 at 12:24 PM

I just have to be a smart ass.

Quote:
Original post by Rattrap
This is all purely academic. So please, do not recommend using someone else's libraries ala boost or just using a constant.

And besides wouldn't you want to do this

template<bool TemplateTricksAreSilly>struct E{    static const float result;};template<bool TemplateTricksAreSilly>const float E<TemplateTricksAreSilly>::result = 2.718281828f;ortemplate<bool TemplateTricksAreSilly = true>struct E{    static const float result;};template<>const float E<true>::result = 2.718281828f;template<>const float E<false>::result = 2.718281828f;

Sorry, in advance.

Sneftel

296

Sneftel

Sep 2, 2005 at 12:29 PM

I don't think it's semantically valid to set TemplateTricksAreSilly=false. [wink]

mattnewport

Sep 2, 2005 at 1:18 PM

You can do it using function calls. In a release build this will be optimized away to loading a floating point constant for each value of limit (I checked):

#include "stdafx.h"using namespace std;template<int n>class factorial{public:    static const int value = n * factorial<n - 1>::value;};template<>class factorial<0>{public:    static const int value = 1;};template<int limit>class e{public:    static double value()    {        return e<limit - 1>::value() + double(1) / double(factorial<limit>::value);    }};template<>class e<0>{public:    static double value()    {        return double(1);    }};int _tmain(int , _TCHAR* []){    cout << factorial<5>::value << endl;    cout << e<0>::value() << endl;    cout << e<1>::value() << endl;    cout << e<2>::value() << endl;    cout << e<3>::value() << endl;    cout << e<4>::value() << endl;    cout << e<5>::value() << endl;    cout << e<6>::value() << endl;    cout << e<7>::value() << endl;    cout << e<8>::value() << endl;    cout << e<9>::value() << endl;    return 0;}

Rattrap

Author

125

Rattrap

Sep 2, 2005 at 1:22 PM

Yeah, I have gotten a version that would work that did involve function calls. Also just for a test, it is actually bad to do the e in order, but it did work if you did 0,1,2,3,4,5,6,7 in a row, just not a value by itself (if a recursive call was needed).

And actually, it is even optimized away to a constant in debug. That is one of the cool things about these "Silly Template Tricks" :) Templates are resolved at compile-time, not at runtime.

mattnewport

Sep 2, 2005 at 1:28 PM

My version isn't optimized to a constant in debug - there are actual function calls to value() in the debug build. The version that doesn't use function calls does resolve to a constant even in debug but as has been pointed out it relies on undefined behaviour and so won't work in general, or indeed in this specific case under VC :-)

nmi

Sep 2, 2005 at 1:39 PM

Look at mathworld:
http://mathworld.wolfram.com/e.html

You can write e as a nested radical:
e = 1 + (1 + (1 + (1 + (...) / 4) / 3) / 2) / 1
(the mathworld version looks more beatiful)

template <unsigned N>struct E{	static double calc(double d)	{		return E<N-1>::calc(1.0 + d / (double)N);	}};template <>struct E<0>{	static double calc(double d)	{		return d;	}};double foo(){	double e = E<55>::calc(1.0);	return e;}/*#include <iostream>#include <iomanip>using namespace std;int main(){	cout << "d: " << setprecision(20) << foo() << endl;	return 0;}*/

The resulting code produced by gcc is:

C:\TEMP>gcc -c -O3 o3_test.ccC:\TEMP>objdump -DCzh o3_test.oo3_test.o:     file format pe-i386Sections:Idx Name          Size      VMA       LMA       File off  Algn  0 .text         00000010  00000000  00000000  000000b4  2**4                  CONTENTS, ALLOC, LOAD, RELOC, READONLY, CODE  1 .data         00000000  00000000  00000000  00000000  2**4                  ALLOC, LOAD, DATA  2 .bss          00000000  00000000  00000000  00000000  2**4                  ALLOC  3 .rdata        00000010  00000000  00000000  000000c4  2**4                  CONTENTS, ALLOC, LOAD, READONLY, DATADisassembly of section .text:00000000 <foo()>:   0:   55                      push   %ebp   1:   89 e5                   mov    %esp,%ebp   3:   5d                      pop    %ebp   4:   dd 05 00 00 00 00       fldl   0x0   a:   c3                      ret   b:   90                      nop   c:   90                      nop   d:   90                      nop   e:   90                      nop   f:   90                      nopDisassembly of section .rdata:00000000 <.rdata>:   0:   69 57 14 8b 0a bf 05    imul   $0x5bf0a8b,0x14(%edi),%edx   7:   40                      inc    %eax   8:   00 00                   add    %al,(%eax)   a:   00 00                   add    %al,(%eax)   c:   00 00                   add    %al,(%eax)   e:   00 00                   add    %al,(%eax)C:\TEMP>gcc -vReading specs from C:/MinGW/bin/../lib/gcc/mingw32/3.4.2/specsConfigured with: ../gcc/configure --with-gcc --with-gnu-ld --with-gnu-as --host=mingw32 --target=mingw32 --prefix=/mingw --enable-threads --disable-nls --enable-languages=c,c++,f77,ada,objc,java --disable-win32-registry --disable-shared --enable-sjlj-exceptions --enable-libgcj --disable-java-awt --without-x --enable-java-gc=boehm --disable-libgcj-debug --enable-interpreter --enable-hash-synchronization --enable-libstdcxx-debugThread model: win32gcc version 3.4.2 (mingw-special)

As you can see, no additional code is generated.
Also no factorials are calculated, since that is not necessary.

BTW a similar formula for Pi is:
pi/2 = 1 + (1/3)*(1 + (2/5)*(1 + (3/7)*(1 + (4/8)*(1 + ...

Rattrap

Author

125

Rattrap

Sep 2, 2005 at 1:39 PM

Hmmm. That's interesting. In the version I used, I am not showing any assmebly calls other than to cout to print it.

#define	_USE_MATH_DEFINES#include <iostream>#include <cmath>template <unsigned int Value>struct Factorial{	enum	{		result = Value * Factorial<(Value - 1)>::result	};};template <>struct Factorial <0>{	enum	{		result = 1	};};template <unsigned int Value, typename T>struct Invert{};template <unsigned int Value>struct Invert <Value, float>{	static const float result;};template <unsigned int Value>const float Invert <Value, float>::result = (1.0f / Value);template <unsigned int Value>struct Invert <Value, double>{	static const double result;};template <unsigned int Value>const double Invert <Value, double>::result = (1.0 / Value);template<unsigned int Limit>struct E{public:	static const double result;	static inline double f()	{		return Invert<Factorial<Limit>::result, double>::result + E<(Limit > 0) ? (Limit - 1) : 0>::f();	}};template<unsigned int Limit>const double E<Limit>::result = f();template <>struct E <0>{	//static const double result;	static inline double f()	{		return 1.0;	}};#define	QUICK_CHANGE	12int main(void){	std::cout << Factorial<QUICK_CHANGE>::result << std::endl;	std::cout.precision(36);	std::cout << Invert<Factorial<QUICK_CHANGE>::result, float>::result << std::endl;	std::cout << Invert<Factorial<QUICK_CHANGE>::result, double>::result << std::endl;	std::cout << "-----------------------------------" << std::endl;	double result = E<12>::result;	std::cout << result << std::endl;	std::cout << E<QUICK_CHANGE>::result << std::endl;	std::cout << M_E << std::endl;	return 0;}

[edit]
This is compiled in VS2003, and 2005 beta 2
[/edit]

Rattrap

Author

125

Rattrap

Sep 2, 2005 at 1:43 PM

I have seen 2 different methods of computing e

The one I used was

e = sum (1 / x!)
x = 0 -> infinity

and another

e = lim (1 + (1/N)^N
N->infinity

Since I had found the sum method first, that is what I went with. This mostly started out as a personal exercise in metalanguage templates.

mattnewport

Sep 2, 2005 at 1:46 PM

Your version has no extra function calls even in debug because you initialize the static data member with the function call whereas my version actually makes the function call when the value is needed in cout.

Polymorphic OOP

Sep 2, 2005 at 2:22 PM

If you want to not have problems with the initialization order then you'll have to use a compile-time representation of, for example, a floating point number -- a template which has a sign, mantissa, and exponent as template parameters. Perform all of your calculations at compile-time with metafunctions on instantatiations of that template to result in a type which represents the value entirely at compile-time. Finally, you initialize the runtime float, double, or long double directly from its compile-time representation. Making a compile-time floating point type is no simple task, especially if you wish to account for infinities, posititive and negative zero representations, proper promotion when working with other compile-time arithmetic types, and have the exact same results as the runtime counter-part.

Because of this, I recommend checking out boost's metamath library which is currently in development and can be downloaded from the vault for your tinkering pleasure. It contains a compile-time double implementation and is capable of all of these things. As well, they are currently working on compile-time trigonometric functions with it and other math functions. Unfortunately, the Boost vault mysteriously got hosed recently, losing just about everything, meaning you'll have to wait for people to re-upload. Once everything is back in order, you will be able to find the metamath library here, in expaler's directory.

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.