Long Num speed

[fermineutron]

For a while now i have been "playing" with a little C program to
compute the factorial of large numbers. Currently it takes aboy 1
second per 1000 multiplications , that is 25000P1000 will take about a
second. It will be longer for 50000P1000 as expected, since more digits
will be in the answer. Now, on the Num Analyses forum/Group there is a
post reaporting that that person wrot java code that computed 1000000!
in about a second. That is about 10000 times faste than I would expect
my code to do it. So the two possiblilities are:
1) I am doing something terribly wrong
2) The othr person is lying

At the moment i am inclined to believe that its number 1.

I am posing my code below, I would like to hear your opinions about
why it is slow and how i can improove its speed.

I know that there are public BIGNUM libraries which are already
optimized for such calculations, but I dont want to use them, bcause i
want to approach this problem on a lower level. I am mostly interested
to find out how to get this code perform faster or what alternative
algorythms i should consider. The factorial calculation is just a test
program.

=============== ====start paste========== =============== =

#include<stdio. h>
#include<stdlib .h>
#include<math.h >

#define al 1024*20
#define base 1000
typedef long int IntegerArrayTyp e;

struct AEI{
IntegerArrayTyp e data[al];
long int digits;
};

void pack(IntegerArr ayType i, struct AEI *N1);
void Amult(struct AEI * A, struct AEI * B, struct AEI * C);
void AEIprintf(struc t AEI * N1);
void AEIfprintf(FILE * fp, struct AEI * N1);
int main(void)
{

struct AEI *N1, *MO, *Ans;
long i = 0, j = 0, ii, NUM, iii;
FILE *ff;

N1=malloc(sizeo f(struct AEI));
MO=malloc(sizeo f(struct AEI));
Ans=malloc(size of(struct AEI));
while (i < al){
N1->data[i] = 0;
MO->data[i] = 0;
Ans->data[i]=0;
i++;
}

printf("Enter integer to Factorialize: ");
scanf("%ld", &NUM);

pack(1, N1);
pack(1, Ans);
ff = fopen("Results. txt", "w");
printf("you entered: %ld", NUM);

i=1;
while(i < NUM ){

iii=0;
while(iii<NUM && iii<1000){
ii = 1;
while (ii < al)
{
MO->data[ii] = 0;
ii++;
}
pack((i+iii), MO);
Amult(N1, MO, N1);
iii++;
}
i+=iii;
Amult(Ans, N1, Ans);
printf("\nProgr ess is: %d",i);
pack(1, N1);
}
if(ff!=NULL){
fprintf(ff,"\n% d\n",i-1);
AEIfprintf(ff, Ans);
}
fclose(ff);

printf("\nProgr ess: 100\%");

return 0;
}
void AEIprintf(struc t AEI *N1){

float fieldLength;
double temp;
char format1[8];
long j, FL0;
j = N1->digits-1;
FL0=(long)log10 ((float)base);
fieldLength = (float)log10((f loat)base);
temp = modf(fieldLengt h, &fieldLength );
format1[0] = '%';
format1[1] = '0';
format1[2] = fieldLength + 48;
format1[3] = 'd';
format1[4] = 0x00;

printf("%*d", FL0, N1->data[j]);
j--;

while (j >= 0)
{
printf(format1, N1->data[j]);

j--;
}

return;
}
void AEIfprintf(FILE * fp, struct AEI *N1){
long j = N1->digits-1;

double fieldLength, temp;
char format0[8], format1[8];

fieldLength = (int)log10(base );
temp = modf(fieldLengt h, &fieldLength );

format0[0] = '%';
format0[1] = fieldLength + 48;
format0[2] = 'd';
format0[3] = 0x00;
format1[0] = '%';
format1[1] = '0';
format1[2] = fieldLength + 48;
format1[3] = 'd';
format1[4] = 0x00;

fprintf(fp,form at0, N1->data[j]);
j--;

while (j >= 0){
fprintf(fp, format1, N1->data[j]);
j--;
}
return;
}

void pack(IntegerArr ayType i, struct AEI * N1)
{
long t = 1, i1, j = 0;

while (t == 1){
i1 = i % base;
N1->data[j] = i1;
i = (i - i1) / base;
j++;
if (i == 0)
t = 0;
}
N1->digits=j;
return;
}


void Amult(struct AEI * A, struct AEI * B, struct AEI * C){
/*C = A * B; */
long i, ii,d, result, carry=0, digits=0;
struct AEI *Ans;
Ans=malloc(size of(struct AEI));
i=0;
d= (A->digits+B->digits-1);
while(i<d){
Ans->data[i]=carry;
carry=0;
ii=0;
while(ii<=i){
if(B->data[ii]!=0){
Ans->data[i]+=A->data[i-ii]*B->data[ii];
carry+=Ans->data[i]/base;
Ans->data[i]=Ans->data[i]%base;
}
ii++;
}
carry+=Ans->data[i]/base;
Ans->data[i]=Ans->data[i]%base;

i++;
}
if(carry!=0){
d++;
Ans->data[i]=carry;
}

C->digits=d;
i=0;
while(i<d){
C->data[i]=Ans->data[i];
i++;
}
return;
}

=============== =====end paste========== =============== ===

I tried to indent the code with spaces instead of tabs, but if some
parts end up not properly indented, I hope no one will hold it against
me.

Thanks ahead

[Ian Collins]

fermineutron wrote:

For a while now i have been "playing" with a little C program to
compute the factorial of large numbers. Currently it takes aboy 1
second per 1000 multiplications , that is 25000P1000 will take about a
second. It will be longer for 50000P1000 as expected, since more digits
will be in the answer. Now, on the Num Analyses forum/Group there is a
post reaporting that that person wrot java code that computed 1000000!
in about a second. That is about 10000 times faste than I would expect
my code to do it. So the two possiblilities are:
1) I am doing something terribly wrong
2) The othr person is lying

At the moment i am inclined to believe that its number 1.

I am posing my code below, I would like to hear your opinions about
why it is slow and how i can improove its speed.

Before anyone reviews the code, have you profiled it?

If not, why not? If you have, where were the bottlenecks?

--
Ian Collins.

[Frederick Gotham]

fermineutron:

Now, on the Num Analyses forum/Group there is a
post reaporting that that person wrot java code that computed 1000000!
in about a second. That is about 10000 times faste than I would expect
my code to do it. So the two possiblilities are:
1) I am doing something terribly wrong
2) The othr person is lying

The burden of proof is on the Java dude.

Probability suggests that he's lying, because only a very small proportion of
proficient programmers waste their time on mickey-mouse hold-my-hand
languages such as Java. Also, Java is slow.

Either that or the algorithm's something stupid like:

char const *Func(void)
{
return "23495072395732 584712579750932 750923750932759 2387509";
}

Ask the Java dude if you can see his code. If he refuses, assume that he's a
liar, then egg his house.

--

Frederick Gotham

[fermineutron]

Ian Collins wrote:

Before anyone reviews the code, have you profiled it?

If not, why not? If you have, where were the bottlenecks?

I profilled it, but there were no obvious bottlenecks which i would not
anticipate to be there by design.

here is the profiler output

http://igorpetrusky.awardspace.com/Temp/RunStats.html

I was thinking that maybe there is some other algorythm that is better
than mine for the long int arithemetic?

Factorial calculation is just a driver for the multiplication function.

[Barry Schwarz]

On 1 Nov 2006 14:28:20 -0800, "fermineutr on" <fr**********@y ahoo.com>
wrote:

>For a while now i have been "playing" with a little C program to
compute the factorial of large numbers. Currently it takes aboy 1
second per 1000 multiplications , that is 25000P1000 will take about a
second. It will be longer for 50000P1000 as expected, since more digits
will be in the answer. Now, on the Num Analyses forum/Group there is a
post reaporting that that person wrot java code that computed 1000000!
in about a second. That is about 10000 times faste than I would expect
my code to do it. So the two possiblilities are:
1) I am doing something terribly wrong
2) The othr person is lying

At the moment i am inclined to believe that its number 1.

I am posing my code below, I would like to hear your opinions about
why it is slow and how i can improove its speed.

I know that there are public BIGNUM libraries which are already
optimized for such calculations, but I dont want to use them, bcause i
want to approach this problem on a lower level. I am mostly interested
to find out how to get this code perform faster or what alternative
algorythms i should consider. The factorial calculation is just a test
program.

I don't know if my comment in Amult addresses the question you voiced
but you need to pay attention to your compiler diagnostics.

>
============== =====start paste========== =============== =

#include<stdio .h>
#include<stdli b.h>
#include<math. h>

#define al 1024*20
#define base 1000
typedef long int IntegerArrayTyp e;

struct AEI{
IntegerArrayTyp e data[al];
long int digits;
};

void pack(IntegerArr ayType i, struct AEI *N1);
void Amult(struct AEI * A, struct AEI * B, struct AEI * C);
void AEIprintf(struc t AEI * N1);
void AEIfprintf(FILE * fp, struct AEI * N1);
int main(void)
{

struct AEI *N1, *MO, *Ans;
long i = 0, j = 0, ii, NUM, iii;
FILE *ff;

N1=malloc(sizeo f(struct AEI));
MO=malloc(sizeo f(struct AEI));
Ans=malloc(size of(struct AEI));
while (i < al){
N1->data[i] = 0;
MO->data[i] = 0;
Ans->data[i]=0;
i++;
}

printf("Enter integer to Factorialize: ");

Please learn to indent consistently. This is like trying to roller
skate on cobblestones.

snip rest of function until

>
printf("\nProgr ess: 100\%");

The correct way to include a '%' in a printf string is %%, not \%.
Didn't your compiler complain about an undefined escape sequence?

>
return 0;
}
void AEIprintf(struc t AEI *N1){

float fieldLength;
double temp;
char format1[8];
long j, FL0;
j = N1->digits-1;
FL0=(long)log10 ((float)base);
fieldLength = (float)log10((f loat)base);
temp = modf(fieldLengt h, &fieldLength );

This is a constraint violation. The second argument should be a
double*, not a float*. During execution it will invoke undefined
behavior. Didn't your compiler complain about a type mismatch?

>

format1[0] = '%';
format1[1] = '0';
format1[2] = fieldLength + 48;

What is 48. You probably mean '0' and should use that.

> format1[3] = 'd';
format1[4] = 0x00;

While the hex value is correct, '\0' is more in keeping with common C
idiom.

>
printf("%*d", FL0, N1->data[j]);
j--;

while (j >= 0)
{
printf(format1, N1->data[j]);

j--;
}

return;
}

snip

>void Amult(struct AEI * A, struct AEI * B, struct AEI * C){
/*C = A * B; */
long i, ii,d, result, carry=0, digits=0;
struct AEI *Ans;
Ans=malloc(size of(struct AEI));
i=0;
d= (A->digits+B->digits-1);
while(i<d){
Ans->data[i]=carry;
carry=0;
ii=0;
while(ii<=i){
if(B->data[ii]!=0){
Ans->data[i]+=A->data[i-ii]*B->data[ii];
carry+=Ans->data[i]/base;
Ans->data[i]=Ans->data[i]%base;
}
ii++;
}
carry+=Ans->data[i]/base;
Ans->data[i]=Ans->data[i]%base;

i++;
}
if(carry!=0){
d++;
Ans->data[i]=carry;
}

C->digits=d;
i=0;
while(i<d){
C->data[i]=Ans->data[i];

You could have simply said *C = *Ans and let the compiler do this
possibly more efficiently.

Is there a reason you did not do the work in C directly (and avoid the
time spent in malloc and copying *Ans to *C)?

i++;

Why do you go to so much trouble to avoid for loops?
for (i = 0; i < d; i++)
C->...;

}
return;

Where do you free Ans? You are causing repeated memory leaks.

>}

============== ======end paste========== =============== ===

I tried to indent the code with spaces instead of tabs, but if some
parts end up not properly indented, I hope no one will hold it against
me.

Thanks ahead

Remove del for email

[Richard Heathfield]

Barry Schwarz said:

On 1 Nov 2006 14:28:20 -0800, "fermineutr on" <fr**********@y ahoo.com>
wrote:

<snip>

>>
float fieldLength;
double temp;
char format1[8];
long j, FL0;
j = N1->digits-1;
FL0=(long)log 10((float)base) ;
fieldLength = (float)log10((f loat)base);
temp = modf(fieldLengt h, &fieldLength );

This is a constraint violation. The second argument should be a
double*, not a float*. During execution it will invoke undefined
behavior. Didn't your compiler complain about a type mismatch?

In my experience, fermineutron pays no heed to correctness issues, and is
concerned only with speed. He has a track record of ignoring corrections. I
think you're wasting your time, Barry.

--
Richard Heathfield
"Usenet is a strange place" - dmr 29/7/1999
http://www.cpax.org.uk
email: rjh at above domain (but drop the www, obviously)

[Samuel Stearley]

Hi
The math software in my ti89 calculator with a 12 MHz 68k CPU can
calculate all 613 digits of 299! in 1 second.
It then takes 5 seconds to convert that to a string for display. But
for giggles I once wrote my own display routine that does the integer
to string conversion in 1/3 of a second.

Frederick Gotham wrote:

fermineutron:

Now, on the Num Analyses forum/Group there is a
post reaporting that that person wrot java code that computed 1000000!
in about a second. That is about 10000 times faste than I would expect
my code to do it. So the two possiblilities are:
1) I am doing something terribly wrong
2) The othr person is lying

The burden of proof is on the Java dude.

Probability suggests that he's lying, because only a very small proportion of
proficient programmers waste their time on mickey-mouse hold-my-hand
languages such as Java. Also, Java is slow.

Either that or the algorithm's something stupid like:

char const *Func(void)
{
return "23495072395732 584712579750932 750923750932759 2387509";
}

Ask the Java dude if you can see his code. If he refuses, assume that he's a
liar, then egg his house.

--

Frederick Gotham

[websnarf]

fermineutron wrote:

Ian Collins wrote:

Before anyone reviews the code, have you profiled it?

If not, why not? If you have, where were the bottlenecks?

Its a factorial calculation. The bottleneck is in the big integer
multiply, which itself should have a bottleneck in the platform's
multiply.

I profilled it, but there were no obvious bottlenecks which i would not
anticipate to be there by design.

here is the profiler output

http://igorpetrusky.awardspace.com/Temp/RunStats.html

I was thinking that maybe there is some other algorythm that is better
than mine for the long int arithemetic?

Probably so. Consider that you are doing nothing more than computing
the straight product of the numbers using no arithmetic short cuts at
all.

So here's what comes off the top of my head: Ask yourself the
following problem. How many factors of 2 are there in 1000000! ?
Certainly every other number is even. But every 4th number has 2
factors of 2, and every 8th number has 3 factors of two in it. So the
answer is:

f(2) = floor(1000000/2) + floor(1000000/4) + floor(1000000/8) + ...

Similarly we can figure out the number of factors of 3s, 5s, 7s, 11s,
and all the primes less than 1000000, as f(3), f(5), etc. Then the
result you are looking for is:

1000000! = pow(2,f(2))*pow (3,f(3))*pow(5, f(5))*...

Now, the question is -- what makes us think this will be any faster?
Well, the pow() function can be computed with successive squaring
tricks. Squaring faster than straight multiplying because (a*q^r + b)
^ 2 = a^2*q^(2*r) + b^2 + 2*a*b*(q^r). And the resulting big number
multiplies that you have to perform here can be accelerated using any
number of big number multiply acceleration tricks (Karatsuba,
Toom-Cook, or DFTs.).

I don't know how much faster, if any, doing things this way would be.
If you find a faster way in the literature, I would be interested in
know it.

--
Paul Hsieh
http://www.pobox.com/~qed/
http://bstring.sf.net/

[websnarf]

fermineutron wrote:

Ian Collins wrote:

Before anyone reviews the code, have you profiled it?

If not, why not? If you have, where were the bottlenecks?

Its a factorial calculation. The bottleneck is in the big integer
multiply, which itself should have a bottleneck in the platform's
multiply.

I profilled it, but there were no obvious bottlenecks which i would not
anticipate to be there by design.

here is the profiler output

http://igorpetrusky.awardspace.com/Temp/RunStats.html

I was thinking that maybe there is some other algorythm that is better
than mine for the long int arithemetic?

Probably so. Consider that you are doing nothing more than computing
the straight product of the numbers using no arithmetic short cuts at
all.

So here's what comes off the top of my head: Ask yourself the
following problem. How many factors of 2 are there in 1000000! ?
Certainly every other number is even. But every 4th number has 2
factors of 2, and every 8th number has 3 factors of two in it. So the
answer is:

f(2) = floor(1000000/2) + floor(1000000/4) + floor(1000000/8) + ...

Similarly we can figure out the number of factors of 3s, 5s, 7s, 11s,
and all the primes less than 1000000, as f(3), f(5), etc. Then the
result you are looking for is:

1000000! = pow(2,f(2))*pow (3,f(3))*pow(5, f(5))*...

Now, the question is -- what makes us think this will be any faster?
Well, the pow() function can be computed with successive squaring
tricks. Squaring faster than straight multiplying because (a*q^r + b)
^ 2 = a^2*q^(2*r) + b^2 + 2*a*b*(q^r). And the resulting big number
multiplies that you have to perform here can be accelerated using any
number of big number multiply acceleration tricks (Karatsuba,
Toom-Cook, or DFTs.).

I don't know how much faster, if any, doing things this way would be.
If you find a faster way in the literature, I would be interested in
knowing it.

--
Paul Hsieh
http://www.pobox.com/~qed/
http://bstring.sf.net/

[Chris Dollin]

Frederick Gotham wrote:

Probability suggests that he's lying, because only a very small proportion of
proficient programmers waste their time on mickey-mouse hold-my-hand
languages such as Java. Also, Java is slow.

Can we drop the insults, please?

--
Chris "unhashedup hashed up hashing" Dollin
"I'm still here and I'm holding the answers" - Karnataka, /Love and Affection/