Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I want to find the product of the first 20 Fibonacci numbers (the fibonorial). I want implementations using While, Do, and For loops.

This is my code using Do:

Do[Sum[Print[a = a*fib[n]], {n, 1, 20}], {a, 1}]

But I can't get it working with While or For.

Edit

My fib code is as follows:

Clear[fib, x, a]
fib[1] = 1; fib[2] = 1;
fib[n_] := fib[n] = fib[n - 1] + fib[n - 2]
fib[10]
share|improve this question
    
Could you post the code that you had trouble with? –  jVincent Feb 15 '13 at 10:08
    
@ jvincent hey vincent, i actually solved it. i failed to set an initial condition. however, i am looking to improve my code, as it gives a result of every sum value, rather then just the last value which I am looking for. is there anyway to do that? \For[n = 1; a = 1, n < 21, n++, Print[a = a*fib[n]]] \a = 1; n = 1; While[n < 21, Print[a = a*fib[n]]; n++] –  julesverne Feb 15 '13 at 10:16
    
In the future please format the code using two ` or four indents. mathematica.stackexchange.com/faq –  Ajasja Feb 15 '13 at 10:21
    
You seem to be misunderstanding what Print does. It prints the results to a new line, but doesn't actually return it. If you remove Print and run either of your examples, the result will be stored in a at the end of the loop. –  jVincent Feb 15 '13 at 10:26
    
wow, thank you so much! –  julesverne Feb 15 '13 at 10:30

2 Answers 2

up vote 3 down vote accepted

According to your specific needs, which seems to be more syntax and not performance related:

While loop:

Clear[a, n]
n = 1; a = 1;
While[n <= 20, Print[a = a*fib[n]]; n++]

For loop:

Clear[a, n]
For[n = 1; a = 1, n <= 20, n++, Print[a = a*fib[n]]]

I guess you know that there are other ways for compute Fibonacci. Roman Maeder's book, for instance, Computer Science with Mathematia or Michael Trott's Guidebook to Programming in Mathematica.

Edit: Approaches to calculate product of Fibonacci Numbers:

Times @@ Array[Fibonacci, 20]

Times @@ Fibonacci[Range@20]

a = 0; Times @@ NestWhileList[a + # &, 1, (a = #1) < Fibonacci@19 &, 2]

fibProduct[n_Integer?Positive] :=
    Module[{f1 = 1, f2 = 0, res = 1},
    Do[{f1, f2} = {f1 + f2, f1}; res *= f1, {n - 1}];
    res]

fibProduct[n_Integer?NonPositive] := Times @@ ((-1)^(n - 1) Fibonacci[Range@-n])

Edit 2: the not for the faint of heart approach:

The .tm file:

:Begin:
:Function:      fibProduct
:Pattern:       FibProduct[n_]
:Arguments:     { n }
:ArgumentTypes: { Integer }
:ReturnType:    Integer
:End:

The .cpp file:

#include <algorithm>
#include <iostream>
#include <numeric>
#include <iterator>
#include <vector>


using namespace std;

template<class F, class X, class S>
X foldl( F && f, X x, const S & s)
{
    return std::accumulate(std::begin(s), std::end(s), std::move(x), std::forward<F>(f));
}

unsigned long long fibProduct( int n )
{
    std::function<int(int)> fibonacci = [&fibonacci](int n) -> int { return (n < 2) ? 
        n : fibonacci(n - 1) + fibonacci(n - 2); };
    vector<int> vec(n);
    int i = 0;
    generate( begin(vec), end(vec), [&i, &fibonacci]() { return fibonacci(i += 1); } );

    return foldl( std::multiplies<int>(), 1, vec );

}

Attention: Don't use n > 10....

Edit 3: the not for the faint of memoized heart approach:

Add these includes:

#include <unordered_map>
#include <map>

Change this in .cpp file:

template <typename ReturnType, typename... Args>
std::function<ReturnType (Args...)>
memoize(ReturnType (*func) (Args...))
{
    auto cache = std::make_shared<std::map<std::tuple<Args...>, ReturnType>>();
    return ([=](Args... args) mutable  {
        std::tuple<Args...> t(args...);
        if (cache->find(t) == cache->end())
            (*cache)[t] = func(args...);
        return (*cache)[t];
    });
}


template <typename F_ret, typename...  F_args>
std::function<F_ret (F_args...)> memoized_recursion(F_ret (*func)(F_args...))
{
    typedef std::function<F_ret (F_args...)> FunctionType;
    static std::unordered_map<decltype(func), FunctionType> functor_map;

    if(functor_map.find(func) == functor_map.end())
        functor_map[func] = memoize(func);

    return functor_map[func];
}

unsigned long fibonacci(unsigned n)
{
    return (n < 2) ? n :
    memoized_recursion(fibonacci)(n - 1) +
    memoized_recursion(fibonacci)(n - 2);
}

unsigned long long fibProduct( int n )
{
    vector<int> vec(n);
    int i = 0;
    generate( begin(vec), end(vec), [&i]() { return fibonacci(i += 1); } );
    return foldl( std::multiplies<int>(), 1, vec );
}

The rest is just proper compilation, linking and link installation in Mathematica...

Btw. you've to link explicitly against libstdc++.6.0.9.dylib (or another version number), since there is all the basic_string code etc. (MacOSX)

Addendum:

Has someone dealt so far, presumably yes, with big Integer issues with MathLink? Any tips for me?

Edit 4: the not for the faint of memoized heart MultiPrecision approach:

The .tm file:

:Begin:
:Function:      fibProductMP
:Pattern:       FibProductMP[n_]
:Arguments:     { n }
:ArgumentTypes: { Integer }
:ReturnType:    Manual
:End:

The .cpp file:

Add this include:

#include <boost/multiprecision/integer.hpp>

Maybe, if you'd like to avoid long namespace names:

using boost::multiprecision::cpp_int;

template<typename T>
inline T Identity(const std::multiplies<T>&) { return T(1); }

void fibProductMP( int n )
{
    vector<int> vec(n);
    int i = 0;
    generate( begin(vec), end(vec), [&i]() { return fibonacci(i += 1); } );
    cpp_int result = foldl( std::multiplies<cpp_int>(), Identity<cpp_int
        (multiplies<cpp_int>()), vec );
    MLNewPacket(stdlink);
    MLPutString(stdlink, result.str().c_str() );
}

Mathematica side:

link = Install["wherever the binary is", LinkMode -> Launch];

FibProductMP[30]

(* 607373569868916007005878071331449502263924414704952629297115029592606\
    043656028160000000 *)

What we've achieved is a multiprecision FibProduct function which is cool, decent C++11 code, exception safe and memoized. Now I'd be tempted to say victory!

P.S.: My life wouldn't be what it is, if there would not still exist one caveat. I return the result as a const char*; so a ToExpression is needed to reuse the result for further calculations.

Just can't get enough (until 47 ;) ... Gosh...could solve that as well, but I've to say, enough is enough :)

FibProductMP[#] & /@ Range[47] // TableForm

enter image description here

That was real fun...

Thank you.

share|improve this answer
    
Yes, it seemed as i didn't specify a value of a, in which it was not computing my loops correctly for While or For. Thank you very much for the prompt response! –  julesverne Feb 15 '13 at 10:43
    
I credit you with reviving my interest in C++! I abandoned it years ago when I got tired of spending a whole day updating dependencies every time I wanted to do 5 minutes worth of coding, but your code above just plain works (up to n = 40 at least!). –  Reb.Cabin Feb 15 '13 at 17:15
    
Thank you :) With C++11 out, C++ feels like a whole new language. Id didn't wanted to exaggerate things, but normally I'd leave out the "1" for foldl and insert: template<typename T> inline T Identity(const std::multiplies<T>&) {return T(1);} foldl(std::multiplies<int>(), Identity<int>(std::multiplies), vec ); That would be, mathematically, clearer. –  Stefan Feb 15 '13 at 17:51

For/Do/While Loops are generally to be avoided in Mathematica.

Instead this can be solved with a single line

Times @@ (fib /@ Range[20])

Range just generates a list of numbers from 1 to 20. The /@ is short for Map, which applies fib to every element of the list. The @@ is short for Apply, which gives you Times[...] instead of List[...].

share|improve this answer
2  
One might also use Fibonacci[n] (especially if you're lazy as I am...) –  Pinguin Dirk Feb 15 '13 at 10:44
    
Being lazy and loop averse too, I would suggest Times @@ Round[GoldenRatio^Range[20]/Sqrt[5]] –  KennyColnago Feb 16 '13 at 0:38

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.