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