1use std::f64;
7
8const SQRT_TWO: f64 = 1.4142135623730950488016887242096980785696718753769;
10const SQRT_TWO_PI: f64 = 2.5066282746310005024157652848110452530069867406099;
11const LN_TWO_PI: f64 = 1.8378770664093454835606594728112352797227949472756;
12const TWO_PI: f64 = 6.283185307179586476925286766559005768394338798750;
13const SQRT_PI_OVER_TWO: f64 = 1.253314137315500251207882642405522626503493370305;
14const SQRT_THREE: f64 = 1.732050807568877293527446341505872366942805253810;
15const SQRT_ONE_OVER_THREE: f64 = 0.577350269189625764509148780501957455647601751270;
16const TWO_PI_OVER_SQRT_TWENTY_SEVEN: f64 = 1.209199576156145233729385505094770488189377498728;
17const SQRT_THREE_OVER_THIRD_ROOT_TWO_PI: f64 = 0.938643487427383566075051356115075878414688769574;
18const PI_OVER_SIX: f64 = 0.523598775598298873077107230546583814032861566563;
19
20const DBL_EPSILON: f64 = f64::EPSILON;
21const DBL_MIN: f64 = f64::MIN_POSITIVE;
22const DBL_MAX: f64 = f64::MAX;
23
24const ETA: f64 = -13.0;
25const VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC: f64 = -f64::MAX;
26const VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM: f64 = f64::MAX;
27
28fn d_int(x: f64) -> f64 {
33 if x > 0.0 {
34 x.floor()
35 } else {
36 -(-x).floor()
37 }
38}
39
40fn smoothened_exponential_of_negative_square(y: f64) -> f64 {
41 let y_tilde = d_int(y * 16.0) / 16.0;
42 let del = (y - y_tilde) * (y + y_tilde);
43 (-y_tilde * y_tilde).exp() * (-del).exp()
44}
45
46fn smoothened_exponential_of_positive_square(x: f64) -> f64 {
47 let x_tilde = d_int(x * 16.0) / 16.0;
48 let del = (x - x_tilde) * (x + x_tilde);
49 (x_tilde * x_tilde).exp() * (del).exp()
50}
51
52const CODY_A: [f64; 5] = [
53 3.1611237438705656,
54 113.864154151050156,
55 377.485237685302021,
56 3209.37758913846947,
57 0.185777706184603153,
58];
59const CODY_B: [f64; 4] = [
60 23.6012909523441209,
61 244.024637934444173,
62 1282.61652607737228,
63 2844.23683343917062,
64];
65const CODY_C: [f64; 9] = [
66 0.564188496988670089,
67 8.88314979438837594,
68 66.1191906371416295,
69 298.635138197400131,
70 881.95222124176909,
71 1712.04761263407058,
72 2051.07837782607147,
73 1230.33935479799725,
74 2.15311535474403846E-8,
75];
76const CODY_D: [f64; 8] = [
77 15.7449261107098347,
78 117.693950891312499,
79 537.181101862009858,
80 1621.38957456669019,
81 3290.79923573345963,
82 4362.61909014324716,
83 3439.36767414372164,
84 1230.33935480374942,
85];
86const CODY_P: [f64; 6] = [
87 0.305326634961232344,
88 0.360344899949804439,
89 0.125781726111229246,
90 0.0160837851487422766,
91 6.58749161529837803E-4,
92 0.0163153871373020978,
93];
94const CODY_Q: [f64; 5] = [
95 2.56852019228982242,
96 1.87295284992346047,
97 0.527905102951428412,
98 0.0605183413124413191,
99 0.00233520497626869185,
100];
101
102fn cody_ab(z: f64) -> f64 {
103 ((((CODY_A[4] * z + CODY_A[0]) * z + CODY_A[1]) * z + CODY_A[2]) * z + CODY_A[3])
104 / ((((z + CODY_B[0]) * z + CODY_B[1]) * z + CODY_B[2]) * z + CODY_B[3])
105}
106
107fn cody_cd(y: f64) -> f64 {
108 ((((((((CODY_C[8] * y + CODY_C[0]) * y + CODY_C[1]) * y + CODY_C[2]) * y + CODY_C[3]) * y
109 + CODY_C[4])
110 * y
111 + CODY_C[5])
112 * y
113 + CODY_C[6])
114 * y
115 + CODY_C[7])
116 / ((((((((y + CODY_D[0]) * y + CODY_D[1]) * y + CODY_D[2]) * y + CODY_D[3]) * y + CODY_D[4])
117 * y
118 + CODY_D[5])
119 * y
120 + CODY_D[6])
121 * y
122 + CODY_D[7])
123}
124
125fn cody_pq(z: f64) -> f64 {
126 z * (((((CODY_P[5] * z + CODY_P[0]) * z + CODY_P[1]) * z + CODY_P[2]) * z + CODY_P[3]) * z
127 + CODY_P[4])
128 / (((((z + CODY_Q[0]) * z + CODY_Q[1]) * z + CODY_Q[2]) * z + CODY_Q[3]) * z + CODY_Q[4])
129}
130
131const ONE_OVER_SQRT_PI: f64 = 0.56418958354775628695;
132const THRESHOLD: f64 = 0.46875;
133const XNEG: f64 = -26.6287357137514;
134const XBIG: f64 = 26.543;
135const XHUGE: f64 = 6.71E7;
136const XMAX: f64 = 2.53E307;
137
138pub fn erfc_cody(x: f64) -> f64 {
139 let y = x.abs();
140 if y <= THRESHOLD {
141 return 1.0 - x * cody_ab(y * y);
142 }
143 let erfc_abs_x = if y >= XBIG {
144 0.0
145 } else {
146 let val = if y <= 4.0 {
147 cody_cd(y)
148 } else {
149 (ONE_OVER_SQRT_PI - cody_pq(1.0 / (y * y))) / y
150 };
151 val * smoothened_exponential_of_negative_square(y)
152 };
153 if x < 0.0 {
154 2.0 - erfc_abs_x
155 } else {
156 erfc_abs_x
157 }
158}
159
160pub fn erf_cody(x: f64) -> f64 {
161 let y = x.abs();
162 if y <= THRESHOLD {
163 return x * cody_ab(y * y);
164 }
165 let erfc_abs_x = if y >= XBIG {
166 0.0
167 } else {
168 let val = if y <= 4.0 {
169 cody_cd(y)
170 } else {
171 (ONE_OVER_SQRT_PI - cody_pq(1.0 / (y * y))) / y
172 };
173 val * smoothened_exponential_of_negative_square(y)
174 };
175 if x < 0.0 {
176 erfc_abs_x - 1.0
177 } else {
178 1.0 - erfc_abs_x
179 }
180}
181
182pub fn erfcx_cody(x: f64) -> f64 {
183 let y = x.abs();
184 if y <= THRESHOLD {
185 let z = y * y;
186 return z.exp() * (1.0 - x * cody_ab(z));
187 }
188 if x < XNEG {
189 return f64::MAX;
190 }
191 let result = if y <= 4.0 {
192 cody_cd(y)
193 } else {
194 (ONE_OVER_SQRT_PI - cody_pq(1.0 / (y * y))) / y
195 };
196 if x < 0.0 {
197 let expx2 = smoothened_exponential_of_positive_square(x);
198 (expx2 + expx2) - result
199 } else {
200 result
201 }
202}
203
204pub fn norm_pdf(x: f64) -> f64 {
209 (1.0 / SQRT_TWO_PI) * (-0.5 * x * x).exp()
210}
211
212pub fn norm_cdf(z: f64) -> f64 {
213 if z <= -10.0 {
214 let mut sum = 1.0;
215 if z >= -1.0 / f64::EPSILON.sqrt() {
216 let zsqr = z * z;
217 let mut i = 1.0;
218 let mut g = 1.0;
219 let mut a = f64::MAX;
220 let mut lasta;
221 loop {
222 lasta = a;
223 let x = (4.0 * i - 3.0) / zsqr;
224 let y = x * ((4.0 * i - 1.0) / zsqr);
225 a = g * (x - y);
226 sum -= a;
227 g *= y;
228 i += 1.0;
229 a = a.abs();
230 if !(lasta > a && a >= (sum * f64::EPSILON).abs()) {
231 break;
232 }
233 }
234 }
235 return -norm_pdf(z) * sum / z;
236 }
237 0.5 * erfc_cody(-z * (1.0 / SQRT_TWO))
238}
239
240fn inverse_norm_cdf_for_low_probabilities(p: f64) -> f64 {
245 let r = (-p.ln()).sqrt();
246 if r < 6.7 {
247 if r < 3.41 {
248 if r < 2.05 {
249 (3.691562302945566191 + r * (4.7170590600740689449E1 + r * (6.5451292110261454609E1 + r * (-7.4594687726045926821E1 + r * (-8.3383894003636969722E1 - 1.3054072340494093704E1 * r))))) / (1.0 + r * (2.0837211328697753726E1 + r * (7.1813812182579255459E1 + r * (5.9270122556046077717E1 + r * (9.2216887978737432303 + 1.8295174852053530579E-4 * r)))))
250 } else {
251 (3.2340179116317970288 + r * (1.449177828689122096E1 + r * (6.8397370256591532878E-1 + r * (-1.81254427791789183E1 + r * (-1.005916339568646151E1 - 1.2013147879435525574E0 * r))))) / (1.0 + r * (8.8820931773304337525 + r * (1.4656370665176799712E1 + r * (7.1369811056109768745 + r * (8.4884892199149255469E-1 + 1.0957576098829595323E-5 * r)))))
252 }
253 } else {
254 (3.1252235780087584807 + r * (9.9483724317036560676 + r * (-5.1633929115525534628 + r * (-1.1070534689309368061E1 + r * (-2.8699061335882526744 - 1.5414319494013597492E-1 * r))))) / (1.0 + r * (7.076769154309171622 + r * (8.1086341122361532407 + r * (2.0307076064309043613 + r * (1.0897972234131828901E-1 + 1.3565983564441297634E-7 * r)))))
255 }
256 } else {
257 if r < 12.9 {
258 (2.6161264950897283681 + r * (2.250881388987032271 + r * (-3.688196041019692267 + r * (-2.9644251353150605663 + r * (-4.7595169546783216436E-1 - 1.612303318390145052E-2 * r))))) / (1.0 + r * (3.2517455169035921495 + r * (2.1282030272153188194 + r * (3.3663746405626400164E-1 + r * (1.1400087282177594359E-2 + 3.0848093570966787291E-9 * r)))))
259 } else {
260 (2.3226849047872302955 + r * (-4.2799650734502094297E-2 + r * (-2.5894451568465728432 + r * (-8.6385181219213758847E-1 + r * (-6.5127593753781672404E-2 - 1.0566357727202585402E-3 * r))))) / (1.0 + r * (1.9361316119254412206 + r * (6.1320841329197493341E-1 + r * (4.6054974512474443189E-2 + r * (7.471447992167225483E-4 + 2.3135343206304887818E-11 * r)))))
261 }
262 }
263}
264
265const U_MAX: f64 = 0.3413447460685429;
266
267fn inverse_norm_cdfm05_for_midrange_probabilities(u: f64) -> f64 {
268 let s = U_MAX * U_MAX - u * u;
269 u * ((2.92958954698308805 + s * (5.0260572167303103E1 + s * (3.01870541922933937E2 + s * (7.4997781456657924E2 + s * (6.90489242061408612E2 + s * (1.34233243502653864E2 - 7.58939881401259242 * s)))))) / (1.0 + s * (1.8918538074574598E1 + s * (1.29404120448755281E2 + s * (3.86821208540417453E2 + s * (4.79123914509756757E2 + 1.79227008508102628E2 * s))))))
270}
271
272pub fn inverse_norm_cdf(p: f64) -> f64 {
273 let u = p - 0.5;
274 if u.abs() < U_MAX {
275 return inverse_norm_cdfm05_for_midrange_probabilities(u);
276 }
277 if u > 0.0 {
278 -inverse_norm_cdf_for_low_probabilities(1.0 - p)
279 } else {
280 inverse_norm_cdf_for_low_probabilities(p)
281 }
282}
283
284pub fn erfinv(e: f64) -> f64 {
285 if e.abs() < (2.0 * U_MAX) {
286 return inverse_norm_cdfm05_for_midrange_probabilities(0.5 * e) * (1.0 / SQRT_TWO);
287 }
288 let val = if e < 0.0 {
289 inverse_norm_cdf_for_low_probabilities(0.5 * e + 0.5)
290 } else {
291 -inverse_norm_cdf_for_low_probabilities(-0.5 * e + 0.5)
292 };
293 val * (1.0 / SQRT_TWO)
294}
295
296const MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE: f64 = -(1.0 - 1.4901161193847656e-08); const MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE: f64 = 2.0 / (f64::EPSILON * f64::EPSILON);
302
303fn is_zero(x: f64) -> bool {
304 x.abs() < f64::MIN_POSITIVE
305}
306
307fn rational_cubic_interpolation(
308 x: f64,
309 x_l: f64,
310 x_r: f64,
311 y_l: f64,
312 y_r: f64,
313 d_l: f64,
314 d_r: f64,
315 r: f64,
316) -> f64 {
317 let h = x_r - x_l;
318 if h.abs() <= 0.0 {
319 return 0.5 * (y_l + y_r);
320 }
321 let t = (x - x_l) / h;
322 if !(r >= MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE) {
323 let omt = 1.0 - t;
324 let t2 = t * t;
325 let omt2 = omt * omt;
326 return (y_r * t2 * t + (r * y_r - h * d_r) * t2 * omt + (r * y_l + h * d_l) * t * omt2
327 + y_l * omt2 * omt)
328 / (1.0 + (r - 3.0) * t * omt);
329 }
330 y_r * t + y_l * (1.0 - t)
331}
332
333fn rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
334 x_l: f64,
335 x_r: f64,
336 y_l: f64,
337 y_r: f64,
338 d_l: f64,
339 d_r: f64,
340 second_derivative_l: f64,
341) -> f64 {
342 let h = x_r - x_l;
343 let numerator = 0.5 * h * second_derivative_l + (d_r - d_l);
344 let denominator = (y_r - y_l) / h - d_l;
345 if is_zero(denominator) {
346 if numerator > 0.0 {
347 MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE
348 } else {
349 MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE
350 }
351 } else {
352 numerator / denominator
353 }
354}
355
356fn rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
357 x_l: f64,
358 x_r: f64,
359 y_l: f64,
360 y_r: f64,
361 d_l: f64,
362 d_r: f64,
363 second_derivative_r: f64,
364) -> f64 {
365 let h = x_r - x_l;
366 let numerator = 0.5 * h * second_derivative_r + (d_r - d_l);
367 let denominator = d_r - (y_r - y_l) / h;
368 if is_zero(denominator) {
369 if numerator > 0.0 {
370 MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE
371 } else {
372 MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE
373 }
374 } else {
375 numerator / denominator
376 }
377}
378
379fn minimum_rational_cubic_control_parameter(
380 d_l: f64,
381 d_r: f64,
382 s: f64,
383 prefer_shape_preservation: bool,
384) -> f64 {
385 let monotonic = d_l * s >= 0.0 && d_r * s >= 0.0;
386 let convex = d_l <= s && s <= d_r;
387 let concave = d_l >= s && s >= d_r;
388 if !monotonic && !convex && !concave {
389 return MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE;
390 }
391 let d_r_m_d_l = d_r - d_l;
392 let d_r_m_s = d_r - s;
393 let s_m_d_l = s - d_l;
394 let mut r1 = f64::MIN;
395 let mut r2 = r1;
396 if monotonic {
397 if !is_zero(s) {
398 r1 = (d_r + d_l) / s;
399 } else if prefer_shape_preservation {
400 r1 = MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE;
401 }
402 }
403 if convex || concave {
404 if !(is_zero(s_m_d_l) || is_zero(d_r_m_s)) {
405 r2 = (d_r_m_d_l / d_r_m_s).abs().max((d_r_m_d_l / s_m_d_l).abs());
406 } else if prefer_shape_preservation {
407 r2 = MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE;
408 }
409 } else if monotonic && prefer_shape_preservation {
410 r2 = MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE;
411 }
412 MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE.max(r1.max(r2))
413}
414
415fn convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
416 x_l: f64,
417 x_r: f64,
418 y_l: f64,
419 y_r: f64,
420 d_l: f64,
421 d_r: f64,
422 second_derivative_l: f64,
423 prefer_shape_preservation: bool,
424) -> f64 {
425 let r = rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
426 x_l,
427 x_r,
428 y_l,
429 y_r,
430 d_l,
431 d_r,
432 second_derivative_l,
433 );
434 let r_min = minimum_rational_cubic_control_parameter(
435 d_l,
436 d_r,
437 (y_r - y_l) / (x_r - x_l),
438 prefer_shape_preservation,
439 );
440 r.max(r_min)
441}
442
443fn convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
444 x_l: f64,
445 x_r: f64,
446 y_l: f64,
447 y_r: f64,
448 d_l: f64,
449 d_r: f64,
450 second_derivative_r: f64,
451 prefer_shape_preservation: bool,
452) -> f64 {
453 let r = rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
454 x_l,
455 x_r,
456 y_l,
457 y_r,
458 d_l,
459 d_r,
460 second_derivative_r,
461 );
462 let r_min = minimum_rational_cubic_control_parameter(
463 d_l,
464 d_r,
465 (y_r - y_l) / (x_r - x_l),
466 prefer_shape_preservation,
467 );
468 r.max(r_min)
469}
470
471fn householder3_factor(nu: f64, h2: f64, h3: f64) -> f64 {
476 (1.0 + 0.5 * h2 * nu) / (1.0 + nu * (h2 + h3 * nu * (1.0 / 6.0)))
477}
478
479fn householder4_factor(nu: f64, h2: f64, h3: f64, h4: f64) -> f64 {
480 (1.0 + nu * (h2 + nu * h3 * (1.0 / 6.0)))
481 / (1.0 + nu * (1.5 * h2 + nu * (h2 * h2 * 0.25 + h3 * (1.0 / 3.0) + nu * h4 * (1.0 / 24.0))))
482}
483
484fn householder_factor_3(nu: f64, h2: f64, h3: f64) -> f64 {
485 householder3_factor(nu, h2, h3)
486}
487
488fn householder_factor_4(nu: f64, h2: f64, h3: f64, h4: f64) -> f64 {
489 householder4_factor(nu, h2, h3, h4)
490}
491
492fn normalised_intrinsic(theta_x: f64) -> f64 {
493 if theta_x <= 0.0 {
494 0.0
495 } else {
496 (0.5 * theta_x).sinh() * 2.0
497 }
498}
499
500fn square(x: f64) -> f64 {
501 x * x
502}
503
504const SIXTEENTH_ROOT_DBL_EPSILON: f64 = 0.10511205190671433; const TAU: f64 = 2.0 * SIXTEENTH_ROOT_DBL_EPSILON;
506
507fn asymptotic_expansion_of_scaled_normalised_black(h: f64, t: f64) -> f64 {
508 let e = square(t / h);
509 let r = (h + t) * (h - t);
510 let q = square(h / r);
511
512 let a0 = 2.0;
513 let a1 = -6.0 - 2.0 * e;
514 let a2 = 30.0 + e * (60.0 + 6.0 * e);
515 let a3 = -2.1E2 + e * (-1.05E3 + e * (-6.3E2 - 30.0 * e));
516 let a4 = 1.89E3 + e * (1.764E4 + e * (2.646E4 + e * (7.56E3 + 2.1E2 * e)));
517 let a5 = -2.079E4 + e * (-3.1185E5 + e * (-8.7318E5 + e * (-6.237E5 + e * (-1.0395E5 - 1.89E3 * e))));
518 let a6 = 2.7027E5
519 + e * (5.94594E6
520 + e * (2.675673E7 + e * (3.567564E7 + e * (1.486485E7 + e * (1.62162E6 + 2.079E4 * e)))));
521 let a7 = -4.05405E6
522 + e * (-1.2297285E8
523 + e * (-8.1162081E8
524 + e * (-1.73918745E9
525 + e * (-1.35270135E9 + e * (-3.6891855E8 + e * (-2.837835E7 - 2.7027E5 * e))))));
526 let a8 = 6.891885E7
527 + e * (2.756754E9
528 + e * (2.50864614E10
529 + e * (7.88431644E10
530 + e * (9.85539555E10
531 + e * (5.01729228E10
532 + e * (9.648639E9 + e * (5.513508E8 + 4.05405E6 * e)))))));
533 let a9 = -1.30945815E9
534 + e * (-6.678236565E10
535 + e * (-8.013883878E11
536 + e * (-3.4726830138E12
537 + e * (-6.3665855253E12
538 + e * (-5.2090245207E12
539 + e * (-1.8699062382E12
540 + e * (-2.671294626E11
541 + e * (-1.178512335E10 - 6.891885E7 * e))))))));
542 let a10 = 2.749862115E10
543 + e * (1.7415793395E12
544 + e * (2.664616389435E13
545 + e * (1.52263793682E14
546 + e * (3.848890340295E14
547 + e * (4.618668408354E14
548 + e * (2.664616389435E14
549 + e * (7.10564370516E13
550 + e * (7.83710702775E12 + e * (2.749862115E11 + 1.30945815E9 * e)))))))));
551 let a11 = -6.3246828645E11
552 + e * (-4.870005805665E13
553 + e * (-9.2530110307635E14
554 + e * (-6.74147946527055E15
555 + e * (-2.24715982175685E16
556 + e * (-3.71802806872497E16
557 + e * (-3.14602375045959E16
558 + e * (-1.34829589305411E16
559 + e * (-2.77590330922905E15
560 + e * (-2.4350029028325E14
561 + e * (-6.95715115095E12 - 2.749862115E10 * e))))))))));
562 let a12 = 1.581170716125E13
563 + e * (1.454677058835E15
564 + e * (3.36030400590885E16
565 + e * (3.04027505296515E17
566 + e * (1.29211689751018875E18
567 + e * (2.81916414002223E18
568 + e * (3.289024830025935E18
569 + e * (2.067387036016302E18
570 + e * (6.8406188691715875E17
571 + e * (1.12010133530295E17
572 + e * (8.0007238235925E15
573 + e * (1.89740485935E14 + 6.3246828645E11 * e)))))))))));
574 let a13 = -4.2691609335375E14
575 + e * (-4.624924344665625E16
576 + e * (-1.2764791191277125E18
577 + e * (-1.40412703104048375E19
578 + e * (-7.41067044160255312E19
579 + e * (-2.06151377739125569E20
580 + e * (-3.17155965752500875E20
581 + e * (-2.74868503652167425E20
582 + e * (-1.33392067948845956E20
583 + e * (-3.51031757760120938E19
584 + e * (-4.6804234368016125E18
585 + e * (-2.774954606799375E17
586 + e * (-5.54990921359875E15 - 1.581170716125E13 * e))))))))))));
587 let a14 = 1.238056670725875E16
588 + e * (1.5599514051146025E18
589 + e * (5.06984206662245812E19
590 + e * (6.66322100184665925E20
591 + e * (4.27556680951827302E21
592 + e * (1.47701398874267613E22
593 + e * (2.89721974714909549E22
594 + e * (3.31110828245610914E22
595 + e * (2.2155209831140142E22
596 + e * (8.55113361903654604E21
597 + e * (1.83238577550783129E21
598 + e * (2.02793682664898325E20
599 + e * (1.01396841332449162E19
600 + e * (1.733279339016225E17 + 4.2691609335375E14 * e)))))))))))));
601 let a15 = -3.8379756792502125E17
602 + e * (-5.56506473491280812E19
603 + e * (-2.10359446979704147E21
604 + e * (-3.25556286992399275E22
605 + e * (-2.49593153360839444E23
606 + e * (-1.04829124411552567E24
607 + e * (-2.55352995361474201E24
608 + e * (-3.72085793241005264E24
609 + e * (-3.28310994036181115E24
610 + e * (-1.74715207352587611E24
611 + e * (-5.49104937393846778E23
612 + e * (-9.76668860977197826E22
613 + e * (-9.11557603578717971E21
614 + e * (-3.89554531443896569E20
615 + e * (-5.75696351887531875E18 - 1.238056670725875E16 * e))))))))))))));
616 let a16 = 1.26653197415257012E19
617 + e * (2.09399953059891594E21
618 + e * (9.10889795810528434E22
619 + e * (1.63960163245895118E24
620 + e * (1.48019591819210871E25
621 + e * (7.42789224401858187E25
622 + e * (2.19979885688242617E26
623 + e * (3.98058840769200926E26
624 + e * (4.47816195865351041E26
625 + e * (3.1425697955463231E26
626 + e * (1.36178024473674001E26
627 + e * (3.55247020366106089E25
628 + e * (5.32870530549159134E24
629 + e * (4.25081904711579936E23
630 + e * (1.57049964794918696E22
631 + e * (2.0264511586441122E20 + 3.8379756792502125E17 * e)))))))))))))));
632
633 let mut omega = 0.0;
634 let thresholds = [
635 12.347, 12.958, 13.729, 14.718, 16.016, 17.769, 20.221, 23.816, 29.419, 38.93, 57.171,
636 99.347,
637 ];
638 let val_to_check = -h - t + TAU + 0.5;
639 let idx = thresholds.iter().position(|&thr| thr > val_to_check).unwrap_or(thresholds.len());
640
641 if idx <= 0 { omega = q * (a16 + omega); }
642 if idx <= 1 { omega = q * (a15 + omega); }
643 if idx <= 2 { omega = q * (a14 + omega); }
644 if idx <= 3 { omega = q * (a13 + omega); }
645 if idx <= 4 { omega = q * (a12 + omega); }
646 if idx <= 5 { omega = q * (a11 + omega); }
647 if idx <= 6 { omega = q * (a10 + omega); }
648 if idx <= 7 { omega = q * (a9 + omega); }
649 if idx <= 8 { omega = q * (a8 + omega); }
650 if idx <= 9 { omega = q * (a7 + omega); }
651 if idx <= 10 { omega = q * (a6 + omega); }
652 if idx <= 11 { omega = q * (a5 + omega); }
653
654 omega = a0 + q * (a1 + q * (a2 + q * (a3 + q * (a4 + omega))));
655
656 (t / r) * omega
657}
658
659fn yprime_tail_expansion_rational_function_part(w: f64) -> f64 {
660 w * (-2.9999999999994663866
661 + w * (-1.7556263323542206288E2
662 + w * (-3.4735035445495633334E3
663 + w * (-2.7805745693864308643E4
664 + w * (-8.3836021460741980839E4 - 6.6818249032616849037E4 * w)))))
665 / (1.0
666 + w * (6.3520877744831739102E1
667 + w * (1.4404389037604337538E3
668 + w * (1.4562545638507033944E4
669 + w * (6.6886794165651675684E4
670 + w * (1.2569970380923908488E5 + 6.9286518679803751694E4 * w))))))
671}
672
673fn yprime(h: f64) -> f64 {
674 if h < -4.0 {
675 let w = 1.0 / (h * h);
676 return w * (1.0 + yprime_tail_expansion_rational_function_part(w));
677 }
678 if h <= -0.46875 {
679 return (1.0000000000594317229
680 - h * (6.1911449879694112749E-1
681 - h * (2.2180844736576013957E-1
682 - h * (4.5650900351352987865E-2
683 - h * (5.545521007735379052E-3
684 - h * (3.0717392274913902347E-4
685 - h * (4.2766597835908713583E-8 + 8.4592436406580605619E-10 * h)))))))
686 / (1.0
687 - h * (1.8724286369589162071
688 - h * (1.5685497236077651429
689 - h * (7.6576489836589035112E-1
690 - h * (2.3677701403094640361E-1
691 - h * (4.6762548903194957675E-2
692 - h * (5.5290453576936595892E-3 - 3.0822020417927147113E-4 * h)))))));
693 }
694 1.0 + h * SQRT_PI_OVER_TWO * erfcx_cody(-(1.0 / SQRT_TWO) * h)
695}
696
697fn small_t_expansion_of_scaled_normalised_black(h: f64, t: f64) -> f64 {
698 let a = yprime(h);
699 let h2 = h * h;
700 let t2 = t * t;
701 let b0 = 2.0 * a;
702 let b1 = (-1.0 + a * (3.0 + h2)) / 3.0;
703 let b2 = (-7.0 - h2 + a * (15.0 + h2 * (10.0 + h2))) / 60.0;
704 let b3 = (-57.0 + (-18.0 - h2) * h2 + a * (105.0 + h2 * (105.0 + h2 * (21.0 + h2)))) / 2520.0;
705 let b4 = (-561.0 + h2 * (-285.0 + (-33.0 - h2) * h2) + a * (945.0 + h2 * (1260.0 + h2 * (378.0 + h2 * (36.0 + h2))))) / 181440.0;
706 let b5 = (-6555.0 + h2 * (-4680.0 + h2 * (-840.0 + (-52.0 - h2) * h2)) + a * (10395.0 + h2 * (17325.0 + h2 * (6930.0 + h2 * (990.0 + h2 * (55.0 + h2)))))) / 19958400.0;
707 let b6 = (-89055.0 + h2 * (-82845.0 + h2 * (-20370.0 + h2 * (-1926.0 + (-75.0 - h2) * h2))) + a * (135135.0 + h2 * (270270.0 + h2 * (135135.0 + h2 * (25740.0 + h2 * (2145.0 + h2 * (78.0 + h2))))))) / 3113510400.0;
708 t * (b0 + t2 * (b1 + t2 * (b2 + t2 * (b3 + t2 * (b4 + t2 * (b5 + b6 * t2))))))
709}
710
711fn normalised_black_with_optimal_use_of_codys_functions(theta_x: f64, s: f64) -> f64 {
712 let codys_threshold = 0.46875;
713 let h = theta_x / s;
714 let t = 0.5 * s;
715 let q1 = -(1.0 / SQRT_TWO) * (h + t);
716 let q2 = -(1.0 / SQRT_TWO) * (h - t);
717 let two_b = if q1 < codys_threshold {
718 if q2 < codys_threshold {
719 (0.5 * theta_x).exp() * erfc_cody(q1) - (-0.5 * theta_x).exp() * erfc_cody(q2)
720 } else {
721 (0.5 * theta_x).exp() * erfc_cody(q1) - (-0.5 * (h * h + t * t)).exp() * erfcx_cody(q2)
722 }
723 } else {
724 if q2 < codys_threshold {
725 (-0.5 * (h * h + t * t)).exp() * erfcx_cody(q1) - (-0.5 * theta_x).exp() * erfc_cody(q2)
726 } else {
727 (-0.5 * (h * h + t * t)).exp() * (erfcx_cody(q1) - erfcx_cody(q2))
728 }
729 };
730 (0.5 * two_b).max(0.0)
731}
732
733fn normalised_vega(x: f64, s: f64) -> f64 {
734 let h = x / s;
735 let t = 0.5 * s;
736 (1.0 / SQRT_TWO_PI) * (-0.5 * (h * h + t * t)).exp()
737}
738
739fn inv_normalised_vega(x: f64, s: f64) -> f64 {
740 let h = x / s;
741 let t = 0.5 * s;
742 SQRT_TWO_PI * (0.5 * (h * h + t * t)).exp()
743}
744
745fn ln_normalised_vega(x: f64, s: f64) -> f64 {
746 let h = x / s;
747 let t = 0.5 * s;
748 -(LN_TWO_PI * 0.5) - 0.5 * (h * h + t * t)
749}
750
751fn is_region_i(theta_x: f64, s: f64) -> bool {
752 theta_x < s * ETA && s * (0.5 * s - (TAU + 0.5 + ETA)) + theta_x < 0.0
753}
754
755fn is_region_ii(theta_x: f64, s: f64) -> bool {
756 s * (s - (2.0 * TAU)) - theta_x / ETA < 0.0
757}
758
759fn normalised_black_internal(theta_x: f64, s: f64) -> f64 {
760 if is_region_i(theta_x, s) {
761 return asymptotic_expansion_of_scaled_normalised_black(theta_x / s, 0.5 * s)
762 * normalised_vega(theta_x, s);
763 }
764 if is_region_ii(theta_x, s) {
765 return small_t_expansion_of_scaled_normalised_black(theta_x / s, 0.5 * s)
766 * normalised_vega(theta_x, s);
767 }
768 normalised_black_with_optimal_use_of_codys_functions(theta_x, s)
769}
770
771fn scaled_normalised_black_and_ln_vega_internal(theta_x: f64, s: f64) -> (f64, f64) {
772 if is_region_i(theta_x, s) {
773 return (
774 asymptotic_expansion_of_scaled_normalised_black(theta_x / s, 0.5 * s),
775 ln_normalised_vega(theta_x, s),
776 );
777 }
778 if is_region_ii(theta_x, s) {
779 return (
780 small_t_expansion_of_scaled_normalised_black(theta_x / s, 0.5 * s),
781 ln_normalised_vega(theta_x, s),
782 );
783 }
784 let ln_vega = ln_normalised_vega(theta_x, s);
785 (
786 normalised_black_with_optimal_use_of_codys_functions(theta_x, s) * (-ln_vega).exp(),
787 ln_vega,
788 )
789}
790
791fn compute_f_lower_map_and_first_two_derivatives(x: f64, s: f64) -> (f64, f64, f64) {
792 let ax = x.abs();
793 let z = SQRT_ONE_OVER_THREE * ax / s;
794 let y = z * z;
795 let s2 = s * s;
796 let phi_z = 0.5 * erfc_cody((1.0 / SQRT_TWO) * z); let pdf_z = norm_pdf(z);
798 let fpp = PI_OVER_SIX * y / (s2 * s) * phi_z
799 * (8.0 * SQRT_THREE * s * ax + (3.0 * s2 * (s2 - 8.0) - 8.0 * x * x) * phi_z / pdf_z)
800 * (2.0 * y + 0.25 * s2).exp();
801 let phi_z2 = phi_z * phi_z;
802 let fp = TWO_PI * y * phi_z2 * (y + 0.125 * s2).exp();
803 let f = TWO_PI_OVER_SQRT_TWENTY_SEVEN * ax * (phi_z2 * phi_z);
804 (f, fp, fpp)
805}
806
807fn inverse_f_lower_map(x: f64, f: f64) -> f64 {
808 (x / (SQRT_THREE
809 * inverse_norm_cdf(
810 SQRT_THREE_OVER_THIRD_ROOT_TWO_PI * f.cbrt() / x.abs().cbrt(),
811 )))
812 .abs()
813}
814
815fn compute_f_upper_map_and_first_two_derivatives(x: f64, s: f64) -> (f64, f64, f64) {
816 let f = 0.5 * erfc_cody((0.5 / SQRT_TWO) * s); let w = square(x / s);
818 let fp = -0.5 * (0.5 * w).exp();
819 let fpp = SQRT_PI_OVER_TWO * (w + 0.125 * s * s).exp() * w / s;
820 (f, fp, fpp)
821}
822
823fn inverse_f_upper_map(f: f64) -> f64 {
824 -2.0 * inverse_norm_cdf(f)
825}
826
827fn one_minus_erfcx(x: f64) -> f64 {
828 if x < -1.0 / 5.0 || x > 1.0 / 3.0 {
829 return 1.0 - erfcx_cody(x);
830 }
831 x * (1.128379167095512573896
832 - x * (1.0000000000000002
833 + x * (1.1514967181784756
834 + x * (5.7689001208873741E-1
835 + x * (1.4069188744609651E-1 + 1.4069285713634565E-2 * x))))
836 / (1.0
837 + x * (1.9037494962421563
838 + x * (1.5089908593742723
839 + x * (6.2486081658640257E-1
840 + x * (1.358008134514386E-1 + 1.2463320728346347E-2 * x))))))
841}
842
843fn implied_normalised_volatility_atm(beta: f64) -> f64 {
844 let beta_max = 0.6826894921370859;
845 if beta <= beta_max {
846 let r = beta_max * beta_max - beta * beta;
847 return beta
848 * ((2.92958954698308816
849 + r * (1.4014698674754995E1
850 + r * (2.44918990556468762E1
851 + r * (1.90763928424894996E1
852 + r * (6.43250149461895996
853 + r * (7.52328633671821543E-1 + 1.38781536163865582E-2 * r))))))
854 / (1.0
855 + r * (5.22443271807813073
856 + r * (1.02258209975070629E1
857 + r * (9.28187483709036392
858 + r * (3.9095549184069553
859 + r * (6.61214199809055912E-1 + 2.89411828874884851E-2 * r)))))));
860 }
861 -2.0 * inverse_norm_cdf_for_low_probabilities(0.5 * (1.0 - beta))
862}
863
864fn b_l_over_b_max(s_c: f64) -> f64 {
865 if s_c >= 2.6267851073127395 {
866 if s_c >= 7.348469228349534 {
867 let s_c_val = s_c;
868 return (1.4500072297240603183E-3 + s_c_val * (-1.5116692485011195757E-3 + s_c_val * (7.1682178310936334831E-2 + s_c_val * (3.921610857820463493E-2 + s_c_val * (2.9342405658628443931E-2 + s_c_val * (5.1832526171631521426E-3 + 1.6930208078421474854E-3 * s_c_val)))))) / (1.0 + s_c_val * (1.6176313502305414664 + s_c_val * (1.6823159175281531664 + s_c_val * (8.4878307567372222113E-1 + s_c_val * (3.7543742137375791321E-1 + s_c_val * (7.126137099644302999E-2 + 1.6116992546788676159E-2 * s_c_val))))));
869 }
870 return (-9.3325115354837883291E-5 + s_c * (5.3118033972794648837E-4 + s_c * (7.4114855448345002595E-2 + s_c * (7.4039658186822817454E-2 + s_c * (3.9225177407687604785E-2 + s_c * (1.0022913378254090083E-2 + 1.7012579407246055469E-3 * s_c)))))) / (1.0 + s_c * (2.2217238132228132256 + s_c * (2.3441816707087403282 + s_c * (1.3912323646271141826 + s_c * (5.3231258443501838354E-1 + s_c * (1.1744005919716101572E-1 + 1.6195405895930935811E-2 * s_c))))));
871 }
872 if s_c >= 0.7099295739719539 {
873 return (1.9795737927598581235E-9 + s_c * (-2.7081288564685588037E-8 + s_c * (7.5610142272549044609E-2 + s_c * (6.917130174466834016E-2 + s_c * (2.9537058950963019803E-2 + s_c * (6.5849252702302307774E-3 + 6.9711400639834715731E-4 * s_c)))))) / (1.0 + s_c * (2.1941448525586579756 + s_c * (2.1297103549995181357 + s_c * (1.1571483187179784072 + s_c * (3.7831622253060456794E-1 + s_c * (7.1714862448829349869E-2 + 6.6361975827861200167E-3 * s_c))))));
874 }
875 let g = (8.0741072372882856924E-2 + s_c * (9.8078911786358897272E-2 + s_c * (3.9760631445677058375E-2 + s_c * (5.9716928459589189876E-3 + s_c * (-6.4036399341479799981E-6 + 4.5425102093616062245E-7 * s_c))))) / (1.0 + s_c * (1.8594977672287664353 + s_c * (1.3658801475711790419 + s_c * (4.6132707108655653215E-1 + 6.1254597049831720643E-2 * s_c))));
876 (s_c * s_c) * (0.07560996640296361767172 + s_c * (s_c * g - 0.09672719281339436290858))
877}
878
879fn b_u_over_b_max(s_c: f64) -> f64 {
880 if s_c > 1.7888543819998317 {
881 if s_c > 6.164414002968976 {
882 let s_c_val = s_c;
883 return (7.91133825948419359E-1 + s_c_val * (1.24653733210880042 + s_c_val * (1.32747426980537386 + s_c_val * (6.95009705717846778E-1 + s_c_val * (3.05965944268228457E-1 + s_c_val * (6.02200363391352887E-2 + 1.29050244454344842E-2 * s_c_val)))))) / (1.0 + s_c_val * (1.58117486714634672 + s_c_val * (1.60144713247629644 + s_c_val * (8.30040185836882436E-1 + s_c_val * (3.53071863813401531E-1 + s_c_val * (6.95901684131758475E-2 + 1.44197580643890011E-2 * s_c_val))))));
884 }
885 return (7.8990640048967596475E-1 + s_c * (1.5993699253596663678 + s_c * (1.6481729039140370242 + s_c * (9.8227188109869200166E-1 + s_c * (3.6313557966186936883E-1 + s_c * (7.8277036261179606301E-2 + 9.3404307364538726214E-3 * s_c)))))) / (1.0 + s_c * (2.0247407005640401446 + s_c * (2.0087454279103740489 + s_c * (1.1627561803056961973 + s_c * (4.2004672123723823581E-1 + s_c * (8.9130862793887234546E-2 + 1.0436767768858021717E-2 * s_c))))));
886 }
887 if s_c >= 0.7745966692414833 {
888 return (7.8990944435755287611E-1 + s_c * (-1.2655410534988972886 + s_c * (-2.8803040699221003256 + s_c * (-2.6936198689113258727 + s_c * (-1.1213067281643205754 + s_c * (-2.1277793801691629892E-1 + 5.1486445905299802703E-6 * s_c)))))) / (1.0 + s_c * (-1.6021222722060444448 + s_c * (-3.7242680976480704555 + s_c * (-3.2083117718907365085 + s_c * (-1.2922333835930958583 - 2.3762328334050001161E-1 * s_c)))));
889 }
890 let g = (-6.063099881233561706E-2 + s_c * (-8.1011946637120604985E-2 + s_c * (-4.2505564862438753828E-2 + s_c * (-8.9880000946868691788E-3 + s_c * (-7.5603072110443268356E-6 + 4.3879556621540147458E-7 * s_c))))) / (1.0 + s_c * (1.8400371530721828756 + s_c * (1.5709283443886143691 + s_c * (6.8913245453611400484E-1 + 1.4703173061720980923E-1 * s_c))));
891 0.7899085945560627246288 + (s_c * s_c) * (0.0614616805805147403487 + s_c * g)
892}
893
894pub fn lets_be_rational(beta: f64, theta_x: f64, n_iterations: i32) -> f64 {
895 if beta <= 0.0 {
896 return if beta == 0.0 { 0.0 } else { VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC };
897 }
898 let b_max = (0.5 * theta_x).exp();
899 if beta >= b_max {
900 return VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM;
901 }
902 if theta_x == 0.0 {
903 return implied_normalised_volatility_atm(beta);
904 }
905
906 let mut iterations = 0;
907 let mut f;
908 let mut s;
909 let mut ds = -DBL_MAX;
910 let mut s_left = DBL_MIN;
911 let mut s_right = DBL_MAX;
912
913 let sqrt_ax = (-theta_x).sqrt();
914 let s_c = SQRT_TWO * sqrt_ax;
915 let ome = one_minus_erfcx(sqrt_ax);
916 let b_c = 0.5 * b_max * ome;
917
918 if beta < b_c {
919 let s_l = s_c - SQRT_PI_OVER_TWO * ome;
920 let b_l = b_l_over_b_max(s_c) * b_max;
921 if beta < b_l {
922 let (f_lower_map_l, d_f_lower_map_l_d_beta, d2_f_lower_map_l_d_beta2) =
923 compute_f_lower_map_and_first_two_derivatives(theta_x, s_l);
924 let r_ll = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
925 0.0,
926 b_l,
927 0.0,
928 f_lower_map_l,
929 1.0,
930 d_f_lower_map_l_d_beta,
931 d2_f_lower_map_l_d_beta2,
932 true,
933 );
934 f = rational_cubic_interpolation(beta, 0.0, b_l, 0.0, f_lower_map_l, 1.0, d_f_lower_map_l_d_beta, r_ll);
935 if !(f > 0.0) {
936 let t = beta / b_l;
937 f = (f_lower_map_l * t + b_l * (1.0 - t)) * t;
938 }
939 s = inverse_f_lower_map(theta_x, f);
940 s_right = s_l;
941
942 let ln_beta = beta.ln();
943 while iterations < n_iterations && ds.abs() > DBL_EPSILON * s {
944 let (bx, ln_vega) = scaled_normalised_black_and_ln_vega_internal(theta_x, s);
945 let ln_b = bx.ln() + ln_vega;
946 let bpob = 1.0 / bx;
947 let b = ln_b.exp();
948 if b > beta && s < s_right {
949 s_right = s;
950 } else if b < beta && s > s_left {
951 s_left = s;
952 }
953
954 let h = theta_x / s;
955 let x2_over_s3 = h * h / s;
956 let b_h2 = x2_over_s3 - s / 4.0;
957 let nu = (ln_beta - ln_b) * ln_b / ln_beta / bpob;
958 let lambda = 1.0 / ln_b;
959 let otl = 1.0 + 2.0 * lambda;
960 let h2 = b_h2 - bpob * otl;
961 let c = 3.0 * (x2_over_s3 / s);
962 let b_h3 = b_h2 * b_h2 - c - 0.25;
963 let mu = 6.0 * lambda * (1.0 + lambda);
964 let h3 = b_h3 + (bpob * bpob) * (2.0 + mu) - (b_h2 * bpob) * 3.0 * otl;
965
966 if theta_x < -190.0 {
967 let b_h4 = b_h2 * (b_h3 - 0.5) - (b_h2 - 2.0 / s) * 2.0 * c;
968 let bpob2 = bpob * bpob;
969 let bppob = b_h2 * bpob;
970 ds = nu * householder_factor_4(
971 nu,
972 h2,
973 h3,
974 b_h4 - bpob * (bpob2 * (6.0 + lambda * (22.0 + lambda * (36.0 + lambda * 24.0))) - bppob * (12.0 + 6.0 * mu)) - bppob * b_h2 * 3.0 * otl - b_h3 * bpob * 4.0 * otl,
975 );
976 } else {
977 ds = nu * householder_factor_3(nu, h2, h3);
978 }
979 s += ds;
980 if s > s_right { s = s_right; } else if s < s_left { s = s_left; }
981 iterations += 1;
982 }
983 return s;
984 } else {
985 let inv_v_c = SQRT_TWO_PI / b_max;
986 let inv_v_l = inv_normalised_vega(theta_x, s_l);
987 let r_lm = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
988 b_l, b_c, s_l, s_c, inv_v_l, inv_v_c, 0.0, false,
989 );
990 s = rational_cubic_interpolation(beta, b_l, b_c, s_l, s_c, inv_v_l, inv_v_c, r_lm);
991 s_left = s_l;
992 s_right = s_c;
993 }
994 } else {
995 let s_u = s_c + SQRT_PI_OVER_TWO * (2.0 - ome);
996 let b_u = b_u_over_b_max(s_c) * b_max;
997 if beta <= b_u {
998 let inv_v_c = SQRT_TWO_PI / b_max;
999 let inv_v_u = inv_normalised_vega(theta_x, s_u);
1000 let r_um = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
1001 b_c, b_u, s_c, s_u, inv_v_c, inv_v_u, 0.0, false,
1002 );
1003 s = rational_cubic_interpolation(beta, b_c, b_u, s_c, s_u, inv_v_c, inv_v_u, r_um);
1004 s_left = s_c;
1005 s_right = s_u;
1006 } else {
1007 let (f_upper_map_h, d_f_upper_map_h_d_beta, d2_f_upper_map_h_d_beta2) =
1008 compute_f_upper_map_and_first_two_derivatives(theta_x, s_u);
1009 let r_uu = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
1010 b_u, b_max, f_upper_map_h, 0.0, d_f_upper_map_h_d_beta, -0.5, d2_f_upper_map_h_d_beta2, true,
1011 );
1012 f = rational_cubic_interpolation(beta, b_u, b_max, f_upper_map_h, 0.0, d_f_upper_map_h_d_beta, -0.5, r_uu);
1013 if f <= 0.0 {
1014 let h = b_max - b_u;
1015 let t = (beta - b_u) / h;
1016 f = (f_upper_map_h * (1.0 - t) + 0.5 * h * t) * (1.0 - t);
1017 }
1018 s = inverse_f_upper_map(f);
1019 s_left = s_u;
1020 if beta > 0.5 * b_max {
1021 let beta_bar = b_max - beta;
1022 while iterations < n_iterations && ds.abs() > DBL_EPSILON * s {
1023 let h = theta_x / s;
1024 let t = s / 2.0;
1025 let gp = (2.0 / SQRT_TWO_PI) / (erfcx_cody((t + h) * (1.0 / SQRT_TWO)) + erfcx_cody((t - h) * (1.0 / SQRT_TWO)));
1026 let b_bar = normalised_vega(theta_x, s) / gp;
1027 if b_bar < beta_bar && s < s_right {
1028 s_right = s;
1029 } else if b_bar > beta_bar && s > s_left {
1030 s_left = s;
1031 }
1032
1033 let g = (beta_bar / b_bar).ln();
1034 let x2_over_s3 = (h * h) / s;
1035 let b_h2 = x2_over_s3 - s / 4.0;
1036 let c = 3.0 * (x2_over_s3 / s);
1037 let b_h3 = b_h2 * b_h2 - c - 0.25;
1038 let nu = -g / gp;
1039 let h2 = b_h2 + gp;
1040 let h3 = b_h3 + gp * (2.0 * gp + 3.0 * b_h2);
1041
1042 if theta_x < -580.0 {
1043 let b_h4 = b_h2 * (b_h3 - 0.5) - (b_h2 - 2.0 / s) * 2.0 * c;
1044 ds = nu * householder_factor_4(nu, h2, h3, b_h4 + gp * (6.0 * gp * (gp + 2.0 * b_h2) + 3.0 * b_h2 * b_h2 + 4.0 * b_h3));
1045 } else {
1046 ds = nu * householder_factor_3(nu, h2, h3);
1047 }
1048 s += ds;
1049 if s > s_right { s = s_right; } else if s < s_left { s = s_left; }
1050 iterations += 1;
1051 }
1052 return s;
1053 }
1054 }
1055 }
1056
1057 while iterations < n_iterations && ds.abs() > DBL_EPSILON * s {
1058 let b = normalised_black_internal(theta_x, s);
1059 let inv_bp = inv_normalised_vega(theta_x, s);
1060 let nu = (beta - b) * inv_bp;
1061 let h = theta_x / s;
1062 let x2_over_s3 = (h * h) / s;
1063 let h2 = x2_over_s3 - s * 0.25;
1064 let h3 = h2 * h2 - 3.0 * (x2_over_s3 / s) - 0.25;
1065 if b > beta && s < s_right {
1066 s_right = s;
1067 } else if b < beta && s > s_left {
1068 s_left = s;
1069 }
1070 ds = nu * householder_factor_3(nu, h2, h3);
1071 s += ds;
1072 if s > s_right { s = s_right; } else if s < s_left { s = s_left; }
1073 iterations += 1;
1074 }
1075 s
1076}
1077
1078pub fn normalised_black(x: f64, s: f64, theta: f64) -> f64 {
1079 if x == 0.0 {
1080 return erf_cody((0.5 / SQRT_TWO) * s);
1081 }
1082 let theta_x = if theta < 0.0 { -x } else { x };
1083 normalised_intrinsic(theta_x) + if s <= 0.0 { 0.0 } else { normalised_black_internal(-x.abs(), s) }
1084}
1085
1086pub fn black(f: f64, k: f64, sigma: f64, t: f64, theta: f64) -> f64 {
1087 let s = sigma * t.sqrt();
1088 if k == f {
1089 return f * erf_cody((0.5 / SQRT_TWO) * s);
1090 }
1091 let mu = if theta < 0.0 {
1092 (k - f).max(0.0)
1093 } else {
1094 (f - k).max(0.0)
1095 };
1096 mu + if s <= 0.0 {
1097 0.0
1098 } else {
1099 (f.sqrt() * k.sqrt()) * normalised_black_internal(-(f / k).ln().abs(), s)
1100 }
1101}
1102
1103pub fn implied_black_volatility(price: f64, f: f64, k: f64, t: f64, theta: f64) -> f64 {
1104 if price >= (if theta < 0.0 { k } else { f }) {
1105 return VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM;
1106 }
1107 let mu = if theta < 0.0 { k - f } else { f - k };
1108 let adjusted_price = if mu > 0.0 { price - mu } else { price };
1109 lets_be_rational(
1110 adjusted_price / (f.sqrt() * k.sqrt()),
1111 -(f / k).ln().abs(),
1112 2,
1113 ) / t.sqrt()
1114}
1115
1116pub fn normalised_implied_black_volatility(beta: f64, x: f64, theta: f64) -> f64 {
1117 let theta_x = if theta < 0.0 { -x } else { x };
1118 lets_be_rational(beta - normalised_intrinsic(theta_x), -x.abs(), 2)
1119}