Binary Square Root Calculation (large Numbers)

Share your advanced PureBasic knowledge/code with the community.
Helle
Enthusiast
Enthusiast
Posts: 178
Joined: Wed Apr 12, 2006 7:59 pm
Location: Germany
Contact:

Binary Square Root Calculation (large Numbers)

Post by Helle »

Just I make tests with older math-codes from me and found this nice (I think) code:

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
User avatar
idle
Always Here
Always Here
Posts: 5915
Joined: Fri Sep 21, 2007 5:52 am
Location: New Zealand

Re: Binary Square Root Calculation (large Numbers)

Post by idle »

interesting, thanks
Post Reply