Calculating with big numbers

[Arnold]

Hi Charles,

this is really fantastic. (if I did this correctly).
I used your functions of BigNumbers2.o2bas for Karatsuba's method:

Code: [Select]

....
....

namespace

  function MatchLength(string *n1,*n2)
  ====================================
  int l1=len n1
  int l2=len n2
  if l1>l2 then n2=stringbytes(l1-l2,48)+n2 : return
  if l2>l1 then n1=stringbytes(l2-l1,48)+n1 : return
  end function


function multiply(string a,b) as string
'http://www.rosettacode.org/wiki/Arbitrary-precision_integers_(included)#Pure_C
'https://en.wikipedia.org/wiki/Multiplication_algorithm#Karatsuba_multiplication
'https://en.wikipedia.org/wiki/Karatsuba_algorithm#Pseudocode

   string high1,high2,low1,low2
   string z0,z1,z2, t1a, t1b, r1,r2,r3, result
   int le,le1,le2,t

   le1=len(a) : le2=len(b)
   if le1 < 64 or le2 < 64 then
     return bignum::mul(a,b)
   end if

   'calculates the size of the numbers
   le=le1
   if le2>le1 then le=le2
   int t=mod(le,2)
   int m=le+t
   int mby2 = m\2 

   'split the digit sequences about the middle
   'high1, low1 = split_at(num1, m2)
   'high2, low2 = split_at(num2, m2)

   high1=left(a,le1-mby2) : low1=mid(a,-mby2)
   if len(high1)=0 then high1="0"   
   high2=left(b,le2-mby2) : low2=mid(b,-mby2)
   if len(high2)=0 then high2="0"

    ' 3 calls made to numbers approximately half the size
    ' z0 = karatsuba(low1,low2)
    ' z1 = karatsuba((low1+high1),(low2+high2))
    ' z2 = karatsuba(high1,high2)

    z2 = multiply(high1,high2)
    z0 = multiply(low1,low2)

      MatchLength(high1,low1)
      t1a=bignum::add(high1,low1)
      MatchLength(high2,low2)
      t1b=bignum::add(high2,low2)
    z1 = multiply(t1a,t1b)

    ' return (z2*10^(2*m2))+((z1-z2-z0)*10^(m2))+(z0)           

      t1 = bignum::sub(z1,z2)
      t2 = bignum::sub(t1,z0)

      r1=z2 & stringbytes(m,48)     
      r2=t2 & stringbytes(mby2,48)
      MatchLength(r1,r2)
      r3 = bignum::add(r1, r2)
      MatchLength(r3,z0)     
    result = bignum::add(r3, z0)

   return bignum::zerotrim result
end function


'TEST

'Use Windows functions
! GetTickCount lib "kernel32.dll" () as sys

function str_exp(int b, int n) as string
  string result, a

  result="1"
  a=str(b)

  if n=0 then return result
 
  while n > 0
    if n and 1 then 'mod(n,2)=1 then
      result = multiply(result, a)     
    end if
    n = n >> 1 '\= 2
    if n=0 then exit while   
    a = multiply(a, a)
  wend

  return result
end function


for z=3 to 4'4
for x=2 to 3
for y=1 to 2
   z1=str_exp(x,y)
   printl z1
   z2=str_exp(z,z1)
   printl z2
   
   t1=GetTickCount()
   printl "Calculatating 5^" z "^" x "^" y " which is 5^" z2 & cr

   z3=str_exp(5,z2)

   t2=(GetTickCount()-t1)/1000
   printl "Length for result of 5^" z "^" x "^" y " is " len(z3)
   printl "First 20 digits: " left(z3,20)
   printl "Last 20 digits:  " mid(z3,-20)
   printl "This took " str(t2,3) " seconds"
   printl
next y
next x
next z

printl "Enter ... " : waitkey

/*
GCC compiled:
Length of 5^4^3^2 is 183231
First 20 digits:  62060698786608744707
Last 20 digits:   92256259918212890625
This took about 300 seconds on my 12 years old notebook
*/

I used this to solve the Rosettacode challenge 5**262144. The execution time is 4 to 5 times faster than the compiled gcc code (which applied only school-grade multiplication).

Most probably the code can still be optimized. But the result is already really cool.

Roland

[Charles Pegge]

Here is the base100million multiplier. It is slightly more complex than base10000 due to the eax register overflow, and I lapse into a short section of assembly code. It could be adapted quite easily for use with BigNum strings, and provide a 50..60X speed improvement over base10 multiplication.

Code: [Select]


uses console

function bhmMul(dword *a,*b,*c,q)
=================================
dword d,i,j,k,ca,cy,oo,oq,m100=1E8
dword aa at @a
dword bb at @b
dword cc at @c
for i=1 to q*2
  aa=0 : @aa+=4 'clearing accum
next
oo=q*4-4 'shifting offset
oq=oo    'fixed offset
for i=1 to q
  @cc=@c+oo     'multiplier
  @bb=@b+oq     'value
  @aa=@a+q*4+oo 'accum
  cyy=0         'clear accum carry
  cy=0          'clear multiply carry
  for j=1 to q
    'k=bb*cc+cy : cy=0
    'if k>=m100 then
    '  cy=k\m100
    '  k=edx 'remainder after cpu division
    'end if
    addr edi,bb
    mov eax,[edi]
    addr edi,cc
    mul [edi]
    mov ecx,m100
    div ecx
    mov d,eax   'output carry
    add edx,cy  'input carry
    mov eax,edx 'remainder == modulus
    cmp eax,ecx
    (
      jb exit
      inc d
      sub eax,ecx
    )
    mov k,eax
    cy=d
    k=k+aa+cyy : cyy=0 'accum add with carry
    if k>=m100 then k-=m100 : cyy=1
    aa=k
    @aa-=4
    @bb-=4
  next
  aa=cyy+cy 'final carries
  oo-=4 'shift left for next line mul
next
end function


function bhmCompare(dword *b,*c,i) as int
=========================================
dword bb at @b
dword cc at @c
while i--
  if bb>cc then return 1
  if bb<cc then return -1
  @bb+=4
  @cc+=4
wend
end function



indexbase 0
int nm[32],nt[32],nr[32] 'result,target,root

/*
int nr[32],n1[32],n2[32]
'n1=14140000
'n2=14140000
'n1=22300002
'n2=22300002
'n1=87654321
'n2=2000000
bhmmul(nr,n1,n2,16)
for i=0 to 31
  print i tab mid("0000000"+str(nr[i]),-8) cr
next
print cr
*/


'square roots
=============

'square root of 2
'using binary step approximator
indexbase 0
for i=0 to 31 : nr[i]=0 : next 'clear nr
nt[0]=02000000   'target value representing 2
int b,c,i
int n at @nr 'strider for square root digits
b=0x8000000 '67,108,864 shifting bit setter
i=16
while i
  n or= b 'set bit
  bhmMul nm,nr,nr,16
  c=bhmCompare(nt,nm,16)
  if c<0 then 'too large
    n xor= b 'clear bit
  end if
  shr b 'next bit down
  if b=0 then
    i--
    b=0x4000000
    @n+=4 'next 100M down
  end if
wend

'display digits in groups of eight
print "square root of 2" cr
for i=0 to 15
  print i tab mid("0000000"+str(nr[i]),-8) cr
next
print cr
'  1.4142135623730950488016887242096980785696718753769480731766797379907324784621


'square root of 5 --> phi
'using binary step approximator
indexbase 0
for i=0 to 31 : nr[i]=0 : next 'clear nr
nt[0]=05000000 'target value representing 5
int b,c,i
int n at @nr 'strider for square root digits
b=0x4000000 '67,108,864 shifting bit setter
i=16
while i
  n or= b 'set bit
  bhmMul nm,nr,nr,16
  c=bhmCompare(nt,nm,32)
  if c<0 then 'too large
    n xor= b 'clear bit
  end if
  shr b 'next bit down
  if b=0 then
    i--
    b=0x4000000
    @n+=4 'next 10k down
  end if
wend

'display digits in groups of eight
print "square root of 5" cr
for i=0 to 15
  print i tab mid("0000000"+str(nr[i]),-8) cr
next
'2.2360679774997896964091736687312762354406183596115257242708972454105209256378



'convert to little phi
nr[0]-=10000000 'subtract 1
'shift right 'divide by 2
addr esi,nr
mov ecx,16
clc
pushf
(
 xor edx,edx
 popf
 (
  jnc exit
  mov edx,50000000
 )
 shr dword [esi] 'shift right carry
 pushf
 add [esi],edx
 add esi,4
 dec ecx
 jg repeat
)
popf


'display digits in groups of eight
print "phi" cr
for i=0 to 15
  print i tab mid("0000000"+str(nr[i]),-8) cr
next
'61803398874989484820458683436563811772030917980576286213544862270526046281890

'print 0x4000000 cr '67,108,864
'print hex 67108864


waitkey
 


[Arnold]

Hi Charles,

your bhmMul routine is of course very fast. I tried this again with the Rosetta challenge:

Code: [Select]

...
,,,

  function ZeroTrim(string a) as string
  =====================================
  byte aa at strptr a
  int la=len a
  int c=1
  while la
    if aa>48 then return mid a,c
    la-- : c++ : @aa++
  wend
  end function

sub ZeroIntTrim(int a[], int *count)
  int aa at @a
  int *cnt = @count
                                     
  int c=1
  while aa[c]=0
     c+=1
  wend
  c-=1
  if c>0 then
    for x=1 to count
       diff=x+c           
       aa[x]=aa[diff]
    next
  end if
  cnt = count-c 
end sub

function big2string(int n[],count) as string
   string s
   
   for x=1 to count
      s += mid("0000000"+str(n[x]),-8)
   next
   return s     
end function


' Test
'Use Windows function
! GetTickCount lib "kernel32.dll" () as sys

printl "Calculting 5^4^3^2 which is 5^262144"
t1=GetTickCount()

int result[25000], n[25000] 'result,target,root
int cn=1                    'count

n={5}
for y = 1 to 18
   bhmMul(result, n,n,cn)
   cn=cn+cn
   for x=1 to cn
     n[x]=result[x]
   next 
   zeroIntTrim(n,cn)
next

string s=zeroTrim big2string(n,cn)

t2=(GetTickCount()-t1)/1000
printl "Length for result of 5^4^3^2 is " len(s)
printl "First 20 digits: " left(s,20)
printl "Last 20 digits:  " mid(s,-20)
printl "This took " str(t2,3) " seconds"
printl

printl "Enter ... " : waitkey

'Length of 5^4^3^2 is 183231
'First 20 digits:  62060698786608744707
'Last 20 digits:   92256259918212890625

I assume my code is not very clever, but the program beats the execution time of the provided C-code easily by a factor of 10 to 14.

Nevertheless I like your routines of BigNumbers2.o2bas very much. All basic calculus are provided, I do not need fixed precision and it gives me a feeling of school-grade calculation. I will use them to evaluate some formulas walking in the footsteps of some mathematical heros, without the need of using libraries like GMP or LibTomMath.

Unfortunately I was not very good in mathematics during my school time. But I suspect my teachers were not very qualified neither. And many findings have been discoverd only in recent years.

I wondered about the use of counting with big numbers. Searching in Internet there are surprisingly several application areas, not only dealing with public-key cryptograpy. And I like the feeling that I can be more precise than my windows calculator.

Roland

[JRS]

Maybe Charles work will lead to faster bitcoin mining. 

[Arnold]

This is a small app to calculate Pi to several places, derived from a pascal program which I found at RosettaCode.org. It applies a dynamic array to store the digits.
I learned that there is a man in Japan who could memorize and recite 100000 digits of Pi in about 16.5 hours. But I will not need so many digits.

Roland

Code: OxygenBasic

1. 'http://rosettacode.org/wiki/Pi#Pascal
2. 'http://www.pi314.net/eng/goutte.php
3.  
4. /* A spigot algorithm for the digits of pi,
5.    Stanley Rabinowitz and Stan Wagon
6. */
7.  
8. include "$/inc/console.inc"
9.  
10. function Pi_Spigot(int n) as string
11.    int  count, j, k, q, nines, predigit
12.    int x,i
13.    int idx
14.  
15.    count = 10*n \ 3
16.  
17.    'Create dynamic array
18.   sys aa=getmemory count * sizeof(uint)
19.    uint a at aa
20.  
21.    string s = stringbytes(n,0)  'nuls  
22.   byte ss at strptr s
23.  
24.    for j = 1 to count
25.       a[j] = 2                  'Start with 2s
26.   next j
27.    nines = 0
28.    predigit = 0                 'First predigit is a 0
29.  
30.    for j = 1 to n
31.       if mod(j,1000) = 0 then print j & cr  
32.    
33.       'Only calculate as far as needed
34.      '+16 for security digits ~5 decimals
35.      i = (n-j)*10 \ 3 +16
36.       if i > count then i = count
37.       q = 0
38.       do                        'Work backwards
39.         x  = 10*a[i] + q*i
40.          q = x \ (2*i - 1)
41.          a[i] = x - q*(2*i - 1) 'mod(x,2*i - 1)
42.         i--
43.          if i<= 0 then exit do
44.       end do
45.  
46.       a[1] = mod(q,10)
47.       q = q \ 10    
48.       if q = 9 then
49.         nines += 1
50.       else
51.         if q = 10 then
52.           idx++ : ss[idx]=predigit+1
53.           for k = 1 to nines : idx++ : ss[idx]=0 : next k   'zeros
54.          predigit = 0
55.           nines = 0
56.         else      
57.           if j > 1 then idx++ : ss[idx]=predigit
58.           predigit = q
59.           if nines <> 0 then
60.             for k = 1 to nines : idx++ : ss[idx]=9 : next k 'nine
61.            nines = 0
62.           end if
63.         end if
64.       end if
65.    next j
66.    idx++ : ss[idx]=predigit
67.  
68.    'convert to ascii
69.   for x=1 to n : ss[x]+=48 : next x
70.  
71.    'release dynamic array
72.   freememory aa
73.  
74.    return s
75. end function
76.  
77. cls
78. print "This will calculate pi and save the result as a text file"
79. printl "How many digits? "
80. int num=val(rtrim ltrim input())
81.  
82. string s=Pi_Spigot(num)
83. print s
84.  
85. putfile("pi.txt",s)  
86.  
87. printl
88. printl "Saved as pi.txt"
89. printl "Enter ... " : waitkey
90.  

.

[Charles Pegge]

Hi Roland,

Many thanks for demonstrating this Pi algorithm. I'm studying it closely, as with the Keratsuba multiplication.

I can get your Karatsuba demo to run twice the speed by using base100 and digit lookup tables inside the bignum::mul function, along with some assembler on the inner loop.

Karatsuba reminds me of merge-sort which breaks the deadly exponential computation of Bubble-sort

[Arnold]

Hi Charles,

you will make my old machine run really fast. But this is also a reason why I still prefer my old notebook to my new one, I can see the gain of speed at once visually. And therefore I know that Oxygenbasic in fact does a very good job.

There is also the BBP formula which is a special spigot algorithm. Details and demos can be found at David H.Bailey's website:
http://www.davidhbailey.com/
http://www.experimentalmath.info/bbp-codes/

I assume the use of the BBP algorithm is more effective in some cases and will save memory and execution time. But unfortunately I have not yet found out how this will work correctly with Oxygenbasic.

Roland

[Arnold]

This is a nice little routine to compute the mathematical constant Euler's number e to a given amount of digits. As e is of great importantance in mathematics it could be useful sometimes to have access to it. Although nowadays one can paste almost everything from Internet.

Code: OxygenBasic

1. ' http://faculty.cs.tamu.edu/klappi/csce411-s15/csce411-set1.pdf
2. ' compute N digits of Euler's number e, using Spigot algorithm
3.  
4. include "$/inc/console.inc"
5.  
6. declare function e_spigot(int N) as string
7.  
8. SetConsoleTitle "Euler's number e"
9. string s = e_Spigot(150)
10. printl s
11.  
12. putfile("e.txt",s)
13. printl
14. printl "Saved as e.txt"
15.  
16. printl "Enter ... " : waitkey
17. ------------------------------
18.  
19. function e_Spigot(int N) as string
20.    int q, idx
21.  
22.    'int A[N]
23.   sys aa=getmemory N * sizeof(int)
24.    int A at aa
25.  
26.    string s = stringbytes(N,0)      'nuls  
27.   byte ss at strptr s
28.    idx++ : ss[idx]=2
29.    
30.    for j = 1 to N
31.       A[j] = 1                      'set all digits of nonintegral part to 1.
32.   next j                          
33.  
34.    for i = 1 to N-2
35.       if mod(i,1000) = 0 then printl i
36.       q = 0
37.       j=N
38.       while j > 0
39.          A[j] = 10 * A[j] + q
40.          q = A[j] \ (j + 1)         'amount that needs to be carried
41. '         A[j] = mod(A[j],(j + 1))   'keep only remainder
42.         A[j] = edx
43.          j--
44.       wend
45.       idx++ : ss[idx] = q
46.    next i
47.    
48.    'convert to ascii
49.   for x=1 to n : ss[x]+=48 : next x
50.  
51.    'release dynamic array
52.   freememory aa
53.  
54.    return s  
55. end function
56.  

[Charles Pegge]

Thanks Roland,

Very similar to the Pi algorithm

I am also looking for a spigot algorithm to approximate square roots. The series for roots cannot be formulated in the same way as Pi and Euler, since it involves recursive division. It cannot be split into discrete digits. But I think it can be done by refining the binary approximation procedure, by eliminating repeated calculations at the upper end, as the digits resolve from left to right.

[Charles Pegge]

Hi Roland,

While studying the left-to-right addition in the Pi algorithm, I found that it gave incorrect values at the end of 200 and 400 digit series. This was due to digits remaining in the nines buffer.

The end of the spigot function should look something like:

Code: [Select]

...
  next j
  idx++ : ss[idx]=predigit 'final predigit
  if nines then 'remaining nines
    for k = 1 to nines
      idx++ : ss[idx]=9
    next k
    nines = 0
  end if
  'to ascii
  for j=1 to len s
    ss[j]+=48
  next
  return s
end function

.

[Arnold]

Hi Charles,

your code of PiAlgorithm2.o2bas looks so clear, I wished I could do this in some way too. And you are right, the last digit produced with my PiSpigot.o2bas is different in some cases. (200, 400 is wrong, 201 and 401 would work). I would never had found this difference. I testet for 100000 digits and some other numbers and compared with http://www.pibel.de/ but I could not find a difference with these counts, so I assumed the code would be ok. But PiAlgorithm2.o2bas should have fixed this remaining problem now.

Roland

Edit: the attached text file includes in fact 100001 digits, I always forget the preceeding 3.

.