### Author Topic: Calculating with big numbers  (Read 989 times)

0 Members and 1 Guest are viewing this topic.

#### Arnold

• Hero Member
• Posts: 510
##### 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 functionfunction 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 resultend function'TEST'Use Windows functions! GetTickCount lib "kernel32.dll" () as sysfunction 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 resultend functionfor z=3 to 4'4for x=2 to 3for 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"   printlnext ynext xnext zprintl "Enter ... " : waitkey/*GCC compiled:Length of 5^4^3^2 is 183231First 20 digits:  62060698786608744707Last 20 digits:   92256259918212890625This 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: 3365
##### 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 consolefunction bhmMul(dword *a,*b,*c,q)=================================dword d,i,j,k,ca,cy,oo,oq,m100=1E8dword aa at @adword bb at @bdword cc at @cfor i=1 to q*2  aa=0 : @aa+=4 'clearing accumnextoo=q*4-4 'shifting offsetoq=oo    'fixed offsetfor 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 mulnextend functionfunction bhmCompare(dword *b,*c,i) as int=========================================dword bb at @bdword cc at @cwhile i--  if bb>cc then return 1  if bb<cc then return -1  @bb+=4  @cc+=4wendend functionindexbase 0int 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=2000000bhmmul(nr,n1,n2,16)for i=0 to 31  print i tab mid("0000000"+str(nr[i]),-8) crnextprint cr*/'square roots============='square root of 2'using binary step approximatorindexbase 0for i=0 to 31 : nr[i]=0 : next 'clear nrnt[0]=02000000   'target value representing 2int b,c,iint n at @nr 'strider for square root digitsb=0x8000000 '67,108,864 shifting bit setteri=16while 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 ifwend'display digits in groups of eightprint "square root of 2" crfor i=0 to 15  print i tab mid("0000000"+str(nr[i]),-8) crnextprint cr'  1.4142135623730950488016887242096980785696718753769480731766797379907324784621'square root of 5 --> phi'using binary step approximatorindexbase 0for i=0 to 31 : nr[i]=0 : next 'clear nrnt[0]=05000000 'target value representing 5int b,c,iint n at @nr 'strider for square root digitsb=0x4000000 '67,108,864 shifting bit setteri=16while 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 ifwend'display digits in groups of eightprint "square root of 5" crfor i=0 to 15  print i tab mid("0000000"+str(nr[i]),-8) crnext'2.2360679774997896964091736687312762354406183596115257242708972454105209256378'convert to little phinr[0]-=10000000 'subtract 1'shift right 'divide by 2addr esi,nrmov ecx,16clcpushf( 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 eightprint "phi" crfor i=0 to 15  print i tab mid("0000000"+str(nr[i]),-8) crnext'61803398874989484820458683436563811772030917980576286213544862270526046281890'print 0x4000000 cr '67,108,864'print hex 67108864waitkey `

#### Arnold

• Hero Member
• Posts: 510
##### 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 functionsub 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 subfunction 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 sysprintl "Calculting 5^4^3^2 which is 5^262144"t1=GetTickCount()int result[25000], n[25000] 'result,target,rootint cn=1                    'countn={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)nextstring s=zeroTrim big2string(n,cn)t2=(GetTickCount()-t1)/1000printl "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"printlprintl "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: 2843
##### 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

• Hero Member
• Posts: 510
##### 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
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: 3365
##### 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

• Hero Member
• Posts: 510
##### 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

• Hero Member
• Posts: 510
##### 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: 3365
##### 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: 3365
##### 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 send function`
« Last Edit: November 10, 2017, 12:29:49 AM by Charles Pegge »

#### Arnold

• Hero Member
• Posts: 510
##### 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 »