1use std::f64;
8
9const SQRT_TWO: f64 = 1.4142135623730950488016887242096980785696718753769;
11const SQRT_TWO_PI: f64 = 2.5066282746310005024157652848110452530069867406099;
12const LN_TWO_PI: f64 = 1.8378770664093454835606594728112352797227949472756;
13const TWO_PI: f64 = 6.283185307179586476925286766559005768394338798750;
14const SQRT_PI_OVER_TWO: f64 = 1.253314137315500251207882642405522626503493370305;
15const SQRT_THREE: f64 = 1.732050807568877293527446341505872366942805253810;
16const SQRT_ONE_OVER_THREE: f64 = 0.577350269189625764509148780501957455647601751270;
17const TWO_PI_OVER_SQRT_TWENTY_SEVEN: f64 = 1.209199576156145233729385505094770488189377498728;
18const SQRT_THREE_OVER_THIRD_ROOT_TWO_PI: f64 = 0.938643487427383566075051356115075878414688769574;
19const PI_OVER_SIX: f64 = 0.523598775598298873077107230546583814032861566563;
20
21const DBL_EPSILON: f64 = f64::EPSILON;
22const DBL_MIN: f64 = f64::MIN_POSITIVE;
23const DBL_MAX: f64 = f64::MAX;
24
25const ETA: f64 = -13.0;
26const VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC: f64 = -f64::MAX;
27const VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM: f64 = f64::MAX;
28
29fn d_int(x: f64) -> f64 {
34 if x > 0.0 {
35 x.floor()
36 } else {
37 -(-x).floor()
38 }
39}
40
41fn smoothened_exponential_of_negative_square(y: f64) -> f64 {
42 let y_tilde = d_int(y * 16.0) / 16.0;
43 let del = (y - y_tilde) * (y + y_tilde);
44 (-y_tilde * y_tilde).exp() * (-del).exp()
45}
46
47fn smoothened_exponential_of_positive_square(x: f64) -> f64 {
48 let x_tilde = d_int(x * 16.0) / 16.0;
49 let del = (x - x_tilde) * (x + x_tilde);
50 (x_tilde * x_tilde).exp() * (del).exp()
51}
52
53const CODY_A: [f64; 5] = [
54 3.1611237438705656,
55 113.864154151050156,
56 377.485237685302021,
57 3209.37758913846947,
58 0.185777706184603153,
59];
60const CODY_B: [f64; 4] = [
61 23.6012909523441209,
62 244.024637934444173,
63 1282.61652607737228,
64 2844.23683343917062,
65];
66const CODY_C: [f64; 9] = [
67 0.564188496988670089,
68 8.88314979438837594,
69 66.1191906371416295,
70 298.635138197400131,
71 881.95222124176909,
72 1712.04761263407058,
73 2051.07837782607147,
74 1230.33935479799725,
75 2.15311535474403846E-8,
76];
77const CODY_D: [f64; 8] = [
78 15.7449261107098347,
79 117.693950891312499,
80 537.181101862009858,
81 1621.38957456669019,
82 3290.79923573345963,
83 4362.61909014324716,
84 3439.36767414372164,
85 1230.33935480374942,
86];
87const CODY_P: [f64; 6] = [
88 0.305326634961232344,
89 0.360344899949804439,
90 0.125781726111229246,
91 0.0160837851487422766,
92 6.58749161529837803E-4,
93 0.0163153871373020978,
94];
95const CODY_Q: [f64; 5] = [
96 2.56852019228982242,
97 1.87295284992346047,
98 0.527905102951428412,
99 0.0605183413124413191,
100 0.00233520497626869185,
101];
102
103fn cody_ab(z: f64) -> f64 {
104 ((((CODY_A[4] * z + CODY_A[0]) * z + CODY_A[1]) * z + CODY_A[2]) * z + CODY_A[3])
105 / ((((z + CODY_B[0]) * z + CODY_B[1]) * z + CODY_B[2]) * z + CODY_B[3])
106}
107
108fn cody_cd(y: f64) -> f64 {
109 ((((((((CODY_C[8] * y + CODY_C[0]) * y + CODY_C[1]) * y + CODY_C[2]) * y + CODY_C[3]) * y
110 + CODY_C[4])
111 * y
112 + CODY_C[5])
113 * y
114 + CODY_C[6])
115 * y
116 + CODY_C[7])
117 / ((((((((y + CODY_D[0]) * y + CODY_D[1]) * y + CODY_D[2]) * y + CODY_D[3]) * y + CODY_D[4])
118 * y
119 + CODY_D[5])
120 * y
121 + CODY_D[6])
122 * y
123 + CODY_D[7])
124}
125
126fn cody_pq(z: f64) -> f64 {
127 z * (((((CODY_P[5] * z + CODY_P[0]) * z + CODY_P[1]) * z + CODY_P[2]) * z + CODY_P[3]) * z
128 + CODY_P[4])
129 / (((((z + CODY_Q[0]) * z + CODY_Q[1]) * z + CODY_Q[2]) * z + CODY_Q[3]) * z + CODY_Q[4])
130}
131
132const ONE_OVER_SQRT_PI: f64 = 0.56418958354775628695;
133const THRESHOLD: f64 = 0.46875;
134const XNEG: f64 = -26.6287357137514;
135const XBIG: f64 = 26.543;
136
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
307#[allow(clippy::too_many_arguments)]
308fn rational_cubic_interpolation(
309 x: f64,
310 x_l: f64,
311 x_r: f64,
312 y_l: f64,
313 y_r: f64,
314 d_l: f64,
315 d_r: f64,
316 r: f64,
317) -> f64 {
318 let h = x_r - x_l;
319 if h.abs() <= 0.0 {
320 return 0.5 * (y_l + y_r);
321 }
322 let t = (x - x_l) / h;
323 if r < MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE {
324 let omt = 1.0 - t;
325 let t2 = t * t;
326 let omt2 = omt * omt;
327 return (y_r * t2 * t + (r * y_r - h * d_r) * t2 * omt + (r * y_l + h * d_l) * t * omt2
328 + y_l * omt2 * omt)
329 / (1.0 + (r - 3.0) * t * omt);
330 }
331 y_r * t + y_l * (1.0 - t)
332}
333
334fn rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
335 x_l: f64,
336 x_r: f64,
337 y_l: f64,
338 y_r: f64,
339 d_l: f64,
340 d_r: f64,
341 second_derivative_l: f64,
342) -> f64 {
343 let h = x_r - x_l;
344 let numerator = 0.5 * h * second_derivative_l + (d_r - d_l);
345 let denominator = (y_r - y_l) / h - d_l;
346 if is_zero(denominator) {
347 if numerator > 0.0 {
348 MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE
349 } else {
350 MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE
351 }
352 } else {
353 numerator / denominator
354 }
355}
356
357fn rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
358 x_l: f64,
359 x_r: f64,
360 y_l: f64,
361 y_r: f64,
362 d_l: f64,
363 d_r: f64,
364 second_derivative_r: f64,
365) -> f64 {
366 let h = x_r - x_l;
367 let numerator = 0.5 * h * second_derivative_r + (d_r - d_l);
368 let denominator = d_r - (y_r - y_l) / h;
369 if is_zero(denominator) {
370 if numerator > 0.0 {
371 MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE
372 } else {
373 MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE
374 }
375 } else {
376 numerator / denominator
377 }
378}
379
380fn minimum_rational_cubic_control_parameter(
381 d_l: f64,
382 d_r: f64,
383 s: f64,
384 prefer_shape_preservation: bool,
385) -> f64 {
386 let monotonic = d_l * s >= 0.0 && d_r * s >= 0.0;
387 let convex = d_l <= s && s <= d_r;
388 let concave = d_l >= s && s >= d_r;
389 if !monotonic && !convex && !concave {
390 return MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE;
391 }
392 let d_r_m_d_l = d_r - d_l;
393 let d_r_m_s = d_r - s;
394 let s_m_d_l = s - d_l;
395 let mut r1 = f64::MIN;
396 let mut r2 = r1;
397 if monotonic {
398 if !is_zero(s) {
399 r1 = (d_r + d_l) / s;
400 } else if prefer_shape_preservation {
401 r1 = MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE;
402 }
403 }
404 if convex || concave {
405 if !(is_zero(s_m_d_l) || is_zero(d_r_m_s)) {
406 r2 = (d_r_m_d_l / d_r_m_s).abs().max((d_r_m_d_l / s_m_d_l).abs());
407 } else if prefer_shape_preservation {
408 r2 = MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE;
409 }
410 } else if monotonic && prefer_shape_preservation {
411 r2 = MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE;
412 }
413 MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE.max(r1.max(r2))
414}
415
416#[allow(clippy::too_many_arguments)]
417fn convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
418 x_l: f64,
419 x_r: f64,
420 y_l: f64,
421 y_r: f64,
422 d_l: f64,
423 d_r: f64,
424 second_derivative_l: f64,
425 prefer_shape_preservation: bool,
426) -> f64 {
427 let r = rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
428 x_l,
429 x_r,
430 y_l,
431 y_r,
432 d_l,
433 d_r,
434 second_derivative_l,
435 );
436 let r_min = minimum_rational_cubic_control_parameter(
437 d_l,
438 d_r,
439 (y_r - y_l) / (x_r - x_l),
440 prefer_shape_preservation,
441 );
442 r.max(r_min)
443}
444
445#[allow(clippy::too_many_arguments)]
446fn convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
447 x_l: f64,
448 x_r: f64,
449 y_l: f64,
450 y_r: f64,
451 d_l: f64,
452 d_r: f64,
453 second_derivative_r: f64,
454 prefer_shape_preservation: bool,
455) -> f64 {
456 let r = rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
457 x_l,
458 x_r,
459 y_l,
460 y_r,
461 d_l,
462 d_r,
463 second_derivative_r,
464 );
465 let r_min = minimum_rational_cubic_control_parameter(
466 d_l,
467 d_r,
468 (y_r - y_l) / (x_r - x_l),
469 prefer_shape_preservation,
470 );
471 r.max(r_min)
472}
473
474fn householder3_factor(nu: f64, h2: f64, h3: f64) -> f64 {
479 (1.0 + 0.5 * h2 * nu) / (1.0 + nu * (h2 + h3 * nu * (1.0 / 6.0)))
480}
481
482fn householder4_factor(nu: f64, h2: f64, h3: f64, h4: f64) -> f64 {
483 (1.0 + nu * (h2 + nu * h3 * (1.0 / 6.0)))
484 / (1.0 + nu * (1.5 * h2 + nu * (h2 * h2 * 0.25 + h3 * (1.0 / 3.0) + nu * h4 * (1.0 / 24.0))))
485}
486
487fn householder_factor_3(nu: f64, h2: f64, h3: f64) -> f64 {
488 householder3_factor(nu, h2, h3)
489}
490
491fn householder_factor_4(nu: f64, h2: f64, h3: f64, h4: f64) -> f64 {
492 householder4_factor(nu, h2, h3, h4)
493}
494
495fn normalised_intrinsic(theta_x: f64) -> f64 {
496 if theta_x <= 0.0 {
497 0.0
498 } else {
499 (0.5 * theta_x).sinh() * 2.0
500 }
501}
502
503fn square(x: f64) -> f64 {
504 x * x
505}
506
507const SIXTEENTH_ROOT_DBL_EPSILON: f64 = 0.10511205190671433; const TAU: f64 = 2.0 * SIXTEENTH_ROOT_DBL_EPSILON;
509
510fn asymptotic_expansion_of_scaled_normalised_black(h: f64, t: f64) -> f64 {
511 let e = square(t / h);
512 let r = (h + t) * (h - t);
513 let q = square(h / r);
514
515 let a0 = 2.0;
516 let a1 = -6.0 - 2.0 * e;
517 let a2 = 30.0 + e * (60.0 + 6.0 * e);
518 let a3 = -2.1E2 + e * (-1.05E3 + e * (-6.3E2 - 30.0 * e));
519 let a4 = 1.89E3 + e * (1.764E4 + e * (2.646E4 + e * (7.56E3 + 2.1E2 * e)));
520 let a5 = -2.079E4 + e * (-3.1185E5 + e * (-8.7318E5 + e * (-6.237E5 + e * (-1.0395E5 - 1.89E3 * e))));
521 let a6 = 2.7027E5
522 + e * (5.94594E6
523 + e * (2.675673E7 + e * (3.567564E7 + e * (1.486485E7 + e * (1.62162E6 + 2.079E4 * e)))));
524 let a7 = -4.05405E6
525 + e * (-1.2297285E8
526 + e * (-8.1162081E8
527 + e * (-1.73918745E9
528 + e * (-1.35270135E9 + e * (-3.6891855E8 + e * (-2.837835E7 - 2.7027E5 * e))))));
529 let a8 = 6.891885E7
530 + e * (2.756754E9
531 + e * (2.50864614E10
532 + e * (7.88431644E10
533 + e * (9.85539555E10
534 + e * (5.01729228E10
535 + e * (9.648639E9 + e * (5.513508E8 + 4.05405E6 * e)))))));
536 let a9 = -1.30945815E9
537 + e * (-6.678236565E10
538 + e * (-8.013883878E11
539 + e * (-3.4726830138E12
540 + e * (-6.3665855253E12
541 + e * (-5.2090245207E12
542 + e * (-1.8699062382E12
543 + e * (-2.671294626E11
544 + e * (-1.178512335E10 - 6.891885E7 * e))))))));
545 let a10 = 2.749862115E10
546 + e * (1.7415793395E12
547 + e * (2.664616389435E13
548 + e * (1.52263793682E14
549 + e * (3.848890340295E14
550 + e * (4.618668408354E14
551 + e * (2.664616389435E14
552 + e * (7.10564370516E13
553 + e * (7.83710702775E12 + e * (2.749862115E11 + 1.30945815E9 * e)))))))));
554 let a11 = -6.3246828645E11
555 + e * (-4.870005805665E13
556 + e * (-9.2530110307635E14
557 + e * (-6.74147946527055E15
558 + e * (-2.24715982175685E16
559 + e * (-3.71802806872497E16
560 + e * (-3.14602375045959E16
561 + e * (-1.34829589305411E16
562 + e * (-2.77590330922905E15
563 + e * (-2.4350029028325E14
564 + e * (-6.95715115095E12 - 2.749862115E10 * e))))))))));
565 let a12 = 1.581170716125E13
566 + e * (1.454677058835E15
567 + e * (3.36030400590885E16
568 + e * (3.04027505296515E17
569 + e * (1.29211689751018875E18
570 + e * (2.81916414002223E18
571 + e * (3.289024830025935E18
572 + e * (2.067387036016302E18
573 + e * (6.8406188691715875E17
574 + e * (1.12010133530295E17
575 + e * (8.0007238235925E15
576 + e * (1.89740485935E14 + 6.3246828645E11 * e)))))))))));
577 let a13 = -4.2691609335375E14
578 + e * (-4.624924344665625E16
579 + e * (-1.2764791191277125E18
580 + e * (-1.40412703104048375E19
581 + e * (-7.41067044160255312E19
582 + e * (-2.06151377739125569E20
583 + e * (-3.17155965752500875E20
584 + e * (-2.74868503652167425E20
585 + e * (-1.33392067948845956E20
586 + e * (-3.51031757760120938E19
587 + e * (-4.6804234368016125E18
588 + e * (-2.774954606799375E17
589 + e * (-5.54990921359875E15 - 1.581170716125E13 * e))))))))))));
590 let a14 = 1.238056670725875E16
591 + e * (1.5599514051146025E18
592 + e * (5.06984206662245812E19
593 + e * (6.66322100184665925E20
594 + e * (4.27556680951827302E21
595 + e * (1.47701398874267613E22
596 + e * (2.89721974714909549E22
597 + e * (3.31110828245610914E22
598 + e * (2.2155209831140142E22
599 + e * (8.55113361903654604E21
600 + e * (1.83238577550783129E21
601 + e * (2.02793682664898325E20
602 + e * (1.01396841332449162E19
603 + e * (1.733279339016225E17 + 4.2691609335375E14 * e)))))))))))));
604 let a15 = -3.8379756792502125E17
605 + e * (-5.56506473491280812E19
606 + e * (-2.10359446979704147E21
607 + e * (-3.25556286992399275E22
608 + e * (-2.49593153360839444E23
609 + e * (-1.04829124411552567E24
610 + e * (-2.55352995361474201E24
611 + e * (-3.72085793241005264E24
612 + e * (-3.28310994036181115E24
613 + e * (-1.74715207352587611E24
614 + e * (-5.49104937393846778E23
615 + e * (-9.76668860977197826E22
616 + e * (-9.11557603578717971E21
617 + e * (-3.89554531443896569E20
618 + e * (-5.75696351887531875E18 - 1.238056670725875E16 * e))))))))))))));
619 let a16 = 1.26653197415257012E19
620 + e * (2.09399953059891594E21
621 + e * (9.10889795810528434E22
622 + e * (1.63960163245895118E24
623 + e * (1.48019591819210871E25
624 + e * (7.42789224401858187E25
625 + e * (2.19979885688242617E26
626 + e * (3.98058840769200926E26
627 + e * (4.47816195865351041E26
628 + e * (3.1425697955463231E26
629 + e * (1.36178024473674001E26
630 + e * (3.55247020366106089E25
631 + e * (5.32870530549159134E24
632 + e * (4.25081904711579936E23
633 + e * (1.57049964794918696E22
634 + e * (2.0264511586441122E20 + 3.8379756792502125E17 * e)))))))))))))));
635
636 let mut omega = 0.0;
637 let thresholds = [
638 12.347, 12.958, 13.729, 14.718, 16.016, 17.769, 20.221, 23.816, 29.419, 38.93, 57.171,
639 99.347,
640 ];
641 let val_to_check = -h - t + TAU + 0.5;
642 let idx = thresholds.iter().position(|&thr| thr > val_to_check).unwrap_or(thresholds.len());
643
644 if idx == 0 { omega = q * (a16 + omega); }
645 if idx <= 1 { omega = q * (a15 + omega); }
646 if idx <= 2 { omega = q * (a14 + omega); }
647 if idx <= 3 { omega = q * (a13 + omega); }
648 if idx <= 4 { omega = q * (a12 + omega); }
649 if idx <= 5 { omega = q * (a11 + omega); }
650 if idx <= 6 { omega = q * (a10 + omega); }
651 if idx <= 7 { omega = q * (a9 + omega); }
652 if idx <= 8 { omega = q * (a8 + omega); }
653 if idx <= 9 { omega = q * (a7 + omega); }
654 if idx <= 10 { omega = q * (a6 + omega); }
655 if idx <= 11 { omega = q * (a5 + omega); }
656
657 omega = a0 + q * (a1 + q * (a2 + q * (a3 + q * (a4 + omega))));
658
659 (t / r) * omega
660}
661
662fn yprime_tail_expansion_rational_function_part(w: f64) -> f64 {
663 w * (-2.9999999999994663866
664 + w * (-1.7556263323542206288E2
665 + w * (-3.4735035445495633334E3
666 + w * (-2.7805745693864308643E4
667 + w * (-8.3836021460741980839E4 - 6.6818249032616849037E4 * w)))))
668 / (1.0
669 + w * (6.3520877744831739102E1
670 + w * (1.4404389037604337538E3
671 + w * (1.4562545638507033944E4
672 + w * (6.6886794165651675684E4
673 + w * (1.2569970380923908488E5 + 6.9286518679803751694E4 * w))))))
674}
675
676fn yprime(h: f64) -> f64 {
677 if h < -4.0 {
678 let w = 1.0 / (h * h);
679 return w * (1.0 + yprime_tail_expansion_rational_function_part(w));
680 }
681 if h <= -0.46875 {
682 return (1.0000000000594317229
683 - h * (6.1911449879694112749E-1
684 - h * (2.2180844736576013957E-1
685 - h * (4.5650900351352987865E-2
686 - h * (5.545521007735379052E-3
687 - h * (3.0717392274913902347E-4
688 - h * (4.2766597835908713583E-8 + 8.4592436406580605619E-10 * h)))))))
689 / (1.0
690 - h * (1.8724286369589162071
691 - h * (1.5685497236077651429
692 - h * (7.6576489836589035112E-1
693 - h * (2.3677701403094640361E-1
694 - h * (4.6762548903194957675E-2
695 - h * (5.5290453576936595892E-3 - 3.0822020417927147113E-4 * h)))))));
696 }
697 1.0 + h * SQRT_PI_OVER_TWO * erfcx_cody(-(1.0 / SQRT_TWO) * h)
698}
699
700fn small_t_expansion_of_scaled_normalised_black(h: f64, t: f64) -> f64 {
701 let a = yprime(h);
702 let h2 = h * h;
703 let t2 = t * t;
704 let b0 = 2.0 * a;
705 let b1 = (-1.0 + a * (3.0 + h2)) / 3.0;
706 let b2 = (-7.0 - h2 + a * (15.0 + h2 * (10.0 + h2))) / 60.0;
707 let b3 = (-57.0 + (-18.0 - h2) * h2 + a * (105.0 + h2 * (105.0 + h2 * (21.0 + h2)))) / 2520.0;
708 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;
709 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;
710 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;
711 t * (b0 + t2 * (b1 + t2 * (b2 + t2 * (b3 + t2 * (b4 + t2 * (b5 + b6 * t2))))))
712}
713
714fn normalised_black_with_optimal_use_of_codys_functions(theta_x: f64, s: f64) -> f64 {
715 let codys_threshold = 0.46875;
716 let h = theta_x / s;
717 let t = 0.5 * s;
718 let q1 = -(1.0 / SQRT_TWO) * (h + t);
719 let q2 = -(1.0 / SQRT_TWO) * (h - t);
720 let two_b = if q1 < codys_threshold {
721 if q2 < codys_threshold {
722 (0.5 * theta_x).exp() * erfc_cody(q1) - (-0.5 * theta_x).exp() * erfc_cody(q2)
723 } else {
724 (0.5 * theta_x).exp() * erfc_cody(q1) - (-0.5 * (h * h + t * t)).exp() * erfcx_cody(q2)
725 }
726 } else {
727 if q2 < codys_threshold {
728 (-0.5 * (h * h + t * t)).exp() * erfcx_cody(q1) - (-0.5 * theta_x).exp() * erfc_cody(q2)
729 } else {
730 (-0.5 * (h * h + t * t)).exp() * (erfcx_cody(q1) - erfcx_cody(q2))
731 }
732 };
733 (0.5 * two_b).max(0.0)
734}
735
736fn normalised_vega(x: f64, s: f64) -> f64 {
737 let h = x / s;
738 let t = 0.5 * s;
739 (1.0 / SQRT_TWO_PI) * (-0.5 * (h * h + t * t)).exp()
740}
741
742fn inv_normalised_vega(x: f64, s: f64) -> f64 {
743 let h = x / s;
744 let t = 0.5 * s;
745 SQRT_TWO_PI * (0.5 * (h * h + t * t)).exp()
746}
747
748fn ln_normalised_vega(x: f64, s: f64) -> f64 {
749 let h = x / s;
750 let t = 0.5 * s;
751 -(LN_TWO_PI * 0.5) - 0.5 * (h * h + t * t)
752}
753
754fn is_region_i(theta_x: f64, s: f64) -> bool {
755 theta_x < s * ETA && s * (0.5 * s - (TAU + 0.5 + ETA)) + theta_x < 0.0
756}
757
758fn is_region_ii(theta_x: f64, s: f64) -> bool {
759 s * (s - (2.0 * TAU)) - theta_x / ETA < 0.0
760}
761
762fn normalised_black_internal(theta_x: f64, s: f64) -> f64 {
763 if is_region_i(theta_x, s) {
764 return asymptotic_expansion_of_scaled_normalised_black(theta_x / s, 0.5 * s)
765 * normalised_vega(theta_x, s);
766 }
767 if is_region_ii(theta_x, s) {
768 return small_t_expansion_of_scaled_normalised_black(theta_x / s, 0.5 * s)
769 * normalised_vega(theta_x, s);
770 }
771 normalised_black_with_optimal_use_of_codys_functions(theta_x, s)
772}
773
774fn scaled_normalised_black_and_ln_vega_internal(theta_x: f64, s: f64) -> (f64, f64) {
775 if is_region_i(theta_x, s) {
776 return (
777 asymptotic_expansion_of_scaled_normalised_black(theta_x / s, 0.5 * s),
778 ln_normalised_vega(theta_x, s),
779 );
780 }
781 if is_region_ii(theta_x, s) {
782 return (
783 small_t_expansion_of_scaled_normalised_black(theta_x / s, 0.5 * s),
784 ln_normalised_vega(theta_x, s),
785 );
786 }
787 let ln_vega = ln_normalised_vega(theta_x, s);
788 (
789 normalised_black_with_optimal_use_of_codys_functions(theta_x, s) * (-ln_vega).exp(),
790 ln_vega,
791 )
792}
793
794fn compute_f_lower_map_and_first_two_derivatives(x: f64, s: f64) -> (f64, f64, f64) {
795 let ax = x.abs();
796 let z = SQRT_ONE_OVER_THREE * ax / s;
797 let y = z * z;
798 let s2 = s * s;
799 let phi_z = 0.5 * erfc_cody((1.0 / SQRT_TWO) * z); let pdf_z = norm_pdf(z);
801 let fpp = PI_OVER_SIX * y / (s2 * s) * phi_z
802 * (8.0 * SQRT_THREE * s * ax + (3.0 * s2 * (s2 - 8.0) - 8.0 * x * x) * phi_z / pdf_z)
803 * (2.0 * y + 0.25 * s2).exp();
804 let phi_z2 = phi_z * phi_z;
805 let fp = TWO_PI * y * phi_z2 * (y + 0.125 * s2).exp();
806 let f = TWO_PI_OVER_SQRT_TWENTY_SEVEN * ax * (phi_z2 * phi_z);
807 (f, fp, fpp)
808}
809
810fn inverse_f_lower_map(x: f64, f: f64) -> f64 {
811 (x / (SQRT_THREE
812 * inverse_norm_cdf(
813 SQRT_THREE_OVER_THIRD_ROOT_TWO_PI * f.cbrt() / x.abs().cbrt(),
814 )))
815 .abs()
816}
817
818fn compute_f_upper_map_and_first_two_derivatives(x: f64, s: f64) -> (f64, f64, f64) {
819 let f = 0.5 * erfc_cody((0.5 / SQRT_TWO) * s); let w = square(x / s);
821 let fp = -0.5 * (0.5 * w).exp();
822 let fpp = SQRT_PI_OVER_TWO * (w + 0.125 * s * s).exp() * w / s;
823 (f, fp, fpp)
824}
825
826fn inverse_f_upper_map(f: f64) -> f64 {
827 -2.0 * inverse_norm_cdf(f)
828}
829
830fn one_minus_erfcx(x: f64) -> f64 {
831 if !(-1.0 / 5.0..=1.0 / 3.0).contains(&x) {
832 return 1.0 - erfcx_cody(x);
833 }
834 x * (1.128379167095512573896
835 - x * (1.0000000000000002
836 + x * (1.1514967181784756
837 + x * (5.7689001208873741E-1
838 + x * (1.4069188744609651E-1 + 1.4069285713634565E-2 * x))))
839 / (1.0
840 + x * (1.9037494962421563
841 + x * (1.5089908593742723
842 + x * (6.2486081658640257E-1
843 + x * (1.358008134514386E-1 + 1.2463320728346347E-2 * x))))))
844}
845
846fn implied_normalised_volatility_atm(beta: f64) -> f64 {
847 let beta_max = 0.6826894921370859;
848 if beta <= beta_max {
849 let r = beta_max * beta_max - beta * beta;
850 return beta
851 * ((2.92958954698308816
852 + r * (1.4014698674754995E1
853 + r * (2.44918990556468762E1
854 + r * (1.90763928424894996E1
855 + r * (6.43250149461895996
856 + r * (7.52328633671821543E-1 + 1.38781536163865582E-2 * r))))))
857 / (1.0
858 + r * (5.22443271807813073
859 + r * (1.02258209975070629E1
860 + r * (9.28187483709036392
861 + r * (3.9095549184069553
862 + r * (6.61214199809055912E-1 + 2.89411828874884851E-2 * r)))))));
863 }
864 -2.0 * inverse_norm_cdf_for_low_probabilities(0.5 * (1.0 - beta))
865}
866
867fn b_l_over_b_max(s_c: f64) -> f64 {
868 if s_c >= 2.6267851073127395 {
869 if s_c >= 7.348469228349534 {
870 let s_c_val = s_c;
871 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))))));
872 }
873 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))))));
874 }
875 if s_c >= 0.7099295739719539 {
876 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))))));
877 }
878 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))));
879 (s_c * s_c) * (0.07560996640296361767172 + s_c * (s_c * g - 0.09672719281339436290858))
880}
881
882fn b_u_over_b_max(s_c: f64) -> f64 {
883 if s_c > 1.7888543819998317 {
884 if s_c > 6.164414002968976 {
885 let s_c_val = s_c;
886 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))))));
887 }
888 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))))));
889 }
890 if s_c >= 0.7745966692414833 {
891 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)))));
892 }
893 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))));
894 0.7899085945560627246288 + (s_c * s_c) * (0.0614616805805147403487 + s_c * g)
895}
896
897pub fn lets_be_rational(beta: f64, theta_x: f64, n_iterations: i32) -> f64 {
898 if beta <= 0.0 {
899 return if beta == 0.0 { 0.0 } else { VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC };
900 }
901 let b_max = (0.5 * theta_x).exp();
902 if beta >= b_max {
903 return VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM;
904 }
905 if theta_x == 0.0 {
906 return implied_normalised_volatility_atm(beta);
907 }
908
909 let mut iterations = 0;
910 let mut f;
911 let mut s;
912 let mut ds = -DBL_MAX;
913 let mut s_left = DBL_MIN;
914 let mut s_right = DBL_MAX;
915
916 let sqrt_ax = (-theta_x).sqrt();
917 let s_c = SQRT_TWO * sqrt_ax;
918 let ome = one_minus_erfcx(sqrt_ax);
919 let b_c = 0.5 * b_max * ome;
920
921 if beta < b_c {
922 let s_l = s_c - SQRT_PI_OVER_TWO * ome;
923 let b_l = b_l_over_b_max(s_c) * b_max;
924 if beta < b_l {
925 let (f_lower_map_l, d_f_lower_map_l_d_beta, d2_f_lower_map_l_d_beta2) =
926 compute_f_lower_map_and_first_two_derivatives(theta_x, s_l);
927 let r_ll = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
928 0.0,
929 b_l,
930 0.0,
931 f_lower_map_l,
932 1.0,
933 d_f_lower_map_l_d_beta,
934 d2_f_lower_map_l_d_beta2,
935 true,
936 );
937 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);
938 if f <= 0.0 {
939 let t = beta / b_l;
940 f = (f_lower_map_l * t + b_l * (1.0 - t)) * t;
941 }
942 s = inverse_f_lower_map(theta_x, f);
943 s_right = s_l;
944
945 let ln_beta = beta.ln();
946 while iterations < n_iterations && ds.abs() > DBL_EPSILON * s {
947 let (bx, ln_vega) = scaled_normalised_black_and_ln_vega_internal(theta_x, s);
948 let ln_b = bx.ln() + ln_vega;
949 let bpob = 1.0 / bx;
950 let b = ln_b.exp();
951 if b > beta && s < s_right {
952 s_right = s;
953 } else if b < beta && s > s_left {
954 s_left = s;
955 }
956
957 let h = theta_x / s;
958 let x2_over_s3 = h * h / s;
959 let b_h2 = x2_over_s3 - s / 4.0;
960 let nu = (ln_beta - ln_b) * ln_b / ln_beta / bpob;
961 let lambda = 1.0 / ln_b;
962 let otl = 1.0 + 2.0 * lambda;
963 let h2 = b_h2 - bpob * otl;
964 let c = 3.0 * (x2_over_s3 / s);
965 let b_h3 = b_h2 * b_h2 - c - 0.25;
966 let mu = 6.0 * lambda * (1.0 + lambda);
967 let h3 = b_h3 + (bpob * bpob) * (2.0 + mu) - (b_h2 * bpob) * 3.0 * otl;
968
969 if theta_x < -190.0 {
970 let b_h4 = b_h2 * (b_h3 - 0.5) - (b_h2 - 2.0 / s) * 2.0 * c;
971 let bpob2 = bpob * bpob;
972 let bppob = b_h2 * bpob;
973 ds = nu * householder_factor_4(
974 nu,
975 h2,
976 h3,
977 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,
978 );
979 } else {
980 ds = nu * householder_factor_3(nu, h2, h3);
981 }
982 s += ds;
983 if s > s_right { s = s_right; } else if s < s_left { s = s_left; }
984 iterations += 1;
985 }
986 return s;
987 } else {
988 let inv_v_c = SQRT_TWO_PI / b_max;
989 let inv_v_l = inv_normalised_vega(theta_x, s_l);
990 let r_lm = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
991 b_l, b_c, s_l, s_c, inv_v_l, inv_v_c, 0.0, false,
992 );
993 s = rational_cubic_interpolation(beta, b_l, b_c, s_l, s_c, inv_v_l, inv_v_c, r_lm);
994 s_left = s_l;
995 s_right = s_c;
996 }
997 } else {
998 let s_u = s_c + SQRT_PI_OVER_TWO * (2.0 - ome);
999 let b_u = b_u_over_b_max(s_c) * b_max;
1000 if beta <= b_u {
1001 let inv_v_c = SQRT_TWO_PI / b_max;
1002 let inv_v_u = inv_normalised_vega(theta_x, s_u);
1003 let r_um = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
1004 b_c, b_u, s_c, s_u, inv_v_c, inv_v_u, 0.0, false,
1005 );
1006 s = rational_cubic_interpolation(beta, b_c, b_u, s_c, s_u, inv_v_c, inv_v_u, r_um);
1007 s_left = s_c;
1008 s_right = s_u;
1009 } else {
1010 let (f_upper_map_h, d_f_upper_map_h_d_beta, d2_f_upper_map_h_d_beta2) =
1011 compute_f_upper_map_and_first_two_derivatives(theta_x, s_u);
1012 let r_uu = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
1013 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,
1014 );
1015 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);
1016 if f <= 0.0 {
1017 let h = b_max - b_u;
1018 let t = (beta - b_u) / h;
1019 f = (f_upper_map_h * (1.0 - t) + 0.5 * h * t) * (1.0 - t);
1020 }
1021 s = inverse_f_upper_map(f);
1022 s_left = s_u;
1023 if beta > 0.5 * b_max {
1024 let beta_bar = b_max - beta;
1025 while iterations < n_iterations && ds.abs() > DBL_EPSILON * s {
1026 let h = theta_x / s;
1027 let t = s / 2.0;
1028 let gp = (2.0 / SQRT_TWO_PI) / (erfcx_cody((t + h) * (1.0 / SQRT_TWO)) + erfcx_cody((t - h) * (1.0 / SQRT_TWO)));
1029 let b_bar = normalised_vega(theta_x, s) / gp;
1030 if b_bar < beta_bar && s < s_right {
1031 s_right = s;
1032 } else if b_bar > beta_bar && s > s_left {
1033 s_left = s;
1034 }
1035
1036 let g = (beta_bar / b_bar).ln();
1037 let x2_over_s3 = (h * h) / s;
1038 let b_h2 = x2_over_s3 - s / 4.0;
1039 let c = 3.0 * (x2_over_s3 / s);
1040 let b_h3 = b_h2 * b_h2 - c - 0.25;
1041 let nu = -g / gp;
1042 let h2 = b_h2 + gp;
1043 let h3 = b_h3 + gp * (2.0 * gp + 3.0 * b_h2);
1044
1045 if theta_x < -580.0 {
1046 let b_h4 = b_h2 * (b_h3 - 0.5) - (b_h2 - 2.0 / s) * 2.0 * c;
1047 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));
1048 } else {
1049 ds = nu * householder_factor_3(nu, h2, h3);
1050 }
1051 s += ds;
1052 if s > s_right { s = s_right; } else if s < s_left { s = s_left; }
1053 iterations += 1;
1054 }
1055 return s;
1056 }
1057 }
1058 }
1059
1060 while iterations < n_iterations && ds.abs() > DBL_EPSILON * s {
1061 let b = normalised_black_internal(theta_x, s);
1062 let inv_bp = inv_normalised_vega(theta_x, s);
1063 let nu = (beta - b) * inv_bp;
1064 let h = theta_x / s;
1065 let x2_over_s3 = (h * h) / s;
1066 let h2 = x2_over_s3 - s * 0.25;
1067 let h3 = h2 * h2 - 3.0 * (x2_over_s3 / s) - 0.25;
1068 if b > beta && s < s_right {
1069 s_right = s;
1070 } else if b < beta && s > s_left {
1071 s_left = s;
1072 }
1073 ds = nu * householder_factor_3(nu, h2, h3);
1074 s += ds;
1075 if s > s_right { s = s_right; } else if s < s_left { s = s_left; }
1076 iterations += 1;
1077 }
1078 s
1079}
1080
1081pub fn normalised_black(x: f64, s: f64, theta: f64) -> f64 {
1082 if x == 0.0 {
1083 return erf_cody((0.5 / SQRT_TWO) * s);
1084 }
1085 let theta_x = if theta < 0.0 { -x } else { x };
1086 normalised_intrinsic(theta_x) + if s <= 0.0 { 0.0 } else { normalised_black_internal(-x.abs(), s) }
1087}
1088
1089pub fn black(f: f64, k: f64, sigma: f64, t: f64, theta: f64) -> f64 {
1090 let s = sigma * t.sqrt();
1091 if k == f {
1092 return f * erf_cody((0.5 / SQRT_TWO) * s);
1093 }
1094 let mu = if theta < 0.0 {
1095 (k - f).max(0.0)
1096 } else {
1097 (f - k).max(0.0)
1098 };
1099 mu + if s <= 0.0 {
1100 0.0
1101 } else {
1102 (f.sqrt() * k.sqrt()) * normalised_black_internal(-(f / k).ln().abs(), s)
1103 }
1104}
1105
1106pub fn implied_black_volatility(price: f64, f: f64, k: f64, t: f64, theta: f64) -> f64 {
1107 if price >= (if theta < 0.0 { k } else { f }) {
1108 return VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM;
1109 }
1110 let mu = if theta < 0.0 { k - f } else { f - k };
1111 let adjusted_price = if mu > 0.0 { price - mu } else { price };
1112 lets_be_rational(
1113 adjusted_price / (f.sqrt() * k.sqrt()),
1114 -(f / k).ln().abs(),
1115 2,
1116 ) / t.sqrt()
1117}
1118
1119pub fn normalised_implied_black_volatility(beta: f64, x: f64, theta: f64) -> f64 {
1120 let theta_x = if theta < 0.0 { -x } else { x };
1121 lets_be_rational(beta - normalised_intrinsic(theta_x), -x.abs(), 2)
1122}