1#include <memory.h>
  2
  3#include "zn_arithmetic_private.h"
  4
  5typedef union {
  6  double dv;
  7  int64_t s64v;
  8  int32_t s32v;
  9  uint64_t u64v;
 10  uint32_t u32v;
 11} di_t;
 12
 13/** reduction mod 1, output in torus32 space */
 14EXPORT void dbl_to_tn32_ref(const MOD_Z* module,              //
 15                            int32_t* res, uint64_t res_size,  // res
 16                            const double* a, uint64_t a_size  // a
 17) {
 18  static const double ADD_CST = 0.5 + (double)(INT64_C(3) << (51 - 32));
 19  static const int32_t XOR_CST = (INT32_C(1) << 31);
 20  const uint64_t msize = res_size < a_size ? res_size : a_size;
 21  for (uint64_t i = 0; i < msize; ++i) {
 22    di_t t = {.dv = a[i] + ADD_CST};
 23    res[i] = t.s32v ^ XOR_CST;
 24  }
 25  memset(res + msize, 0, (res_size - msize) * sizeof(int32_t));
 26}
 27
 28/** real centerlift mod 1, output in double space */
 29EXPORT void tn32_to_dbl_ref(const MOD_Z* module,               //
 30                            double* res, uint64_t res_size,    // res
 31                            const int32_t* a, uint64_t a_size  // a
 32) {
 33  static const uint32_t XOR_CST = (UINT32_C(1) << 31);
 34  static const di_t OR_CST = {.dv = (double)(INT64_C(1) << (52 - 32))};
 35  static const double SUB_CST = 0.5 + (double)(INT64_C(1) << (52 - 32));
 36  const uint64_t msize = res_size < a_size ? res_size : a_size;
 37  for (uint64_t i = 0; i < msize; ++i) {
 38    uint32_t ai = a[i] ^ XOR_CST;
 39    di_t t = {.u64v = OR_CST.u64v | (uint64_t)ai};
 40    res[i] = t.dv - SUB_CST;
 41  }
 42  memset(res + msize, 0, (res_size - msize) * sizeof(double));
 43}
 44
 45/** round to the nearest int, output in i32 space */
 46EXPORT void dbl_round_to_i32_ref(const MOD_Z* module,              //
 47                                 int32_t* res, uint64_t res_size,  // res
 48                                 const double* a, uint64_t a_size  // a
 49) {
 50  static const double ADD_CST = (double)((INT64_C(3) << (51)) + (INT64_C(1) << (31)));
 51  static const int32_t XOR_CST = INT32_C(1) << 31;
 52  const uint64_t msize = res_size < a_size ? res_size : a_size;
 53  for (uint64_t i = 0; i < msize; ++i) {
 54    di_t t = {.dv = a[i] + ADD_CST};
 55    res[i] = t.s32v ^ XOR_CST;
 56  }
 57  memset(res + msize, 0, (res_size - msize) * sizeof(int32_t));
 58}
 59
 60/** small int (int32 space) to double */
 61EXPORT void i32_to_dbl_ref(const MOD_Z* module,               //
 62                           double* res, uint64_t res_size,    // res
 63                           const int32_t* a, uint64_t a_size  // a
 64) {
 65  static const uint32_t XOR_CST = (UINT32_C(1) << 31);
 66  static const di_t OR_CST = {.dv = (double)(INT64_C(1) << 52)};
 67  static const double SUB_CST = (double)((INT64_C(1) << 52) + (INT64_C(1) << 31));
 68  const uint64_t msize = res_size < a_size ? res_size : a_size;
 69  for (uint64_t i = 0; i < msize; ++i) {
 70    uint32_t ai = a[i] ^ XOR_CST;
 71    di_t t = {.u64v = OR_CST.u64v | (uint64_t)ai};
 72    res[i] = t.dv - SUB_CST;
 73  }
 74  memset(res + msize, 0, (res_size - msize) * sizeof(double));
 75}
 76
 77/** round to the nearest int, output in int64 space */
 78EXPORT void dbl_round_to_i64_ref(const MOD_Z* module,              //
 79                                 int64_t* res, uint64_t res_size,  // res
 80                                 const double* a, uint64_t a_size  // a
 81) {
 82  static const double ADD_CST = (double)(INT64_C(3) << (51));
 83  static const int64_t AND_CST = (INT64_C(1) << 52) - 1;
 84  static const int64_t SUB_CST = INT64_C(1) << 51;
 85  const uint64_t msize = res_size < a_size ? res_size : a_size;
 86  for (uint64_t i = 0; i < msize; ++i) {
 87    di_t t = {.dv = a[i] + ADD_CST};
 88    res[i] = (t.s64v & AND_CST) - SUB_CST;
 89  }
 90  memset(res + msize, 0, (res_size - msize) * sizeof(int64_t));
 91}
 92
 93/** small int (int64 space) to double */
 94EXPORT void i64_to_dbl_ref(const MOD_Z* module,               //
 95                           double* res, uint64_t res_size,    // res
 96                           const int64_t* a, uint64_t a_size  // a
 97) {
 98  static const uint64_t ADD_CST = UINT64_C(1) << 51;
 99  static const uint64_t AND_CST = (UINT64_C(1) << 52) - 1;
100  static const di_t OR_CST = {.dv = (INT64_C(1) << 52)};
101  static const double SUB_CST = INT64_C(3) << 51;
102  const uint64_t msize = res_size < a_size ? res_size : a_size;
103  for (uint64_t i = 0; i < msize; ++i) {
104    di_t t = {.u64v = ((a[i] + ADD_CST) & AND_CST) | OR_CST.u64v};
105    res[i] = t.dv - SUB_CST;
106  }
107  memset(res + msize, 0, (res_size - msize) * sizeof(double));
108}