Fast Hilbert curve algorithm by http://threadlocalmutex.com/ Ported from C++ https://github.com/rawrunprotected/hilbert_curves (public domain)
335525
function hilbert(x, y) {
    let a = x ^ y;
    let b = 0xFFFF ^ a;
    let c = 0xFFFF ^ (x | y);
    let d = x & (y ^ 0xFFFF);

    let A = a | (b >> 1);
    let B = (a >> 1) ^ a;
    let C = ((c >> 1) ^ (b & (d >> 1))) ^ c;
    let D = ((a & (c >> 1)) ^ (d >> 1)) ^ d;

    a = A; b = B; c = C; d = D;
    A = ((a & (a >> 2)) ^ (b & (b >> 2)));
    B = ((a & (b >> 2)) ^ (b & ((a ^ b) >> 2)));
    C ^= ((a & (c >> 2)) ^ (b & (d >> 2)));
    D ^= ((b & (c >> 2)) ^ ((a ^ b) & (d >> 2)));

    a = A; b = B; c = C; d = D;
    A = ((a & (a >> 4)) ^ (b & (b >> 4)));
    B = ((a & (b >> 4)) ^ (b & ((a ^ b) >> 4)));
    C ^= ((a & (c >> 4)) ^ (b & (d >> 4)));
    D ^= ((b & (c >> 4)) ^ ((a ^ b) & (d >> 4)));

    a = A; b = B; c = C; d = D;
    C ^= ((a & (c >> 8)) ^ (b & (d >> 8)));
    D ^= ((b & (c >> 8)) ^ ((a ^ b) & (d >> 8)));

    a = C ^ (C >> 1);
    b = D ^ (D >> 1);

    let i0 = x ^ y;
    let i1 = b | (0xFFFF ^ (i0 | a));

    i0 = (i0 | (i0 << 8)) & 0x00FF00FF;
    i0 = (i0 | (i0 << 4)) & 0x0F0F0F0F;
    i0 = (i0 | (i0 << 2)) & 0x33333333;
    i0 = (i0 | (i0 << 1)) & 0x55555555;

    i1 = (i1 | (i1 << 8)) & 0x00FF00FF;
    i1 = (i1 | (i1 << 4)) & 0x0F0F0F0F;
    i1 = (i1 | (i1 << 2)) & 0x33333333;
    i1 = (i1 | (i1 << 1)) & 0x55555555;

    return ((i1 << 1) | i0) >>> 0;
}
560250
def pts_on_bezier_curve(P=[(0.0, 0.0)], n_seg=0):
    '''
    '''
    assert isinstance(P, list)
    assert len(P) > 0
    for p in P:
        assert isinstance(p, tuple)
        for i in p:
            assert len(p) > 1
            assert isinstance(i, float)
    assert isinstance(n_seg, int)
    assert n_seg >= 0
    
    return [pt_on_bezier_curve(P, float(i)/n_seg) for i in range(n_seg)] + [P[-1]]
931314
def bitwise_or(self, t):
        
        """
        This implementation combines the approaches used by 'WYSINWYX: what you see is not what you execute'
        paper and 'Signedness-Agnostic Program Analysis: Precise Integer Bounds for Low-Level Code'. The
        first paper provides an sound way to approximate the stride, whereas the second provides a way
        to calculate the or operation using wrapping intervals.
        Note that, even though according Warren's work 'Hacker's delight', one should follow different
        approaches to calculate the minimun/maximum values of an or operations according on the type
        of the operands (signed/unsigned). On the other other hand, by splitting the wrapping-intervals
        at the south pole, we can safely and soundly only use the Warren's functions for unsigned
        integers.
        """
        s = self
        result_interval = list()

        for u in s._ssplit():
            for v in t._ssplit():
                w = u.bits
                # u |w v
                if u.is_integer:
                    s_t = StridedInterval._ntz(v.stride)
                elif v.is_integer:
                    s_t = StridedInterval._ntz(u.stride)
                else:
                    s_t = min(StridedInterval._ntz(u.stride), StridedInterval._ntz(v.stride))

                if u.is_integer and u.lower_bound == 0:
                    new_stride = v.stride
                elif v.is_integer and v.lower_bound == 0:
                    new_stride = u.stride
                else:
                    new_stride = 2 ** s_t
                mask = (1 << s_t) - 1
                r = (u.lower_bound & mask) | (v.lower_bound & mask)
                m = (2 ** w) - 1
                low_bound = WarrenMethods.min_or(u.lower_bound & (~mask & m), u.upper_bound & (~mask & m), v.lower_bound & (~mask & m), v.upper_bound & (~mask & m), w)
                upper_bound = WarrenMethods.max_or(u.lower_bound & (~mask & m), u.upper_bound & (~mask & m), v.lower_bound & (~mask & m), v.upper_bound & (~mask & m), w)
                if low_bound == upper_bound:
                    new_stride = 0

                new_interval = StridedInterval(lower_bound=((low_bound & (~mask & m)) | r), upper_bound=((upper_bound & (~mask & m)) | r), bits=w, stride=new_stride)
                result_interval.append(new_interval)

        return StridedInterval.least_upper_bound(*result_interval).normalize()
1096760
function X9ECParameters(name) {
	var params = algs.curves[name];
	assert.object(params);

	var p = new jsbn(params.p);
	var a = new jsbn(params.a);
	var b = new jsbn(params.b);
	var n = new jsbn(params.n);
	var h = jsbn.ONE;
	var curve = new ec.ECCurveFp(p, a, b);
	var G = curve.decodePointHex(params.G.toString('hex'));

	this.curve = curve;
	this.g = G;
	this.n = n;
	this.h = h;
}
1341556
public int[] makeLut() {
        int numKnots = x.length;
        float[] nx = new float[numKnots + 2];
        float[] ny = new float[numKnots + 2];
        System.arraycopy(x, 0, nx, 1, numKnots);
        System.arraycopy(y, 0, ny, 1, numKnots);
        nx[0] = nx[1];
        ny[0] = ny[1];
        nx[numKnots + 1] = nx[numKnots];
        ny[numKnots + 1] = ny[numKnots];

        int[] table = new int[256];
        for (int i = 0; i < 1024; i++) {
            float f = i / 1024.0f;
            int x = (int) (255 * Curve.Spline(f, nx.length, nx) + 0.5f);
            int y = (int) (255 * Curve.Spline(f, nx.length, ny) + 0.5f);
            x = x > 255 ? 255 : x;
            x = x < 0 ? 0 : x;
            y = y > 255 ? 255 : y;
            y = y < 0 ? 0 : y;
            table[x] = y;
        }
        return table;
    }
552411
def examples(self):
        
        # TODO: adapt for tensor-valued functions

        # Get the points and calculate some statistics on them
        mins = self.domain.min()
        maxs = self.domain.max()
        means = (maxs + mins) / 2.0
        stds = (maxs - mins) / 4.0
        ndim = getattr(self.domain, 'ndim', 0)

        # Zero and One
        yield ('Zero', self.zero())
        yield ('One', self.one())

        # Indicator function in first dimension
        def step_fun(x):
            return (x[0] > means[0])

        yield ('Step', self.element(step_fun))

        # Indicator function on hypercube
        def cube_fun(x):
            result = True
            for points, mean, std in zip(x, means, stds):
                result = np.logical_and(result, points < mean + std)
                result = np.logical_and(result, points > mean - std)
            return result

        yield ('Cube', self.element(cube_fun))

        # Indicator function on a ball
        if ndim > 1:  # Only if ndim > 1, don't duplicate cube
            def ball_fun(x):
                r = sum((xi - mean) ** 2 / std ** 2
                        for xi, mean, std in zip(x, means, stds))
                return r < 1.0

            yield ('Ball', self.element(ball_fun))

        # Gaussian function
        def gaussian_fun(x):
            r2 = sum((xi - mean) ** 2 / (2 * std ** 2)
                     for xi, mean, std in zip(x, means, stds))
            return np.exp(-r2)

        yield ('Gaussian', self.element(gaussian_fun))

        # Gradient in each dimensions
        for axis in range(ndim):
            def gradient_fun(x):
                return (x[axis] - mins[axis]) / (maxs[axis] - mins[axis])

            yield ('Grad {}'.format(axis), self.element(gradient_fun))

        # Gradient in all dimensions
        if ndim > 1:  # Only if ndim > 1, don't duplicate grad 0
            def all_gradient_fun(x):
                return sum((xi - xmin) / (xmax - xmin)
                           for xi, xmin, xmax in zip(x, mins, maxs))

            yield ('Grad all', self.element(all_gradient_fun))
844847
def hilbertrot(n, x, y, rx, ry):
    
    if ry == 0:
        if rx == 1:
            x = n - 1 - x
            y = n - 1 - y
        return y, x
    return x, y
890878
def to_curve_spline(obj):
    '''
    
    '''
    if   is_curve_spline(obj):            return obj
    elif is_tuple(obj) and len(obj) == 2: (crds,opts) = obj
    else:                                 (crds,opts) = (obj,{})
    if pimms.is_matrix(crds) or is_curve_spline(crds): crds = [crds]
    spls = [c for c in crds if is_curve_spline(c)]
    opts = dict(opts)
    if 'weights' not in opts and len(spls) == len(crds):
        if all(c.weights is not None for c in crds):
            opts['weights'] = np.concatenate([c.weights for c in crds])
    if 'order' not in opts and len(spls) > 0:
        opts['order'] = np.min([c.order for c in spls])
    if 'smoothing' not in opts and len(spls) > 0:
        sm = set([c.smoothing for c in spls])
        if len(sm) == 1: opts['smoothing'] = list(sm)[0]
        else: opts['smoothing'] = None
    crds = [x.crds if is_curve_spline(crds) else np.asarray(x) for x in crds]
    crds = [x if x.shape[0] == 2 else x.T for x in crds]
    crds = np.hstack(crds)
    return curve_spline(crds, **opts)
1343269
function() {

            var start = this.start;
            var controlPoint1 = this.controlPoint1;
            var controlPoint2 = this.controlPoint2;
            var end = this.end;

            var x0 = start.x;
            var y0 = start.y;
            var x1 = controlPoint1.x;
            var y1 = controlPoint1.y;
            var x2 = controlPoint2.x;
            var y2 = controlPoint2.y;
            var x3 = end.x;
            var y3 = end.y;

            var points = new Array(); // local extremes
            var tvalues = new Array(); // t values of local extremes
            var bounds = [new Array(), new Array()];

            var a, b, c, t;
            var t1, t2;
            var b2ac, sqrtb2ac;

            for (var i = 0; i < 2; ++i) {

                if (i === 0) {
                    b = 6 * x0 - 12 * x1 + 6 * x2;
                    a = -3 * x0 + 9 * x1 - 9 * x2 + 3 * x3;
                    c = 3 * x1 - 3 * x0;

                } else {
                    b = 6 * y0 - 12 * y1 + 6 * y2;
                    a = -3 * y0 + 9 * y1 - 9 * y2 + 3 * y3;
                    c = 3 * y1 - 3 * y0;
                }

                if (abs(a) < 1e-12) { // Numerical robustness
                    if (abs(b) < 1e-12) { // Numerical robustness
                        continue;
                    }

                    t = -c / b;
                    if ((0 < t) && (t < 1)) tvalues.push(t);

                    continue;
                }

                b2ac = b * b - 4 * c * a;
                sqrtb2ac = sqrt(b2ac);

                if (b2ac < 0) continue;

                t1 = (-b + sqrtb2ac) / (2 * a);
                if ((0 < t1) && (t1 < 1)) tvalues.push(t1);

                t2 = (-b - sqrtb2ac) / (2 * a);
                if ((0 < t2) && (t2 < 1)) tvalues.push(t2);
            }

            var j = tvalues.length;
            var jlen = j;
            var mt;
            var x, y;

            while (j--) {
                t = tvalues[j];
                mt = 1 - t;

                x = (mt * mt * mt * x0) + (3 * mt * mt * t * x1) + (3 * mt * t * t * x2) + (t * t * t * x3);
                bounds[0][j] = x;

                y = (mt * mt * mt * y0) + (3 * mt * mt * t * y1) + (3 * mt * t * t * y2) + (t * t * t * y3);
                bounds[1][j] = y;

                points[j] = { X: x, Y: y };
            }

            tvalues[jlen] = 0;
            tvalues[jlen + 1] = 1;

            points[jlen] = { X: x0, Y: y0 };
            points[jlen + 1] = { X: x3, Y: y3 };

            bounds[0][jlen] = x0;
            bounds[1][jlen] = y0;

            bounds[0][jlen + 1] = x3;
            bounds[1][jlen + 1] = y3;

            tvalues.length = jlen + 2;
            bounds[0].length = jlen + 2;
            bounds[1].length = jlen + 2;
            points.length = jlen + 2;

            var left = min.apply(null, bounds[0]);
            var top = min.apply(null, bounds[1]);
            var right = max.apply(null, bounds[0]);
            var bottom = max.apply(null, bounds[1]);

            return new Rect(left, top, (right - left), (bottom - top));
        }
1371562
function fft2d(re, im) {
    var tre = [],
        tim = [],
        i = 0;
    // x-axis
    for (var y = 0; y < _n; y++) {
        i = y * _n;
        for (var x1 = 0; x1 < _n; x1++) {
            tre[x1] = re[x1 + i];
            tim[x1] = im[x1 + i];
        }
        fft1d(tre, tim);
        for (var x2 = 0; x2 < _n; x2++) {
            re[x2 + i] = tre[x2];
            im[x2 + i] = tim[x2];
        }
    }
    // y-axis
    for (var x = 0; x < _n; x++) {
        for (var y1 = 0; y1 < _n; y1++) {
            i = x + y1 * _n;
            tre[y1] = re[i];
            tim[y1] = im[i];
        }
        fft1d(tre, tim);
        for (var y2 = 0; y2 < _n; y2++) {
            i = x + y2 * _n;
            re[i] = tre[y2];
            im[i] = tim[y2];
        }
    }
}
1254711
function(prp, adata) {
    var i,
        times2 = sjcl.mode.ocb2._times2,
        w = sjcl.bitArray,
        xor = w._xor4,
        checksum = [0,0,0,0],
        delta = prp.encrypt([0,0,0,0]),
        bi;

    delta = xor(delta,times2(times2(delta)));

    for (i=0; i+4<adata.length; i+=4) {
      delta = times2(delta);
      checksum = xor(checksum, prp.encrypt(xor(delta, adata.slice(i,i+4))));
    }

    bi = adata.slice(i);
    if (w.bitLength(bi) < 128) {
      delta = xor(delta,times2(delta));
      bi = w.concat(bi,[0x80000000|0,0,0,0]);
    }
    checksum = xor(checksum, bi);
    return prp.encrypt(xor(times2(xor(delta,times2(delta))), checksum));
  }
282039
def read_turbomole(basis_lines, fname):
    '''
    '''

    skipchars = '*#$'
    basis_lines = [l for l in basis_lines if l and not l[0] in skipchars]

    bs_data = create_skel('component')

    i = 0
    while i < len(basis_lines):
        line = basis_lines[i]
        elementsym = line.split()[0]

        element_Z = lut.element_Z_from_sym(elementsym)
        element_Z = str(element_Z)

        if not element_Z in bs_data['elements']:
            bs_data['elements'][element_Z] = {}

        element_data = bs_data['elements'][element_Z]

        if "ecp" in line.lower():
            if not 'ecp_potentials' in element_data:
                element_data['ecp_potentials'] = []

            i += 1
            line = basis_lines[i]

            lsplt = line.split('=')
            maxam = int(lsplt[2])
            n_elec = int(lsplt[1].split()[0])

            amlist = [maxam]
            amlist.extend(list(range(0, maxam)))

            i += 1
            for shell_am in amlist:
                shell_am2 = lut.amchar_to_int(basis_lines[i][0])[0]
                if shell_am2 != shell_am:
                    raise RuntimeError("AM not in expected order?")

                i += 1

                ecp_shell = {
                    'ecp_type': 'scalar_ecp',
                    'angular_momentum': [shell_am],
                }
                ecp_exponents = []
                ecp_rexponents = []
                ecp_coefficients = []

                while i < len(basis_lines) and basis_lines[i][0].isalpha() is False:
                    lsplt = basis_lines[i].split()
                    ecp_exponents.append(lsplt[2])
                    ecp_rexponents.append(int(lsplt[1]))
                    ecp_coefficients.append(lsplt[0])
                    i += 1

                ecp_shell['r_exponents'] = ecp_rexponents
                ecp_shell['gaussian_exponents'] = ecp_exponents
                ecp_shell['coefficients'] = [ecp_coefficients]
                element_data['ecp_potentials'].append(ecp_shell)

            element_data['ecp_electrons'] = n_elec

        else:
            if not 'electron_shells' in element_data:
                element_data['electron_shells'] = []

            i += 1
            while i < len(basis_lines) and basis_lines[i][0].isalpha() == False:
                lsplt = basis_lines[i].split()
                shell_am = lut.amchar_to_int(lsplt[1])
                nprim = int(lsplt[0])

                if max(shell_am) <= 1:
                    func_type = 'gto'
                else:
                    func_type = 'gto_spherical'

                shell = {
                    'function_type': func_type,
                    'region': '',
                    'angular_momentum': shell_am
                }

                exponents = []
                coefficients = []

                i += 1
                for j in range(nprim):
                    line = basis_lines[i].replace('D', 'E')
                    line = line.replace('d', 'E')
                    lsplt = line.split()
                    exponents.append(lsplt[0])
                    coefficients.append(lsplt[1:])
                    i += 1

                shell['exponents'] = exponents

                # We need to transpose the coefficient matrix
                # (we store a matrix with primitives being the column index and
                # general contraction being the row index)
                shell['coefficients'] = list(map(list, zip(*coefficients)))

                element_data['electron_shells'].append(shell)

    return bs_data
104319
function isEar (data, ear, minX, minY, size) {

  var a = ear.prev.i,
    b = ear.i,
    c = ear.next.i,

    ax = data[a], ay = data[a + 1],
    bx = data[b], by = data[b + 1],
    cx = data[c], cy = data[c + 1],

    abd = ax * by - ay * bx,
    acd = ax * cy - ay * cx,
    cbd = cx * by - cy * bx,
    A = abd - acd - cbd

  if (A <= 0) {
    return false
  } // reflex, can't be an ear

  // now make sure we don't have other points inside the potential ear;
  // the code below is a bit verbose and repetitive but this is done for performance

  var cay = cy - ay,
    acx = ax - cx,
    aby = ay - by,
    bax = bx - ax,
    i, px, py, s, t, k, node

  // if we use z-order curve hashing, iterate through the curve
  if (minX !== undefined) {

    // triangle bbox; min & max are calculated like this for speed
    var minTX = ax < bx ? (ax < cx ? ax : cx) : (bx < cx ? bx : cx),
      minTY = ay < by ? (ay < cy ? ay : cy) : (by < cy ? by : cy),
      maxTX = ax > bx ? (ax > cx ? ax : cx) : (bx > cx ? bx : cx),
      maxTY = ay > by ? (ay > cy ? ay : cy) : (by > cy ? by : cy),

      // z-order range for the current triangle bbox;
      minZ = zOrder(minTX, minTY, minX, minY, size),
      maxZ = zOrder(maxTX, maxTY, minX, minY, size)

    // first look for points inside the triangle in increasing z-order
    node = ear.nextZ

    while (node && node.z <= maxZ) {
      i = node.i
      node = node.nextZ
      if (i === a || i === c) {
        continue
      }

      px = data[i]
      py = data[i + 1]

      s = cay * px + acx * py - acd
      if (s >= 0) {
        t = aby * px + bax * py + abd
        if (t >= 0) {
          k = A - s - t
          if ((k >= 0) && ((s && t) || (s && k) || (t && k))) {
            return false
          }
        }
      }
    }

    // then look for points in decreasing z-order
    node = ear.prevZ

    while (node && node.z >= minZ) {
      i = node.i
      node = node.prevZ
      if (i === a || i === c) {
        continue
      }

      px = data[i]
      py = data[i + 1]

      s = cay * px + acx * py - acd
      if (s >= 0) {
        t = aby * px + bax * py + abd
        if (t >= 0) {
          k = A - s - t
          if ((k >= 0) && ((s && t) || (s && k) || (t && k))) {
            return false
          }
        }
      }
    }

    // if we don't use z-order curve hash, simply iterate through all other points
  } else {
    node = ear.next.next

    while (node !== ear.prev) {
      i = node.i
      node = node.next

      px = data[i]
      py = data[i + 1]

      s = cay * px + acx * py - acd
      if (s >= 0) {
        t = aby * px + bax * py + abd
        if (t >= 0) {
          k = A - s - t
          if ((k >= 0) && ((s && t) || (s && k) || (t && k))) {
            return false
          }
        }
      }
    }
  }

  return true
}
730271
def description(self):
        ''''''
        if 0 < self.algo <= 3:
            return "Reed-Solomon with polynomials in Galois field of characteristic %i (2^%i) with generator=%s, prime poly=%s and first consecutive root=%s." % (self.field_charac, self.c_exp, self.gen_nb, hex(self.prim), self.fcr)
        elif self.algo == 4:
            return "Reed-Solomon with polynomials in Galois field of characteristic %i (2^%i) under US FAA ADSB UAT RS FEC standard with generator=%s, prime poly=%s and first consecutive root=%s." % (self.field_charac, self.c_exp, self.gen_nb, hex(self.prim), self.fcr)
        else:
            return "No description for this ECC algorithm."
941612
async function kdf(hash_algo, X, length, param, stripLeading=false, stripTrailing=false) {
  // Note: X is little endian for Curve25519, big-endian for all others.
  // This is not ideal, but the RFC's are unclear
  // https://tools.ietf.org/html/draft-ietf-openpgp-rfc4880bis-02#appendix-B
  let i;
  if (stripLeading) {
    // Work around old go crypto bug
    for (i = 0; i < X.length && X[i] === 0; i++);
    X = X.subarray(i);
  }
  if (stripTrailing) {
    // Work around old OpenPGP.js bug
    for (i = X.length - 1; i >= 0 && X[i] === 0; i--);
    X = X.subarray(0, i + 1);
  }
  const digest = await hash.digest(hash_algo, util.concatUint8Array([
    new Uint8Array([0, 0, 0, 1]),
    X,
    param
  ]));
  return digest.subarray(0, length);
}
743060
public static KernelFunction autoKernel(Vec dataPoints )
    {
        if(dataPoints.length() < 30)
            return GaussKF.getInstance();
        else if(dataPoints.length() < 1000)
            return EpanechnikovKF.getInstance();
        else//For very large data sets, Uniform is FAST and just as accurate
            return UniformKF.getInstance();
    }
296473
function clippingRectangleRounded($x1, $y1, $w, $h, $rTL, $rTR, $rBR, $rBL)
    {
        $this->save();

        // start: top edge, left end
        $this->addContent(sprintf("\n%.3F %.3F m ", $x1, $y1 - $rTL + $h));

        // line: bottom edge, left end
        $this->addContent(sprintf("\n%.3F %.3F l ", $x1, $y1 + $rBL));

        // curve: bottom-left corner
        $this->ellipse($x1 + $rBL, $y1 + $rBL, $rBL, 0, 0, 8, 180, 270, false, false, false, true);

        // line: right edge, bottom end
        $this->addContent(sprintf("\n%.3F %.3F l ", $x1 + $w - $rBR, $y1));

        // curve: bottom-right corner
        $this->ellipse($x1 + $w - $rBR, $y1 + $rBR, $rBR, 0, 0, 8, 270, 360, false, false, false, true);

        // line: right edge, top end
        $this->addContent(sprintf("\n%.3F %.3F l ", $x1 + $w, $y1 + $h - $rTR));

        // curve: bottom-right corner
        $this->ellipse($x1 + $w - $rTR, $y1 + $h - $rTR, $rTR, 0, 0, 8, 0, 90, false, false, false, true);

        // line: bottom edge, right end
        $this->addContent(sprintf("\n%.3F %.3F l ", $x1 + $rTL, $y1 + $h));

        // curve: top-right corner
        $this->ellipse($x1 + $rTL, $y1 + $h - $rTL, $rTL, 0, 0, 8, 90, 180, false, false, false, true);

        // line: top edge, left end
        $this->addContent(sprintf("\n%.3F %.3F l ", $x1 + $rBL, $y1));

        // Close & clip
        $this->addContent(" W n");
    }