factorial function

[CJC0117]

Can someone tell me why my factorial function gives different output than wolframalpha.com after 22! ? Is double not precise enough?

#include <iostream>
#include <sstream>
#include <iomanip>
#include <limits>

int main()
{
    double factorial(int number);
    std::cout<<std::fixed<<std::setprecision(0)<<factorial(23)<<std::endl;
    std::cout<<std::endl<<"Press ENTER to continue...";
    std::cin.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
    return 0;
}

double factorial(int number)
{
    int n;
    double factorial=number;
    for(n=number;n>1;n--)
    {
        factorial*=(n-1);
    }
    return factorial;
}

[Moschops]

It's not precise enough.

http://en.wikipedia.org/wiki/Double-precision_floating-point_format

22! is roughly 2^70

[CJC0117]

Is there anything more precise? I've tried long double, but it's still not enough to calculate 100!, which is ultimately what I need to do.

EDIT: will I need to use arbitrary-precision arithmetic?

[firedraco]

Use a bignum library.

[htirwin]

You don't need to be using floating point variables for factorials. No matter what a factorial is going to be an integer.

You can use unsigned long long, but still probably not long enough for 100!. Or like firedraco said, you can use bignum library.

I seams that unsigned long long will hold up to 65!

[Duthomhas]

An unsigned long long is 64 bits? The maximum factorial you can get in that is 20!

If your unsigned long long is 128 bits then the maximum you can get in that is 34!

[edit] If you are using floating point numbers you already have only approximations to the correct number after 18!.

[edit 2] I can't resist: http://www.cplusplus.com/forum/lounge/32041/

[CJC0117]

Thanks for all of the replies.

@firedraco, I probably will have to use a bignum library eventually but for now I want to try to find my own solution. The reason I am doing this is actually to solve problem 20 of project euler.

@iseeplusplus, unsigned long long seems to have less precision when I try it.

@ Duoas, thanks for the link, but I don't want to look at that thread quite yet. First I want to see if I am in the right direction with the following source code, which I have spent the last few hours writing. It's ugly and it doesn't seem to work well though.

EDIT: Nevermind, I finally got it to work. Here's my finished code:

#include <iostream>
#include <cmath>
#include <vector>
#include <sstream>
#include <iomanip>
#include <limits>

int main()
{
    int i;
    std::stringstream ss;
    std::string s,product="100";
    std::string multiply(std::string a,std::string b);
    for(i=100;i>1;i--)
    {
        ss<<i-1;
        ss>>s;
        ss.clear();
        product=multiply(product,s);
    }
    std::cout<<product<<std::endl;
    std::cout<<std::endl<<"Press ENTER to continue...";
    std::cin.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
    return 0;
}

std::string multiply(std::string a,std::string b)
{
    int i,j,k,digit1,digit2,carry=0,sum=0;
    std::stringstream ss;
    std::string s,product;
    std::vector<std::vector<int> >rows;
    std::vector<int>row;
    std::vector<int>reverseProduct;
    for(i=(int) b.size()-1,k=0;i>=0;i--,k++)
    {
        for(j=0;j<k;j++)
            row.push_back(0);
        for(j=(int) a.size()-1;j>=0;j--)
        {
            ss<<b[i]<<" "<<a[j];
            ss>>digit1>>digit2;
            ss.clear();
            row.push_back((digit1*digit2+carry)%10);
            carry=floor((digit1*digit2+carry)/10);
        }
        if(carry!=0)
        {
            row.push_back(carry);
            carry=0;
        }
        rows.push_back(row);
        row.clear();
    }
    for(i=0;i<(int) rows.size()-1;i++)
    {
        k=rows[rows.size()-1].size()-(int) rows[i].size();
        for(j=0;j<k;j++)
            rows[i].push_back(0);
    }
    carry=0;
    for(j=0;j<(int) rows[0].size();j++)
    {
        for(i=0;i<(int) rows.size();i++)
            sum+=rows[i][j];
        reverseProduct.push_back((sum+carry)%10);
        carry=floor((sum+carry)/10);
        sum=0;
    }
    if(carry!=0)
        reverseProduct.push_back(carry);
    for(i=(int) reverseProduct.size()-1;i>=0;i--)
        ss<<reverseProduct[i];
    ss>>product;
    return product;
}