with trace
-- HexMath.e
-- A series of routines for performing integer math on sequences of up to (2^32)-1
-- bytes.  Integer math, for purposes of this library, includes only positive
-- numbers (and zero), and the following operations: sum, diff, mult, power, xor,
-- concat, mod, rdiv, rand, gt, lt, eq.
-- All functions are returned with byte[1] as most significant.
-- Parameters passed as sequences are assumed to have the most significant bytes first.
-- The machine.e function int_to_bytes returns least significant bytes first.

-- Written by Michael J. Sabal

include get.e
include machine.e
include misc.e

constant DEBUG = 0
atom DEBUG_LOG
--DEBUG_LOG = open("hexmath.log","w")

constant HALF_MAX = power(2,31)
--constant LOG_2 = log(2)
sequence SQUARES
SQUARES = repeat(0,255)
for i = 1 to 255 do
  SQUARES[i] = i * i
end for

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

global function hex_trim(object a)
  if atom(a) then
    return a
  end if
  for ctr = 1 to length(a) do
    if sequence(a[ctr]) or a[ctr]!=0 then
      return a[ctr..length(a)]
    end if
  end for
  return {0}
end function

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

global function hex_eq(object a, object b)
  atom la,lb
  if atom(a) and atom(b) and a=b then
    return 1
  end if
  if atom(a) then
    a = reverse(int_to_bytes(a))
  end if
  if atom(b) then
    b = reverse(int_to_bytes(b))
  end if
  la = length(a)
  lb = length(b)
  if la < lb then
    a = repeat(0,lb-la)&a
    la = lb
  end if
  if lb < la then
    b = repeat(0,la-lb)&b
    lb = la
  end if
  if compare(a,b)=0 then
    return 1
  end if
  return 0
end function

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

global function hex_lt(object a, object b)
  atom la,lb
  if atom(a) and atom(b) and a<b then
    return 1
  end if
  if atom(a) then
    a = reverse(int_to_bytes(a))
  end if
  if atom(b) then
    b = reverse(int_to_bytes(b))
  end if
  la = length(a)
  lb = length(b)
  if la < lb then
    a = repeat(0,lb-la)&a
    la = lb
  end if
  if lb < la then
    b = repeat(0,la-lb)&b
    lb = la
  end if
  if compare(a,b)<0 then
    return 1
  end if
  return 0
end function

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

global function hex_gt(object a, object b)
  atom la,lb
  if atom(a) and atom(b) and a>b then
    return 1
  end if
  if atom(a) then
    a = reverse(int_to_bytes(a))
  end if
  if atom(b) then
    b = reverse(int_to_bytes(b))
  end if
  la = length(a)
  lb = length(b)
  if la < lb then
    a = repeat(0,lb-la)&a
    la = lb
  end if
  if lb < la then
    b = repeat(0,la-lb)&b
    lb = la
  end if
  if compare(a,b)>0 then
    return 1
  end if
  return 0
end function

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

global function hex_sum(object a, object b)

  atom la, lb, s, tim
  sequence carry, rtn

  if atom(a) and atom(b) and a < HALF_MAX and b < HALF_MAX then
    return reverse(int_to_bytes(a+b))
  end if
  if atom(a) then
    a = reverse(int_to_bytes(a))
  end if
  if atom(b) then
    b = reverse(int_to_bytes(b))
  end if
  la = length(a)
  lb = length(b)
  if la < lb then
    a = repeat(0,lb-la)&a
    la = lb
  end if
  if lb < la then
    b = repeat(0,la-lb)&b
    lb = la
  end if
  carry = repeat(0,la+1)
  rtn = repeat(0,la)
  tim = time()
  for ctr = lb to 1 by -1 do
    s = carry[ctr+1]+a[ctr]+b[ctr]
    while s > #FF do
      s = s - #100
      carry[ctr] += 1
    end while
    rtn[ctr] = s
  end for
  if carry[1] > 0 then
    rtn = carry[1] & rtn
  end if
  if DEBUG then printf(DEBUG_LOG,"Sum time: %f\n", time()-tim) end if

  return rtn

end function

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

global function hex_diff(object a, object b)

  atom la, lb
  sequence rtn

  if atom(a) and atom(b) then
    if a > b then
      return reverse(int_to_bytes(a-b))
    else
      return reverse(int_to_bytes(b-a))
    end if
  end if
  if atom(a) then
    a = reverse(int_to_bytes(a))
  end if
  if atom(b) then
    b = reverse(int_to_bytes(b))
  end if
  la = length(a)
  lb = length(b)
  if la < lb then
    a = repeat(0,lb-la)&a
    la = lb
  end if
  if lb < la then
    b = repeat(0,la-lb)&b
    lb = la
  end if
  if compare(a,b)<0 then
    -- Hex_diff only returns a positive integer
    rtn = a
    a = b
    b = rtn
  end if
  rtn = repeat(0,la)
  for ctr = lb to 1 by -1 do
    if a[ctr] < b[ctr] then
      a[ctr-1] -= 1
      a[ctr] += #100
    end if
    rtn[ctr] = a[ctr] - b[ctr]
  end for

  return rtn

end function

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

global function hex_xor(object a, object b)

  atom la,lb
  sequence rtn

  if atom(a) and atom(b) then
    return xor_bits(a,b)
  end if
  if atom(a) then
    a = reverse(int_to_bytes(a))
  end if
  if atom(b) then
    b = reverse(int_to_bytes(b))
  end if
  la = length(a)
  lb = length(b)
  if la < lb then
    a = repeat(0,lb-la)&a
    la = lb
  end if
  if lb < la then
    b = repeat(0,la-lb)&b
    lb = la
  end if
  rtn = repeat(0,la)
  for ctr = 1 to la do
    rtn[ctr] = xor_bits(a[ctr],b[ctr])
  end for
  return rtn

end function

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

global function hex_concat(object a, object b)
  if atom(a) then
    a = reverse(int_to_bytes(a))
  end if
  if atom(b) then
    b = reverse(int_to_bytes(b))
  end if
  return a&b 
end function

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

global function hex_rdiv_trad(object a, object b)

  -- This is the traditional repeated subtraction algorithm.

  atom la, lb
  sequence rem
  object div

  div = 0
  rem = {0}
  if atom(a) then
    a = reverse(int_to_bytes(a))
  end if
  if atom(b) then
    b = reverse(int_to_bytes(b))
  end if
  la = length(a)
  lb = length(b)
  if la < lb then
    a = repeat(0,lb-la)&a
    la = lb
  end if
  if lb < la then
    b = repeat(0,la-lb)&b
    lb = la
  end if

  while not hex_lt(a,b) do
    a = hex_diff(a,b)
    if atom(div) and div < 2100000000 then
      div += 1
    else  -- Only use expensive hex_sum when necessary
      div = hex_sum(div,1)
    end if
  end while
  if atom(div) then
    div = reverse(int_to_bytes(div))
  end if
  rem = a

  return {div,rem}

end function

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

global function hex_rdiv(object a, object b)

  -- Standard long division algorithm

  sequence accum, reg
  object div
  atom aptr, la, lb

  if atom(a) and atom(b) then
    return {floor(a/b),remainder(a,b)}
  end if
  div = {}
  accum = {0}
  if atom(a) then
    a = reverse(int_to_bytes(a))
  end if
  if atom(b) then
    b = reverse(int_to_bytes(b))
  end if
  a = hex_trim(a)
  b = hex_trim(b)
  la = length(a)
  lb = length(b)
  if lb > la then
    return {{0},a}
  elsif la = lb and hex_lt(a,b) then
    return {{0},a}
  end if
  aptr = 1
  accum = {}
  while aptr <= la do
    accum = accum & a[aptr]
    aptr = aptr + 1
    if length(accum) > lb or hex_lt(b,accum) then
      reg = hex_rdiv_trad(accum,b)
      div = div & hex_trim(reg[1])
      accum = hex_trim(reg[2])
    else
      div = div & #00
    end if
  end while
  return {hex_trim(div),accum}

end function

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

global function hex_mod(object a, object b)
  
  sequence rtn
  rtn = hex_rdiv(a,b)
  return rtn[2]

end function

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

global function hex_square(object a)

  atom la, tim
  sequence accum
  accum = {0}
  la = length(a)
  tim = time()
  for i = 1 to la do
    if a[i]!=0 then
      if i = la then
        accum = hex_sum(accum,{SQUARES[a[i]]})
      else
        if length(accum) <= (2*(la-i)) then
          accum = repeat(0,(2*(la-i))-length(accum)+2) & accum
        end if
        accum[length(accum)-(2*(la-i))] += SQUARES[a[i]] -- We'll resolve carries at the end
      end if
      for j = i+1 to la do
          accum[length(accum)-((la-i)+(la-j))] += 2*a[i]*a[j]
--        accum = hex_sum(accum,{2*a[i]*a[j]}&repeat(0,(la-i)+(la-j)))
      end for
    end if
  end for
  if DEBUG then printf(DEBUG_LOG,"Square time: %f\n", time()-tim) end if
  return hex_trim(accum)

end function

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

global function hex_mult4(object a, object b)

  -- a * b = (floor(a/b)*(b^2)) + (b*remainder(a/b))

  atom la,lb,p,tim
  sequence temp, sqq, squ, prsq, prrem, ret

  if atom(a) then
    a = reverse(int_to_bytes(a))
  end if
  if atom(b) then
    b = reverse(int_to_bytes(b))
  end if
  la = length(a)
  lb = length(b)
  if hex_lt(a,b) then
    -- a must be > b
    temp = a
    a = b
    b = temp
  end if

  tim = time()
  -- Calculate the square of b and pull it out.
  sqq = hex_rdiv(a,b)
  if hex_gt(sqq[1],0) then
    squ = hex_square(b)
  else
    squ = {0}
  end if

  -- [a/b]*b^2
  la = length(sqq[1])
  lb = length(squ)
  prsq = repeat(0,(la+lb))
  for i = lb to 1 by -1 do
    for j = la to 1 by -1 do
      p = squ[i] * sqq[1][j]
      prsq[i+j] += p
    end for
  end for

  -- Calculate the product of b and the remainder
  la = length(sqq[2])
  lb = length(b)
  prrem = repeat(0,(la+lb))
  for i = lb to 1 by -1 do
    for j = la to 1 by -1 do
      p = b[i] * sqq[2][j]
      prrem[i+j] += p
    end for
  end for

  -- Add the two together
  ret = hex_sum(prsq,prrem)
  if DEBUG then printf(DEBUG_LOG,"Mult4 time: %f\n", time()-tim) end if
  return hex_trim(ret)

end function

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

global function hex_mult(object a, object b)

  atom la, lb, p, tim
  sequence rtn

  if atom(a) then
    a = reverse(int_to_bytes(a))
  end if
  if atom(b) then
    b = reverse(int_to_bytes(b))
  end if
  la = length(a)
  lb = length(b)

  tim = time()
  rtn = repeat(0,(la+lb))
  for i = lb to 1 by -1 do
    for j = la to 1 by -1 do
      p = b[i] * a[j]
      rtn[i+j] += p
    end for
  end for

  -- Clean up overflows
  for k = la+lb to 2 by -1 do
    p = floor(rtn[k]/#100)
    rtn[k] = remainder(rtn[k],#100)
    rtn[k-1] += p
  end for
  if rtn[1] > #100 then
    rtn = floor(rtn[1]/#100) & rtn
    rtn[2] = remainder(rtn[2],#100)
  end if
--  if DEBUG then printf(DEBUG_LOG,"Mult time: %f\n", time()-tim) end if

  return hex_trim(rtn)

end function

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

global function hex_power(object a, object b)

  -- See http://www.daimi.au.dk/~ivan/FastExpproject.pdf for a comparison
  -- of algorithms.  This function implements left-to-right binary exp.

  sequence rtn, bb
  atom tim

  tim = time()
  rtn = {1}
  if atom(a) then
    a = reverse(int_to_bytes(a))
  end if
  if atom(b) then
    b = reverse(int_to_bytes(b))
  end if
  bb = {}
  for k = 1 to length(b) do
    bb = bb & reverse(int_to_bits(b[k],8))
  end for
  for i = 1 to length(bb) do
    rtn = hex_trim(hex_square(rtn))
    if bb[i] then
      rtn = hex_trim(hex_mult(rtn,a))
    end if
  end for
  if DEBUG then printf(DEBUG_LOG,"Power time: %f\n", time()-tim) end if
  return rtn

end function

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

global function hex_power_mod(object a, object b, object m)

  -- See http://www.daimi.au.dk/~ivan/FastExpproject.pdf for a comparison
  -- of algorithms.  This function implements left-to-right binary exp.

  -- Adding the modulo check after each round, 128-byte a ^ 65537
  -- mod 128-byte m takes about 5 seconds on a P4. Smaller m = faster
  -- results.

  sequence rtn, bb
  atom tim

  tim = time()
  rtn = {1}
  if atom(a) then
    a = hex_trim(reverse(int_to_bytes(a)))
  end if
  if atom(b) then
    b = hex_trim(reverse(int_to_bytes(b)))
  end if
  if atom(m) then
    m = hex_trim(reverse(int_to_bytes(m)))
  end if
  bb = {}
  for k = 1 to length(b) do
    bb = bb & reverse(int_to_bits(b[k],8))
  end for
  for i = 1 to length(bb) do
    rtn = hex_trim(hex_mult(rtn,rtn))
    if hex_gt(rtn,m) then
      rtn = hex_mod(rtn,m)
    end if
    if bb[i] then
      rtn = hex_trim(hex_mult(rtn,a))
      rtn = hex_mod(rtn,m)
    end if
  end for
  if DEBUG then printf(DEBUG_LOG,"Powermod time: %f\n", time()-tim) end if
  return hex_mod(rtn,m)

end function

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

global function hex_rand(atom num_bytes)
  sequence rtn
  rtn = repeat(0,num_bytes)
  for ctr = 1 to num_bytes do
    rtn[ctr] = rand(256)-1
    if rtn[ctr] > 255 then
      rtn[ctr] = 255
    end if
    if rtn[ctr] < 0 then
      rtn[ctr] = 0
    end if
  end for
  return rtn
end function

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

global function hex_gcd(object a, object b)
  -- Uses the Euclidean algorithm
  -- http://en.wikipedia.org/wiki/Euclidean_algorithm

  object t
  if atom(a) then
    a = reverse(int_to_bytes(a))
  end if
  while not hex_eq(b,0) do
    t = b
    b = hex_mod(a,b)
    a = t
  end while
  return a

end function

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

global function hex_lcm(object a, object b)

  object gcd

  gcd = hex_gcd(a,b)
  return hex_rdiv(hex_mult(a,b),gcd)

end function

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

procedure hex_test()

  atom tim
  sequence a, b, x, m, n

  a = hex_rand(90)
  b = hex_rand(80)
if 0 then
  ? a
  ? b
  tim = time()
  x = hex_rdiv(a,b)
  ? x
  ? length(x[1])
  ? length(x[2])
  ? time()-tim
  m = hex_mult(b,x[1])
  n = hex_sum(m,x[2])
  if hex_eq(n,a) then
    puts(1,"Sanity test PASS\n")
  else
    puts(1,"Sanity test FAIL\n")
    ? a
    ? n
  end if
else
  tim = time()
--  a = {15}
--  x = hex_power(a,{1,0,5})
--  ? x
--  ? a
--  ? length(x)
--  ? hex_square(repeat(#FF,500))
  ? hex_power_mod(hex_rand(128),{1,0,1},hex_rand(32))
  ? time()-tim
  printf(DEBUG_LOG,"Total time: %f\n", time()-tim)
end if

end procedure
--hex_test()
--close(DEBUG_LOG)
--if wait_key() then end if