Calculating with big numbers

[Arnold]

Hi Charles,

this is a programming task of Rosettacode
http://www.rosettacode.org/wiki/Arbitrary-precision_integers_(included)

and I am not very happy with my solution.
Calculating up to 1000 digits there is no problem but then the execution time grows very strongly. To calculate 5^19683 will take as much time on my old notebook as for gcc compiled 5^262144. I could not calculate this power with the Oxygenbasic code and I gave up after 2 hours.
Can I improve this code in some way? I know I use too many garbage collections, but I do not know how I can get rid of them.

Roland

Code: [Select]

'http://www.rosettacode.org/wiki/Arbitrary-precision_integers_(included)
'http://www.rosettacode.org/wiki/Arbitrary-precision_integers_(included)#Pure_C

$ filename "Arbitr_Prec.exe"
'include "$/inc/RTL32.inc"

include "$/inc/console.inc"

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

' return = a * b
' Handling of negatives and zeros is not here, since not needed.
function str_mult(string First, string Second) as string

    string a=First, b=Second, r ' result
   
    int ax = 0, bx = 0, rx = 0, al, bl
    int carry, n
    int an, bn, rn, t 'auxiliary

    al = len(a) : bl = len(b)   
    r = space (al+bl)
       
    ' grade-school method of multiplication
    for ax = al to 1 step -1
        carry = 0
        for bx = bl to 1 step -1
           rx = ax + bx
           an=mid(a,ax,1) : bn=mid(b,bx,1) : rn=mid(r,rx,1)                   
           n=an * bn + rn + carry
           mid r,rx, str(mod(n, 10))           
           carry = n \ 10             
        next bx
        rx--       
        int t=mid(r,rx,1) + carry
        mid r,rx,t
    next ax

    ' omit possible first place - zero
    if left(r,1) = "0" then return mid(r,2)   
    return r
end function


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

  res="1"
  a=str(b)

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

  return res
end function


for z=3 to 3 '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 xxx seconds
*/

output:

2
9
Calculatating 5^3^2^1 which is 5^9

Length for result of 5^3^2^1 is 7
First 20 digits: 1953125
Last 20 digits:  1953125
This took 0 seconds

4
81
Calculatating 5^3^2^2 which is 5^81

Length for result of 5^3^2^2 is 57
First 20 digits: 41359030627651383743
Last 20 digits:  82906055450439453125
This took 0.015 seconds

3
27
Calculatating 5^3^3^1 which is 5^27

Length for result of 5^3^3^1 is 19
First 20 digits: 7450580596923828125
Last 20 digits:  7450580596923828125
This took 0 seconds

9
19683
Calculatating 5^3^3^2 which is 5^19683

Length for result of 5^3^3^2 is 13758
First 20 digits: 67080354005587923286
Last 20 digits:  06702899932861328125
This took 277.9 seconds

Enter ...

[Charles Pegge]

Hi Roland,

Here is your multiplier with most of the stringware replaced with byte overlays. I have verified it against examples/math/mulhuge.o2bas.

It does your large calculation in 1.4 seconds on my notebook, but I have not looked at the rest of the code.

Code: [Select]

function str_mult(string a,b) as string

    string r ' result
   
    int ax = 0, bx = 0, rx = 0, al, bl
    int carry, n
    int an, bn, rn, t 'auxiliary

    al = len(a) : bl = len(b)   
    r = stringbytes (al+bl,48) 'FILL 0000..
    byte aa at strptr a
    byte bb at strptr b
    byte rr at strptr r
    indexbase 1
    macro midval(b,i){ b[i]-48 } 'ASCII TO DIGIT VALUE
       
    ' grade-school method of multiplication
    for ax = al to 1 step -1
        carry = 0
        for bx = bl to 1 step -1
           rx = ax + bx
           an=midval(aa,ax)
           bn=midval(bb,bx)
           rn=midval(rr,rx)                   
           n=an * bn + rn + carry
           rr[rx]=mod(n,10)+48
           carry = n \ 10             
        next bx
        rx--       
        'int t=midval(r,rx) + carry
        rr[rx]+=carry 'FINAL CARRY
    next ax

    ' omit possible first place - zero
    if rr[1]=48 then return mid(r,2)   
    return r
end function

/*
a="4225678345671239269412376809064563457830987094214727890460503258124429509352125309123467889090967803"
b="9034098905607687535809534264324877930963452357175890408065643678402109067842097453251206046098340056"

ans1="38175196118078646234796018047078508876075927016257021775208337789644882612626336290568135391711676360295755214293161291134825045621348587392814424177752560047244220841334894056033458232344500841216968"
ans2=str_mult(a,b)
printl ans1
printl ans2
waitkey
*/



[Arnold]

Hi Charles,

this is really impressive. Your modifications for str_mult reduced the execution time by almost 99% (!). I also get a correct result for 5**262144. This example demonstrates again how important it is to deal with strings efficiently in time-consuming tasks and I will study the code very carefully. I learned about a new keyword: stringbytes; I already missed this functionality sometimes.

It is also remarkable that the function str_mult runs almost as fast as the multiply function in \examples\math\MulHuge.o2bas which is realized with assembler statements. This is another proof of OxygenBasic's speed.

Thank you again for your help. It will help me to understand some more examples where this approach is applied too.

Roland

[Charles Pegge]

Hi Roland,

I'm glad you found it helpful.

There are lots of ways to make it go even faster. For instance, working with base10000 digits instead of base10 digits. In theory, this would give you a 16x boost in speed.

I've also got an algorithm for big number division, something I have always wanted to do. It is based on a series of shifting subtractions, similar to school-teacher punishment long division

Will post it later.

[Charles Pegge]

First, the numbers are converted from ascii to binary by subtracting 48 from each digit. Then a lookup table is prepared with multiples of the divisor x1 x2 x3 .. x9. Then we are ready for the main algorithm.

Comparisons are made using the table values, and the selected table index becomes the value of the current digit. The selected number is then subtracted leaving a possible remainder. We then shift 1 place to the right for the next digit, and repeat the process until the right hand end is reached.

The devil is in the detail

Code: [Select]

  function div(string sa,sb) as string
  ====================================
  string a="0"+zerotrim(sa)
  int la=len a
  string b="0"+zerotrim(sb)
  int lb=len b
  tobin a : tobin b
  int lo=la-lb+1
  string o=nuls lo
  '
  'CREATE LOOKUP TABLE 1xB 2xB ..9xB
  indexbase 1
  string tbl[9]
  int i,j,k
  sbyte *ti,*tj
  tbl[1]=b : j=1
  for i=2 to 9
    tbl[i]=b
    @ti = strptr(tbl[i])+lb-1
    @tj = strptr(tbl[j])+lb-1
    cy=0
    for k=1 to lb 'add digits
      ti=ti+tj+cy
      if ti>=10 then ti-=10 : cy=1 else cy=0
      @ti-- : @tj--
    next
    j++
  next
  '
  sbyte aa at strptr(a)
  sbyte bb at strptr(b)
  sbyte oo at strptr(o)
  sbyte d,carry
  int i,j,k,c,cm
  sys pa=@aa
  for i = 1 to lo 'setting dividend digits
    c=9 'default digit
    for k=1 to 9
      cm=CompareMid(i,a,tbl[k]) 'lookup multiple of b
      if cm<0 then c=k-1 : exit for 'set digit
    next
    if c then 'perform subtraction
      @aa=pa+i+lb-2
      @bb=strptr(tbl[c])+lb-1
      cb=lb
      carry=0
      for j=1 to lb 'subtraction
        d=aa-bb-carry
        if d<0 then d+=10 : carry=1 else carry=0
        aa=d
        @bb-- : @aa--
      next 'next subtraction digit
    end if 'subtraction
    oo=c 'set dividend digit as subraction count
    @oo++
  next  'next dividend digit
  remainder=nullTrim a
  o=nulltrim o
  if not o then o=chr 0
  toasc o
  toasc remainder
  return o
  end function
  '

supporting functions

Code: [Select]

  macro tobin(s,  bb,le) {
  ========================
  int le=len s
  sbyte bb at strptr s
  while le
    bb-=48
    @bb++ : le--
  wend
  }
  macro toasc(s,  bb,le) {
  ========================
  int le=len s
  sbyte bb at strptr s
  while le
    bb+=48
    @bb++ : le--
  wend
  }

  function compareMid(int i, string a,b) as int
  =============================================
  i-- 'convert to base 0
  int la=len(a)-i
  int lb=len b
  byte aa at i+strptr a
  byte bb at strptr b
  if lb>la
    return -1
  else
    while lb
      if aa>bb then return 1
      if aa<bb then return -1
      lb-- : @aa++ : @bb++
    wend
  end if
  end function
  '
  function nullTrim(string a) as string
  =====================================
  byte aa at strptr a
  int la=len a
  int c=1
  while la
    if aa>0 then return mid a,c
    la-- : c++ : @aa++
  wend
  end function
  '
  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


.

[Arnold]

Hey Charles,

you have already finished what I was about to do. But since I will never achieve your accuracy I will focus on your code. Will you work out the functions as a library? There are already many things possible, e.g.:

The idea to explore big numbers arose from Robbek's message:
http://www.oxygenbasic.org/forum/index.php?topic=1231.msg11994#msg11994

and this can be done easily with your functions:

Code: [Select]

...
print "fibonacci" cr
sa=0
sb=1
time1=GetTickCount()
for i=1 to 4000
  sr=add sa,sb
'  print i tab sr cr
  sa=sb
  sb=sr
next
i-=1 ' adjust i
time=(GetTickCount()-time1) / 1000
print "fibo " i " = " sr cr
printl "This took " time " Seconds" cr
...

With a modern pc this will take 0.0078 seconds or less and this is really fast if one takes into account that the code does not play with changing the bases and/or grouping the digits to integers.

I have not yet explored the functions intensively but I noticed you also added some code for negative numbers. Is it already possible to calculate 1-1000, or -100*1 or -2*-3 ? I was about to implement Karatsuba's method for multiplication, but I noticed that I will need to calculate negative interim results in some cirumstances.

Roland

[Charles Pegge]

Hi Roland,

Here is a base10000 multiplier. It works with integer arrays instead of bytes. Unlike straight binary Base10000 is very easy to manage in our decimalised world, and is more efficient than using decimal byte digits.

Code: [Select]

function btkMul(int *a,*b,*c,q)
===============================
dword i,j,k,ca,cy,oo,oq
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>=10000 then
      cy=k\10000
      k=edx 'remainder after cpu division
    end if
    k=k+aa+cyy : cyy=0
    if k>=10000 then k-=10000 : 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


[Arnold]

Hi Charles,

I tried this test:

Code: [Select]

int num1[]={9999,9999,9999,9999}
int num2[]={9999,9999,9999,9999}
int result[]={0,0,0,0,0,0,0,0}

btkMul(result[],num1[],num2[],8)

string res=""
for x=1 to 8 : res &= str(result[x]) & " " : next
print res

I get this result:
9999 9999 9999 9998 9999 9999 9999 9999

but I would expect something like this:
9999 9999 9999 9998 0000 0000 0000 0001

Do I use the function incorrectly?
In theory base 100 000 000 could also be used?

Roland

[Charles Pegge]

Hi Roland,

It should work correctly if you specify 4 instead of 8. It assumes the result buffer is twice the size.

Base100M arithmetic is a little more complicated since the 10^16 product overflows into the edx register, but I think it can be done with a little assembler to ensure unsigned mul and div are used.

[Arnold]

Thank you Charles. I thought I must use the count of the result array.
And if I apply:
int num2[]={0,9999,9999,9999}
then the result is shown also correctly:
0 9999 9999 9998 9999 0 0 1

[Arnold]

Hi Charles,

my implementation of Karatsuba's method worked somehow but it was a lot slower than the normal school-grade multiplication. I wasted too  much time with adding the minus sign with addition / subtraction. Therefore I now started to try the functions of BigNumbers2.o2bas. But there is a problem with addition:

Code: [Select]

printl "Test bignum::add"
n1="56712345646897989899898" : n2="4321078787800789457899799978151512"
printl "n1 = " n1
printl "n2 = " n2
printl "n1 + n2 = " bignum::add(n1,n2)

Output:
Test bignum::add
n1 = 56712345646897989899898
n2 = 4321078787800789457899799978151512
n1 + n2 = ¥}NZ07┴D87857501803546697968051410
Enter ...

I expect as result:
4321078787857501803546697968051410

Is there a limit of accepted digits for addition / subtraction?

Roland

[Charles Pegge]

I've used it to approximate the square roots of 2 and 5.

The square root of 5 is then used to calculate little phi: ( sqr(5) -1 ) / 2, which ties in nicely with the RobbeK link

The binary approximation is very simple with bit setting and clearing. Long division is not required, but a multiplication is required after each bit-setting to test the product against the target value.

Code: [Select]


uses console

function btkMul(int *a,*b,*c,q)
===============================
dword i,j,k,ca,cy,oo,oq
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>=10000 then
      cy=k\10000
      k=edx 'remainder after cpu division
    end if
    k=k+aa+cyy : cyy=0
    if k>=10000 then k-=10000 : 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 btkCompare(int *b,*c,i) as int
=======================================
dword bb at @b
dword cc at @c
while i
  i--
  if bb>cc then return 1
  if bb<cc then return -1
  @bb+=4
  @cc+=4
wend
end function


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

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

'square root of 2
'using binary step approximator
indexbase 0
for i=0 to 15 : nr[i]=0 : next 'clear nr
nt[0]=200   'target value representing 2
nr[0]=1410  'estimate of square root
int b,c,i
int n at @nr 'strider for square root digits
b=8192 '0x2000 'shifting bit setter
i=16
while i
  n or= b 'set bit
  btkMul nm,nr,nr,16
  c=btkCompare(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=8192
    @n+=4 'next 10k down
  end if
wend

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

'square root of 5 --> phi
'using binary step approximator
indexbase 0
for i=0 to 15 : nr[i]=0 : next 'clear nr
nt[0]=500   'target value representing 5
nr[0]=2200  'estimate of square root
int b,c,i
int n at @nr 'strider for square root digits
b=8192 '0x2000 'shifting bit setter
i=16
while i
  n or= b 'set bit
  btkMul nm,nr,nr,16
  c=btkCompare(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=8192
    @n+=4 'next 10k down
  end if
wend


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


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


waitkey
 

little phi : Fibonacci ratio
0       0618
1       0339
2       8874
3       9894
4       8482
5       0458
6       6834
7       3656
8       3811
9       7720
10      3091
11      7980
12      5762
13      8621
14      3544
15      8622


PS: is there a missing tobin call? You are crunching ascii characters instead of numbers

[Arnold]

I refer to BigNumbers2.o2bas which you provided with reply#4. I only added the code after namespace for testing addition. I tried to use tobin / toasc:

Code: [Select]

printl "Test bignum::add"
string n1,n2,result
n1="56712345646897989899898" : n2="4321078787800789457899799978151512"
printl "n1 = " n1
printl "n2 = " n2
bignum::tobin(n1)
bignum::tobin(n2)
result = bignum::add(n1,n2)
bignum::toasc(result)
printl result

but this does not help neither. Should I do something differently?
Edit: I see that you already applied tobin and toasc macro in the function add. So there must be another reason for the wrong result.

[Charles Pegge]

Yes, I see toasc and tobin are both handled internally ok.

The problem is mismatched lengths for n1 and n2

You can either pad with zeros on the left side or use a function to match the lengths for you:

Code: [Select]

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

...

printl "Test bignum::add"
string n1,n2,result
n1="56712345646897989899898"
n2="4321078787800789457899799978151512"
MatchLength(n1,n2)
printl "n1 = " n1
printl "n2 = " n2
printl "sum= " add(n1,n2)



[Arnold]

Hi Charles,

thank you again. I also just noticed that n1 and n2 expect the same length. I will use Matchlength to pad the zeros.

Roland