exchange_api_stefan.c (8108B)
1 /* 2 This file is part of TALER 3 Copyright (C) 2023 Taler Systems SA 4 5 TALER is free software; you can redistribute it and/or modify it under the 6 terms of the GNU General Public License as published by the Free Software 7 Foundation; either version 3, or (at your option) any later version. 8 9 TALER is distributed in the hope that it will be useful, but WITHOUT ANY 10 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 11 A PARTICULAR PURPOSE. See the GNU General Public License for more details. 12 13 You should have received a copy of the GNU General Public License along with 14 TALER; see the file COPYING. If not, see 15 <http://www.gnu.org/licenses/> 16 */ 17 /** 18 * @file lib/exchange_api_stefan.c 19 * @brief calculations on the STEFAN curve 20 * @author Christian Grothoff 21 */ 22 #include "taler/platform.h" 23 #include "taler/taler_json_lib.h" 24 #include <gnunet/gnunet_curl_lib.h> 25 #include "exchange_api_handle.h" 26 #include <math.h> 27 28 29 /** 30 * Determine smallest denomination in @a keys. 31 * 32 * @param keys exchange response to evaluate 33 * @return NULL on error (no denominations) 34 */ 35 static const struct TALER_Amount * 36 get_unit (const struct TALER_EXCHANGE_Keys *keys) 37 { 38 const struct TALER_Amount *min = NULL; 39 40 for (unsigned int i = 0; i<keys->num_denom_keys; i++) 41 { 42 const struct TALER_EXCHANGE_DenomPublicKey *dk 43 = &keys->denom_keys[i]; 44 45 if ( (NULL == min) || 46 (1 == TALER_amount_cmp (min, 47 /* > */ 48 &dk->value)) ) 49 min = &dk->value; 50 } 51 GNUNET_break (NULL != min); 52 return min; 53 } 54 55 56 /** 57 * Convert amount to double for STEFAN curve evaluation. 58 * 59 * @param a input amount 60 * @return (rounded) amount as a double 61 */ 62 static double 63 amount_to_double (const struct TALER_Amount *a) 64 { 65 double d = (double) a->value; 66 67 d += a->fraction / ((double) TALER_AMOUNT_FRAC_BASE); 68 return d; 69 } 70 71 72 /** 73 * Convert double to amount for STEFAN curve evaluation. 74 * 75 * @param dv input amount 76 * @param currency deisred currency 77 * @param[out] rval (rounded) amount as a double 78 */ 79 static void 80 double_to_amount (double dv, 81 const char *currency, 82 struct TALER_Amount *rval) 83 { 84 double rem; 85 86 GNUNET_assert (GNUNET_OK == 87 TALER_amount_set_zero (currency, 88 rval)); 89 rval->value = floorl (dv); 90 rem = dv - ((double) rval->value); 91 if (rem < 0.0) 92 rem = 0.0; 93 rem *= TALER_AMOUNT_FRAC_BASE; 94 rval->fraction = floorl (rem); 95 if (rval->fraction >= TALER_AMOUNT_FRAC_BASE) 96 { 97 /* Strange, multiplication overflowed our range, 98 round up value instead */ 99 rval->fraction = 0; 100 rval->value += 1; 101 } 102 } 103 104 105 enum GNUNET_GenericReturnValue 106 TALER_EXCHANGE_keys_stefan_b2n ( 107 const struct TALER_EXCHANGE_Keys *keys, 108 const struct TALER_Amount *brut, 109 struct TALER_Amount *net) 110 { 111 const struct TALER_Amount *min; 112 double log_d = amount_to_double (&keys->stefan_log); 113 double lin_d = keys->stefan_lin; 114 double abs_d = amount_to_double (&keys->stefan_abs); 115 double bru_d = amount_to_double (brut); 116 double min_d; 117 double fee_d; 118 double net_d; 119 120 if (TALER_amount_is_zero (brut)) 121 { 122 *net = *brut; 123 return GNUNET_NO; 124 } 125 min = get_unit (keys); 126 if (NULL == min) 127 return GNUNET_SYSERR; 128 if (1.0 <= keys->stefan_lin) 129 { 130 /* This cannot work, linear STEFAN fee estimate always 131 exceed any gross amount. */ 132 GNUNET_break_op (0); 133 return GNUNET_SYSERR; 134 } 135 min_d = amount_to_double (min); 136 fee_d = abs_d 137 + log_d * log2 (bru_d / min_d) 138 + lin_d * bru_d; 139 if (fee_d > bru_d) 140 { 141 GNUNET_assert (GNUNET_OK == 142 TALER_amount_set_zero (brut->currency, 143 net)); 144 return GNUNET_NO; 145 } 146 net_d = bru_d - fee_d; 147 double_to_amount (net_d, 148 brut->currency, 149 net); 150 return GNUNET_OK; 151 } 152 153 154 /** 155 * Our function 156 * f(x) := ne + ab + lo * log2(x/mi) + li * x - x 157 * for #newton(). 158 */ 159 static double 160 eval_f (double mi, 161 double ab, 162 double lo, 163 double li, 164 double ne, 165 double x) 166 { 167 return ne + ab + lo * log2 (x / mi) + li * x - x; 168 } 169 170 171 /** 172 * Our function 173 * f'(x) := lo / log(2) / x + li - 1 174 * for #newton(). 175 */ 176 static double 177 eval_fp (double mi, 178 double lo, 179 double li, 180 double ne, 181 double x) 182 { 183 return lo / log (2) / x + li - 1; 184 } 185 186 187 /** 188 * Use Newton's method to find x where f(x)=0. 189 * 190 * @return x where "eval_f(x)==0". 191 */ 192 static double 193 newton (double mi, 194 double ab, 195 double lo, 196 double li, 197 double ne) 198 { 199 const double eps = 0.00000001; /* max error allowed */ 200 double min_ab = ne + ab; /* result cannot be smaller than this! */ 201 /* compute lower bounds by various heuristics */ 202 double min_ab_li = min_ab + min_ab * li; 203 double min_ab_li_lo = min_ab_li + log2 (min_ab_li / mi) * lo; 204 double min_ab_lo = min_ab + log2 (min_ab / mi) * lo; 205 double min_ab_lo_li = min_ab_lo + min_ab_lo * li; 206 /* take global lower bound */ 207 double x_min = GNUNET_MAX (min_ab_lo_li, 208 min_ab_li_lo); 209 double x = x_min; /* use lower bound as starting point */ 210 211 /* Objective: invert 212 ne := br - ab - lo * log2 (br/mi) - li * br 213 to find 'br'. 214 Method: use Newton's method to find root of: 215 f(x) := ne + ab + lo * log2 (x/mi) + li * x - x 216 using also 217 f'(x) := lo / log(2) / x + li - 1 218 */ 219 /* Loop to abort in case of divergence; 220 100 is already very high, 2-4 is normal! */ 221 for (unsigned int i = 0; i<100; i++) 222 { 223 double fx = eval_f (mi, ab, lo, li, ne, x); 224 double fxp = eval_fp (mi, lo, li, ne, x); 225 double x_new = x - fx / fxp; 226 227 if (fabs (x - x_new) <= eps) 228 { 229 GNUNET_log (GNUNET_ERROR_TYPE_INFO, 230 "Needed %u rounds from %f to result BRUT %f => NET: %f\n", 231 i, 232 x_min, 233 x_new, 234 x_new - ab - li * x_new - lo * log2 (x / mi)); 235 return x_new; 236 } 237 if (x_new < x_min) 238 { 239 GNUNET_log (GNUNET_ERROR_TYPE_WARNING, 240 "Divergence, obtained very bad estimate %f after %u rounds!\n", 241 x_new, 242 i); 243 return x_min; 244 } 245 x = x_new; 246 } 247 GNUNET_log (GNUNET_ERROR_TYPE_WARNING, 248 "Slow convergence, returning bad estimate %f!\n", 249 x); 250 return x; 251 } 252 253 254 enum GNUNET_GenericReturnValue 255 TALER_EXCHANGE_keys_stefan_n2b ( 256 const struct TALER_EXCHANGE_Keys *keys, 257 const struct TALER_Amount *net, 258 struct TALER_Amount *brut) 259 { 260 const struct TALER_Amount *min; 261 double lin_d = keys->stefan_lin; 262 double log_d = amount_to_double (&keys->stefan_log); 263 double abs_d = amount_to_double (&keys->stefan_abs); 264 double net_d = amount_to_double (net); 265 double min_d; 266 double brut_d; 267 268 if (TALER_amount_is_zero (net)) 269 { 270 *brut = *net; 271 return GNUNET_NO; 272 } 273 min = get_unit (keys); 274 if (NULL == min) 275 return GNUNET_SYSERR; 276 if (1.0 <= keys->stefan_lin) 277 { 278 /* This cannot work, linear STEFAN fee estimate always 279 exceed any gross amount. */ 280 GNUNET_break_op (0); 281 return GNUNET_SYSERR; 282 } 283 min_d = amount_to_double (min); 284 brut_d = newton (min_d, 285 abs_d, 286 log_d, 287 lin_d, 288 net_d); 289 double_to_amount (brut_d, 290 net->currency, 291 brut); 292 return GNUNET_OK; 293 } 294 295 296 void 297 TALER_EXCHANGE_keys_stefan_round ( 298 const struct TALER_EXCHANGE_Keys *keys, 299 struct TALER_Amount *val) 300 { 301 const struct TALER_Amount *min; 302 uint32_t mod; 303 uint32_t frac; 304 uint32_t lim; 305 306 if (0 == val->fraction) 307 { 308 /* rounding of non-fractions not supported */ 309 return; 310 } 311 min = get_unit (keys); 312 if (NULL == min) 313 return; 314 if (0 == min->fraction) 315 { 316 frac = TALER_AMOUNT_FRAC_BASE; 317 } 318 else 319 { 320 frac = min->fraction; 321 } 322 lim = frac / 2; 323 mod = val->fraction % frac; 324 if (mod < lim) 325 val->fraction -= mod; /* round down */ 326 else 327 val->fraction += frac - mod; /* round up */ 328 }