Register forum user name Search FAQ

Gammon Forum

Notice: Any messages purporting to come from this site telling you that your password has expired, or that you need to verify your details, confirm your email, resolve issues, making threats, or asking for money, are spam. We do not email users with any such messages. If you have lost your password you can obtain a new one by using the password reset link.
 Entire forum ➜ MUSHclient ➜ Lua ➜ Doing arbitrary precision (bignum) arithmetic in Lua
Doing arbitrary precision (bignum) arithmetic in Lua

It is now over 60 days since the last post. This thread is closed.     Refresh page

Pages: 1 2  

Posted by Nick Gammon   Australia  (23,120 posts)  Bio   Forum Administrator
Date Mon 02 Apr 2007 06:16 AM (UTC)

Amended on Thu 05 Apr 2007 11:21 AM (UTC) by Nick Gammon

Message
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

www.gammon.com.au, www.mushclient.com
Top
Posted by Nick Gammon   Australia  (23,120 posts)  Bio   Forum Administrator
Date Reply #1 on Mon 02 Apr 2007 10:52 PM (UTC)
Message
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

www.gammon.com.au, www.mushclient.com
Top
Posted by Nick Gammon   Australia  (23,120 posts)  Bio   Forum Administrator
Date Reply #2 on Wed 18 Apr 2007 12:45 AM (UTC)
Message
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

www.gammon.com.au, www.mushclient.com
Top
Posted by Nick Gammon   Australia  (23,120 posts)  Bio   Forum Administrator
Date Reply #3 on Wed 18 Apr 2007 05:31 AM (UTC)

Amended on Wed 18 Apr 2007 04:48 PM (UTC) by Nick Gammon

Message
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

www.gammon.com.au, www.mushclient.com
Top
Posted by Nick Gammon   Australia  (23,120 posts)  Bio   Forum Administrator
Date Reply #4 on Wed 18 Apr 2007 06:21 AM (UTC)
Message
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

www.gammon.com.au, www.mushclient.com
Top
Posted by Nick Gammon   Australia  (23,120 posts)  Bio   Forum Administrator
Date Reply #5 on Wed 18 Apr 2007 06:53 AM (UTC)
Message
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

www.gammon.com.au, www.mushclient.com
Top
Posted by Nick Gammon   Australia  (23,120 posts)  Bio   Forum Administrator
Date Reply #6 on Wed 18 Apr 2007 07:37 AM (UTC)

Amended on Wed 18 Apr 2007 07:39 AM (UTC) by Nick Gammon

Message
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.
- Nick Gammon

www.gammon.com.au, www.mushclient.com
Top
Posted by David Haley   USA  (3,881 posts)  Bio
Date Reply #7 on Wed 18 Apr 2007 07:39 AM (UTC)
Message
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




:-)
David Haley aka Ksilyan
Head Programmer,
Legends of the Darkstone

http://david.the-haleys.org
Top
Posted by Nick Gammon   Australia  (23,120 posts)  Bio   Forum Administrator
Date Reply #8 on Wed 18 Apr 2007 07:41 AM (UTC)
Message
Yes, but my numbers include the digits 4 through to 9. :P
- Nick Gammon

www.gammon.com.au, www.mushclient.com
Top
Posted by David Haley   USA  (3,881 posts)  Bio
Date Reply #9 on Wed 18 Apr 2007 07:48 AM (UTC)
Message
Foiled again!

[in a dark, ominous voice:]
I'll be back!
David Haley aka Ksilyan
Head Programmer,
Legends of the Darkstone

http://david.the-haleys.org
Top
Posted by Nick Gammon   Australia  (23,120 posts)  Bio   Forum Administrator
Date Reply #10 on Wed 25 Apr 2007 07:18 AM (UTC)

Amended on Wed 25 Apr 2007 07:32 AM (UTC) by Nick Gammon

Message
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

www.gammon.com.au, www.mushclient.com
Top
Posted by Nick Gammon   Australia  (23,120 posts)  Bio   Forum Administrator
Date Reply #11 on Wed 25 Apr 2007 07:53 AM (UTC)

Amended on Thu 26 Apr 2007 03:38 AM (UTC) by Nick Gammon

Message
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

www.gammon.com.au, www.mushclient.com
Top
Posted by Nick Gammon   Australia  (23,120 posts)  Bio   Forum Administrator
Date Reply #12 on Wed 25 Apr 2007 07:58 AM (UTC)

Amended on Thu 26 Apr 2007 03:39 AM (UTC) by Nick Gammon

Message
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

www.gammon.com.au, www.mushclient.com
Top
Posted by Nick Gammon   Australia  (23,120 posts)  Bio   Forum Administrator
Date Reply #13 on Thu 26 Apr 2007 03:46 AM (UTC)
Message
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

www.gammon.com.au, www.mushclient.com
Top
Posted by Nick Gammon   Australia  (23,120 posts)  Bio   Forum Administrator
Date Reply #14 on Thu 26 Apr 2007 04:38 AM (UTC)

Amended on Thu 26 Apr 2007 05:36 AM (UTC) by Nick Gammon

Message
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

www.gammon.com.au, www.mushclient.com
Top

The dates and times for posts above are shown in Universal Co-ordinated Time (UTC).

To show them in your local time you can join the forum, and then set the 'time correction' field in your profile to the number of hours difference between your location and UTC time.

57,412 views.

This is page 1, subject is 2 pages long: 1 2  [Next page]

It is now over 60 days since the last post. This thread is closed.     Refresh page

Go to topic:           Search the forum

[Go to top] top

Information and images on this site are licensed under the Creative Commons Attribution 3.0 Australia License unless stated otherwise.

Home