Doing arbitrary precision (bignum) arithmetic in Lua

Posted by Nick Gammon on Mon 02 Apr 2007 06:16 AM — 22 posts, 80,447 views.

Nick Gammon Australia Forum Administrator Mon 02 Apr 2007 06:16 AM #0

I am pleased to release a copy of the bc (bignum) library.

You can download the DLL, source files, and the GPL document from:

http://www.gammon.com.au/files/mushclient/lua5.1_extras/bc.zip

All you really need is bc.dll, which you can place in the same directory as MUSHclient.exe.

Then you can load the DLL and use it along these lines:


-- test bc library

assert (package.loadlib ("bc.dll","luaopen_bc")) ()

------------------------------------------------------------------------------
print(bc.version)

------------------------------------------------------------------------------

print"Factorials"

function factorial(n,f)
 for i=2,n do f=f*i end
 return f
end

one=bc.number(1)
for i=1,40 do
 print(i,factorial(i,1),factorial(i,one))
end

This test demonstrates doing factorials. First it uses a straight Lua number, secondly a "bc" number. The results show the Lua number quickly become floating point numbers, whereas the bc number retains full accuracy:


bc library for Lua 5.1 / Aug 2006 / based on GNU bc-1.06
Factorials
1 1 1
2 2 2
3 6 6
4 24 24
5 120 120
6 720 720
7 5040 5040
8 40320 40320
9 362880 362880
10 3628800 3628800
11 39916800 39916800
12 479001600 479001600
13 6227020800 6227020800
14 87178291200 87178291200
15 1307674368000 1307674368000
16 20922789888000 20922789888000
17 3.55687428096e+014 355687428096000
18 6.402373705728e+015 6402373705728000
19 1.2164510040883e+017 121645100408832000
20 2.4329020081766e+018 2432902008176640000
21 5.1090942171709e+019 51090942171709440000
22 1.1240007277776e+021 1124000727777607680000
23 2.5852016738885e+022 25852016738884976640000
24 6.2044840173324e+023 620448401733239439360000
25 1.5511210043331e+025 15511210043330985984000000
26 4.0329146112661e+026 403291461126605635584000000
27 1.0888869450418e+028 10888869450418352160768000000
28 3.0488834461171e+029 304888344611713860501504000000
29 8.8417619937397e+030 8841761993739701954543616000000
30 2.6525285981219e+032 265252859812191058636308480000000
31 8.2228386541779e+033 8222838654177922817725562880000000
32 2.6313083693369e+035 263130836933693530167218012160000000
33 8.6833176188119e+036 8683317618811886495518194401280000000
34 2.952327990396e+038 295232799039604140847618609643520000000
35 1.0333147966386e+040 10333147966386144929666651337523200000000
36 3.719933267899e+041 371993326789901217467999448150835200000000
37 1.3763753091226e+043 13763753091226345046315979581580902400000000
38 5.230226174666e+044 523022617466601111760007224100074291200000000
39 2.0397882081197e+046 20397882081197443358640281739902897356800000000
40 8.159152832479e+047 815915283247897734345611269596115894272000000000

Nick Gammon Australia Forum Administrator Mon 02 Apr 2007 10:52 PM #1

Operations supported by the bc library are:

bc.number (n) - make a new "big" number. eg.
```
i = bc.number (2)
i = i ^ 50
print (i)  --> 1267650600228229401496703205376
```
The argument can be a string, eg.

PI = bc.number ("3.141592653589793238462643383279502884197")
bc.digits (n) - do subsequent calculations to n digits after the decimal point. eg.
```
i = bc.sqrt (2)
print (i)  --> 1

bc.digits (40)
i = bc.sqrt (2)
print (i)  --> 1.4142135623730950488016887242096980785696
```
bc.digits returns what the number of digits was previously (so you could save and restore it).

bc.compare (a, b) - returns -1 if a < b, 0 if a == b, and 1 if a > b. eg.


print (bc.compare (1, 2)) --> -1
print (bc.compare (2, 2)) --> 0
print (bc.compare (3, 2)) --> 1

bc.add (a, b) -- adds a + b
bc.sub (a, b) -- subtracts a - b
bc.div (a, b) -- divides a / b
bc.mod (a, b) -- returns a modulus b
bc.mul (a, b) -- multiplies a * b
bc.isneg (a) -- returns true if a is negative, false otherwise
bc.iszero (a) -- returns true if a is zero, false otherwise
bc.pow (a, b) -- raises a to the power b
bc.sqrt (a) -- returns the square root of a
bc.tostring (a) -- converts its argument to a string (eg. for printing)

When you create a "big number" with bc.number, you get a userdata with a metatable associated with it. The metatable supports the following operations:

__add
__div
__eq
__lt
__mul
__pow
__sub
__tostring
__unm

This means you can directly add, subtract, divide etc. big numbers (eg. a = b + c) without having to use functions like bc.add.

Nick Gammon Australia Forum Administrator Wed 18 Apr 2007 12:45 AM #2

This library is now included natively in MUSHclient 4.03 onwards. In other words, you can use it without needing to download or install the DLL.

Nick Gammon Australia Forum Administrator Wed 18 Apr 2007 05:31 AM #3

Out of interest, this algorithm seems to calculate pi correctly - I'm not sure exactly where the author (Luiz Henrique de Figueiredo) got it.

He says:

http://pauillac.inria.fr/algo/bsolve/constant/pi/pi.html

But I can't see it on that page directly.



bc.digits(1000)   -- calculate 1000 digits of pi

function calc_pi ()
 local x = bc.sqrt (2)
 local p = 2 + x
 local y = bc.sqrt (x)
 x = (y + 1 / y) / 2
 p = p * (x + 1) / (y + 1)
 for i = 1, 200 do
  local P = p
  local t = bc.sqrt (x)
  y = (y * t + 1 / t) / (y + 1)
  x = (t + 1 / t) / 2
  p = p * (x + 1) / (y + 1)
  if p == P then 
    break 
  end
 end
 return p
end

print (calc_pi ())

Output
3.141592653589793238462643383279502884197169399375
10582097494459230781640628620899862803482534211706
79821480865132823066470938446095505822317253594081
28481117450284102701938521105559644622948954930381
96442881097566593344612847564823378678316527120190
91456485669234603486104543266482133936072602491412
73724587006606315588174881520920962829254091715364
36789259036001133053054882046652138414695194151160
94330572703657595919530921861173819326117931051185
48074462379962749567351885752724891227938183011949
12983367336244065664308602139494639522473719070217
98609437027705392171762931767523846748184676694051
32000568127145263560827785771342757789609173637178
72146844090122495343014654958537105079227968925892
35420199561121290219608640344181598136297747713099
60518707211349999998372978049951059731732816096318
59502445945534690830264252230825334468503526193118
81710100031378387528865875332083814206171776691473
03598253490428755468731159562863882353787593751957
78185778053217122680661300192787661119590921642019
89

Nick Gammon Australia Forum Administrator Wed 18 Apr 2007 06:21 AM #4

There are times when you might need greater precision than usual. For example, this simple loop gives somewhat surprising results:


for i = -1, 1, 0.1 do 
  print (i) 
end

Output
-1
-0.9
-0.8
-0.7
-0.6
-0.5
-0.4
-0.3
-0.2
-0.1
-1.3877787807814e-016
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1

The value for zero is wrong (it is close to zero, but not exactly zero).

Also it turns out that the final number isn't really 1. If you add "print (i == 1)" to the loop, it always prints false.

However using the bc library you get more precise results:


i = bc.number (-1)
bc.digits (1)

repeat
  print (i)
  i = i + 0.1
until i > bc.number (1)

Output

-1
-0.9
-0.8
-0.7
-0.6
-0.5
-0.4
-0.3
-0.2
-0.1
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0

A test at the end of this loop shows that the final number is indeed exactly 1.

Nick Gammon Australia Forum Administrator Wed 18 Apr 2007 06:53 AM #5

Compare the digits of pi printed above to this page:

http://www.joyofpi.com/pi.html

The first 1000 seem the same to me. :)

Nick Gammon Australia Forum Administrator Wed 18 Apr 2007 07:37 AM #6

And if you want to fill your screen with numbers, try this:


print (bc.number (2) ^ 10000)

I must admit that this takes a while to calculate:


print (bc.number (2) ^ 1000000)

Took me 13 seconds. :)

Still, that filled up the output window with 3,763 lines of numbers.

David Haley USA Wed 18 Apr 2007 07:39 AM #7

Hey, I can do that without using your fancy schmancy library:


for i = 1, 1000000 do
  for j = 1, 1000000 do
    for k = 1, 1000000 do
      print("123")
    end
  end
end

:-)

Nick Gammon Australia Forum Administrator Wed 18 Apr 2007 07:41 AM #8

Yes, but my numbers include the digits 4 through to 9. :P

David Haley USA Wed 18 Apr 2007 07:48 AM #9

Foiled again!

[in a dark, ominous voice:]
I'll be back!

Nick Gammon Australia Forum Administrator Wed 25 Apr 2007 07:18 AM #10

Here's another interesting algorithm. This calculates e to 100 decimal places:


bc.digits (103)

n = bc.number (1)
e = bc.number (1)

for i = 1, 200 do
  local E = e
  n = n * i   -- n is i! (factorial of i)
  e = e + (1 / n)
  if e == E then
    break
  end -- if no change
end -- for

print (e) 

Output
 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274229

Note that the last 3 digits are wrong - to get the first 100 right I had to calculate to 103 digits of precision.

A debugging display shows it stops on the 72nd iteration, which means it got it right on the 71st iteration.

Nick Gammon Australia Forum Administrator Wed 25 Apr 2007 07:53 AM #11

This algorithm calculates sin (x) to a specified precision - the precision is the count of times you want to go through the loop:


bc.digits (100)

function sine (x, precision)
   local val = bc.number (1)
   x = bc.number (x)
   while precision > 0 do
       val = 1 - val * x * x / (2 * precision) / (2 * precision + 1)
       precision = precision - 1
   end -- while
   val = x * val
   return val
end -- function sine

print (sine (1, 35))

Output

0.8414709848078965066525023216302989996225630607983710656727517099919104043912396689486397435430526959

Nick Gammon Australia Forum Administrator Wed 25 Apr 2007 07:58 AM #12

Given the sine, we can calculate the cosine:


pi = bc.number (
"3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706")

function cosine (x, precision)
  return sine (pi / 2 - x, precision)
end -- cosine

print (cosine (1, 35))

Output

0.5403023058681397174009366074429766037323104206179222276700972553811003947744717645179518560871830860

Nick Gammon Australia Forum Administrator Thu 26 Apr 2007 03:46 AM #13

And now we have the cosine we can calculate the tangent:


function tangent (x, precision)
  return (sine (x, precision) / cosine (x, precision))
end -- tangent

print (tangent (1, 35))

Output

1.5574077246549022305069748074583601730872507723815200383839466056988613971517272895550999652022429935

Nick Gammon Australia Forum Administrator Thu 26 Apr 2007 04:38 AM #14

And now to raise e to a power:


bc.digits (100)

function exp (x)

  local x = bc.number (x)  -- argument converted to bignum

  -- if x > 1 halve it and square the result
  if bc.compare (x, 1) == 1 then
    return exp (x / 2) ^ 2
  end -- if x > 1

  local e = bc.number (1)  -- result

  local n = bc.number (1)  -- denominator: i!
  local t = bc.number (1)  -- numerator:   x^i

  for i = 1, 200 do
    local E = e
    n = n * i  -- n is i!  (factorial of i)
    t = t * x  -- t is x^i (x to the power i)
    e = e + (t / n)  -- new value

    -- break if no change
    if e == E then
      break
    end -- if no change
  end -- for

  return e

end -- function exp

print (exp (5))

Output

148.4131591025766034211155800405522796234876675938789890467528451109120648209585760796884094598990189279

Nick Gammon Australia Forum Administrator Thu 26 Apr 2007 05:31 AM #15

And now the inverse operation, taking the natural logarithm:


bc.digits (100)

function ln (x)
  local x = bc.number (x)  -- argument converted to bignum

  -- if x <= 0.5 or >= 2 take the square root and double the result
  if bc.compare (x, 0.5) ~= 1 or bc.compare (x, 2) == 1 then
    return ln (bc.sqrt (x)) * 2
  end -- if x <= 0.5 or >= 2

  local a = bc.number ((x - 1) / (x + 1))
  local e = bc.number (0)
  local t = bc.number (a)

  for i = 1, 200, 2 do
   local E = e
   
   e = e + (t / i)
   t = t * a * a
    
   -- break if no change
    if e == E then
      break
    end -- if no change

  end -- for 

  return e * 2
end -- function ln

print (ln (100))

Output

4.6051701859880913680359829093687284152022029772575459520666558019351452193547049604719944101791965216

David Haley USA Thu 26 Apr 2007 06:45 AM #16

I agree that this is really nifty stuff, and some very clever coding, but I'm curious: what do people use bignum math for in MUSHclient?

Nick Gammon Australia Forum Administrator Thu 26 Apr 2007 07:37 AM #17

They probably don't, I must admit, but it is fun stuff to play with.

I think someone famous said "after the third decimal place, no-one gives a damn".

David Haley USA Thu 26 Apr 2007 07:54 AM #18

Ah, ok, I was just checking. :)

Quote:
I think someone famous said "after the third decimal place, no-one gives a damn".

Hrmph, if only I could say the same. :-) In some of the code I am writing, I unfortunately care very much about some numbers after the third decimal... if I can make the inner loop a few tens of microseconds faster, that could mean speeding up the entire program considerably...

Nick Gammon Australia Forum Administrator Thu 26 Apr 2007 08:08 AM #19

Yes, well, I bought an interesting book recently (Coincidences, Chaos and All That Math Jazz), and it suggested this experiment:

In your favourite spreadsheet, put the number 0.5 in cell A1, and then make A2 be "=A1^2-2", and then "fill down" the same calculation over a few dozen rows. My first 20 look like this:


0.5
-1.75
1.0625
-0.87109375
-1.241195679
-0.459433287
-1.788921055
1.20023854
-0.559427448
-1.687040931
0.846107103
-1.284102771
-0.351080073
-1.876742782
1.52216347
0.316981629
-1.899522647
1.608186286
0.58626313
-1.656295543

Now do the same in column B, but start with a very slightly different number: 0.50001.

This time the results are this:


0.50001
-1.74999
1.062465
-0.871168124
-1.241066099
-0.459754937
-1.788625398
1.199180813
-0.561965379
-1.684194913
0.836512505
-1.300246828
-0.309358186
-1.904297513
1.626349018
0.645011127
-1.583960646
0.508931328
-1.740988903
1.031042361

Notice that by the 19th row, the sign isn't even the same! (compare 0.58626313 to -1.740988903).

This is an interesting example of how retaining the decimal places is absolutely crucial to getting a meaningful answer.

Nick Gammon Australia Forum Administrator Mon 16 Aug 2010 01:41 AM #20

Another interesting algorithm - GCD (Greatest common divisor)


-- See: http://en.wikipedia.org/wiki/Greatest_common_divisor
-- and  http://en.wikipedia.org/wiki/Euclidean_algorithm

function gcd (a, b)
local old_digits = bc.digits (0)  -- we are using integers
local temp

  a = bc.number (a)
  b = bc.number (b)

  while not bc.iszero (b) do  
     temp = b
     b = bc.mod (a, b)
     a = temp
  end -- while
  
  bc.digits (old_digits)  -- put old decimal places back
  return a
end -- gcd

print (gcd (8, 12))   --> 4
print (gcd (42, 56))  --> 14
print (gcd (206, 40)) --> 2
print (gcd (12070, 41820)) --> 170
print (gcd (10362, 12397)) --> 11

In mathematics, the greatest common divisor (GCD) of two or more non-zero integers, is the largest positive integer that divides the numbers without a remainder. For example, the GCD of 8 and 12 is 4.

Nick Gammon Australia Forum Administrator Mon 16 Aug 2010 01:47 AM #21

Given the above, we can also work out the LCM (least common multiple):


-- See: http://en.wikipedia.org/wiki/Least_common_multiple

function lcm (a, b)
  a = bc.number (a)
  b = bc.number (b)
  return (a * b) / gcd (a, b)
end -- lcm

print (lcm (21, 6))  --> 42
print (lcm (4, 6))   --> 12

In arithmetic and number theory, the least common multiple or (LCM) of two rational numbers a and b is the smallest positive rational number that is an integer multiple of both a and b.