Author Topic: Calculating with big numbers  (Read 597 times)

0 Members and 1 Guest are viewing this topic.

Arnold

  • Sr. Member
  • ****
  • Posts: 469
Re: Calculating with big numbers
« Reply #15 on: October 24, 2017, 09:36:02 AM »
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

  • Author
  • *****
  • Posts: 3290
    • Oxygen Basic
Re: Calculating with big numbers
« Reply #16 on: October 25, 2017, 10:12:20 PM »
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

  • Sr. Member
  • ****
  • Posts: 469
Re: Calculating with big numbers
« Reply #17 on: October 30, 2017, 03:37:54 AM »
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

John

  • Hero Member
  • *****
  • Posts: 2813
Re: Calculating with big numbers
« Reply #18 on: October 30, 2017, 05:57:42 PM »
Maybe Charles work will lead to faster bitcoin mining.  :)
« Last Edit: October 30, 2017, 08:49:32 PM by John »

Arnold

  • Sr. Member
  • ****
  • Posts: 469
Re: Calculating with big numbers
« Reply #19 on: November 03, 2017, 02:12:10 AM »
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

  • Author
  • *****
  • Posts: 3290
    • Oxygen Basic
Re: Calculating with big numbers
« Reply #20 on: November 03, 2017, 11:56:49 PM »
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

  • Sr. Member
  • ****
  • Posts: 469
Re: Calculating with big numbers
« Reply #21 on: November 04, 2017, 03:35:42 AM »
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

  • Sr. Member
  • ****
  • Posts: 469
Re: Calculating with big numbers
« Reply #22 on: November 04, 2017, 07:31:44 AM »
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

  • Author
  • *****
  • Posts: 3290
    • Oxygen Basic
Re: Calculating with big numbers
« Reply #23 on: November 07, 2017, 09:35:28 PM »
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

  • Author
  • *****
  • Posts: 3290
    • Oxygen Basic
Re: Calculating with big numbers
« Reply #24 on: November 09, 2017, 06:03:01 PM »
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
« Last Edit: November 10, 2017, 12:29:49 AM by Charles Pegge »

Arnold

  • Sr. Member
  • ****
  • Posts: 469
Re: Calculating with big numbers
« Reply #25 on: November 10, 2017, 08:02:40 AM »
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.
« Last Edit: November 10, 2017, 08:09:53 AM by Arnold »