Sign in to follow this  

Template Recursion

This topic is 4490 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

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<some number, float> 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]

Share this post


Link to post
Share on other sites
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 once


namespace 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:
24
2
2.5
2.6666667461395264
2.7083334922790527

Which seems fine.

If I remove the first three Test::e templates, I get this:
24
0.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]

Share this post


Link to post
Share on other sites

#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.

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites

template<bool TemplateTricksAreSilly>
struct E
{
static const float result;
};

template<>
const float E<true>::result = 2.718281828f;


[grin]

Share this post


Link to post
Share on other sites
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;

or

template<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.

Share this post


Link to post
Share on other sites
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;
}

Share this post


Link to post
Share on other sites
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<x> 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.

Share this post


Link to post
Share on other sites
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 :-)

Share this post


Link to post
Share on other sites
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.cc

C:\TEMP>objdump -DCzh o3_test.o

o3_test.o: file format pe-i386

Sections:
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, DATA
Disassembly 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 nop
Disassembly 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 -v
Reading specs from C:/MinGW/bin/../lib/gcc/mingw32/3.4.2/specs
Configured 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 --e
nable-sjlj-exceptions --enable-libgcj --disable-java-awt --without-x --enable-ja
va-gc=boehm --disable-libgcj-debug --enable-interpreter --enable-hash-synchroniz
ation --enable-libstdcxx-debug
Thread model: win32
gcc 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 + ...

Share this post


Link to post
Share on other sites
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 12

int 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]

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites

This topic is 4490 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this