Project Lovecraft Problem 7

[danbaron]

Problem 7:

Make a thinBasic script that finds the quotient,

q = a \ b (integer division),

and the remainder, r ( 0 <= r < b ), such that,

a = b * q + r.

Where, a and b, are as follows.

a = 903847349374823506590013665374655721156283703993015605482624365303626229042859483018997843762967582541408949179799124552

b = 75638641208056310245

DTB

[danbaron]

Here is one way to do it.

Dan


'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
'file = Lovecraft7.tbasicc
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Uses "console"
Uses "file"
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%neg = -1
%pos = +1

%intsize = 1000

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Type genint
digit(%intsize) As Byte
nsdigits As DWord
sign As Integer
error As Byte
End Type

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Global filestring As String = "Lovecraft7.txt"
Global outfile As DWord

'##########################################################################################

Function TBMain()
Local s As String
Local tf As DWord
Local a, b, q, r, ck As genint

outfile = FILE_Open(filestring, "output")

s = "903847349374823506590013665374655721156283703993015605482624365303626229042859483018997843762967582541408949179799124552"

genintsetvalue(a, s)
genintstringrep(a, s)
writeoutputline("a =")
writeoutputline("")
writeoutputline(s)
writesep()

If geninterror(a) Then
geninterrormessage("bad a: --> script termination.")
Exit Function
End If

s = "75638641208056310245"

genintsetvalue(b, s)
genintstringrep(b, s)
writeoutputline("b =")
writeoutputline("")
writeoutputline(s)
writesep()

If geninterror(b) Then
geninterrormessage("bad b: --> script termination.")
Exit Function
End If

Console_WriteLine("working ..")
Console_WriteLine("")
genintdiv(a, b, q, r)

genintstringrep(q, s)
writeoutputline("q = a \ b, =")
writeoutputline("")
writeoutputline(s)
writesep()

genintstringrep(r, s)
writeoutputline("r = ")
writeoutputline("")
writeoutputline(s)
writesep()

genintmul(b, q, ck)
genintadd(ck, r, ck)

genintstringrep(ck, s)
writeoutputline("CHECK:")
writeoutputline("")
writeoutputline("a = b * q + r, =")
writeoutputline("")
writeoutputline(s)
writesep()

Console_WriteLine("The same output has been written to the file, 'Lovecraft7.txt'.")
Console_WriteLine("")
Console_WriteLine("Done.")
Console_WriteLine("")
Console_WriteLine("Press a key.")

FILE_Close(outfile)

WaitKey
End Function

'##########################################################################################

Function genintsetvalue(ByRef gn As genint, ByRef vs As String)
Local i, ln As DWord
Local v As Byte

ln = Len(vs)

v = Asc(Mid$(vs, 1, 1))

Select Case v
Case 43
genintsetsign(gn, %pos)
Case 45
genintsetsign(gn, %neg)
Case 48 To 57
genintsetsign(gn, %pos)
gn.digit(ln) = v - 48
Case Else
genintseterror(gn)
Return
End Select

For i = 2 To ln
v = Asc(Mid$(vs, i, 1))

If v < 48 Then
genintseterror(gn)
Return
EndIf

If v > 57 Then
genintseterror(gn)
Return
EndIf

gn.digit(ln - i + 1) = v - 48
Next
genintsetnsdigits(gn)
End Function

'--------------------------------------------------------------

Function genintsetdigit(ByRef gn As genint, index As DWord, digit As Byte)
gn.digit(index) = digit
genintsetnsdigits(gn)
End Function

'--------------------------------------------------------------

Function genintsetequal(ByRef a As genint, ByRef b As genint)
'b = a.

Local i As DWord

For i = 1 To %intsize
b.digit(i) = a.digit(i)
Next

b.nsdigits = a.nsdigits
b.sign = a.sign
b.error = a.error

End Function

'--------------------------------------------------------------

Function genintstringrep(ByRef gn As genint, ByRef vs As String)
Local i, nsd As DWord

nsd = gn.nsdigits

If (nsd = 0 Or nsd = 1) Then
If gn.digit(1) = 0 Then
vs = "0"
Return
EndIf
EndIf

If genintnegsign(gn) Then
vs = "-"
Else
vs = "+"
End If

For i = 1 To nsd
vs = vs & Chr$(48 + gn.digit(nsd - i + 1))
Next

End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function genintequalto(ByRef a As genint, ByRef b As genint) As Byte
'|a| = |b|, ?

Local i, asize, bsize As DWord
asize = a.nsdigits
bsize = b.nsdigits
If asize <> bsize Then Return FALSE
For i = 1 To asize
If a.digit(i) <> b.digit(i) Then Return FALSE
Next
Return TRUE
End Function

'--------------------------------------------------------------

Function genintlessthan(ByRef a As genint, ByRef b As genint) As Byte
'|a| < |b|, ?

Local i, j, asize, bsize As DWord
asize = a.nsdigits
bsize = b.nsdigits
If asize < bsize Then Return TRUE
If asize > bsize Then Return FALSE
For i = 1 To asize
j = asize + 1 - i
If a.digit(j) < b.digit(j) Then Return TRUE
If a.digit(j) > b.digit(j) Then Return FALSE
Next
Return FALSE
End Function

'--------------------------------------------------------------

Function genintgreaterthan(ByRef a As genint, ByRef b As genint) As Byte
'|a| > |b|, ?

Local i, j, asize, bsize As DWord
asize = a.nsdigits
bsize = b.nsdigits
If asize > bsize Then Return TRUE
If asize < bsize Then Return FALSE
For i = 1 To asize
j = asize + 1 - i
If a.digit(j) > b.digit(j) Then Return TRUE
If a.digit(j) < b.digit(j) Then Return FALSE
Next
Return FALSE
End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function genintsetnsdigits(ByRef gn As genint)
Local i As DWord

i = 0

Do
If %intsize - i = 0 Then Return 1
If gn.digit(%intsize - i) > 0 Then Exit Do
i += 1
Loop

gn.nsdigits = %intsize - i
If gn.nsdigits = 0 Then gn.nsdigits = 1
End Function

'--------------------------------------------------------------

Function genintzero(ByRef gn As genint)
Local i As DWord
For i = 1 To %intsize
gn.digit(i) = 0
Next
genintsetnsdigits(gn)
End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function genintsetsign(ByRef gn As genint, s As Integer)
gn.sign = s
End Function

'--------------------------------------------------------------

Function genintgetsign(ByRef gn As genint) As Integer
If gn.sign < 0 Then
Return -1
Else
Return 1
EndIf
End Function

'--------------------------------------------------------------

Function genintnegsign(ByRef gn As genint) As Byte
If gn.sign > -1 Then
Return FALSE
Else
Return TRUE
End If
End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function genintseterror(ByRef gn As genint)
gn.error = TRUE
End Function

'--------------------------------------------------------------

Function geninterror(ByRef gn As genint) As Byte
If gn.error Then
Return TRUE
Else
Return FALSE
End If
End Function

'--------------------------------------------------------------

Function geninterrormessage(es As String)
writeoutputline(es)
FILE_Close(outfile)
WaitKey
End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function genintshiftleftandadddigit(ByRef gn As genint, digit As Byte)

genintshiftleft(gn, 1)
genintsetdigit(gn, 1, digit)
genintsetnsdigits(gn)

End Function

'--------------------------------------------------------------

Function genintshiftleft(ByRef gn As genint, nplaces As DWord)
Local sz, i, j As DWord


genintsetnsdigits(gn)
sz = gn.nsdigits

For i = 1 To sz
j = sz - i + 1
gn.digit(j + nplaces) = gn.digit(j)
Next

For i = 1 To nplaces
gn.digit(i) = 0
Next

genintsetnsdigits(gn)
End Function

'--------------------------------------------------------------

Function genintmulbydigit(ByRef a As genint, ByRef b As genint, digit As Byte)
'b = digit * a

Local i As DWord
Local p, v, k As Byte

genintzero(b)

If digit = 10 Then
genintsetequal(a, b)
genintshiftleft(b, 1)
Return
EndIf

k = 0
For i = 1 To a.nsdigits + 1
p = digit * a.digit(i) + k
v = Mod(p, 10)
k = p \ 10
b.digit(i) = v
Next

genintsetnsdigits(b)
End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function genintadd(ByRef a As genint, ByRef b As genint, ByRef c As genint)
'c = a + b.

Local asign, bsign As Integer
Local asignpos, asignneg, bsignpos, bsignneg

asign = genintgetsign(a)
If asign = %pos Then
asignpos = TRUE
asignneg = FALSE
Else
asignpos = FALSE
asignneg = TRUE
EndIf

bsign = genintgetsign(b)
If bsign = %pos Then
bsignpos = TRUE
bsignneg = FALSE
Else
bsignpos = FALSE
bsignneg = TRUE
EndIf

If (asignpos And bsignpos) Then
genintadd1(a, b, c)
genintsetnsdigits(c)
Return
EndIf

If (asignpos And bsignneg) Then
genintsetsign(b, %pos)
genintsub(a, b, c)
genintsetsign(b, %neg)
genintsetnsdigits(c)
Return
EndIf

If (asignneg And bsignpos) Then
genintsetsign(a, %pos)
genintsub(b, a, c)
genintsetsign(a, %neg)
genintsetnsdigits(c)
Return
EndIf

If (asignneg And bsignneg) Then
genintadd1(a, b, c)
genintsetsign(c, %neg)
genintsetnsdigits(c)
Return
EndIf

End Function

'--------------------------------------------------------------

Function genintadd1(ByRef a As genint, ByRef b As genint, ByRef c As genint)
'Work function for genintadd.

Local i As DWord
Local s, v, k As Byte
Local temp As genint

Local asize, bsize, maxsize As DWord

asize = a.nsdigits
bsize = b.nsdigits
maxsize = asize
If bsize > maxsize Then maxsize = bsize

genintzero(temp)

k = 0
For i = 1 To maxsize + 1
s = a.digit(i) + b.digit(i) + k
v = Mod(s, 10)
k = s \ 10
temp.digit(i) = v
Next

For i = 1 To %intsize
c.digit(i) = temp.digit(i)
Next

End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function genintsub(ByRef a As genint, ByRef b As genint, ByRef c As genint)
'c = a - b.

Local asign, bsign As Integer
Local asignpos, asignneg, bsignpos, bsignneg

asign = genintgetsign(a)
If asign = %pos Then
asignpos = TRUE
asignneg = FALSE
Else
asignpos = FALSE
asignneg = TRUE
EndIf

bsign = genintgetsign(b)
If bsign = %pos Then
bsignpos = TRUE
bsignneg = FALSE
Else
bsignpos = FALSE
bsignneg = TRUE
EndIf

If (asignpos And bsignpos) Then
If genintlessthan(a, b) Then
genintsub1(b, a, c)
genintsetsign(c, %neg)
genintsetnsdigits(c)
Return
Else
genintsub1(a, b, c)
genintsetnsdigits(c)
Return
End If
EndIf

If (asignpos And bsignneg) Then
genintsetsign(b, %pos)
genintadd(a, b, c)
genintsetsign(b, %neg)
genintsetnsdigits(c)
Return
EndIf

If (asignneg And bsignpos) Then
genintsetsign(a, %pos)
genintadd(a, b, c)
genintsetsign(a, %neg)
genintsetsign(c, %neg)
genintsetnsdigits(c)
Return
EndIf

If (asignneg And bsignneg) Then
If genintlessthan(a, b) Then
genintsub1(b, a, c)
genintsetnsdigits(c)
Return
Else
genintsub1(a, b, c)
genintsetsign(c, %neg)
genintsetnsdigits(c)
Return
End If
EndIf

End Function

'--------------------------------------------------------------

Function genintsub1(ByRef a As genint, ByRef b As genint, ByRef c As genint)
'Work function for genintsub.

Local bsize As DWord
Local i As DWord
Local d, v, bor As Integer
Local temp As genint

bsize = b.nsdigits

genintzero(temp)

For i = 1 To bsize + 1
d = a.digit(i) - b.digit(i) - bor

If d < 0 Then
d += 10
bor = 1
Else
bor = 0
End If

temp.digit(i) = d
Next

For i = 1 To %intsize
c.digit(i) = temp.digit(i)
Next

End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function genintmul(ByRef a As genint, ByRef b As genint, ByRef c As genint)
'c = a * b.

Local i, j As DWord
Local p, v, k As Byte
Local asize, bsize, tsize As DWord
Local temp As genint

genintzero(c)
genintzero(temp)

asize = a.nsdigits
bsize = b.nsdigits
tsize = asize + bsize
If tsize > %intsize Then
genintseterror(c)
Return
EndIf

For i = 1 To asize
k = 0
For j = 1 To bsize + 1
p = a.digit(i) * b.digit(j) + k
v = Mod(p, 10)
k = p \ 10
temp.digit(j) = v
Next

genintsetnsdigits(temp)
genintshiftleft(temp, i - 1)

genintadd(c, temp, c)
genintzero(temp)
Next

If ((genintgetsign(a) * genintgetsign(b)) < 0) Then genintsetsign(c, %neg)

End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function genintdiv(ByRef a As genint, ByRef b As genint, ByRef c As genint, ByRef d As genint)
'c = a \ b. (integer division)
'd = the remainder of a \ b.

Local i, qindex As DWord
Local asign, bsign As Integer
Local dr, dd, qu, diff As genint
Local testdr As genint

If genintgreaterthan(b, a) Then
genintzero(c)
genintsetequal(b, d)
If ((genintgetsign(a) * genintgetsign(b)) < 0) Then genintsetsign(d, %neg)
Return
EndIf

asign = a.sign
bsign = b.sign

a.sign = %pos
b.sign = %pos

genintsetequal(b, dr)

For i = 1 To dr.nsdigits
dd.digit(i) = a.digit(a.nsdigits - dr.nsdigits + i)
Next

genintsetnsdigits(dd)
qindex = a.nsdigits - dd.nsdigits + 1

Do

While genintgreaterthan(dr, dd)
qindex -= 1
If qindex = 0 Then Exit Do
genintshiftleftandadddigit(dd, a.digit(qindex))
Wend

i = 1

Do
i += 1
genintmulbydigit(dr, testdr, i)
If genintgreaterthan(testdr, dd) Then Exit Do
Loop

genintsetdigit(qu, qindex, i - 1)
genintzero(testdr)
genintmulbydigit(dr, testdr, i - 1)
genintsub(dd, testdr, diff)
genintsetequal(diff, dd)

Loop

a.sign = asign
b.sign = bsign
genintsetequal(qu, c)
genintsetequal(dd, d)

If asign * bsign < 0 Then
genintsetsign(c, %neg)
genintsetsign(d, %neg)
EndIf

End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function genintsquareroot(ByRef a As genint, ByRef b As genint)
' b = a ^ (1/2)

Local bot, top, two, test, square, r As genint
Local s As String

s = "2"
genintsetvalue(two, s)
s = "0"
genintsetvalue(bot, s)
genintsetequal(a, top)

Do

genintadd(bot, top, test)
genintdiv(test, two, test, r)
genintmul(test, test, square)
If genintequalto(square, a) Then
genintsetequal(test, b)
Return
EndIf

If genintlessthan(square, a) Then
genintsetequal(test, bot)
EndIf

If genintgreaterthan(square, a) Then
genintsetequal(test, top)
EndIf

genintstringrep(square, s)
Console_WriteLine(s)

Loop

End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function writesep()
writeoutputline("")
writeoutputline("-------------------------")
writeoutputline("")
End Function

'--------------------------------------------------------------

Function writeoutputline(s As String)
Console_WriteLine(s)
FILE_LinePrint(outfile, s)
End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

[Johannes]

Saw this thread today. At work. I'm glad I don't have thinBasic at work because I would have done this on the boss' dime.

Dan, first of all, thank you. I'm into the eighth month of the current project (three more months to go :() and your task is the first thing not from my own mind that I had to shift my brain out of first gear for. Hell, out of idle. But straight into overdrive. I cheated a little though: I reused some of my variable-precision stuff I had lying around. That was fixed-point with oodles of decimals and the change to integer was minimal.

Second, I have no idea what you're doing in your program. I purposely did not look at your solution until I had mine working and your stuff baffles me. At a guess you're working entirely decimal in an array while I went binary. Well, 25 years of assembly coding the hard way will do that to you.:D I will definitely give your program a close look as I spotted a square root function in there. I love square roots.

Anyway, see attached. Not as extended as your program, and it doesn't have the nice output or the final multiply-and-add to check the results. I could have gone for full functionality of the big integer routines, but that wasn't the task you set. ;)

At least we got the same result. Given the two radically different approaches that's more of an assurance than Windows Calculator.



'----------------------------------------------------------------------------
' !Lovecraft VII
'
' Project Lovecraft Problem 7
'
' Quoth Dan Baron:
'
' Make a thinBasic script that finds the quotient,
' q = a \ b (Integer division),
' and the remainder, r ( 0 <= r < b ), such that,
' a = b * q + r.
' Where, a And b, are As follows.
'
' a = 90384734937482350659001366537465572115628370399301
' 56054826243653036262290428594830189978437629675825 41408949179799124552
'
' b = 75638641208056310245
'
' /Quoth
'
' http://www.thinbasic.com/community/showthread.php?10494-Project-Lovecraft-Problem-7
'
' Challenge accepted.
'
' 9 December 2010
'
'----------------------------------------------------------------------------
Uses "CONSOLE"
Dim a,eh,b,bee,q,cue,r,are As String
String n = String$(4,0)
Local i,j As Long
Local k,l As Quad
Local x As Quad = 65536*65536
a = "90384734937482350659001366537465572115628370399301 56054826243653036262290428594830189978437629675825 41408949179799124552"
b = "75638641208056310245"
'a = "10000000000000000000000000000000000000000000000000000000000000000"
'b = "100000000000000000000000000000000"
' Print a and b.
PrintL "Dividend :",a
PrintL "Divisor :",b
PrintL
' First convert the numbers to a more manageable format.
eh = varNum(a)
bee = varNum(b)
' Test print to see if all is well.
'PrintL varHex(eh)
'PrintL varHex(bee)
'PrintL
' Divide those puppies.
varDiv(eh,bee,cue,are)
' Test print the results.
'PrintL varHex(cue)
'PrintL varHex(are)
'PrintL
' Convert to decimal.
q = varStr(cue)
r = varStr(are)
' Print the actual results.
PrintL "Quotient :",q
PrintL "Remainder:",r
PrintL
WaitKey

' Convert ASCII number to internal.
'
Function varNum(a As String) As String
Local b As String = n
For i=1 To Len(a)
k=InStr("0123456789",Mid$(a,i,1))
If k>0 Then
k-=1
For j=Len(b)/4 To 1 Step -1
l=10*CVDWD(Mid$(b,4*j-3,4))+k
k=Int(l/x)
l=Mod(l,x)
Mid$(b,4*j-3,4)=MKDWD$(l)
Next j
If k>0 Then b=MKDWD$(k)+b
EndIf
Next i
Function = b
End Function
' Convert internal number to ASCII.
'
Function varStr(a As String) As String
Local b As String
Do
k=0
For i=1 To Len(a)/4
l=k*x+CVDWD(Mid$(a,4*i-3,4))
k=Mod(l,10)
l=Int(l/10)
Mid$(a,4*i-3,4)=MKDWD$(l)
Next i
b=TStr$(k)+b
If LEFT$(a,4)=n Then a=Mid$(a,5)
Loop Until a=""
Function = b
End Function
' Print hexadecimal representation of internal number.
'
Function varHex(a As String) As String
Local b As String
For i=1 To Len(a)/4
b+=Hex$(CVDWD(Mid$(a,4*i-3,4)),8)
Next i
Function = b
End Function
' Integer divide with remainder two internal numbers.
'
Sub varDiv(a As String,b As String,ByRef c As String,ByRef d As String)
Do While LEFT$(a,4)=n
a=Mid$(a,5)
Loop
Do While LEFT$(b,4)=n
b=Mid$(b,5)
Loop
If b="" Then
PrintL "Now you've done it. Vortex imminent. Division by zero."
End
EndIf
If Len(a)<Len(b) Or (Len(a)=Len(b) And a<b) Then
c=n
d=a
Else
' Straightforward binary division.
' Shift out as many bits as possible beforehand.
c=n
d=n+n+LEFT$(a,Len(b)-4)
a=Mid$(a,Len(b)-3)
b=n+b
' Do for number of remaining bits.
For j=1 To 8*Len(a)
' Shift MSB out of dividend.
k=0
For i=Len(a)/4 To 1 Step -1
l=2*CVDWD(Mid$(a,4*i-3,4))+k
k=Int(l/x)
l=Mod(l,x)
Mid$(a,4*i-3,4)=MKDWD$(l)
Next i
' Shift into LSB of remainder.
For i=Len(d)/4 To 1 Step -1
l=2*CVDWD(Mid$(d,4*i-3,4))+k
k=Int(l/x)
l=Mod(l,x)
Mid$(d,4*i-3,4)=MKDWD$(l)
Next i
' Is subtraction possible?
If varGOE(d,b) Then
' Subtract divisor from remainder.
k=0
For i=Len(d)/4 To 1 Step -1
l=CVDWD(Mid$(d,4*i-3,4))-CVDWD(Mid$(b,4*i-3,4))-k
If l<0 Then
k=1
l+=x
Else
k=0
EndIf
Mid$(d,4*i-3,4)=MKDWD$(l)
Next i
k=1
Else
k=0
EndIf
' Shift result bit (k) into quotient.
For i=Len(c)/4 To 1 Step -1
l=2*CVDWD(Mid$(c,4*i-3,4))+k
k=Int(l/x)
l=Mod(l,x)
Mid$(c,4*i-3,4)=MKDWD$(l)
Next i
If k>0 Then c=MKDWD$(k)+c
If RIGHT$(a,4)=n Then a=LEFT$(a,Len(a)-4)
Next j
Do While LEFT$(d,4)=n And Len(d)>4
d=Mid$(d,5)
Loop
EndIf
End Sub
' Compare if one internal number is greater or equal than the other.
'
Function varGOE(a As String,b As String) As Boolean
Local f As Boolean = TRUE
i=1
do
k=CVDWD(Mid$(a,4*i-3,4))
l=CVDWD(Mid$(b,4*i-3,4))
If k>l Then
i=Len(a)/4
Else
f=(k=l)
EndIf
i+=1
Loop Until i>Len(a)/4 Or Not f
Function = f
End Function

[danbaron]

Hi Johannes.

You're new here. So, first of all, we moved to a new forum program about a month ago. I've noticed some differences with this one.

If I am correct, this forum will not display a string that is longer than 50 characters. That is the reason that in Post #1, "a", has two spaces in it.

I tried to edit Post #2. I hadn't seen this thread since before the move. The code is no longer formatted as code. Also, I like my posts to be displayed in a
fixed width font. But, I couldn't save my changes to the post. The previous forum had a 20,000 character maximum. This one has a 10,000 character maximum. Post
#2 is approximately 18,000 characters.

I checked, and I agree with you, we both got the same answers. Yours is shorter and faster. And, you did it all in one day. So, I surrender.

To be honest, I don't understand your code, and, I don't understand much of mine. When I make a program, I always immediately forget about it. And, if I look at
it later, at first, I wonder who made it.

Usually, it's a terrible thankless task to have to work on someone else's code, isn't it?

All of my schooling and work have been in engineering. I've never done anything in assembly language. And, to me, the learning curve, seems steep; for instance,
to make a compiler (but, Charles Pegge (another member here), might disagree).

You are correct that everything I did was in decimal. As far as I am able to determine now, you did everything with strings (I never have used any of those
functions, CVDWD, MKDWD$, etc., maybe, I should), and, I did everything with arrays of bytes. I guess, my way is wasted effort. I guess, if we wanted to not
waste any RAM, we would both would have used base 256.

Concerning, doing your own stuff at work, I am the same way. My experience is that as long as you get your work-work done, no one cares.

--------------------------------------------------------

On another topic, I have been trying to work with a new language,

http://racket-lang.org/

It is a descendant of Scheme, which, is a descendant of Lisp.

It is free, and entirely self-contained - the IDE, compiler, documentation.

If you download and install it, it will go into, "c:\program files\racket".

The IDE is called, DrRacket.

From within DrRacket, you can run programs, or try stuff out in the "REPL", the read-eval-print loop.

If you so choose, download Racket, start DrRacket, and click on "Help Desk", from within the Help menu.

Your browser will open the index to the documentation it installed on your hard drive. You will have access to the guide and reference, among other things. In
the beginning, you only need the guide and reference.

Then, you can start learning it. The guide starts very gently, from zero. And, to me, that's a good thing. Because, Racket is like trying to learn Chinese, compared
to other languages I have been exposed to. But, I have the feeling that the language has advantages, which make it worth trying out.

REPL, is very useful, while you are going through the documentation.

My opinion so far, is that the entire Racket production, is made very well.

:p
Dan

[Johannes]

To be honest, I don't understand your code, and, I don't understand much of mine. When I make a program, I always immediately forget about it. And, if I look at it later, at first, I wonder who made it.

Thank God for that, I thought I was the only one. ;) I have written code that I did not recognise a month later, and took me more than a week to figure out again why it worked. (I could see it worked, but couldn't stand not knowing anymore why.)



Usually, it's a terrible thankless task to have to work on someone else's code, isn't it?

It's my job. :( Hired gun in software development. And to make matters worse, office automation. Order and work flow. At the very best you could call it a state machine, but it usually isn't.



All of my schooling and work have been in engineering. I've never done anything in assembly language. And, to me, the learning curve, seems steep; for instance, to make a compiler (but, Charles Pegge (another member here), might disagree).

I started assembly on the C64, many a moon ago. Back then it was either the crappy BASIC that was in the C64 or assembly. So assembly it was. I then switched to the Acorn Archimedes (and later RiscPC) and that machine was able to blend (proper) BASIC with assembly much like thinBasic can. (A bit more integrated on the Acorn but not to the degree that I couldn't work around it in thinBasic.) When speed is not of the utmost importance I'm not doing assembly. Although the learning curve is not an issue anymore for me, developing and testing is cumbersome. And let's face it, structured programming has it advantages.



You are correct that everything I did was in decimal. As far as I am able to determine now, you did everything with strings (I never have used any of those functions, CVDWD, MKDWD$, etc., maybe, I should), and, I did everything with arrays of bytes. I guess, my way is wasted effort. I guess, if we wanted to not waste any RAM, we would both would have used base 256.

The strings are the veil that obscures the true meaning. ;) It's just a way of implementing a dynamic array of doublewords. Because of the Acorn I'm used to working with doublewords ("words" on the Acorn because it's a 32-bit machine) or bytes. "Halfwords" on the Acorn are too much trouble so I never use them natively. Where you use base-ten digits I use base-four-gig (2^32) "words", witness the definition of the quad x in the beginning. Little if any memory wastage and of course very efficient should I ever want to implement this in assembler.

And I might. Big integer-arithmetic can come in handy. I might even go for total variable precision where you define a (maximum) number of decimals and the integer part is allowed to grow.

So much to do, so little time... :(

[ErosOlmi]

Sorry for posting limitation to 10000 chars but is is the default option of more than some thousand options this new forum software has. So it is not easy to fix all as we would like.

I've rised max post size up to 40000 chars per post.
Hope it is enough otherwise let me know.

Eros

[danbaron]

Eros fixed the post size limitation.

I know that both Python and Ruby have big integer arithmetic.

But, neither has big floating point operations.

For instance,

x = 1.568957384758712758796476859380397089386970894321096879498867463545219069584736475980796822666453345525526648879809799993665774776884277752099798099968546372647563475869204996812400978695847362109

Anyone who thinks he is not going to lose accuracy in almost any computer language implementation when he calculates sin(x), is dreaming.

In fact, in almost every computer language implementation, accuracy loss begins with x itself (floating point numbers are limited to 8 or 10 bytes). I think the limitation is due to computer architecture. Accuracy is sacrificed for speed.

But, I think Mathematica can do it. It can register x completely as it appears above. And then, it will not lose any accuracy when calculating its sine. I think it normally operates so that it does not lose accuracy during calculations.

I don't know how often it is really useful to be able to add, subtract, multiply, divide, and exponentiate giant integers. And, I don't know how often it is really useful to be able to perform floating point operations on numbers with a huge number of significant digits.

It's hard for me to think of uses for giant integers. I think the upper limit on the estimate of the number of particles in the universe, is 10^87.

http://www.strangehorizons.com/2001/20010402/biggest_numbers.shtml

I think similarly about extremely precise floating point numbers, mostly, no one needs them. I think the most precise calculations are done in astrophysics. The most precise calculation I remember reading about, was accurate to 1 part in 10^120. In the latest book by Hawking and Mlodinow, they claim that they have calculated the number of universes, to be in the neighborhood of 10^500.

Of course, the human mind can always contemplate numbers which are bigger or more precise than exist in physics. So, such numbers will always be useful for exercising the human imagination. And so, for instance, people have calculated pi to at least one million decimal places. But, I can't imagine that this precision will ever be required in a physical calculation.

So, basically, I think we are talking about being able to manipulate on the computer, numbers with an arbitrary number of significant digits, without losing any accuracy. I admit that being able to do so, is a neat, "parlor trick", which we can use to amaze and dumbfound our friends. And, implementing the code to do it, is fun, because, it's hard, but possible.

:p
Dan

[zak]

Hawking and Mlodinow, they claim that they have calculated the number of universes, to be in the neighborhood of 10^500
too big number, but my intuition said that the number of universes must be true infinite, and this is the only mechanical representation i can envisage to represent or to symbolize the shadows of the existence itself, anyway may be there is in one of the universes the same thinbasic forum with the same dan and and johannes and zak with this thread but after this "this" there is another sentence, and at the end of this message there is no dot.
regarding pi i found an article "Lady Pi" www.cadaeic.net/ladypi.htm (http://www.cadaeic.net/ladypi.htm) ,
the author said that he can see a lady in the hexagonal spiral of pi in binary base. look at the figure in that link i can't see the lady, who can please post arrows where that lady is. but on the infinite line of pi there is somewhere any possible combination, the web site http://www.cadaeic.net/ have a treasure of articles .

[Johannes]

x = 1.5689573847587127587964768593803970893869708943210968794[...]
Anyone who thinks he is not going to lose accuracy in almost any computer language implementation when he calculates sin(x), is dreaming.
But that's because you have a grasp of what it takes to store "a number". For the layman a number is a number is a number. They don't, can't, see the difference between a doubleword and an extended-precision float.

In fact, in almost every computer language implementation, accuracy loss begins with x itself (floating point numbers are limited to 8 or 10 bytes). I think the limitation is due to computer architecture. Accuracy is sacrificed for speed.
On x86 hardware you have single, double and extended precision using 4, 8 and 10 bytes respectively. Some of the bits are used for the exponent and so they end up with 24, 53 and 64 bits for the actual number. You can store log10(2) decimals in a bit so you end up with 7, 15 and 19 full significant decimals. With 12 significant decimals you can define PI so that you can describe any circle that fits inside the Earth to the millimetre (1/25th of an inch).

I don't know how often it is really useful to be able to add, subtract, multiply, divide, and exponentiate giant integers. And, I don't know how often it is really useful to be able to perform floating point operations on numbers with a huge number of significant digits.
As far as I know such calculations are only useful in a theoretical sense. Calculating the <n>th decimal of a certain irrational number, for whatever reason. For real-life purposes you can almost always make do with a certain number of significant digits and that's where floating-point comes in.

In the latest book by Hawking and Mlodinow, they claim that they have calculated the number of universes, to be in the neighborhood of 10^500.
And for that you don't even need calculations with a large precision. The number is so tenuous that they don't even supply a significant digit, only the magnitude. If it turns out to be 10^600 (I'd like to see someone prove that conclusively either way :)) I'd be more impressed by Hawking than I already am.

And so, for instance, people have calculated pi to at least one million decimal places. But, I can't imagine that this precision will ever be required in a physical calculation.
Been there, done that. :) With 34 decimals (35 significant digits) in PI you can describe any circle that fits inside the observable universe with a precision of the radius (diameter?) of a hydrogen atom. The latest IEEE 754 standard defines 128-bit floating point values which have 34 significant decimals. I think that should be enough for mere mortals.

And, implementing the code to do it, is fun, because, it's hard, but possible.
Been there, definitely done that. :) It's like climbing that mountain: you wrote the code because you could.

[Johannes]

anyway may be there is in one of the universes the same thinbasic forum with the same dan and and johannes and zak with this thread but after this "this" there is another sentence, and at the end of this message there is no dot.
You're talking about "branching" universes now, where "nearby" universes will be very much alike. Although in that other universe the other Johannes will disagree with me.


regarding pi i found an article "Lady Pi" www.cadaeic.net/ladypi.htm (http://www.cadaeic.net/ladypi.htm), the author said that he can see a lady in the hexagonal spiral of pi in binary base. look at the figure in that link i can't see the lady, who can please post arrows where that lady is. but on the infinite line of pi there is somewhere any possible combination, the web site http://www.cadaeic.net/ have a treasure of articles
Can't see it either, not even with the different colours. I can supply you with 3.32 million bits of pi should you feel the need to explore the picture further.

If it is true that PI is irrational and that it never repeats, then by necessity you will find anything that can be represented binary (or decimally) somewhere in PI. Somewhere in all those bits that represent PI you will find all of Shakespeare's works, cover to cover, in reverse, in Klingon.

[Charles Pegge]

It is fascinating thar you can calculate any binary/hexadecimal digit of pi without going through all the previous digits. Some deep Mathematical principle at work here!

http://pi.nersc.gov/

Charles

[John Spikowski]

Wow. Something cast in stone as pi having new revelations is amazing stuff.

My effort in this area this weekend was of the Dutch Apple variety.

[Charles Pegge]

Formula for the fair distribution of Pie between 8 guests.

pie/8 = 1-1/3+1/5-1/7+1/9-1/11 ...

where pie=2*pi :)

Charles

[danbaron]

I like mountains, too. I wish I could go and climb some, at least in the Alps.

My favorite mountain, is K2. I believe it is the most dangerous mountain in the world. And, it is hundreds of miles from anywhere. (If you want to get killed, then, try to climb K2. And, if you want to guarantee your success (??), then, do it alone.)

But, believe it or not, unless I am mistaken, the mountain with the greatest vertical distance (from the base to the top), is in Alaska, Mt. McKinley (Denali). If I remember correctly, its vertical is approximately 17,000 feet, and Mt. Everest's is approximately 13,000 feet (the bottom of Mt. Everest (at the base camp) is at approximately 16,000 feet above sea level).

If I remember correctly, the top of Mt. Everest is 29,028 feet above sea level. There are at least two bad things about being that high. The first is that (from memory), the oxygen content of the air is approximately 1/3 of what it is at sea level. The second, is that the air pressure on your body is very low (I don't remember how low). The low air pressure can cause lots of potentially fatal physical problems. As I recall, your lungs can fill with blood, you can have a stroke, etc.

By the way, as I remember, most people who are killed on Mt. Everest, die on the way down.

And, I think that climbing Mt. Everest, or any of the hardest mountains, is actually a relatively short trip. Say, you start from the base camp of Mt. Everest at 16,000 feet, and you climb to the top, at 29,000 feet. There is a relatively easy ridge that some people take. Say, the average angle of that ridge is 45 degrees. In that case, how far is your trip to the top? It would be approximately,

13,000 / sin(pi/4) = 18385 feet, = 3.48 miles.

In that case, the entire trip, from the base camp, to the top, and back to the base camp, would be approximately 7 miles. And, people pay lots and lots of money to make that 7 mile trip.

----------------------------------

Concerning pi, I don't know where Charles finds all of this stuff.

Maybe it's true, but the idea that there is a formula to calculate the nth digit of pi, for any value of n, seems unbelievable to me. In that case, we could make a thinBasic function that takes one parameter, n, and returns the nth digit of pi. I'd like to see it. (Incidentally, since n can be any finite digit, it would be impossible for any computer to have a data type that could accommodate its entire range, yes?)

It seems to me there must be something which the summary on Charles' link, is leaving out. At the moment, I don't have the inclination to read the references.

I think common sense would indicate that if there is such a formula, then, obviously the digits of pi are not random. I think that by definition, a random sequence is a sequence for which there is no formula for finding the nth digit.

But, even if there was no such formula, I think you could still make the case that the pi sequence is not random, just because, it is calculable. (In a sense, even before this new formula was discovered, there already were many formulas for calculating the nth digit of pi, for any value of n. They just all contained first calculating the n-1 digits.)

I'm no expert. I know in mathematics, the "water", can be extremely deep. For the hobbyist, confusion is always one inch away. I think it's the same way in advanced physics. Most of us are not in a situation in which we can devote our lives to becoming masters - modern wizards.

:p

[zak]

i have found the missing paper digits.pdf refered to in the site mentioned by Charles, using web.archive.org
i attached here, they said in the paper digits.pdf "We demonstrate this technique by computing the ten billionth hexadecimal digit of pi", the paper is extremely technical, yes it will be good to have a "top secret" thinbasic function for that "with no source code allowed".
i have searched for my name in http://pi.nersc.gov/ ,birthday date ... .the site depends on an already available database of pi digits.

[Charles Pegge]

This is an algorithm published by David Bailey for calculating hexadecimal pi digits
(derived from the Bailey–Borwein–Plouffe formula).

Would anyone like to port it?

Charles

C code


/*
This program implements the BBP algorithm to generate a few hexadecimal
digits beginning immediately after a given position id, or in other words
beginning at position id + 1. On most systems using IEEE 64-bit floating-
point arithmetic, this code works correctly so long as d is less than
approximately 1.18 x 10^7. If 80-bit arithmetic can be employed, this limit
is significantly higher. Whatever arithmetic is used, results for a given
position id can be checked by repeating with id-1 or id+1, and verifying
that the hex digits perfectly overlap with an offset of one, except possibly
for a few trailing digits. The resulting fractions are typically accurate
to at least 11 decimal digits, and to at least 9 hex digits.
*/

/* David H. Bailey 2006-09-08 */

#include <stdio.h>
#include <math.h>

main()
{
double pid, s1, s2, s3, s4;
double series (int m, int n);
void ihex (double x, int m, char c[]);
int id = 1000000;
#define NHX 16
char chx[NHX];

/* id is the digit position. Digits generated follow immediately after id. */

s1 = series (1, id);
s2 = series (4, id);
s3 = series (5, id);
s4 = series (6, id);
pid = 4. * s1 - 2. * s2 - s3 - s4;
pid = pid - (int) pid + 1.;
ihex (pid, NHX, chx);
printf (" position = %i\n fraction = %.15f \n hex digits = %10.10s\n",
id, pid, chx);
}

void ihex (double x, int nhx, char chx[])

/* This returns, in chx, the first nhx hex digits of the fraction of x. */

{
int i;
double y;
char hx[] = "0123456789ABCDEF";

y = fabs (x);

for (i = 0; i < nhx; i++){
y = 16. * (y - floor (y));
chx[i] = hx[(int) y];
}
}

double series (int m, int id)

/* This routine evaluates the series sum_k 16^(id-k)/(8*k+m)
using the modular exponentiation technique. */

{
int k;
double ak, eps, p, s, t;
double expm (double x, double y);
#define eps 1e-17

s = 0.;

/* Sum the series up to id. */

for (k = 0; k < id; k++){
ak = 8 * k + m;
p = id - k;
t = expm (p, ak);
s = s + t / ak;
s = s - (int) s;
}

/* Compute a few terms where k >= id. */

for (k = id; k <= id + 100; k++){
ak = 8 * k + m;
t = pow (16., (double) (id - k)) / ak;
if (t < eps) break;
s = s + t;
s = s - (int) s;
}
return s;
}

double expm (double p, double ak)

/* expm = 16^p mod ak. This routine uses the left-to-right binary
exponentiation scheme. */

{
int i, j;
double p1, pt, r;
#define ntp 25
static double tp[ntp];
static int tp1 = 0;

/* If this is the first call to expm, fill the power of two table tp. */

if (tp1 == 0) {
tp1 = 1;
tp[0] = 1.;

for (i = 1; i < ntp; i++) tp[i] = 2. * tp[i-1];
}

if (ak == 1.) return 0.;

/* Find the greatest power of two less than or equal to p. */

for (i = 0; i < ntp; i++) if (tp[i] > p) break;

pt = tp[i-1];
p1 = p;
r = 1.;

/* Perform binary exponentiation algorithm modulo ak. */

for (j = 1; j <= i; j++){
if (p1 >= pt){
r = 16. * r;
r = r - (int) (r / ak) * ak;
p1 = p1 - pt;
}
pt = 0.5 * pt;
if (pt >= 1.){
r = r * r;
r = r - (int) (r / ak) * ak;
}
}

return r;
}


http://en.wikipedia.org/wiki/Pi
http://www.experimentalmath.info/bbp-codes/
http://www.experimentalmath.info/bbp-codes/bbp-alg.pdf

[zak]

thank you Charles very much for posting the algorithm, i have compiled it with vc6 changing the int id = 10000000; it takes about 2 minutes to produce :hex digits = 7AF5863EFF
to check the program for small range also i changed int id = 5;
it produced A8885A308D , and going to http://www.wolframalpha.com/input/?i=pi+in+hexadecimal
it gives
3.243f6a8885a308d313198a2e03707344a4093822299f31d0082efa98ec4e6c89452821e638d01377be5466cf34e90c6cc0ac29b7c97c50dd3f84d5b5b5470..._16
so it gives a correct answer, this is like a dream.
poor william james who said:
"If we take the geometrical relations, the thousandth decimal of pi sleeps there, though no one may ever try to compute it.”
William James, The Meaning of Truth
http://www.authorama.com/meaning-of-truth-9.html

[danbaron]

I changed the program a little bit, so that it will print the elapsed time.

I compiled and ran it three times, using Pelles C.

________id_____seconds
___100,000_______0.660
_1,000,000_______6.707
10,000,000______72.231

You can see that the time increases faster than linearly.

Here is a summary of the BBP algorithm.

http://en.wikipedia.org/wiki/BBP_algorithm

According to the link,

"The algorithm is the fastest way to compute the nth digit (or a few digits in a neighborhood of the nth), but π-computing algorithms using large data types remain faster when the goal is to compute all the digits from 1 to n."


/*
This program implements the BBP algorithm to generate a few hexadecimal
digits beginning immediately after a given position id, or in other words
beginning at position id + 1. On most systems using IEEE 64-bit floating-
point arithmetic, this code works correctly so long as d is less than
approximately 1.18 x 10^7. If 80-bit arithmetic can be employed, this limit
is significantly higher. Whatever arithmetic is used, results for a given
position id can be checked by repeating with id-1 or id+1, and verifying
that the hex digits perfectly overlap with an offset of one, except possibly
for a few trailing digits. The resulting fractions are typically accurate
to at least 11 decimal digits, and to at least 9 hex digits.
*/

/* David H. Bailey 2006-09-08 */

#include <stdio.h>
#include <math.h>
#include <time.h>

int main()
{
clock_t c1, c2;
double et;
double pid, s1, s2, s3, s4;
double series (int m, int n);
void ihex (double x, int m, char c[]);
int id = 100000;
#define NHX 16
char chx[NHX];

/* id is the digit position. Digits generated follow immediately after id. */

c1 = clock();
s1 = series (1, id);
s2 = series (4, id);
s3 = series (5, id);
s4 = series (6, id);
pid = 4. * s1 - 2. * s2 - s3 - s4;
pid = pid - (int) pid + 1.;
ihex (pid, NHX, chx);
c2 = clock();
et = (double)(c2 - c1)/ CLOCKS_PER_SEC;
printf (" position = %i\n fraction = %.15f \n hex digits = %10.10s\n", id, pid, chx);
printf(" elapsed time = %7.3f seconds\n", et);
}

void ihex (double x, int nhx, char chx[])

/* This returns, in chx, the first nhx hex digits of the fraction of x. */

{
int i;
double y;
char hx[] = "0123456789ABCDEF";

y = fabs (x);

for (i = 0; i < nhx; i++){
y = 16. * (y - floor (y));
chx[i] = hx[(int) y];
}
}

double series (int m, int id)

/* This routine evaluates the series sum_k 16^(id-k)/(8*k+m)
using the modular exponentiation technique. */

{
int k;
double ak, p, s, t;
double expm (double x, double y);
#define eps 1e-17

s = 0.;

/* Sum the series up to id. */

for (k = 0; k < id; k++){
ak = 8 * k + m;
p = id - k;
t = expm (p, ak);
s = s + t / ak;
s = s - (int) s;
}

/* Compute a few terms where k >= id. */

for (k = id; k <= id + 100; k++){
ak = 8 * k + m;
t = pow (16., (double) (id - k)) / ak;
if (t < eps) break;
s = s + t;
s = s - (int) s;
}
return s;
}

double expm (double p, double ak)

/* expm = 16^p mod ak. This routine uses the left-to-right binary
exponentiation scheme. */

{
int i, j;
double p1, pt, r;
#define ntp 25
static double tp[ntp];
static int tp1 = 0;

/* If this is the first call to expm, fill the power of two table tp. */

if (tp1 == 0) {
tp1 = 1;
tp[0] = 1.;

for (i = 1; i < ntp; i++) tp[i] = 2. * tp[i-1];
}

if (ak == 1.) return 0.;

/* Find the greatest power of two less than or equal to p. */

for (i = 0; i < ntp; i++) if (tp[i] > p) break;

pt = tp[i-1];
p1 = p;
r = 1.;

/* Perform binary exponentiation algorithm modulo ak. */

for (j = 1; j <= i; j++){
if (p1 >= pt){
r = 16. * r;
r = r - (int) (r / ak) * ak;
p1 = p1 - pt;
}
pt = 0.5 * pt;
if (pt >= 1.){
r = r * r;
r = r - (int) (r / ak) * ak;
}
}

return r;
}

[zak]

Thanks Dan for adding time, using your ranges, with vc6 my results like this:
for 100000 -> 0.515 seconds
for 1000000 -> 6.265 seconds
for 10000000 -> 71.500 seconds
in the code it is said that If 80-bit arithmetic can be employed, this limit (10^7) is significantly higher, i remember that there is a 80 bit calc in thinbasic, and may be in powerbasic.
amazing this world, here http://www.youtube.com/watch?v=CSGK8IEXvXU
"Marc Umile" write pi from memory to 1000 digits. and in an article "On Pi Day" http://www.cnn.com/2010/TECH/03/12/pi.day.math/ it is said that he typed out 15,314 digits from memory in 2007.

[danbaron]

I think that both thinBasic and PowerBasic have the 80 bit floating point data type.

I remember reading about Andrew Wiles, the mathematician who proved Fermat's Last Theorem. His wife is also a mathematician. Apparently, they both independently memorized pi to 1000 decimal places. Maybe that was what initially attracted them to each other. Anyway, one day at the university, when they were dating, they were outside reciting the first thousand digits of pi, when someone else walked up in the middle of their recitation, and was able to correctly finish the sequence with them!

Here is something along the same line.

http://www.telegraph.co.uk/news/newstopics/howaboutthat/1940420/The-woman-who-can-remember-everything.html

And, here is the most amazing feat of a human brain that I am aware of.

http://en.wikipedia.org/wiki/Alexis_Lemaire

[zak]

it seems getting nth digit of Pi in decimal is much harder,
the mathematicians are aliens , there are already algorithms for hexa and bin and 4 and 5 but they forgot themselves in that they buy goods from the supermarket in decimal and not hexadecimal.

[peter]

cool Zak! lol

[Johannes]

This is an algorithm published by David Bailey for calculating hexadecimal pi digits (derived from the Bailey–Borwein–Plouffe formula).

Would anyone like to port it?
Challenge accepted. :cool:

Mind you, I didn't port the C program but wrote the program using the Wikipedia description of the algorithm.

It is not as fast as the C program (duh) but after tomorrow (ahem, looking at the time that should be "today") I've got the rest of the year off. I might give floating-point assembler a shot and see what I can come up with. But although my program lacks in speed I think it's more elegant and readable than the C example given.

This is why I hate C, in any form, and I was right to solemnly vow, all those years ago, to never, ever code a single line in anything called 'C'. (I made the vow a little more specific, but I don't think you'd appreciate the NC17 langauge. If you want to know my precise vow, send me a PM, with proof of age.)

Not to rain on Mr Bailey's parade, but just look at that C mess. Sometimes I think people who code in C are not clever enough, or too clever for their own good. Take for instance his running divisor (variable ak) which is constantly evaluated as 8*k+m. I don't care how efficient your C compiler is, my way of doing the same by repeatedly adding 8 (variable a) is shorter and quicker in any language or compiler.

If you want to make the program quicker, change the data type from Double to Single. (Perhaps changing the Quads to Longs will help as well. But then you'd be limited to skipping only 2 billion hex digits. ;)) If you want more hex digits change it to Ext.

Enjoy,

Johannes


'----------------------------------------------------------------------------
' !PI Hex Digit
'
' Calculate nth hexadecimal digit in PI
'
' http://en.wikipedia.org/wiki/BBP_algorithm
'
' 22 december 2010
'----------------------------------------------------------------------------

Long n = 0 ' Number of hexadecimal digits to skip
Double hd,s,t,u ' Overall precision calculation (Single, Double, Ext)
Quad i,j,k,m,p,q,r,x
String digits
HiResTimer_Init
HiResTimer_Get
BBPpart(1) : hd+=4*s
BBPpart(4) : hd-=2*s
BBPpart(5) : hd-=s
BBPpart(6) : hd-=s
x=HiResTimer_Delta
Select Case SizeOf(hd) ' Determine data type for number of whole hexadecimal digits
Case = 4
i=6
Case = 8
i=13
Case = 10
i=16
End Select
For j=1 To i-1 ' Ignore last hexadecimal digit (may be incorrect)
hd*=16
digits+=Hex$(Int(hd),1)
Next i
MsgBox 0,"Hexadecimal PI digits after #"+TStr$(n)+": "+digits+$CRLF+$CRLF+"Calculation time: "+TStr$(x/1000000)+" seconds"

Sub BBPpart(a As Quad) ' calculate BBP term
s=0 ' running summation
k=0
Do While k<n ' skips when n=0
p=n-k ' exponent for 16
q=1 ' running 16^(n-k) mod 8k+1
r=Mod(16,a) ' mod is allowed everywhere when multiplying
While p ' binary exponentation
If (p And 1) Then
q=Mod(q*r,a) ' multiply with running factor if bit is on
EndIf
r=Mod(r*r,a) ' running exponentation factor
p=Int(p/2) ' process all bits
Wend
s+=Frac(q/a) ' divide by 8k+1
a+=8 ' running 8k+1
k+=1 ' next k
Loop
s=Frac(s) ' discard integer part
u=1 ' running 16^(n-k), k is implied in routine
Do
t=u/a ' next summation element
s+=t ' running summation
u/=16 ' running power of 16
a+=8 ' running 8k+1
Loop Until t=0 ' get best precision
s=Frac(s) ' discard integer part
End Function

[danbaron]

_____________________I ran Johannes' program.

(I didn't change anything to make it faster.)

________________________Here are the results.

_______________________________n______seconds
_________________________100,000_______24.903
_______________________1,000,000______297.235
______________________10,000,000_____3521.433

___________________________________________:p
__________________________________________Dan

[Johannes]

Dan, thanks for the speed comparison. As I said, it's lacking in that department. :p

When I woke up this morning I realised that I have made a mistake in the calculation of the second part (the very last loop in the BBPpart subroutine). I am ending the loop when variable t reaches zero, but that should be when t becomes smaller than (from the top of my head) 1/16^number_of_hex_digits_in_data_type. Perhaps one or two higher. For single precision this means that I'm running that loop now up to 32 times instead of 7 or 8, for double precision (which Dan timed) up to 256 times instead of 14 or 15, and for extended it's a whopping 4096 times instead of 17 or 18.

Will update the program once I get back from work (where I am now).

[Johannes]

The big time consumer is the "skipping" of the leading hex digits so the speed increase described above will be negligible for the millionth or ten millionth hex digit. I have also been thinking about improving that skipping loop (dwords instead of quads, quitting when variable r becomes 0 or 1) and will introduce those changes as well.

And the more I think about it the more I am convinced that pure assembler might see a huge speed improvement.

I think I'm going to take the afternoon off. :D

[zak]

amazing that my results now using thinbasic is speedier significantly than Dan results, but about the same when using c compiler as listed above.
my results using thinbasic
for 100,000 14 seconds
for 1000,000 166 seconds

this thread turns out to be very usefull either for educational purposes or for benchmarking
thanks for all

[danbaron]

Maybe your version of thinBasic is newer than mine, Zak. Mine is 1.8.0. I don't know why else there would be such a speed difference, since our times for the C version were similar. But, I'm not sure, the great time differences for the thinBasic version seem strange to me.

Johannes, I don't know exactly what you mean when you mention, "floating-point assembler", and, "pure assembler". I am ignorant in this area. I know there are various, "assemblers", which will translate assembly code into an "exe" file. I know that PowerBasic allows you to include assembly code in a source code file, which it then compiles. If I'm correct, in thinBasic, you can use the Oxygen module, which will compile thinBasic code, just before it executes.

[zak]

yes my thinbasic version is newer it is the latest 1.8.6, but may be other factors playing here. such as software/ hardware related and how every programming languages deals with.

[Johannes]

Johannes, I don't know exactly what you mean when you mention, "floating-point assembler", and, "pure assembler". I am ignorant in this area. I know there are various, "assemblers", which will translate assembly code into an "exe" file. I know that PowerBasic allows you to include assembly code in a source code file, which it then compiles. If I'm correct, in thinBasic, you can use the Oxygen module, which will compile thinBasic code, just before it executes.
A misnomer, really. When I say "floating-point assembler" I really mean "regular assembler using floating-point coprocessor commands". Naturally I would do this using Oxygen as it provides the most simple interface to assember/machine code for me.

I just changed my program (remote desktop, I love it) to change the final loop from the threshold for t=0 to t<1/(16^number of hex digits+1) and for double precision that gave an eightfold speed increase. (I skipped zero digits, to be able to time just that part.) However, that amounted to only 0.0018 seconds absolute.

I also did a speed test for 100,000 digits and my computer ran the program in 20.2 seconds, slightly faster than Dan's machine. The change described above had therefore no significant effect at all. Also, for any substantially large n it doesn't really matter which precision (single, double, ext) you use because the skipping of the leading hex digits will take up the bulk of the time.

Couldn't get the afternoon off so the other optimisations will have to wait until this evening. I'll do some tests to check the differences between the FP precisions with large values for n. I'm also thinking about doing the 16^n mod (8k+m) in assembler since that is the bottleneck. If I "limit" n to a DWord (instead of a Quad) I have great hopes for kicking dust in that C program's (inter)face. ;)

[zak]

porting johannes code to freebasic giving the following results:
for n=100,000 :
35EA16CA2B7C
1.24029894123305 seconds
for n=1000,000 :
6C65E51B0046
14.78548740882252 seconds
for n=10,000,000:
7AF58A7C0000
171.4573891846558 seconds

so we can say freebasic is speedy but slower than c (in this realm), compare with c results in page 2
here the freebasic code is almost the same as thinbasic with very small changes

'----------------------------------------------------------------------------
' !PI Hex Digit
'
' Calculate nth hexadecimal digit in PI
'
' http://en.wikipedia.org/wiki/BBP_algorithm
'
' 22 december 2010
'----------------------------------------------------------------------------
Dim As Double start,finish
Dim Shared As Long n
n = 100000 ' Number of hexadecimal digits to skip
Dim Shared As Double hd,s,t,u ' Overall precision calculation (Single, Double, Ext)
Dim Shared As LongInt i,j,k,m,p,q,r,x
Dim As String digits
Sub BBPpart(a As LongInt) ' calculate BBP term
s=0 ' running summation
k=0
Do While k<n ' skips when n=0
p=n-k ' exponent for 16
q=1 ' running 16^(n-k) mod 8k+1
'r=Mod(16,a) ' mod is allowed everywhere when multiplying
r=16 Mod a
While p ' binary exponentation
If (p And 1) Then
'q=Mod(q*r,a) ' multiply with running factor if bit is on
q= q*r Mod a
EndIf
'r=Mod(r*r,a) ' running exponentation factor
r=r*r Mod a
p=Int(p/2) ' process all bits
Wend
s+=Frac(q/a) ' divide by 8k+1
a+=8 ' running 8k+1
k+=1 ' next k
Loop
s=Frac(s) ' discard integer part
u=1 ' running 16^(n-k), k is implied in routine
Do
t=u/a ' next summation element
s+=t ' running summation
u/=16 ' running power of 16
a+=8 ' running 8k+1
Loop Until t=0 ' get best precision
s=Frac(s) ' discard integer part
End Sub
start = Timer
BBPpart(1) : hd+=4*s
BBPpart(4) : hd-=2*s
BBPpart(5) : hd-=s
BBPpart(6) : hd-=s
Select Case SizeOf(hd) ' Determine data type for number of whole hexadecimal digits
Case 4
i=6
Case 8
i=13
Case 10
i=16
End Select
For j=1 To i-1 ' Ignore last hexadecimal digit (may be incorrect)
hd*=16
digits+=Hex(Int(hd),1)
Next j
Print digits
finish =Timer
Print (finish - start)
Sleep
End

[Johannes]

porting johannes code to freebasic giving the following results:
for n=100,000 :
35EA16CA2B7C
1.24029894123305 seconds
I'm impressed by the speed but I have a problem with Freebasic's results. And after checking with my own full PI digit calculator (its output is decimal but the internal result is of course binary and I dumped that internal result as hex to disc) I have a problem with my own program's results as well. Hex digits 100,001 through 100,016 are 35EA16C406363A30B and both Freebasic and thinBasic go wrong after the 'C', and in the same way. Looks like double precision is only enough for 7 hex digits (28 bits). With extended precision the next two hex digits are also correct, but after that 9th hex digit it stops. :(

But I have some interesting timings. The code as I uploaded before, with n=0 (so I can test the speeds of the different floating point precisions) runs on my machine as follows.

Extended precision (80 bits data, 64 bits precision): 0,203 seconds
Double precision (64 bits data, 53 bits precision): 0,00208 seconds
Single precision (32 bits data, 24 bits precision): 0,000090 seconds

Then I corrected that final calculation loop, to stop calculation once the precision mentioned above is crossed (variable v).


'----------------------------------------------------------------------------
' !PI Hex Digit
'
' Calculate nth hexadecimal digit in PI
'
' http://en.wikipedia.org/wiki/BBP_algorithm
'
' 22, 23 december 2010
'----------------------------------------------------------------------------

Quad n = 100000 ' Number of hexadecimal digits to skip
Ext hd,s,t,u,v ' Overall precision calculation (Single, Double, Ext)
Quad i,j,k,m,p,q,r
String digits

Select Case SizeOf(hd) ' Determine data type for number of whole hexadecimal digits
Case = 4
j=6
Case = 8
j=13
Case = 10
j=16
End Select
v=1/(16^(j+1))

HiResTimer_Init
HiResTimer_Get
BBPpart(1) : hd+=4*s
BBPpart(4) : hd-=2*s
BBPpart(5) : hd-=s
BBPpart(6) : hd-=s
Quad x=HiResTimer_Delta

For i=1 To j-2 ' Ignore last two hexadecimal digits (may be incorrect) -- May be? Are! :-)
hd*=16
digits+=Hex$(Int(hd),1)
Next j

MsgBox 0,"Hexadecimal PI digits after #"+TStr$(n)+": "+digits+$CRLF+$CRLF+"Calculation time: "+TStr$(x/1000000)+" seconds"

Sub BBPpart(a As Quad) ' calculate BBP term
s=0 ' running summation
k=0
Do While k<n ' skips when n=0
p=n-k ' exponent for 16
q=1 ' running 16^(n-k) mod 8k+1
r=Mod(16,a) ' mod is allowed everywhere when multiplying
While p ' binary exponentation
If (p And 1) Then
q=Mod(q*r,a) ' multiply with running factor if bit is on
EndIf
r=Mod(r*r,a) ' running exponentation factor
p=Int(p/2) ' process all bits
Wend
s+=Frac(q/a) ' divide by 8k+1
a+=8 ' running 8k+1
k+=1 ' next k
Loop
s=Frac(s) ' discard integer part
u=1 ' running 16^(n-k), k is implied in routine
Do
t=u/a ' next summation element
s+=t ' running summation
u/=16 ' running power of 16
a+=8 ' running 8k+1
Loop Until t<v ' get best precision
s=Frac(s) ' discard integer part
End Function


Extended: 0,000114 seconds
Double: 0,000118 seconds
Single: 0,000118 seconds

I'm saying: they run in equal time this way. So extended precision would be the best choice because it gives you 8 or maybe 9 hex digits you can trust in the same time the others take for less hex digits.

But now for my next trick: speeding up the main digit skipping loop...

[Johannes]

I changed the program a little bit, so that it will print the elapsed time.

I compiled and ran it three times, using Pelles C.

________id_____seconds
___100,000_______0.660
_1,000,000_______6.707
10,000,000______72.231

You can see that the time increases faster than linearly.

I've converted the main "digit skipping" routine into assembler. I really mis my sixteen all-purpose registers in the ARM. But you work with the tools you have, and on the x86 among those tools is a 64-by-32-bit integer division with remainder. As Neo said: "Whoa."

____id (n)_____seconds
___100,000_______0.321
_1,000,000_______3.435
10,000,000______35.595

As you can see my times are roughly half those of Dan's tests, showing the same more-than-linear increase. That's most likely because of the eponentation routine which is not linear (n log2 n?) and takes up the bulk of the time. For a proper comparison Dan would have to time the program on his computer but I'm guessing a full assembler version would run about 1.5 times as fast as the C version.



'----------------------------------------------------------------------------
' !PI Hex Digit
'
' Calculate nth hexadecimal digit in PI
'
' http://en.wikipedia.org/wiki/BBP_algorithm
'
' 22, 23 december 2010
'----------------------------------------------------------------------------

Uses "Oxygen"

DWord n = 10000000
' Number of hexadecimal digits to skip. Do NOT enter a value larger than
' 536870903. I'd have to look at my assembler code very carefully to see
' what would happen If you did enter a value larger than that. The results,
' if any were given, would probably be meaningless.

Ext hd,s,t,u,v
DWord a,i,k,m,p,q
String digits
PROCassemble

v=1/(16^17)

HiResTimer_Init
HiResTimer_Get
a=1 : BBPpart : hd+=4*s
a=4 : BBPpart : hd-=2*s
a=5 : BBPpart : hd-=s
a=6 : BBPpart : hd-=s
Quad x=HiResTimer_Delta

For i=1 To 8 ' More then 9 hexadecimal digits will *definitely* give incorrect results.
hd*=16
digits+=Hex$(Int(hd),1)
Next i

MsgBox 0,"Hexadecimal PI digits after #"+TStr$(n)+": "+digits+$CRLF+$CRLF+"Calculation time: "+TStr$(x/1000000)+" seconds"

Sub BBPpart()
s=0
k=0
Do While k<n
O2_EXEC
s+=Frac(q/a)
a+=8
k+=1
Loop
s=Frac(s)
u=1
Do
t=u/a
s+=t
u/=16
a+=8
Loop Until t<v
s=Frac(s)
End Function

' Copy of the above subroutine, with comments. Above subroutine compacted for speed.
' Main "skipping" routine is still in basic.
'Sub BBPpart(a As Quad) ' calculate BBP term
's=0 ' running summation
'k=0
'Do While k<n ' skips when n=0
' p=n-k ' exponent for 16
' q=1 ' running 16^(n-k) mod 8k+1
' r=Mod(16,a) ' mod is allowed everywhere when multiplying
' While p ' binary exponentation
' If (p And 1) Then
' q=Mod(q*r,a) ' multiply with running factor if bit is on
' EndIf
' r=Mod(r*r,a) ' running exponentation factor
' p=Int(p/2) ' process all bits
' Wend
' s+=Frac(q/a) ' divide by 8k+1
' a+=8 ' running 8k+1
' k+=1 ' next k
'Loop
's=Frac(s) ' discard integer part
'u=1 ' running 16^(n-k), k is implied in routine
'Do
' t=u/a ' next summation element
' s+=t ' running summation
' u/=16 ' running power of 16
' a+=8 ' running 8k+1
'Loop Until t<v ' get best precision
's=Frac(s) ' discard integer part
'End Function

Sub PROCassemble()
String src = "

.powermod
pushad ; save all registers

; --------------------- p=n-k
mov ebp,[#n] ; n
mov edi,[#k] ; k
sub ebp,edi ; n-k (never zero: k<n)

; --------------------- q=1
mov edi,1

; --------------------- r=Mod(16,a)
mov ebx,[#a]
mov edx,0
mov eax,16 ; 16
div ebx
mov ecx,edx ; 16 mod a

; --------------------- While p
.powermod1

; --------------------- If (p And 1) Then q=Mod(q*r,a)
mov eax,1
and eax,ebp ; lowest bit set?
jz powermod2
mov eax,edi ; q
mul ecx ; q * r
div ebx
mov edi,edx ; q * r mod a
.powermod2

; --------------------- r=Mod(r*r,a)
mov eax,ecx ; r
mul ecx ; r * r
div ebx
mov ecx,edx ; r * r mod a

; --------------------- p=Int(p/2)
shr ebp,1 ; logical shift right

; --------------------- Wend
cmp ebp,0 ; repeat until all '1' bits processed
jne powermod1

mov [#q],edi ; return q

popad ; restore all registers
ret ; return

"

O2_ASMO src
If Len(O2_ERROR) Then
MsgBox 0, O2_ERROR
Stop
EndIf
End Sub

[danbaron]

(Obviously, I guess, all of these times, are from using my (desktop, and only) machine, which is nothing special. The operating system is Windows 7; I think it is the, "junk", edition (version? - whatever Microsoft calls it).)

:p



Pelles C

n = 100,000 time = 0.660 seconds
n = 1,000,000 time = 6.707 seconds
n = 10,000,000 time = 72.231 seconds

Newest version by Johannes (run from thinAir, using thinBasic 1.8.0)

n = 100,000 time = 1.729 seconds
n = 1,000,000 time = 13.473 seconds
n = 10,000,000 time = 140.402 seconds

Newest version by Johannes (run from thinAir, using thinBasic 1.8.6.0)

n = 100,000 time = 1.274 seconds
n = 1,000,000 time = 11.800 seconds
n = 10,000,000 time = 122.654 seconds

Newest version by Johannes (run by double-clicking on file, using thinBasic 1.8.6.0)

n = 100,000 time = 1.150 seconds
n = 1,000,000 time = 12.215 seconds
n = 10,000,000 time = 124.430 seconds

gcc main.c -O3 -o main

n = 100,000 time = 0.441 seconds
n = 1,000,000 time = 5.304 seconds
n = 10,000,000 time = 61.266 seconds

Bloodshed Dev-C++

n = 100,000 time = 0.655 seconds
n = 1,000,000 time = 7.261 seconds
n = 10,000,000 time = 82.921 seconds

[zak]

i have astonished with the new enhancement of speed using oxygen with thinbasic , the results using Johannes recent version using oxygen:
for n=100,000 time 0.837
for n=1000,000 time 8.846
for n=10,000,000 time 92.98

but i prefer the other programming languages, i can't understand assembly. also there are another speedy language "qb64" it is an enhanced quick basic 45 and can run the the old code, and it is like bcx converting basic code to c then to exe, i will try it.
making the pi nth digit in basic language will make the algorithm available for a wide variety of people, it is like that the algorithm was in prison available only for the mathematicians, and now liberated to see the sun at the christmas of 2010.

[Johannes]

Newest version by Johannes (run from thinAir, using thinBasic 1.8.6.0)

n = 100,000 time = 1.274 seconds
n = 1,000,000 time = 11.800 seconds
n = 10,000,000 time = 122.654 seconds

gcc main.c -O3 -o main

n = 100,000 time = 0.441 seconds
n = 1,000,000 time = 5.304 seconds
n = 10,000,000 time = 61.266 seconds


That pisses me off. ;)

But I should have expected this; even in the first test the interpreted thinBasic driving loop is executed one hundred thousand times. Well, time to see how the floating-point opcodes work. I have to get my program at least twice (for n=100,000 three times) as fast to beat GCC. Should have known that the open source compiler would be the best one.

But first: clearing the (fresh) snow from the driveway so I can go out and stock up.

[Johannes]

but i prefer the other programming languages, i can't understand assembly.
Zak, I only did the bit in assembler because there speed is of the utmost importance. If you're doing something 10 million times any speed improvement in the main loop will have great advantages.

I understand you not having any use for assembly (starting from scratch in assembly is very hard these days because the Windows environment is so complicated) but I feel that way about C. Whenever I see C code I wonder why they didn't use assembler. (Although I refuse to program in C I can read it. I once needed to write a CRC routine in assembler and the only example I had was in C. I learned to read C then and decided to never program in it.) The only advantage C has over assembler is evaluated calculations, print statements and file I/O, but Oxygen takes care of that very nicely as well.

As with Acorn's BASIC V, for me thinBasic is the ideal scripting/programming tool as it provides gigantic functionality (*) in a proper structured programming language, and it has the option of doing stuff in assembler, should I need the speed.

(*) I recently rewrote in thinBasic a random Opera e-mail signature program I had on the Acorn. thinBasic's string handling is the best I've seen so far. I ended up with one-third of the code I needed in BASIC V. I'll try to remember to put that little program on the forums as well.

Oh, right. Driveway and snow. :p

[Johannes]

Never again. :)

That stack system for the floating-point registers is near-impossible to work with. I'm certain that my assembler code could be optimised but I'm declaring victory and I'm getting the hell out.

After some hard, hard work and a lot of debugging the execution times on my system are as follows.

__number_of______pure______partly_______fully____speed
_____digits_____BASIC___assembler___assembler___factor
____100,000_____6.762_______0.321_______0.074_____4.34
__1,000,000____80.050_______3.435_______0.916_____3.75
_10,000,000___992.752______35.595______10.927_____3.26
100,000,000_______n/a_____384.095_____128.562_____2.99

As you can see there is a significant speed increase when the surrounding program logic is also written in assembler. Also, for higher values of n the speed increase becomes less significant. This is logical because the bulk of the execution time is taken up by the skipping of the digits, which is the part that is already written in assembler in the second column, and cannot be sped up anymore.

Anyway, I'll admit there's a reason to program in a higher-level language but Oxygen Basic would be my choice then, certainly not C. And for a final comment: the speed increase from pure thinBasic to pure assembler for 1,000,000 digits is only a factor of 87. Truly a magnificent achievement from the thinBasic developers. My hat is off to them.

Merry Christmas.


'----------------------------------------------------------------------------
' !PI Hex Digit
'
' Calculate nth hexadecimal digit in PI
'
' http://en.wikipedia.org/wiki/BBP_algorithm
'
' 22-25 december 2010
'----------------------------------------------------------------------------

Uses "Oxygen"

DWord n = 1000000

DWord a,i,q
String d = "00000000"
Word c16=16
If n>536870903 Then MsgBox 0,"I'm afraid I can't do that.":Stop
PROCassemble

HiResTimer_Init
HiResTimer_Get
O2_EXEC
Quad x=HiResTimer_Delta

MsgBox 0,"Hexadecimal PI digits after #"+TStr$(n)+": "+d+$CRLF+$CRLF+"Calculation time: "+TStr$(x/1000000)+" seconds"



' Original thinBasic code.

' v=1/(16^17)

' a=1 : BBPpart : hd+=4*s
' a=4 : BBPpart : hd-=2*s
' a=5 : BBPpart : hd-=s
' a=6 : BBPpart : hd-=s

' For i=1 To 8 ' More then 9 hexadecimal digits will *definitely* give incorrect results.
' hd*=16
' digits+=Hex$(Int(hd),1)
' Next i

' Sub BBPpart(a As Quad) ' calculate BBP term
' s=0 ' running summation
' k=0
' Do While k<n ' skips when n=0
' p=n-k ' exponent for 16
' q=1 ' running 16^(n-k) mod 8k+1
' r=Mod(16,a) ' mod is allowed everywhere when multiplying
' While p ' binary exponentation
' If (p And 1) Then
' q=Mod(q*r,a) ' multiply with running factor if bit is on
' EndIf
' r=Mod(r*r,a) ' running exponentation factor
' p=Int(p/2) ' process all bits
' Wend
' s+=Frac(q/a) ' divide by 8k+1
' a+=8 ' running 8k+1
' k+=1 ' next k
' Loop
' s=Frac(s) ' discard integer part
' u=1 ' running 16^(n-k), k is implied in routine
' Do
' t=u/a ' next summation element
' s+=t ' running summation
' u/=16 ' running power of 16
' a+=8 ' running 8k+1
' Loop Until t<v ' get best precision
' s=Frac(s) ' discard integer part
' End Function



Sub PROCassemble()
String src = "

;--------------------------------------------------------------------------------
;
; Calculate hexadecimal PI digits after nth hexadecimal digit
;
.BBP
pushad

mov ebp,[#n] ; n

; --------------------- a=1 : BBPpart : hd+=4*s
fld1 ; make sure hd never becomes negative FP stack: hd
mov ebx,1
call BBPpart ; FP stack: s hd
faddp ; FP stack: hd
fadd st(0),st(0) ; 2 * s FP stack: 2hd
fadd st(0),st(0) ; 4 * s FP stack: 4hd

; --------------------- a=4 : BBPpart : hd-=2*s
mov ebx,4
call BBPpart ; FP stack: s hd
fadd st(0),st(0) ; 2 * s FP stack: 2s hd
fsubp ; - 2 * s FP stack: hd-2s

; --------------------- a=5 : BBPpart : hd-=s
mov ebx,5
call BBPpart ; FP stack: s hd
fsubp ; - s FP stack: hd-s

; --------------------- a=6 : BBPpart : hd-=s
mov ebx,6
call BBPpart ; FP stack: s hd
fsubp ; - s FP stack: hd-s

; --------------------- For i=1 To 8
mov ecx,[#d]
mov edx,0

fld st(0) ; copy of hd FP stack: hd hd
fisttp dword [#i] ; pop & store int(hd) FP stack: hd
fild dword [#i] ; int(hd) FP stack: int(hd) hd
fsubp ; frac(hd) FP stack: frac(hd)
.hexdigits

; --------------------- hd*=16
fimul word [#c16] ; * 16 FP stack: 16hd

; --------------------- digits+=Hex$(Int(hd),1)
fld st(0) ; copy of hd FP stack: hd hd
fisttp dword [#i] ; pop & store int(hd) FP stack: hd
fild dword [#i] ; int(hd) FP stack: int(hd) hd
fsubp ; frac(hd) FP stack: frac(hd)
mov eax,[#i] ; convert to ASCII
add eax,48
cmp eax,58
jl hexdigits1
add eax,7
.hexdigits1
mov [ecx+edx],al ; store character

; --------------------- Next i
inc edx
cmp edx,8
jne hexdigits

popad
ret



;--------------------------------------------------------------------------------
;
; Calculate term of BBP function
;
.BBPpart
pushad ; save all registers

; --------------------- s=0
fldz ; 0 FP stack: s=0

; --------------------- If n>0 Then
cmp ebp,0 ; skip any hex digits at all?
je calcdigits

; --------------------- k=0
mov esi,0 ; 0

; --------------------- Do
.skipdigits
call powermod ; q = 16^(n-k) mod a

; --------------------- s+=q/a
mov [#q],edi
mov [#a],ebx
fild dword [#q] ; q FP stack: q s
fidiv dword [#a] ; q / a FP stack: q/a s
faddp ; s += q / a FP stack: s+q/a
fld st(0) ; copy of s FP stack: s s
fisttp dword [#i] ; pop & store int(s) FP stack: s
fild dword [#i] ; int(s) FP stack: int(s) s
fsubp ; frac(s) FP stack: frac(s)

; --------------------- a+=8
add ebx,8 ; a + 8

; --------------------- k+=1
inc esi ; k + 1
cmp esi,ebp ; k<n?

; --------------------- Loop Until k=n
jl skipdigits

; --------------------- EndIf
.calcdigits

; --------------------- u=1
mov [#a],ebx
fld1 ; 1 FP stack: u=1 s
mov ecx,16 ; maximum precision is 64 bits = 16 nibbles

; --------------------- Do
.digitcalc

; --------------------- t=u/a
fld st(0) ; u FP stack: u u s
fidiv dword [#a] ; u / a FP stack: t=u/a u s

; --------------------- s+=t
faddp st(2),st(0) ; FP stack: u s+t

; --------------------- u/=16
fidiv word [#c16] ; u / 16 FP stack: u/16 s

; --------------------- a+=8
add [#a],8 ; a + 8

; --------------------- Loop Until t<v
dec ecx ; use counter instead
jnz digitcalc

; --------------------- s=Frac(s)
fcomp ; remove u from stack FP stack: s
fld st(0) ; copy of s FP stack: s s
fisttp dword [#i] ; pop & store int(s) FP stack: s
fild dword [#i] ; int(s) FP stack: int(s)
fsubp ; frac(s) FP stack: frac(s)

popad ; restore all registers
ret ; return



;--------------------------------------------------------------------------------
;
; Calculate 16^(n-k) mod a
;
.powermod
push eax ; save registers
push ecx
push edx
push ebp

; --------------------- p=n-k
sub ebp,esi ; n-k

; --------------------- q=1
mov edi,1

; --------------------- r=Mod(16,a)
mov edx,0
mov eax,16 ; 16
div ebx
mov ecx,edx ; 16 mod a

; --------------------- While p
.powermod1

; --------------------- If (p And 1) Then q=Mod(q*r,a)
mov eax,1
and eax,ebp ; lowest bit set?
jz powermod2
mov eax,edi ; q
mul ecx ; q * r
div ebx
mov edi,edx ; q * r mod a
.powermod2

; --------------------- r=Mod(r*r,a)
mov eax,ecx ; r
mul ecx ; r * r
div ebx
mov ecx,edx ; r * r mod a

; --------------------- p=Int(p/2)
shr ebp,1 ; logical shift right

; --------------------- Wend
cmp ebp,0 ; repeat until all '1' bits processed
jne powermod1

pop ebp ; restore registers
pop edx
pop ecx
pop eax
ret ; return

"

O2_ASMO src
If Len(O2_ERROR) Then
MsgBox 0, O2_ERROR
Stop
EndIf
End Sub