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

[fermineutron]

>Barry Schwarz wrote:
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 do not work in C directly to allow the caller to use the same
argument for one of the multipliers as the result. Doing calculations
in C directly will screw up the math.

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

OOPS, it did not copy over. but i did have it in the code.

[Ancient_Hacker]

fermineutron wrote:

Ancient_Hacker wrote:

You asked this question before, and I gave you the answer.

You're doing one digit at a time. Doing it with 32 bits at a time will
be at least 400,000,000 times faster. Why,? because with 32 bits you
can store numbers up to about 4 billion, AND you don't need to do a
divide and a mod each time.

I am doing 3 digits at the time

Good! Then the speedup will only be by a factor of 400,000 times.

I would approach it this way:

First go to using base 2. Shouldnt be hard, just change all your 10's
to 2's.
Then get the code going as well as before, only using base 2.

Then you can try stepping up the number of digits from 3 up to whatever
will fit in your datatype. Better stop a few bits short the first time
to avoid problems with the sign bit.

Then notice that all your divides and mod's by 2^digits can be replaced
by shifts and masks. Although a good compiler will do this for you.

So by using most of the datatype width, and removing all the divs and
mods, things should get super speedy.

[fermineutron]

Ancient_Hacker wrote:

First go to using base 2. Shouldnt be hard, just change all your 10's
to 2's.
Then get the code going as well as before, only using base 2.

Then you can try stepping up the number of digits from 3 up to whatever
will fit in your datatype. Better stop a few bits short the first time
to avoid problems with the sign bit.

Then notice that all your divides and mod's by 2^digits can be replaced
by shifts and masks. Although a good compiler will do this for you.

So by using most of the datatype width, and removing all the divs and
mods, things should get super speedy.

Hmm, i forgot about super speedy shifts, going to rework the code now.

I know there are better agorithms for computing factorial then straight
up multiplication, but it is the bignum arithemetic that i am mostly
curious about so factoriial calculations is just a driver/benchmark.

[Arthur J. O'Dwyer]

On Thu, 2 Nov 2006, Ancient_Hacker wrote:

fermineutron wrote:

>Ancient_Hack er wrote:

>>You asked this question before, and I gave you the answer.

You're doing one digit at a time. Doing it with 32 bits at a time will
be at least 400,000,000 times faster. Why,? because with 32 bits you
can store numbers up to about 4 billion, AND you don't need to do a
divide and a mod each time.

I am doing 3 digits at the time

Good! Then the speedup will only be by a factor of 400,000 times.

I would approach it this way:

First go to using base 2. Shouldnt be hard, just change all your 10's
to 2's.
Then get the code going as well as before, only using base 2.

Then you can try stepping up the number of digits from 3 up to whatever
will fit in your datatype. Better stop a few bits short the first time
to avoid problems with the sign bit.

Then notice that all your divides and mod's by 2^digits can be replaced
by shifts and masks. Although a good compiler will do this for you.

So by using most of the datatype width, and removing all the divs and
mods, things should get super speedy.

Also, there are almost certainly faster ways to multiply than O(n^2).

However, when fermineutron goes to base 2^k, he will need to invest
some time (or Google searching) in a new "AEIfprintf " routine. This
shouldn't be a big deal --- his existing routine is totally brain-dead
and should be shot on sight --- but it might take a while to get all
the new bugs out, especially since he's working in an unfamiliar language.

Followups set to comp.programmin g, since this thread has lost any
connection to the C programming language.

-Arthur

[Keith Thompson]

"fermineutr on" <fr**********@y ahoo.comwrites:

Ancient_Hacker wrote:

[...]

>Then notice that all your divides and mod's by 2^digits can be replaced
by shifts and masks. Although a good compiler will do this for you.

Hmm, i forgot about super speedy shifts, going to rework the code now.

[...]

I predict that he will ignore the "Although a good compiler will do
this for you" statement, and just manually replace all his
multiplications and divisions by shift without first checking whether
it's necessary. (I could be wrong, of course.)

--
Keith Thompson (The_Other_Keit h) ks***@mib.org <http://www.ghoti.net/~kst>
San Diego Supercomputer Center <* <http://users.sdsc.edu/~kst>
We must do something. This is something. Therefore, we must do this.

[Keith Thompson]

"fermineutr on" <fr**********@y ahoo.comwrites:

>>Barry Schwarz wrote:
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 do not work in C directly to allow the caller to use the same
argument for one of the multipliers as the result. Doing calculations
in C directly will screw up the math.

For those tuning in late, "C" is the name of a variable. I suggest
that, like "a", "C" is a really bad variable name.

--
Keith Thompson (The_Other_Keit h) ks***@mib.org <http://www.ghoti.net/~kst>
San Diego Supercomputer Center <* <http://users.sdsc.edu/~kst>
We must do something. This is something. Therefore, we must do this.

[John Smith]

fermineutron wrote:

For a while now i have been "playing" with a little C program to
compute the factorial of large numbers.

What's the largest factorial you've been able to compute with
your program?

[fermineutron]

John Smith wrote:

What's the largest factorial you've been able to compute with
your program?

with old version the larges one i did in real time is 150000

after i swiched to base 1024 i have not tried yet, still stuck on
printing function, i seem to be unable to come up with algoryth to go
from base "base" to base 10 in a long array
Standard aproach such as
NumInBase10=Num InBaseB*(B/10)^N
will not work since N can be several thousand easy.

[fermineutron]

Here is the function i have for printing to the screen an array encoded
integer, but it does not work, would someone please fix it. I dont
usually ask people to write code for me completely but i fried my brain
trying to figure this one out today. Dont know whats wrong with me
today, but i just seem to be unable to figure this one out today.

here is the code
=============== ======Start paste========== ===========
void AEIprintf(struc t AEI *p) {
int fieldLength = (int)log10(base );
long fieldbase=pow(1 0,fieldLength);
char format[10];
size_t j = p->digits;
long dpf, dpf1;
sprintf(format, "%%0%dd", fieldLength);
dpf = 0;
dpf1 = (p->data[j] + dpf)/fieldbase;
dpf = (p->data[j] + dpf)%fieldbase;
printf("%*d", fieldLength, dpf1);
for (; j 0; j--) {
dpf1 = (p->data[j] + dpf)/fieldbase;
dpf = (p->data[j] + dpf)%fieldbase;
printf(format, dpf1);
}
return;
}

=============== ======End Paste========== ===========

[fermineutron]

can such a printing operation be done at all? I think there may be a
way if to 1st reencode the array into a decimal base and then print,
but honestly speaking i have no idea.