• Jump To … +
    caesar.coffee commitment.coffee format.coffee hash.coffee key.coffee kts.coffee message.coffee opse.coffee ots.coffee searchable.coffee tree.coffee
  • opse.coffee

  • ¶

    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.

    1. key is a random symmetric cipher key, in the form of caesar.key.createRandom().
    2. 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