Ton

1
/*
2
    This file is part of TON Blockchain Library.
3

4
    TON Blockchain Library is free software: you can redistribute it and/or modify
5
    it under the terms of the GNU Lesser General Public License as published by
6
    the Free Software Foundation, either version 2 of the License, or
7
    (at your option) any later version.
8

9
    TON Blockchain Library is distributed in the hope that it will be useful,
10
    but WITHOUT ANY WARRANTY; without even the implied warranty of
11
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12
    GNU Lesser General Public License for more details.
13

14
    You should have received a copy of the GNU Lesser General Public License
15
    along with TON Blockchain Library.  If not, see <http://www.gnu.org/licenses/>.
16

17
    Copyright 2017-2020 Telegram Systems LLP
18
*/
19

20
#include "LossSender.h"
21

22
#include "td/utils/logging.h"
23

24
#if TON_HAVE_GSL
25
#include <gsl/gsl_cdf.h>
26
#endif
27

28
#include <cmath>
29

30
namespace ton {
31
namespace rldp2 {
32
namespace {
33
// works for 1e-x, where x in {1..10}
34
double ndtri_fast(double p) {
35
  if (p < 2e-10) {
36
    return 6.361340902404;
37
  }
38
  if (p < 2e-9) {
39
    return 5.997807015008;
40
  }
41
  if (p < 2e-8) {
42
    return 5.612001244175;
43
  }
44
  if (p < 2e-7) {
45
    return 5.199337582193;
46
  }
47
  if (p < 2e-6) {
48
    return 4.753424308823;
49
  }
50
  if (p < 2e-5) {
51
    return 4.264890793923;
52
  }
53
  if (p < 2e-4) {
54
    return 3.719016485456;
55
  }
56
  if (p < 2e-3) {
57
    return 3.090232306168;
58
  }
59
  if (p < 2e-2) {
60
    return 2.326347874041;
61
  }
62
  return 1.281551565545;
63
}
64
}  // namespace
65

66
LossSender::LossSender(double loss, double p) : loss_(loss), p_(p) {
67
  v_.resize(2);
68
  v_[0] = 1;
69
  res_.push_back(0);
70
  S_ = ndtri_fast(p_);
71
  sigma_ = p * (1 - p);
72
  //LOG(ERROR) << S_ << " " << ndtri(1 - p_);
73
  //CHECK(fabs(S_ - ndtri(1 - p_)) < 1e-6);
74
}
75

76
int LossSender::send_n(int n) {
77
  if (n < 50) {
78
    return send_n_exact(n);
79
  }
80
  return send_n_approx_nbd(n);
81
}
82

83
int LossSender::send_n_exact(int n) {
84
  while ((int)res_.size() <= n) {
85
    step();
86
  }
87
  return res_[n];
88
}
89

90
int LossSender::send_n_approx_norm(int n) {
91
  double a = (1 - loss_) * (1 - loss_);
92
  double b = loss_ * (loss_ - 1) * (2 * n + S_ * S_);
93
  double c = loss_ * loss_ * n * n + S_ * S_ * n * loss_ * (loss_ - 1);
94
  double x = ((-b + sqrt(b * b - 4 * a * c)) / (2 * a));
95
  return (int)(x + n + 1);
96
}
97

98
int LossSender::send_n_approx_nbd(int n) {
99
#if TON_HAVE_GSL
100
  auto mean = n * loss_ / (1 - loss_);
101
  auto var = sqrt(mean / (1 - loss_));
102
  auto min_k = static_cast<int>(mean + var);
103
  auto max_k = min_k + static_cast<int>(var + 1) * 15;
104
  while (min_k + 1 < max_k) {
105
    int k = (min_k + max_k) / 2;
106
    if (gsl_cdf_negative_binomial_P(k, 1 - loss_, n) > 1 - p_) {
107
      max_k = k;
108
    } else {
109
      min_k = k;
110
    }
111
  }
112
  return max_k + n;
113
#endif
114
  return send_n_approx_norm(n);
115
}
116

117
int LossSender::send_n_approx_pd(int n) {
118
#if TON_HAVE_GSL
119
  for (int k = 0;; k++) {
120
    if (gsl_cdf_poisson_P(k, (n + k) * loss_) > 1 - p_) {
121
      return k + n;
122
    }
123
  }
124
#endif
125
  return send_n_approx_norm(n);
126
}
127
bool LossSender::has_good_approx() {
128
#if TON_HAVE_GSL
129
  return true;
130
#else
131
  return false;
132
#endif
133
}
134

135
void LossSender::step() {
136
  n_++;
137
  v_.push_back(0);
138
  v_[n_] = v_[n_ - 1];
139
  for (int j = n_; j >= 0; j--) {
140
    v_[j + 1] += v_[j] * loss_;
141
    v_[j] *= (1 - loss_);
142
  }
143

144
  while (res_i_ < n_ && v_[res_i_] < 1 - p_) {
145
    res_i_++;
146
  }
147
  auto left_ = n_ - res_i_;
148
  if ((int)res_.size() == left_) {
149
    res_.push_back(n_);
150
  }
151
}
152

153
}  // namespace rldp2
154
}  // namespace ton
155