Big Numbers Math

[woodybrison]

I've created a method for handling huge numbers, perhaps it might be useful to someone else. Maybe someone can find an error in it...

I was up against a problem: in a set of 1 million elements, how many ways can they be paired? The first el could be paired with 999,999 other choices. For each of those, the 2nd el could be paired with 999,997 other choices, so 999,999 x 999,997 x 999,995 x ... x 1. There would be a half million pairs:

     499,999
      ∏ 999999-2i  = ?
     i=0

Wolfram/Alpha was unable to compute this.

I wondered if Free Pascal could do it.

Free Pascal's largest type (that I know of) is Extended, with a max value of 1.1x10^4932.

I set up a type with separate mantissa (an Extended) and exponent (a uint64).

type bignum = record mant :Extended; expn :uint64; end;

This might be more useful if the exponent be a signed int64 so it could handle very small numbers also.

To multiply two such numbers, you add the exponents and multiply the mantissas, and if the mantissas' product is larger than some limit you divide it by ten and add one to the exponent.

The program below computes the above product. I ran a test, shown at the end of the program - that worked; and I followed the output, it looks correct. The answer it computes is:

     499,999
      ∏ 999999-2i  ≈  8.12014 x 10^2782852
     i=0

This seems to be correct. It is smaller than 1M^500K, which is 10^3M.

The program took 1 second to finish on a 2 GHz Pentium.

Code: Pascal  [Select][+][-]

1. Program bigmath;
2. // handles very large numbers
3.  
4. type
5.   bignum = record
6.     mant :Extended;
7.     expn :uint64;
8.     end;
9.  
10. var
11. //  AA, BB, prod :bignum;  // test
12.   i :uint32;
13.   accum, multiplier :bignum;
14.  
15. //====================================
16. procedure adjust( var X :bignum );
17.   begin
18.   while X.mant >= 10 do
19.     begin
20.     X.mant /= 10;
21.     inc( X.expn );
22.     end;
23.   end;
24.  
25. //====================================
26. function multiply( A,B :bignum ) :bignum;
27.   var temp :bignum;
28.   begin
29.   temp.mant := 1;
30.   temp.expn := 0;
31.   temp.mant := A.mant * B.mant;
32.   temp.expn := A.expn + B.expn;
33.   adjust( temp );
34.   exit( temp );
35.   end;
36.  
37. //=============== MAIN ===============
38. begin
39. accum.mant := 1;  // accumulates the product
40. accum.expn := 0; //  set it to 1 to begin with
41.  
42. for i := 0 to 499999 do
43.   begin
44.   multiplier.mant := 999999-2*i;
45.   multiplier.expn := 0;
46.   adjust( multiplier );
47.  
48.   accum := multiply( accum, multiplier );
49.  
50.   // tracking progress:
51.   if (i < 50) or (i > 499949) or (i mod 500 = 0) then
52.     writeln( i:7, ' acc=', accum.mant:0:5, 'x10^', accum.expn,
53.              ' multi=', multiplier.mant:0:5, 'x10^', multiplier.expn );
54.   end;
55.  
56. // displaying the final answer:
57. writeln( '        acc=', accum.mant:0:5, 'x10^', accum.expn,
58.          ' multi=', multiplier.mant:0:5, 'x10^', multiplier.expn, ' <====' );
59. end.
60.  
61. // test: 999,999 * 999,997 = 999996000003
62. //
63. // AA.mant := 999999;
64. // AA.expn := 0;
65. // BB.mant := 999997;
66. // BB.expn := 0;
67. // prod := multiply( AA, BB );
68. // write( AA.mant:0:0, 'x10^', AA.expn,
69. //  ' x ', BB.mant:0:0, 'x10^', BB.expn,
70. //  ' = ', prod.mant:0:19, 'x10^', prod.expn );
71. //
72. // test result
73. // 999999x10^0 x 999997x10^0 = 9.9999600000300024000x10^11
74. // note the round-off error creeping up
75.  

[AlexTP]

Why not to use forum CODE tag for the code block?

[woodybrison]

I plead ignorance. Now I've put the tags in. I had to snuffle around a bit to find out how. My only motive here is the help the universe.

[Thaddy]

Well, the notation seems a good idea, but note there are many good biginteger - and so by extension through scaling also floats - for Object Pascal. Even one or two in the standard distribution. Their range is usually in principle infinite at the cost of being exponentially very slow.

[iLya2IK]

Hello! Here is my approach

Code: Pascal  [Select][+][-]

1. Program project1;
2. // handles very large numbers
3.  
4. type
5.   bignum = record
6.     mant :Extended;
7.     expn :uint64;
8.     end;
9.  
10. var
11. //  AA, BB, prod :bignum;  // test
12.   i :uint32;
13.   accum, multiplier :bignum;
14.  
15. //====================================
16. procedure multiply( var A :bignum; B :bignum );
17.   begin
18.   A.mant *= B.mant;
19.   A.expn += B.expn;
20.   if (A.mant >= 10.0) then
21.   begin
22.     A.mant/=10.0;
23.     Inc(A.expn);
24.   end;
25.   end;
26.  
27. const MAX_ORDER = 5;
28. const dividers : array [0..MAX_ORDER] of extended   = (100000, 10000, 1000, 100, 10, 1);
29. const exps     : array [0..MAX_ORDER] of integer    = (5,      4,     3,    2,   1,  0);
30. const multip   : array [0..MAX_ORDER+1] of integer  = (999999, 99999, 9999, 999, 99, 9, 0);
31. //=============== MAIN ===============
32. var
33.   m, mn:integer;
34. begin
35. accum.mant := 1;  // accumulates the product
36. accum.expn := 0; //  set it to 1 to begin with
37.  
38. mn := multip[0];
39. for i := 0 to MAX_ORDER do
40. begin
41.   m  := mn;
42.   mn := multip[i + 1];
43.   while (m > mn) do
44.   begin
45.     multiplier.mant := Extended(m) / dividers[i];
46.     multiplier.expn := exps[i];
47.     multiply( accum, multiplier );
48.  
49.     m := m - 2;
50.   end;
51.   writeln( i:7, ' acc=', accum.mant:0:20, 'x10^', accum.expn,
52.              ' multi=', multiplier.mant:0:20, 'x10^', multiplier.expn );
53. end;
54.  
55. // displaying the final answer:
56. writeln( '        acc=', accum.mant:0:20, 'x10^', accum.expn,
57.          ' multi=', multiplier.mant:0:20, 'x10^', multiplier.expn, ' <====' );
58. end.
59.                
60.  

The result is:
(for original code)
acc=8.12013661044774699965x10^2782852 multi=1.00000000000000000000x10^0 <====
(for modified code)
acc=8.12013661044774706643x10^2782852 multi=1.00000000000000000000x10^0 <====