initial commit

scalprod.spyx

@@ -0,0 +1,55 @@
+#
+# scalprod.spyx: faster implementation of ciphertext expand
+#
+# Copyright (c) 2012 Jean-Sebastien Coron <jean-sebastien.coron@uni.lu>
+# and Mehdi Tibouchi <mehdi.tibouchi@normalesup.org>
+#
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License version 2 as published
+# by the Free Software Foundation.
+
+from sage.libs.gmp.all cimport mpz_t, mpz_export
+from sage.rings.integer cimport Integer
+from sage.misc.misc import cputime
+from sage.all import ZZ
+from libc.stdlib cimport calloc, free
+
+cdef extern from "directscal.c":
+    unsigned long directScal(unsigned long kap, mpz_t c, mpz_t y)
+    unsigned getGMP_NUMB_BITS()
+
+def directProd(c, y, kap):
+    cdef Integer yz = <Integer>y
+    cdef Integer cz = <Integer>c
+    return ZZ(directScal(kap, cz.value, yz.value)) >> (getGMP_NUMB_BITS() - 32)
+
+def partialProd(c, y, kap):
+    cdef int w = sizeof(int) * 4  # 16 bits
+    cdef int nw = (kap + 1) // w
+    assert nw * w == kap + 1
+
+    cdef Integer yz = <Integer>y
+    cdef size_t ny
+    cdef unsigned int* yw = <unsigned int*> calloc(nw, sizeof(int))
+
+    cdef int order = -1  # least significant word first
+    cdef int endian = 0
+    cdef size_t nails = sizeof(int) * 4
+    mpz_export(yw, &ny, order, sizeof(int), endian, nails, yz.value)
+    assert ny <= nw
+
+    cdef Integer cz = <Integer>c
+    cdef size_t nc
+    cdef unsigned int* cw = <unsigned int*> calloc(nw, sizeof(int))
+    mpz_export(cw, &nc, order, sizeof(int), endian, nails, cz.value)
+    assert nc <= nw
+
+    cdef unsigned int val = 0
+    for i in range(nw - 2):
+        val += cw[i] * ((yw[nw-1-i] << w) + yw[nw-2-i]) + ((cw[i] * yw[nw-3-i]) >> w)
+
+    val += cw[nw-2] * ((yw[1] << w) + yw[0]) + cw[nw-1] * (yw[0] << w)
+
+    free(cw)
+    free(yw)
+    return ZZ(val)