'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
*/
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
*/
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
'
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
...
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
...
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
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
printl "Test bignum::add"
n1="56712345646897989899898" : n2="4321078787800789457899799978151512"
printl "n1 = " n1
printl "n2 = " n2
printl "n1 + n2 = " bignum::add(n1,n2)
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
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
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)
....
....
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
*/
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
...
,,,
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
...
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