crypto = require 'crypto'
{BigDecimal} = require 'bigdecimal'
{BigInteger} = require 'bigdecimal'
Order Preserving Symmetric Encryption (OPSE)
Encrypts numeric values in such a way that their original value is confidential, but inequalities can be performed on numbers encrypted under the same key.
crypto = require 'crypto'
{BigDecimal} = require 'bigdecimal'
{BigInteger} = require 'bigdecimal'
Encrypt a number.
key
is a random symmetric cipher key, in the form of
caesar.key.createRandom()
.num
is a 16 bit number to encrypt.exports.encrypt = (key, num) ->
top = 131072
cursor = 65536
bottom = 0
while true # Narrow search down.
k = rhyper key, 65536, 65536, cursor
if Math.abs(num - k) < 20 then break
if num > k then bottom = cursor
if num < k then top = cursor
if cursor % 4 is 1
cursor = Math.floor((top + bottom) / 2)
else if cursor % 4 is 3
cursor = Math.ceil((top + bottom) / 2)
else
cursor = (top + bottom) / 2
ocursor = cursor + 20 # Get as accurate as possible.
if k > num then cursor -= 80
if cursor < 0 then cursor = 0
bcursor = cursor
okay = false
while not okay and cursor <= ocursor
cursor = cursor + 2
k = rhyper key, 65536, 65536, cursor
if k is num then okay = true
if not okay
while bcursor <= ocursor
bcursor = bcursor + 2
k = rhyper key, 65536, 65536, bcursor
if k > num
bcursor -= 1
break
cursor = bcursor
cursor
The following five functions have been ported from the R Project for Statistical Computing and appropriately modified. Don’t ask me how or why they work.
afc = (i) ->
al =
0: '0' # ln(0!)=ln(1)
1: '0' # ln(1!)=ln(1)
2: '0.69314718055994530941723212145817' # ln(2)
3: '1.79175946922805500081247735838070' # ln(6)
4: '3.17805383034794561964694160129705' # ln(24)
5: '4.78749174278204599424770093452324'
6: '6.57925121201010099506017829290394'
7: '8.52516136106541430016553103634712'
8: '10.60460290274525022841722740072165'
if i < 0 then throw 'i less than 0. Should not happen.'
if i <= 8 then return al[i]
pi = 0.5 * Math.log(2 * Math.PI)
(i+0.5) * Math.log(i) - i + (1/12) / i - (1/360) / i / i / i + pi
imax2 = (x, y) -> if x < y then y else x
imin2 = (x, y) -> if x < y then x else y
class PRNG
constructor: (@coin) ->
@cipher = crypto.createCipher 'aes-256-ctr', @coin
@blank = new Buffer 16
@blank.fill 0
draw: ->
@cipher.write @blank
out = @cipher.read()
numer = new BigInteger out.toString('hex'), 16
numer = new BigDecimal numer.toString()
denom = new BigDecimal '340282366920938463463374607431768211456'
parseFloat numer.divide(denom, 100, BigDecimal.ROUND_HALF_UP).toString()
rhyper = (coin, nn1, nn2, kk) ->
prng = new PRNG coin
con = 57.56462733
deltal = 0.0078
deltau = 0.0034
ks = -1
n1s = -1
n2s = -1
scale = 1e25
if nn1 is Infinity or nn2 is Infinity or kk is Infinity then throw 'NaN'
nn1 = Math.floor(nn1 + 0.5)
nn2 = Math.floor(nn2 + 0.5)
kk = Math.floor(kk + 0.5)
if nn1 < 0 or nn2 < 0 or kk < 0 or kk > (nn1 + nn2) then throw 'NaN'
if nn1 isnt n1s or nn2 isnt n2s
setup1 = true
setup2 = true
else if kk isnt ks
setup1 = false
setup2 = true
else
setup1 = false
setup2 = false
if setup1
n1s = nn1
n2s = nn2
tn = nn1 + nn2
if nn1 <= nn2
n1 = nn1
n2 = nn2
else
n1 = nn2
n2 = nn1
if setup2
ks = kk
tk = kk + kk
if (kk + kk) >= tn
k = tn - kk
else
k = kk
if setup1 or setup2
m = (k + 1) * (n1 + 1) / (tn + 2)
minjx = imax2 0, (k - n2)
maxjx = imin2 n1, k
Generate random variate — three basic cases.
if minjx is maxjx # I : Degenerate distribution.
ix = maxjx
else if (m - minjx) < 10 # II : Inverse transformation.
if setup1 or setup2
if k < n2
w = Math.exp(con + afc(n2) + afc(n1 + n2 - k) - afc(n2 - k) -
afc(n1 + n2))
else
w = Math.exp(con + afc(n1) + afc(k) - afc(k-n2) - afc(n1+n2))
l10 = (w, minjx) -> # L10
p = w
ix = minjx
u = prng.draw() * 1e25
[p, ix, u]
[p, ix, u] = l10 w, minjx
while u > p # L20
u -= p
p *= (n1 - ix) * (k - ix)
++ix
p = p / ix / (n2 - k + ix)
if ix > maxjx then [p, ix, u] = l10 w, minjx
else # III : h2pe
if setup1 or setup2
s = Math.sqrt((tn - k) * k * n1 * n2 / (tn - 1) / tn / tn)
d = Math.floor((1.5 * s) + 0.5)
xl = m - d + 0.5
xr = m + d + 0.5
a = afc(m) + afc(n1 - m) + afc(k - m) + afc(n2 - k + m)
kl = Math.exp(a - afc(Math.floor(xl)) - afc(Math.floor(n1 - xl)) -
afc(Math.floor(k - xl)) - afc(Math.floor(n2 - k + xl)))
kr = Math.exp(a - afc(Math.floor(xr - 1)) -
afc(Math.floor(n1 - xr + 1)) - afc(Math.floor(k - xr + 1)) -
afc(Math.floor(n2 - k + xr - 1)))
lamdl = -Math.log(xl * (n2 - k + xl) / (n1 - xl + 1) / (k - xl + 1))
lamdr = -Math.log((n1 - xr + 1) * (k - xr + 1) / xr / (n2 - k + xr))
p1 = d + d
p2 = p1 + kl / lamdl
p3 = p2 + kr / lamdr
while true # L30
u = prng.draw() * p3
v = prng.draw()
if u < p1
ix = Math.floor(xl + u)
else if (u <= p2)
ix = Math.floor(xl + Math.log(v) / lamdl)
if ix < minjx then continue
v = v * (u - p1) * lamdl
else
ix = Math.floor(xr - Math.log(v) / lamdr)
if ix > maxjx then continue
v = v * (u - p2) * lamdr
Acceptance/Rejection Test
if m < 100 or ix <= 50
f = 1
if m < ix
i = m + 1
while i <= ix
f = f * (n1 - i + 1) * (k - i + 1) / (n2 - k + i) / i
++i
else if m > ix
i = ix + 1
while i <= m
f = f * i * (n2 - k + i) / (n1 - i + 1) / (k - i + 1)
++i
if v <= f then break
else
Squeeze using upper and lower bounds.
y = ix
y1 = y + 1
ym = y - m
yn = n1 - y + 1
yk = k - y + 1
nk = n2 - k + y1
r = -ym / y1
s = ym / yn
t = ym / yk
e = -ym / nk
g = yn * yk / (y1 * nk) - 1
dg = 1
if g < 0 then dg = 1 + g
gu = g * (1 + g * (-0.5 + g / 3))
gl = gu - .25 * (g * g * g * g) / dg
xm = m + 0.5
xn = n1 - m + 0.5
xk = k - m + 0.5
nm = n2 - k + xm
ub = y * gu - m * gl + deltau +
xm * r * (1 + r * (-0.5 + r / 3)) +
xn * s * (1 + s * (-0.5 + s / 3)) +
xk * t * (1 + t * (-0.5 + t / 3)) +
nm * e * (1 + e * (-0.5 + e / 3))
Test against upper bound.
alv = Math.log v
if alv > ub then continue
Test against lower bound.
dr = xm * (r * r * r * r)
if r < 0 then dr /= (1 + r)
ds = xn * (s * s * s * s)
if s < 0 then ds /= (1 + s)
dt = xk * (t * t * t * t)
if t < 0 then dt /= (1 + t)
de = nm * (e * e * e * e)
if e < 0 then de /= (1 + e)
cand = ub - 0.25 * (dr+ds+dt+de) + (y + m) * (gl - gu) - deltal
if alv < cand then break
else
Stirling’s formula to machine accuracy
cand = (a - afc(ix) - afc(n1-ix) - afc(k-ix) - afc(n2-k+ix))
if alv <= cand then break else continue
Return appropriate variate.
if (kk + kk) >= tn
if nn1 > nn2
ix = kk - nn2 + ix
else ix = nn1 - ix
else
if nn1 > nn2
ix = kk - ix
ix