
'uses minwin
uses console

type hpmtype
  byte b[64]
  int sg,pp
end type

dim hpmtype aa,bb,cc,rr
'
dim int a[1024] '256 bytes
dim int expon
dim int negat
'
'byte offsets
'A     0  .. 63
'B     64 .. 127
'R    128 .. 255
'T    256 .. 319


def hpmInit
===========
  addr edi,a
  add  edi,60
  mov  cl,16
  mov  ch,16
  clc
end def

def hpmAdd
==========
  hpmBitShiftAddExpon
  hpmInit
  pushf
  (
    mov eax,[edi] 'A
    popf
    adc eax,[edi+64] 'B
    pushf
    mov [edi+192],eax 'R
    sub edi,4
    dec cl
    jnz repeat
  )
  popf
  adc [edi+192],0 'final carry
end def

def hpmSub
==========
  hpmBitShiftSubExpon
  hpmInit
  pushf
  (
    mov eax,[edi] 'A
    popf
    sbb eax,[edi+64] 'B
    pushf
    mov [edi+192],eax 'R
    sub edi,4
    dec cl
    jnz repeat
  )
  popf
  sbb [edi+192],0 'final carry
end def

def hpmCompare
==============
  addr edi,a 'dwords left to right
  mov cl,16
  xor edx,edx
  (
    mov eax,[edi] 'A
    cmp eax,[edi+64] 'B
    jb fwd ncmp
    add edi,4
    dec cl
    jnz repeat
  )
  inc edx
  ncmp:
end def

def hpmBitShiftLeft
===================
  clc
  pushf
  (
    popf
    rcl dword [edi],1
    pushf
    sub edi,4
    dec cl
    jnz repeat
  )
  popf
end def

def hpmBitShiftRight
====================
  clc
  pushf
  (
    popf
    rcr dword [edi],1
    pushf
    add edi,4
    dec cl
    jnz repeat
  )
  popf
end def

def hpmBitShiftAddExpon
=======================
  mov edx,0
  mov ecx,512
  addr edi,a
  add edi,60
  push edi
   (
    mov al,1
    and al,[edi]
    jnz exit
    mov al,1
    and al,[edi+64]
    jnz exit
    inc edx
    push ecx
    push edi
    hpmBitShiftRightA
    hpmBitShiftRightB
    pop edi
    pop ecx
    dec ecx
    jnz repeat
  )
  mov expon,edx
  pop edi
end def

def hpmBitShiftSubExpon
=======================
hpmBitShiftAddExpon
end def

def hpmBitShiftDivExpon
=======================
  mov edx,0
  mov ecx,512
  addr edi,a
  push edi
   (
    mov al,0x80
    and al,[edi+3]
    jnz exit
    dec edx
    push ecx
    push edi
    hpmBitShiftLeftA
    pop edi
    pop ecx
    dec ecx
    jnz repeat
  )
  addr edi,a
  add edi,64
  mov ecx,512
   (
    mov al,0x80
    and al,[edi+3]
    jnz exit
    inc edx
    push ecx
    push edi
    hpmBitShiftLeftB
    pop edi
    pop ecx
    dec ecx
    jnz repeat
  )
  mov expon,edx
  pop edi
end def

def hpmBitShiftMulExpon
=======================
  mov edx,0
  mov ecx,512
  addr edi,a
  add edi,60
  push edi
   (
    mov al,1
    and al,[edi]
    jnz exit
    inc edx
    push ecx
    push edi
    hpmBitShiftRightA
    pop edi
    pop ecx
    dec ecx
    jnz repeat
  )
  addr edi,a
  add edi,124
  mov ecx,512
   (
    mov al,1
    and al,[edi]
    jnz exit
    inc edx
    push ecx
    push edi
    hpmBitShiftRightB
    pop edi
    pop ecx
    dec ecx
    jnz repeat
  )
  mov expon,edx
  pop edi
end def

def hpmBitShiftLeftA
====================
  addr edi,a
  add edi,60
  mov cl,16
  hpmBitShiftLeft
end def

def hpmBitShiftLeftB
====================
  addr edi,a
  add edi,124
  mov cl,16
  hpmBitShiftLeft
end def

def hpmBitShiftLeftR
====================
  addr edi,a
  add edi,252
  mov cl,32
  hpmBitShiftLeft
end def

def hpmBitShiftRightA
=====================
  addr edi,a
  mov cl,16
  hpmBitShiftRight
end def

def hpmBitShiftRightB
=====================
  addr edi,a
  add edi,64
  mov cl,16
  hpmBitShiftRight
end def

def hpmBitShiftRightR
=====================
  addr edi,a
  add edi,128
  mov cl,32
  hpmBitShiftRight
end def

def hpmBitRotateLeft
====================
push edi
mov eax,ecx
dec eax
shl eax,2
add edi,eax 'lowest dword
push edi
clc
pushf
(
  popf
  rcl [edi]
  pushf
  sub edi,4
  dec ecx
  jnz repeat
)
popf
pop edi
(
  jnc exit
  or byte [edi],0x1 'first bit
)
pop edi
end def

def hpmBitRotateRight
=====================
push edi 'highest dword
clc
pushf
(
  popf
  rcr [edi]
  pushf
  add edi,4
  dec ecx
  jnz repeat
)
popf
pop edi
(
  jnc exit
  or byte [edi+3],0x80 'final bit
)
end def

def hpmBitRotateRightB
======================
  addr edi,a
  add edi,64
  mov ecx,16
  hpmBitRotateRight
end def

def hpmAddAR
============
  hpmInit
  pushf
  (
    mov eax,[edi] 'A
    popf
    adc [edi+192],eax 'R+=A
    pushf
    sub edi,4
    dec cl
    jnz repeat
  )
  popf
  adc [edi+192],0 'final carry
end def

def hpmBitRotateRightR
======================
  addr edi,a
  add edi,128
  mov ecx,32
  hpmBitRotateRight
end def

def hpmMulBits
==============
hpmClearR
addr edi,a
mov esi, 512 'bits in multiplier b
(
  (
    addr edi,a
    test byte [edi+124],1 'lowest bit of B
    jz exit
    hpmAddAR
  )
  hpmBitRotateRightR
  hpmBitRotateRightB
  dec esi
  jnz repeat  
)
'swap upper and lower parts of R
addr edi,a
add edi,128
mov cl,16
(
  mov eax,[edi]
  mov edx,[edi+64]
  mov [edi],edx
  mov [edi+64],eax
  add edi,4
  dec cl
  jnz repeat
)
addr edi,a
end def

def hpmMul
==========
int ccl,cch
hpmClearR
hpmBitShiftMulExpon
hpminit
mov ecx,16 'dword displacement
addr esi,a
add esi,64 'address offset of muliplier
mov dword cch,16
(
  xor edx,edx
  mov dword ccl,16
  push edi
  (
    mov eax,[edi] 'digit
    add eax,edx 'carry from prior mul
    mul [esi]
    add [edi+ecx*4+192],eax 'add to accum
    adc [edi+ecx*4+188],0   'carry bit
    sub edi,4             'next binary digit
    dec dword ccl
    jnz repeat
  )
  pop edi 'restore position
  dec ecx
  sub esi,4
  dec dword cch
  jnz repeat
)
end def

def hpmClearR
=============
addr edi,a
add edi,128
mov ecx,32
xor edx,edx
(
  mov [edi],edx
  add edi,4
  dec ecx
  jnz repeat
)
end def

def hpmClearABR
===============
addr edi,a
mov ecx,128 '512 bytes
xor edx,edx
(
  mov [edi],edx
  add edi,4
  dec ecx
  jnz repeat
)
end def

def hpmCopy64
mov cl,16
(
  mov eax,[esi]
  mov [edi],eax
  add edi,4
  add esi,4
  dec cl
  jnz repeat
)
end def

def hpmLoadA
============
addr edi,a
addr esi,%1
hpmCopy64
end def

def hpmLoadB
============
addr edi,a
add edi,64
addr esi,%1
hpmCopy64
end def

def hpmSave
===========
addr esi,a
add esi,192
addr edi,%1
hpmCopy64
end def

def hpmCmpZeroR
===============
addr edi,a
add edi,128
mov cl,32
(
  cmp dword [edi],0
  jnz exit
  add edi,4
  dec cl
  jnz repeat
)
addr edi,a
end def

def hpmCmpR
===========
addr edi,a
add edi,128
mov eax,square
(
  cmp [edi],eax
  jnz exit
  add edi,4
  mov cl,31
  (
    cmp dword [edi],0
    jnz exit
    add edi,4
    dec cl
    jnz repeat
  )
)
addr edi,a
end def

def hpmSqrt2
============
hpmClearABR
addr edi,a
int sqriter=512
int square=0x80000000
hpmSqrt
end def

def hpmSqrt5
============
hpmClearABR
addr edi,a
int sqriter=512
int square=0x50000000
hpmSqrt
end def

def hpmSqrt
===========
mov dword [edi],    0x80000000 'init A
mov dword [edi+64], 0x80000000 'init B
mov dword [edi+256],0x80000000 'init T
mov esi,0
mov ecx,edi
(
  push ecx
  hpmMulBits
  hpmCmpR
  pop ecx
  'jz exit 
  (
    jb exit
    mov eax,[edi+256]
    xor [ecx],eax 'clear prior bit in A
    xor [ecx+64],eax 'clear prior bit in B
  )
  ror dword [edi+256],1
  (
    cmp [edi+256],0x80000000 'roll over
    jnz exit
    add ecx,4 'next dword for bit setting
  )
  mov eax,[edi+256]
  or [ecx],eax 'set bit in A
  or [ecx+64],eax 'set bit in B  
  dec sqriter
  jle exit
  repeat
)
end def

def hpmShowA
============
scope
int i
for i=16 to 16
  print mid(hex(aa[i],8),-8) cr
next
print cr
end scope
end def

def hpmShowB
============
scope
int i
for i=16 to 16
  print mid(hex(bb[i],8),-8) cr
next
print cr
end scope
end def

def hpmShowR
============
scope
int i
for i=32 to 32
  print mid(hex(rr[i],8),-8) cr
next
print cr
end scope
end def

indexbase 1
dword nDivBitsA=512
dword nDivBitsB=512
dword nDivisor32=16
dword aa at @a
dword bb at @a+64
dword rr at @a+128

'hpmAdd
'hpmSub
'hpmMul
'hpmMulBits
'hpmSqrt2
'hpmCmp
'hpmClearR
'hpmClearABR
'hpmLoadA aa
'hpmLoadB bb
'hpmSave rr

'int nDivisor32=1
'hpmDiv

extern lib kernel32.dll
int GetTickCount()
end extern

/*
'TIME TESTS
===========
int i,t1,t2
t1=GetTickCount
for i=1 to 1000 'empty loop 15 nanoseconds
 'hpmAdd '250 nanoseconds (bitshiftaddexpon 260 microseconds)
 'hpmMulBits '406 microseconds
 'hpmMul '265 microseconds
 'hpmDiv '14 milliseconds
next
'hpmSqrt2 '234 milliseconds
t2=GetTickCount
print t2-t1
waitkey
*/

/*
'CALC TESTS
===========
aa=0xe0000000
bb=0xe0000000
'
'aa[16]=0x4444
'bb[16]=0x3000
'
'rr[32]=0x7
'hpmBitRotateRightB
'hpmAdd
'hpmSub
'hpmMul
'hpmMulBits

hpmSqrt5
hpmshowA

print cr cr
print mid( hex(aa,8),-8)  cr
print hex(sqr(5)*0x40000000) cr
'
print aa/0x40000000 cr
print sqr(5) cr
print cr
hpmclearR
hpmmulbits
hpmshowR
*/

def hpmCompareDiv32
======================
  addr edi,a 'bytes left to right
  mov ecx,nDivisor32
  (
    mov eax,[edi+64]
    cmp [edi],eax
    jnz exit
    add edi,4
    dec ecx
    jg repeat
  )
end def

def hpmSubDiv32
===============
  addr edi,a 'bytes right to left
  mov ecx,nDivisor32
  dec ecx 'base 0
  shl ecx,2 'mult by 4
  add edi,ecx
  clc
  pushf
  (
    mov eax,[edi+64] 'B
    popf
    sbb [edi],eax 'A-=B
    pushf
    sub edi,4
    sub ecx,4
    jge repeat
  )
  popf
end def


def hpmDivBits
==============
hpmClearR
hpmInit
mov esi,nBitsA
sub esi,nBitsB
(
  hpmCompareDiv32
  (
   jb exit
   hpmSubDiv32
   addr edi,a
   or byte [edi+252],1 'set lowest bit in R
  )
  dec esi
  jl exit
  hpmBitShiftLeftA
  hpmBitShiftLeftR
  repeat
)
end def

def hpmRol
==========
addr edi,%1
mov esi,%2
(
  dec esi
  jl exit
  mov ecx,16
  hpmBitRotateLeft
  repeat
)
end def

def hpmRor
==========
addr edi,%1
mov esi,%2
(
  dec esi
  jl exit
  mov ecx,16
  hpmBitRotateRight
  repeat
)
end def

hpmClearABR
aa[16]=987654321
bb[16]=10
int nBitsA=31
int nBitsB=4
inc nbitsA
inc nbitsB
int i
nDivisor32=1
hpmRor aa nBitsA
hpmRor bb nBitsB
for i=1 to 10
 hpmDivBits
 hpmRol aa nBitsB 'remainder digit
 print a[16] 'lowest to highest digit
 aa=rr[32] 'repeat divide result
 rr[32]=0
next
print cr cr
hpmRol aa nBitsB
hpmRol bb nBitsB

'hpmShowA
'hpmShowB
'hpmShowR
print aa[16] cr cr
print bb[16] cr cr
print rr[32] cr cr

waitkey

