
uses console

function bhmMul(dword *a,*b,*c,q)
=================================
dword d,i,j,k,ca,cy,cyy,oo,oq,m100=1E8
dword aa at @a
dword bb at @b
dword cc at @c
j=1
if q>1 then
  j=q-1
  @aa=@aa+j*4-4
end if
for i=j to q+4
  aa=0 : @aa+=4 'partial clearing accum
next
oo=q*4-4    'shifting offset
oq=oo       'fixed offset
sys bq=@a+j*4-8
sys pq=@a+17*4 'max resolution + 4
for i=1 to q
  @cc=@c+oo  'multiplier
  @bb=@b+oq  'value
  @aa=@a+q*4+oo 'accum
  cyy=0         'clear accum carry
  cy=0          'clear multiply carry
  for j=1 to q
    'if @aa<bq then exit for
    if @aa<=pq 'ignore digits beyond resolution
      'k=bb*cc+cy : cy=0
      'if k>=m100 then
      '  cy=k\m100
      '  k=edx 'remainder after cpu division
      'end if
      addr edi,bb
      mov eax,[edi]
      addr edi,cc
      mul [edi]
      mov ecx,m100
      div ecx
      mov d,eax   'mul carry
      add edx,cy  'input carry
       mov eax,edx 'remainder
      cmp eax,ecx
      (
        jb exit
        inc d
        sub eax,ecx
      )
      mov k,eax
      cy=d
      k=k+aa+cyy : cyy=0 'accum add with carry
      if k>=m100 then k-=m100 : cyy=1
      aa=k
    end if '<=pq
    @aa-=4
    @bb-=4
  next
  aa=cyy+cy 'final carries
  oo-=4 'shift left for next line mul
next
end function


function bhmCompare(dword *b,*c,i,*j=0) as int
==============================================
dword bb at @b
dword cc at @c
while i--
  if bb>cc then return 1
  if bb<cc then return -1
  @bb+=4
  @cc+=4
  j++
wend
end function

function bhmSub(dword *b,*c,i) as int
=====================================
int cy
int bb at @b
int cc at @c+(i*4)-4
while i--
  bb=bb-cc-cy
  if bb<0 then bb+=1e8 : cy=1 else cy=0
  @bb-=4
  @cc-=4
wend
bb-=cy 'final carry
end function

function bhmRotateLeft(dword*n,i)
=================================
dword nn at @n+(i*4)-4
dword ee at @nn
dword i,cy
while i--
  addr edi,nn
  mov eax,[edi]
  shl eax
  add eax,cy
  mov dword cy,0
  cmp eax,1e8
  (
   jle exit
   sub eax,1e8
   mov cy,1
  )
  mov [edi],eax
  @nn-=4
wend
ee+=cy 'final carry
end function

function bhmRotateRight(dword*n,i)
==================================
dword nn at @n
dword d,i,cy
while i--
  addr edi,nn
  mov eax,[edi]
  mov edx,0
  shr eax
  (
   jnc exit
   mov edx,5e7
  )
  add eax,cy
  mov dword cy,edx
  mov [edi],eax
  @nn+=4
wend
n+=cy 'final carry
end function

function bhmClear(dword *n,i)
=============================
dword nn at @n
while i--
  nn=0
  @nn+=4
wend
end function

function bhmDiv(dword *r,*n,*d,i,j)
===================================
dword bi=i*32 'to bits
dword bj=j*32
dword bc=bi-bj+1
int c
bhmClear(r,i)
while bc--
  c=bhmCompare(n,d,j)
  if c>=0 then
    c=1
    bhmSub(n,d,j)
  else
    c=0
  end if
  r or= c
  if bc then
    bhmRotateLeft(r,i)
    bhmRotateLeft(n,i)
  end if
wend
end function


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

/*
int nr[32],n1[32],n2[32]
'n1=14140000
'n2=14140000
'n1=22300002
'n2=22300002
'n1=87654321
'n2=2000000
bhmmul(nr,n1,n2,16)
for i=0 to 31
  print i tab mid("0000000"+str(nr[i]),-8) cr
next
print cr
*/


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

'square root of 2
'using binary step approximator
indexbase 0
int i
for i=0 to 31 : nr[i]=0 : next 'clear nr
nt[0]=02000000   'target value representing 2
int b,c,i,j,m
int n at @nr 'strider for square root digits
b=0x4000000 '67,108,864 shifting bit setter
i=16
j=1
while i
  n or= b 'set bit
  'm=nm[j-1]
  bhmMul nm,nr,nr,j
  c=bhmCompare(nt,nm,j+1,0)
  if c<=0 then 'too large
    n xor= b  'clear bit
    'nm[j-1]=m 'undo carry
  end if
  shr b 'next bit down
  if b=0 then
    i--
    b=0x4000000
    @n+=4 'next 100M down
    j++
  end if
wend

'display digits in groups of eight
print "square root of 2" cr
for i=0 to 15
  print i tab mid("0000000"+str(nr[i]),-8)
  print   tab mid("0000000"+str(nm[i]),-8) cr
next
print cr
'  1.4142135623730950488016887242096980785696718753769480731766797379907324784621


'square root of 5 --> phi
'using binary step approximator
indexbase 0
for i=0 to 31 : nr[i]=0 : next 'clear nr
nt[0]=05000000 'target value representing 5
int b,c,i
int n at @nr 'strider for square root digits
b=0x4000000 '67,108,864 shifting bit setter
i=16
while i
  n or= b 'set bit
  bhmMul nm,nr,nr,16
  c=bhmCompare(nt,nm,32)
  if c<0 then 'too large
    n xor= b 'clear bit
  end if
  shr b 'next bit down
  if b=0 then
    i--
    b=0x4000000
    @n+=4 'next 10k down
  end if
wend

'display digits in groups of eight
print "square root of 5" cr
for i=0 to 15
  print i tab mid("0000000"+str(nr[i]),-8) cr
next
'2.2360679774997896964091736687312762354406183596115257242708972454105209256378



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


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


'print 0x4000000 cr '67,108,864
'print hex 67108864


waitkey
 

 
 


