Binary Square Root Calculation (large Numbers)
Posted: Sun Dec 06, 2009 3:55 pm
Just I make tests with older math-codes from me and found this nice (I think) code:
Have Fun!
Helle
Code: Select all
;- binary calculation of square root of (pos.) Numbers, without optimization
;- "Helle" Klaus Helbing, 06.12.2009, PB4.40
;for procedures, or another way
Global P$ ;product-string
Global S$ ;sum-string
Procedure.s Mul2(A$)
p_A = @A$ ;pointer of A$
l_A = Len(A$) ;length of A$ (or SSE-routine ;-) )
If PeekB(p_A) > $34 ;greater then 4 -> carry
l_P = l_A
Else
l_P = l_A - 1
EndIf
P$ = Space(l_P + 1)
While l_A
l_A - 1
A = PeekB(p_A + l_A) & $F ;or -$30
A + A + C
If A > 9
C = 1 : A - 10
Else
C = 0
EndIf
PokeB(@P$ + l_P, A + $30)
l_P - 1
Wend
If C
PokeB(@P$ + l_P, $31) ;(last) set carry (cipher 1)
EndIf
EndProcedure
Procedure.s AddStr(A$, B$) ;A$ is >= B$!
p_A = @A$ ;pointer
l_A = Len(A$) ;length
p_B = @B$ ;pointer
l_B = Len(B$) ;length
S$ = Space(l_A + 1) ;sum, with 1 place more, remove space later (if exist)
p_S = @S$
While l_A
l_A - 1
A = PeekB(p_A + l_A) & $F ;or -$30
If l_B > 0
l_B - 1
B = PeekB(p_B + l_B) & $F ;or -$30
Else
B = 0
EndIf
A + B + C
If A > 9
C = 1 : A - 10
Else
C = 0
EndIf
PokeB(p_S + l_A + 1, A + $30)
Wend
If C
PokeB(p_S + l_A, $31) ;(last) set carry (cipher 1)
Else
S$ = LTrim(S$)
EndIf
EndProcedure
Global Freq.q ;4 variables for time-measuring
Global Start.q
Global Ende.q
Global Time.d
;no test of plausibility!
EingabeO$ = InputRequester("Calculation Square Root - Input radicand", "e.g. 12345 or 12345.6789 or 0.00012345", "")
NachKommaO = Val(InputRequester("Number of places after decimal point", "pos. integer-value (e.g. 1000) :", ""))
Eingabe$ = EingabeO$ ;save the original input-string (for output)
NachKomma = Len(Eingabe$) << 3 ;if places after decimal point too short
If NachKomma < NachKommaO
NachKomma = NachKommaO
EndIf
;SetThreadAffinityMask_(GetCurrentThread_(), 1) ;for test, set the thread to Core0 and read the time from Core0
QueryPerformanceFrequency_(@Freq)
QueryPerformanceCounter_(@Start)
L = Len(Eingabe$)
For i = 1 To L
If Mid(Eingabe$, i, 1) = "." ;test for decimal point
Eingabe$ = Mid(Eingabe$, 1, i - 1) + Mid(Eingabe$, i + 1, L - i) ;remove decimal point, if exist
Korr = L - i
L - 1
If Korr & 1 ;odd number of after-decimal-point-places
Eingabe$ + "0" ;make even with "0"
Korr + 1
EndIf
Korr >> 1
For j = i To L
If Mid(Eingabe$, j, 1) <> "0" ;test for zero
Break
EndIf
NullNK + 1
Next
Break
EndIf
Next
NullNK >> 1
If j = L + 1
NullNK = 0 ;without interest, only zeros after decimal point
EndIf
X0 = NachKomma << 2
X1 = (NachKomma << 2) - 1
X2 = (NachKomma << 2) - 2
X4 = (NachKomma << 2) - 4
EingabeVoll$ = Eingabe$
For i = 1 To NachKomma
EingabeVoll$ + "00"
Next
Dec$ = EingabeVoll$
DecBin$ = Space(Len(Dec$) << 2) ;calculate the length of the binary-string: 1x dec = 4x bin (avoid risks)
l_DecBin = Len(DecBin$)
p_DecBin = @DecBin$ + l_DecBin
;Phase 1: Convert the decimal-input-string into a binary-string
Repeat
Mul2(Dec$) ;no division with 2 (time!)
Mul2(P$) ;multiply with 5 with place-adaptation
AddStr(P$, Dec$) ;now equal like a multiplication with 5
p_DecBin - 1
If Mid(S$, Len(S$), 1) = "5"
PokeB(p_DecBin, $31)
Else
PokeB(p_DecBin, $30)
EndIf
Dec$ = Mid(S$, 1, Len(S$) - 1)
If Dec$ = "1"
p_DecBin - 1
PokeB(p_DecBin, $31)
Break
ElseIf Dec$ = "0"
p_DecBin - 1
PokeB(p_DecBin, $30)
Break
EndIf
ForEver
DecBin$ = LTrim(DecBin$)
l_DecBin = Len(DecBin$) ;length of binary-string
k = 4 - L_DecBin % 4
For i = 1 To k
Vornullen$ + "0"
l_DecBin + 1
Next
DecBin$ = Vornullen$ + DecBin$
DecBin$ = LTrim(DecBin$)
p_DecBin = @DecBin$
Binaer$ = Space(X0)
Binaer1$ = Space(X0)
Rest$ = Space(X0)
p_Binaer = @Binaer$
l_Binaer = Len(Binaer$)
p_Binaer1 = @Binaer1$
l_Binaer1 = Len(Binaer1$)
p_Rest = @Rest$
l_Rest = Len(Rest$)
For i = 0 To X4 Step 4 ;fill the strings with "0"
PokeL(p_Binaer + i, $30303030)
PokeL(p_Rest + i, $30303030)
Next
;Phase 2: Calculate binary the square root
Iter = (l_DecBin >> 1) - 1 ;number of iterations
For j = 0 To Iter
For i = 2 To X2 Step 2 ;equal shift left 2 of Rest$
PokeW(p_Rest + i - 2, PeekW(p_Rest + i))
Next
PokeW(p_Rest + l_Rest - 2, PeekW(p_DecBin + ZE))
ZE + 2
l_A = j + 3
For i = 0 To l_A ;copy Binaer$ in Binaer1$
PokeB(p_Binaer1 + l_Binaer1 - l_A + i, PeekB(p_Binaer + l_Binaer -l_A + i))
Next
For i = 2 To X2 Step 2 ;equal shift left 2 of Binaer1$
PokeW(p_Binaer1 + i - 2, PeekW(p_Binaer1 + i))
Next
PokeW(p_Binaer1 + i - 2, $3030) ;set last 2 places to "0"
l_A + 1 ;length
p_RestX = p_Rest + l_Rest - l_A
p_Binaer1X = p_Binaer1 + l_Binaer1 - l_A
Signum = p_Rest + l_Rest - j - 4
If PeekB(Signum) = $30 ;now binary subtraction
PokeB(p_Binaer1 + l_Binaer1 - 1, $31) ;set last place to "1"
Carry = 0
Komp2 = 1 ;for the 2-complement
While l_A ;subtraction = addition of the 2-complement
l_A - 1
A = PeekB(p_RestX + l_A) & $F ;or -$30
If l_A > 0
B = PeekB(p_Binaer1X + l_A) & $F ;or -$30
B = (~B) & 1
B + Komp2
Select B
Case 0, 1
Komp2 = 0
Case 2
B = 0 : Komp2 = 1
Case 3
B = 1 : Komp2 = 1
EndSelect
Else
B = 1
EndIf
A + B + Carry
Select A
Case 0, 1
Carry = 0
Case 2
A = 0 : Carry = 1
Case 3
A = 1 : Carry = 1
EndSelect
PokeB(p_RestX + l_A, A + $30)
Wend
Else ;now binary addition
PokeW(p_Binaer1 + l_Binaer1 - 2, $3131) ;set last 2 places to "1"
Carry = 0
While l_A
l_A - 1
A = PeekB(p_RestX + l_A) & $F ;or -$30
If l_A > 0
B = PeekB(p_Binaer1X + l_A) & $F ;or -$30
Else
B = 0
EndIf
A + B + Carry
Select A
Case 0, 1
Carry = 0
Case 2
A = 0 : Carry = 1
Case 3
A = 1 : Carry = 1
EndSelect
PokeB(p_RestX + l_A, A + $30)
Wend
EndIf
For i = 1 To X1 ;equal shift left 1 of Binaer$
PokeB(p_Binaer + i - 1, PeekB(p_Binaer + i))
Next
PokeB(p_Binaer + i - 1, $30) ;set last place to zero
If PeekB(Signum) = $30
PokeB(p_Binaer + l_Binaer - 1, $31)
EndIf
Next
;Phase 3: Convert the binary-string into the output-decimal-string
If Mid(Binaer$, Len(Binaer$), 1) = "0" ;set LSB from hand
S$ = "0"
Else
S$ = "1"
EndIf
P$ = "1" ;start-value for power-of-2-calculation
l_Bin = Len(Binaer$) - 1
p_Bin = @Binaer$
While L_Bin
Mul2(P$) ;calculate the powers of 2
l_Bin - 1
If PeekB(p_Bin + l_Bin) > $30
AddStr(P$, S$) ;add, if the bit is set
EndIf
Wend
;Output, set decimal point
Vor$ = Mid(S$, 1, Len(S$) - NachKomma - Korr)
If Vor$ = "" ;without this is the output e.g. ".123"
Vor$ = "0"
Else
NullNK = 0 ;a (non-zero) value in front of the decimal point exist
EndIf
Erg$ = Vor$ + "."
For i = 1 To NullNK
Erg$ + "0"
Next
Nach$ = Mid(S$, Len(S$) - NachKomma + 1 - Korr, NachKommaO - NullNK)
Erg$ + Nach$
QueryPerformanceCounter_(@Ende)
Time = (Ende - Start) / Freq
MessageRequester("Square Root of " + EingabeO$ + " with " + Str(NachKommaO) + " places after decimal point in " + StrD(Time, 3)+" s", Erg$)
Have Fun!
Helle