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;
135
136
137pub fn erfc_cody(x: f64) -> f64 {
138 let y = x.abs();
139 if y <= THRESHOLD {
140 return 1.0 - x * cody_ab(y * y);
141 }
142 let erfc_abs_x = if y >= XBIG {
143 0.0
144 } else {
145 let val = if y <= 4.0 {
146 cody_cd(y)
147 } else {
148 (ONE_OVER_SQRT_PI - cody_pq(1.0 / (y * y))) / y
149 };
150 val * smoothened_exponential_of_negative_square(y)
151 };
152 if x < 0.0 {
153 2.0 - erfc_abs_x
154 } else {
155 erfc_abs_x
156 }
157}
158
159pub fn erf_cody(x: f64) -> f64 {
160 let y = x.abs();
161 if y <= THRESHOLD {
162 return x * cody_ab(y * y);
163 }
164 let erfc_abs_x = if y >= XBIG {
165 0.0
166 } else {
167 let val = if y <= 4.0 {
168 cody_cd(y)
169 } else {
170 (ONE_OVER_SQRT_PI - cody_pq(1.0 / (y * y))) / y
171 };
172 val * smoothened_exponential_of_negative_square(y)
173 };
174 if x < 0.0 {
175 erfc_abs_x - 1.0
176 } else {
177 1.0 - erfc_abs_x
178 }
179}
180
181pub fn erfcx_cody(x: f64) -> f64 {
182 let y = x.abs();
183 if y <= THRESHOLD {
184 let z = y * y;
185 return z.exp() * (1.0 - x * cody_ab(z));
186 }
187 if x < XNEG {
188 return f64::MAX;
189 }
190 let result = if y <= 4.0 {
191 cody_cd(y)
192 } else {
193 (ONE_OVER_SQRT_PI - cody_pq(1.0 / (y * y))) / y
194 };
195 if x < 0.0 {
196 let expx2 = smoothened_exponential_of_positive_square(x);
197 (expx2 + expx2) - result
198 } else {
199 result
200 }
201}
202
203pub fn norm_pdf(x: f64) -> f64 {
208 (1.0 / SQRT_TWO_PI) * (-0.5 * x * x).exp()
209}
210
211pub fn norm_cdf(z: f64) -> f64 {
212 if z <= -10.0 {
213 let mut sum = 1.0;
214 if z >= -1.0 / f64::EPSILON.sqrt() {
215 let zsqr = z * z;
216 let mut i = 1.0;
217 let mut g = 1.0;
218 let mut a = f64::MAX;
219 let mut lasta;
220 loop {
221 lasta = a;
222 let x = (4.0 * i - 3.0) / zsqr;
223 let y = x * ((4.0 * i - 1.0) / zsqr);
224 a = g * (x - y);
225 sum -= a;
226 g *= y;
227 i += 1.0;
228 a = a.abs();
229 if !(lasta > a && a >= (sum * f64::EPSILON).abs()) {
230 break;
231 }
232 }
233 }
234 return -norm_pdf(z) * sum / z;
235 }
236 0.5 * erfc_cody(-z * (1.0 / SQRT_TWO))
237}
238
239fn inverse_norm_cdf_for_low_probabilities(p: f64) -> f64 {
244 let r = (-p.ln()).sqrt();
245 if r < 6.7 {
246 if r < 3.41 {
247 if r < 2.05 {
248 (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)))))
249 } else {
250 (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)))))
251 }
252 } else {
253 (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)))))
254 }
255 } else {
256 if r < 12.9 {
257 (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)))))
258 } else {
259 (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)))))
260 }
261 }
262}
263
264const U_MAX: f64 = 0.3413447460685429;
265
266fn inverse_norm_cdfm05_for_midrange_probabilities(u: f64) -> f64 {
267 let s = U_MAX * U_MAX - u * u;
268 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))))))
269}
270
271pub fn inverse_norm_cdf(p: f64) -> f64 {
272 let u = p - 0.5;
273 if u.abs() < U_MAX {
274 return inverse_norm_cdfm05_for_midrange_probabilities(u);
275 }
276 if u > 0.0 {
277 -inverse_norm_cdf_for_low_probabilities(1.0 - p)
278 } else {
279 inverse_norm_cdf_for_low_probabilities(p)
280 }
281}
282
283pub fn erfinv(e: f64) -> f64 {
284 if e.abs() < (2.0 * U_MAX) {
285 return inverse_norm_cdfm05_for_midrange_probabilities(0.5 * e) * (1.0 / SQRT_TWO);
286 }
287 let val = if e < 0.0 {
288 inverse_norm_cdf_for_low_probabilities(0.5 * e + 0.5)
289 } else {
290 -inverse_norm_cdf_for_low_probabilities(-0.5 * e + 0.5)
291 };
292 val * (1.0 / SQRT_TWO)
293}
294
295const 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);
301
302fn is_zero(x: f64) -> bool {
303 x.abs() < f64::MIN_POSITIVE
304}
305
306fn rational_cubic_interpolation(
307 x: f64,
308 x_l: f64,
309 x_r: f64,
310 y_l: f64,
311 y_r: f64,
312 d_l: f64,
313 d_r: f64,
314 r: f64,
315) -> f64 {
316 let h = x_r - x_l;
317 if h.abs() <= 0.0 {
318 return 0.5 * (y_l + y_r);
319 }
320 let t = (x - x_l) / h;
321 if !(r >= MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE) {
322 let omt = 1.0 - t;
323 let t2 = t * t;
324 let omt2 = omt * omt;
325 return (y_r * t2 * t + (r * y_r - h * d_r) * t2 * omt + (r * y_l + h * d_l) * t * omt2
326 + y_l * omt2 * omt)
327 / (1.0 + (r - 3.0) * t * omt);
328 }
329 y_r * t + y_l * (1.0 - t)
330}
331
332fn rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
333 x_l: f64,
334 x_r: f64,
335 y_l: f64,
336 y_r: f64,
337 d_l: f64,
338 d_r: f64,
339 second_derivative_l: f64,
340) -> f64 {
341 let h = x_r - x_l;
342 let numerator = 0.5 * h * second_derivative_l + (d_r - d_l);
343 let denominator = (y_r - y_l) / h - d_l;
344 if is_zero(denominator) {
345 if numerator > 0.0 {
346 MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE
347 } else {
348 MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE
349 }
350 } else {
351 numerator / denominator
352 }
353}
354
355fn rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
356 x_l: f64,
357 x_r: f64,
358 y_l: f64,
359 y_r: f64,
360 d_l: f64,
361 d_r: f64,
362 second_derivative_r: f64,
363) -> f64 {
364 let h = x_r - x_l;
365 let numerator = 0.5 * h * second_derivative_r + (d_r - d_l);
366 let denominator = d_r - (y_r - y_l) / h;
367 if is_zero(denominator) {
368 if numerator > 0.0 {
369 MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE
370 } else {
371 MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE
372 }
373 } else {
374 numerator / denominator
375 }
376}
377
378fn minimum_rational_cubic_control_parameter(
379 d_l: f64,
380 d_r: f64,
381 s: f64,
382 prefer_shape_preservation: bool,
383) -> f64 {
384 let monotonic = d_l * s >= 0.0 && d_r * s >= 0.0;
385 let convex = d_l <= s && s <= d_r;
386 let concave = d_l >= s && s >= d_r;
387 if !monotonic && !convex && !concave {
388 return MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE;
389 }
390 let d_r_m_d_l = d_r - d_l;
391 let d_r_m_s = d_r - s;
392 let s_m_d_l = s - d_l;
393 let mut r1 = f64::MIN;
394 let mut r2 = r1;
395 if monotonic {
396 if !is_zero(s) {
397 r1 = (d_r + d_l) / s;
398 } else if prefer_shape_preservation {
399 r1 = MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE;
400 }
401 }
402 if convex || concave {
403 if !(is_zero(s_m_d_l) || is_zero(d_r_m_s)) {
404 r2 = (d_r_m_d_l / d_r_m_s).abs().max((d_r_m_d_l / s_m_d_l).abs());
405 } else if prefer_shape_preservation {
406 r2 = MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE;
407 }
408 } else if monotonic && prefer_shape_preservation {
409 r2 = MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE;
410 }
411 MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE.max(r1.max(r2))
412}
413
414fn convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
415 x_l: f64,
416 x_r: f64,
417 y_l: f64,
418 y_r: f64,
419 d_l: f64,
420 d_r: f64,
421 second_derivative_l: f64,
422 prefer_shape_preservation: bool,
423) -> f64 {
424 let r = rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
425 x_l,
426 x_r,
427 y_l,
428 y_r,
429 d_l,
430 d_r,
431 second_derivative_l,
432 );
433 let r_min = minimum_rational_cubic_control_parameter(
434 d_l,
435 d_r,
436 (y_r - y_l) / (x_r - x_l),
437 prefer_shape_preservation,
438 );
439 r.max(r_min)
440}
441
442fn convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
443 x_l: f64,
444 x_r: f64,
445 y_l: f64,
446 y_r: f64,
447 d_l: f64,
448 d_r: f64,
449 second_derivative_r: f64,
450 prefer_shape_preservation: bool,
451) -> f64 {
452 let r = rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
453 x_l,
454 x_r,
455 y_l,
456 y_r,
457 d_l,
458 d_r,
459 second_derivative_r,
460 );
461 let r_min = minimum_rational_cubic_control_parameter(
462 d_l,
463 d_r,
464 (y_r - y_l) / (x_r - x_l),
465 prefer_shape_preservation,
466 );
467 r.max(r_min)
468}
469
470fn householder3_factor(nu: f64, h2: f64, h3: f64) -> f64 {
475 (1.0 + 0.5 * h2 * nu) / (1.0 + nu * (h2 + h3 * nu * (1.0 / 6.0)))
476}
477
478fn householder4_factor(nu: f64, h2: f64, h3: f64, h4: f64) -> f64 {
479 (1.0 + nu * (h2 + nu * h3 * (1.0 / 6.0)))
480 / (1.0 + nu * (1.5 * h2 + nu * (h2 * h2 * 0.25 + h3 * (1.0 / 3.0) + nu * h4 * (1.0 / 24.0))))
481}
482
483fn householder_factor_3(nu: f64, h2: f64, h3: f64) -> f64 {
484 householder3_factor(nu, h2, h3)
485}
486
487fn householder_factor_4(nu: f64, h2: f64, h3: f64, h4: f64) -> f64 {
488 householder4_factor(nu, h2, h3, h4)
489}
490
491fn normalised_intrinsic(theta_x: f64) -> f64 {
492 if theta_x <= 0.0 {
493 0.0
494 } else {
495 (0.5 * theta_x).sinh() * 2.0
496 }
497}
498
499fn square(x: f64) -> f64 {
500 x * x
501}
502
503const SIXTEENTH_ROOT_DBL_EPSILON: f64 = 0.10511205190671433; const TAU: f64 = 2.0 * SIXTEENTH_ROOT_DBL_EPSILON;
505
506fn asymptotic_expansion_of_scaled_normalised_black(h: f64, t: f64) -> f64 {
507 let e = square(t / h);
508 let r = (h + t) * (h - t);
509 let q = square(h / r);
510
511 let a0 = 2.0;
512 let a1 = -6.0 - 2.0 * e;
513 let a2 = 30.0 + e * (60.0 + 6.0 * e);
514 let a3 = -2.1E2 + e * (-1.05E3 + e * (-6.3E2 - 30.0 * e));
515 let a4 = 1.89E3 + e * (1.764E4 + e * (2.646E4 + e * (7.56E3 + 2.1E2 * e)));
516 let a5 = -2.079E4 + e * (-3.1185E5 + e * (-8.7318E5 + e * (-6.237E5 + e * (-1.0395E5 - 1.89E3 * e))));
517 let a6 = 2.7027E5
518 + e * (5.94594E6
519 + e * (2.675673E7 + e * (3.567564E7 + e * (1.486485E7 + e * (1.62162E6 + 2.079E4 * e)))));
520 let a7 = -4.05405E6
521 + e * (-1.2297285E8
522 + e * (-8.1162081E8
523 + e * (-1.73918745E9
524 + e * (-1.35270135E9 + e * (-3.6891855E8 + e * (-2.837835E7 - 2.7027E5 * e))))));
525 let a8 = 6.891885E7
526 + e * (2.756754E9
527 + e * (2.50864614E10
528 + e * (7.88431644E10
529 + e * (9.85539555E10
530 + e * (5.01729228E10
531 + e * (9.648639E9 + e * (5.513508E8 + 4.05405E6 * e)))))));
532 let a9 = -1.30945815E9
533 + e * (-6.678236565E10
534 + e * (-8.013883878E11
535 + e * (-3.4726830138E12
536 + e * (-6.3665855253E12
537 + e * (-5.2090245207E12
538 + e * (-1.8699062382E12
539 + e * (-2.671294626E11
540 + e * (-1.178512335E10 - 6.891885E7 * e))))))));
541 let a10 = 2.749862115E10
542 + e * (1.7415793395E12
543 + e * (2.664616389435E13
544 + e * (1.52263793682E14
545 + e * (3.848890340295E14
546 + e * (4.618668408354E14
547 + e * (2.664616389435E14
548 + e * (7.10564370516E13
549 + e * (7.83710702775E12 + e * (2.749862115E11 + 1.30945815E9 * e)))))))));
550 let a11 = -6.3246828645E11
551 + e * (-4.870005805665E13
552 + e * (-9.2530110307635E14
553 + e * (-6.74147946527055E15
554 + e * (-2.24715982175685E16
555 + e * (-3.71802806872497E16
556 + e * (-3.14602375045959E16
557 + e * (-1.34829589305411E16
558 + e * (-2.77590330922905E15
559 + e * (-2.4350029028325E14
560 + e * (-6.95715115095E12 - 2.749862115E10 * e))))))))));
561 let a12 = 1.581170716125E13
562 + e * (1.454677058835E15
563 + e * (3.36030400590885E16
564 + e * (3.04027505296515E17
565 + e * (1.29211689751018875E18
566 + e * (2.81916414002223E18
567 + e * (3.289024830025935E18
568 + e * (2.067387036016302E18
569 + e * (6.8406188691715875E17
570 + e * (1.12010133530295E17
571 + e * (8.0007238235925E15
572 + e * (1.89740485935E14 + 6.3246828645E11 * e)))))))))));
573 let a13 = -4.2691609335375E14
574 + e * (-4.624924344665625E16
575 + e * (-1.2764791191277125E18
576 + e * (-1.40412703104048375E19
577 + e * (-7.41067044160255312E19
578 + e * (-2.06151377739125569E20
579 + e * (-3.17155965752500875E20
580 + e * (-2.74868503652167425E20
581 + e * (-1.33392067948845956E20
582 + e * (-3.51031757760120938E19
583 + e * (-4.6804234368016125E18
584 + e * (-2.774954606799375E17
585 + e * (-5.54990921359875E15 - 1.581170716125E13 * e))))))))))));
586 let a14 = 1.238056670725875E16
587 + e * (1.5599514051146025E18
588 + e * (5.06984206662245812E19
589 + e * (6.66322100184665925E20
590 + e * (4.27556680951827302E21
591 + e * (1.47701398874267613E22
592 + e * (2.89721974714909549E22
593 + e * (3.31110828245610914E22
594 + e * (2.2155209831140142E22
595 + e * (8.55113361903654604E21
596 + e * (1.83238577550783129E21
597 + e * (2.02793682664898325E20
598 + e * (1.01396841332449162E19
599 + e * (1.733279339016225E17 + 4.2691609335375E14 * e)))))))))))));
600 let a15 = -3.8379756792502125E17
601 + e * (-5.56506473491280812E19
602 + e * (-2.10359446979704147E21
603 + e * (-3.25556286992399275E22
604 + e * (-2.49593153360839444E23
605 + e * (-1.04829124411552567E24
606 + e * (-2.55352995361474201E24
607 + e * (-3.72085793241005264E24
608 + e * (-3.28310994036181115E24
609 + e * (-1.74715207352587611E24
610 + e * (-5.49104937393846778E23
611 + e * (-9.76668860977197826E22
612 + e * (-9.11557603578717971E21
613 + e * (-3.89554531443896569E20
614 + e * (-5.75696351887531875E18 - 1.238056670725875E16 * e))))))))))))));
615 let a16 = 1.26653197415257012E19
616 + e * (2.09399953059891594E21
617 + e * (9.10889795810528434E22
618 + e * (1.63960163245895118E24
619 + e * (1.48019591819210871E25
620 + e * (7.42789224401858187E25
621 + e * (2.19979885688242617E26
622 + e * (3.98058840769200926E26
623 + e * (4.47816195865351041E26
624 + e * (3.1425697955463231E26
625 + e * (1.36178024473674001E26
626 + e * (3.55247020366106089E25
627 + e * (5.32870530549159134E24
628 + e * (4.25081904711579936E23
629 + e * (1.57049964794918696E22
630 + e * (2.0264511586441122E20 + 3.8379756792502125E17 * e)))))))))))))));
631
632 let mut omega = 0.0;
633 let thresholds = [
634 12.347, 12.958, 13.729, 14.718, 16.016, 17.769, 20.221, 23.816, 29.419, 38.93, 57.171,
635 99.347,
636 ];
637 let val_to_check = -h - t + TAU + 0.5;
638 let idx = thresholds.iter().position(|&thr| thr > val_to_check).unwrap_or(thresholds.len());
639
640 if idx <= 0 { omega = q * (a16 + omega); }
641 if idx <= 1 { omega = q * (a15 + omega); }
642 if idx <= 2 { omega = q * (a14 + omega); }
643 if idx <= 3 { omega = q * (a13 + omega); }
644 if idx <= 4 { omega = q * (a12 + omega); }
645 if idx <= 5 { omega = q * (a11 + omega); }
646 if idx <= 6 { omega = q * (a10 + omega); }
647 if idx <= 7 { omega = q * (a9 + omega); }
648 if idx <= 8 { omega = q * (a8 + omega); }
649 if idx <= 9 { omega = q * (a7 + omega); }
650 if idx <= 10 { omega = q * (a6 + omega); }
651 if idx <= 11 { omega = q * (a5 + omega); }
652
653 omega = a0 + q * (a1 + q * (a2 + q * (a3 + q * (a4 + omega))));
654
655 (t / r) * omega
656}
657
658fn yprime_tail_expansion_rational_function_part(w: f64) -> f64 {
659 w * (-2.9999999999994663866
660 + w * (-1.7556263323542206288E2
661 + w * (-3.4735035445495633334E3
662 + w * (-2.7805745693864308643E4
663 + w * (-8.3836021460741980839E4 - 6.6818249032616849037E4 * w)))))
664 / (1.0
665 + w * (6.3520877744831739102E1
666 + w * (1.4404389037604337538E3
667 + w * (1.4562545638507033944E4
668 + w * (6.6886794165651675684E4
669 + w * (1.2569970380923908488E5 + 6.9286518679803751694E4 * w))))))
670}
671
672fn yprime(h: f64) -> f64 {
673 if h < -4.0 {
674 let w = 1.0 / (h * h);
675 return w * (1.0 + yprime_tail_expansion_rational_function_part(w));
676 }
677 if h <= -0.46875 {
678 return (1.0000000000594317229
679 - h * (6.1911449879694112749E-1
680 - h * (2.2180844736576013957E-1
681 - h * (4.5650900351352987865E-2
682 - h * (5.545521007735379052E-3
683 - h * (3.0717392274913902347E-4
684 - h * (4.2766597835908713583E-8 + 8.4592436406580605619E-10 * h)))))))
685 / (1.0
686 - h * (1.8724286369589162071
687 - h * (1.5685497236077651429
688 - h * (7.6576489836589035112E-1
689 - h * (2.3677701403094640361E-1
690 - h * (4.6762548903194957675E-2
691 - h * (5.5290453576936595892E-3 - 3.0822020417927147113E-4 * h)))))));
692 }
693 1.0 + h * SQRT_PI_OVER_TWO * erfcx_cody(-(1.0 / SQRT_TWO) * h)
694}
695
696fn small_t_expansion_of_scaled_normalised_black(h: f64, t: f64) -> f64 {
697 let a = yprime(h);
698 let h2 = h * h;
699 let t2 = t * t;
700 let b0 = 2.0 * a;
701 let b1 = (-1.0 + a * (3.0 + h2)) / 3.0;
702 let b2 = (-7.0 - h2 + a * (15.0 + h2 * (10.0 + h2))) / 60.0;
703 let b3 = (-57.0 + (-18.0 - h2) * h2 + a * (105.0 + h2 * (105.0 + h2 * (21.0 + h2)))) / 2520.0;
704 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;
705 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;
706 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;
707 t * (b0 + t2 * (b1 + t2 * (b2 + t2 * (b3 + t2 * (b4 + t2 * (b5 + b6 * t2))))))
708}
709
710fn normalised_black_with_optimal_use_of_codys_functions(theta_x: f64, s: f64) -> f64 {
711 let codys_threshold = 0.46875;
712 let h = theta_x / s;
713 let t = 0.5 * s;
714 let q1 = -(1.0 / SQRT_TWO) * (h + t);
715 let q2 = -(1.0 / SQRT_TWO) * (h - t);
716 let two_b = if q1 < codys_threshold {
717 if q2 < codys_threshold {
718 (0.5 * theta_x).exp() * erfc_cody(q1) - (-0.5 * theta_x).exp() * erfc_cody(q2)
719 } else {
720 (0.5 * theta_x).exp() * erfc_cody(q1) - (-0.5 * (h * h + t * t)).exp() * erfcx_cody(q2)
721 }
722 } else {
723 if q2 < codys_threshold {
724 (-0.5 * (h * h + t * t)).exp() * erfcx_cody(q1) - (-0.5 * theta_x).exp() * erfc_cody(q2)
725 } else {
726 (-0.5 * (h * h + t * t)).exp() * (erfcx_cody(q1) - erfcx_cody(q2))
727 }
728 };
729 (0.5 * two_b).max(0.0)
730}
731
732fn normalised_vega(x: f64, s: f64) -> f64 {
733 let h = x / s;
734 let t = 0.5 * s;
735 (1.0 / SQRT_TWO_PI) * (-0.5 * (h * h + t * t)).exp()
736}
737
738fn inv_normalised_vega(x: f64, s: f64) -> f64 {
739 let h = x / s;
740 let t = 0.5 * s;
741 SQRT_TWO_PI * (0.5 * (h * h + t * t)).exp()
742}
743
744fn ln_normalised_vega(x: f64, s: f64) -> f64 {
745 let h = x / s;
746 let t = 0.5 * s;
747 -(LN_TWO_PI * 0.5) - 0.5 * (h * h + t * t)
748}
749
750fn is_region_i(theta_x: f64, s: f64) -> bool {
751 theta_x < s * ETA && s * (0.5 * s - (TAU + 0.5 + ETA)) + theta_x < 0.0
752}
753
754fn is_region_ii(theta_x: f64, s: f64) -> bool {
755 s * (s - (2.0 * TAU)) - theta_x / ETA < 0.0
756}
757
758fn normalised_black_internal(theta_x: f64, s: f64) -> f64 {
759 if is_region_i(theta_x, s) {
760 return asymptotic_expansion_of_scaled_normalised_black(theta_x / s, 0.5 * s)
761 * normalised_vega(theta_x, s);
762 }
763 if is_region_ii(theta_x, s) {
764 return small_t_expansion_of_scaled_normalised_black(theta_x / s, 0.5 * s)
765 * normalised_vega(theta_x, s);
766 }
767 normalised_black_with_optimal_use_of_codys_functions(theta_x, s)
768}
769
770fn scaled_normalised_black_and_ln_vega_internal(theta_x: f64, s: f64) -> (f64, f64) {
771 if is_region_i(theta_x, s) {
772 return (
773 asymptotic_expansion_of_scaled_normalised_black(theta_x / s, 0.5 * s),
774 ln_normalised_vega(theta_x, s),
775 );
776 }
777 if is_region_ii(theta_x, s) {
778 return (
779 small_t_expansion_of_scaled_normalised_black(theta_x / s, 0.5 * s),
780 ln_normalised_vega(theta_x, s),
781 );
782 }
783 let ln_vega = ln_normalised_vega(theta_x, s);
784 (
785 normalised_black_with_optimal_use_of_codys_functions(theta_x, s) * (-ln_vega).exp(),
786 ln_vega,
787 )
788}
789
790fn compute_f_lower_map_and_first_two_derivatives(x: f64, s: f64) -> (f64, f64, f64) {
791 let ax = x.abs();
792 let z = SQRT_ONE_OVER_THREE * ax / s;
793 let y = z * z;
794 let s2 = s * s;
795 let phi_z = 0.5 * erfc_cody((1.0 / SQRT_TWO) * z); let pdf_z = norm_pdf(z);
797 let fpp = PI_OVER_SIX * y / (s2 * s) * phi_z
798 * (8.0 * SQRT_THREE * s * ax + (3.0 * s2 * (s2 - 8.0) - 8.0 * x * x) * phi_z / pdf_z)
799 * (2.0 * y + 0.25 * s2).exp();
800 let phi_z2 = phi_z * phi_z;
801 let fp = TWO_PI * y * phi_z2 * (y + 0.125 * s2).exp();
802 let f = TWO_PI_OVER_SQRT_TWENTY_SEVEN * ax * (phi_z2 * phi_z);
803 (f, fp, fpp)
804}
805
806fn inverse_f_lower_map(x: f64, f: f64) -> f64 {
807 (x / (SQRT_THREE
808 * inverse_norm_cdf(
809 SQRT_THREE_OVER_THIRD_ROOT_TWO_PI * f.cbrt() / x.abs().cbrt(),
810 )))
811 .abs()
812}
813
814fn compute_f_upper_map_and_first_two_derivatives(x: f64, s: f64) -> (f64, f64, f64) {
815 let f = 0.5 * erfc_cody((0.5 / SQRT_TWO) * s); let w = square(x / s);
817 let fp = -0.5 * (0.5 * w).exp();
818 let fpp = SQRT_PI_OVER_TWO * (w + 0.125 * s * s).exp() * w / s;
819 (f, fp, fpp)
820}
821
822fn inverse_f_upper_map(f: f64) -> f64 {
823 -2.0 * inverse_norm_cdf(f)
824}
825
826fn one_minus_erfcx(x: f64) -> f64 {
827 if x < -1.0 / 5.0 || x > 1.0 / 3.0 {
828 return 1.0 - erfcx_cody(x);
829 }
830 x * (1.128379167095512573896
831 - x * (1.0000000000000002
832 + x * (1.1514967181784756
833 + x * (5.7689001208873741E-1
834 + x * (1.4069188744609651E-1 + 1.4069285713634565E-2 * x))))
835 / (1.0
836 + x * (1.9037494962421563
837 + x * (1.5089908593742723
838 + x * (6.2486081658640257E-1
839 + x * (1.358008134514386E-1 + 1.2463320728346347E-2 * x))))))
840}
841
842fn implied_normalised_volatility_atm(beta: f64) -> f64 {
843 let beta_max = 0.6826894921370859;
844 if beta <= beta_max {
845 let r = beta_max * beta_max - beta * beta;
846 return beta
847 * ((2.92958954698308816
848 + r * (1.4014698674754995E1
849 + r * (2.44918990556468762E1
850 + r * (1.90763928424894996E1
851 + r * (6.43250149461895996
852 + r * (7.52328633671821543E-1 + 1.38781536163865582E-2 * r))))))
853 / (1.0
854 + r * (5.22443271807813073
855 + r * (1.02258209975070629E1
856 + r * (9.28187483709036392
857 + r * (3.9095549184069553
858 + r * (6.61214199809055912E-1 + 2.89411828874884851E-2 * r)))))));
859 }
860 -2.0 * inverse_norm_cdf_for_low_probabilities(0.5 * (1.0 - beta))
861}
862
863fn b_l_over_b_max(s_c: f64) -> f64 {
864 if s_c >= 2.6267851073127395 {
865 if s_c >= 7.348469228349534 {
866 let s_c_val = s_c;
867 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))))));
868 }
869 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))))));
870 }
871 if s_c >= 0.7099295739719539 {
872 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))))));
873 }
874 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))));
875 (s_c * s_c) * (0.07560996640296361767172 + s_c * (s_c * g - 0.09672719281339436290858))
876}
877
878fn b_u_over_b_max(s_c: f64) -> f64 {
879 if s_c > 1.7888543819998317 {
880 if s_c > 6.164414002968976 {
881 let s_c_val = s_c;
882 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))))));
883 }
884 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))))));
885 }
886 if s_c >= 0.7745966692414833 {
887 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)))));
888 }
889 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))));
890 0.7899085945560627246288 + (s_c * s_c) * (0.0614616805805147403487 + s_c * g)
891}
892
893pub fn lets_be_rational(beta: f64, theta_x: f64, n_iterations: i32) -> f64 {
894 if beta <= 0.0 {
895 return if beta == 0.0 { 0.0 } else { VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC };
896 }
897 let b_max = (0.5 * theta_x).exp();
898 if beta >= b_max {
899 return VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM;
900 }
901 if theta_x == 0.0 {
902 return implied_normalised_volatility_atm(beta);
903 }
904
905 let mut iterations = 0;
906 let mut f;
907 let mut s;
908 let mut ds = -DBL_MAX;
909 let mut s_left = DBL_MIN;
910 let mut s_right = DBL_MAX;
911
912 let sqrt_ax = (-theta_x).sqrt();
913 let s_c = SQRT_TWO * sqrt_ax;
914 let ome = one_minus_erfcx(sqrt_ax);
915 let b_c = 0.5 * b_max * ome;
916
917 if beta < b_c {
918 let s_l = s_c - SQRT_PI_OVER_TWO * ome;
919 let b_l = b_l_over_b_max(s_c) * b_max;
920 if beta < b_l {
921 let (f_lower_map_l, d_f_lower_map_l_d_beta, d2_f_lower_map_l_d_beta2) =
922 compute_f_lower_map_and_first_two_derivatives(theta_x, s_l);
923 let r_ll = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
924 0.0,
925 b_l,
926 0.0,
927 f_lower_map_l,
928 1.0,
929 d_f_lower_map_l_d_beta,
930 d2_f_lower_map_l_d_beta2,
931 true,
932 );
933 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);
934 if !(f > 0.0) {
935 let t = beta / b_l;
936 f = (f_lower_map_l * t + b_l * (1.0 - t)) * t;
937 }
938 s = inverse_f_lower_map(theta_x, f);
939 s_right = s_l;
940
941 let ln_beta = beta.ln();
942 while iterations < n_iterations && ds.abs() > DBL_EPSILON * s {
943 let (bx, ln_vega) = scaled_normalised_black_and_ln_vega_internal(theta_x, s);
944 let ln_b = bx.ln() + ln_vega;
945 let bpob = 1.0 / bx;
946 let b = ln_b.exp();
947 if b > beta && s < s_right {
948 s_right = s;
949 } else if b < beta && s > s_left {
950 s_left = s;
951 }
952
953 let h = theta_x / s;
954 let x2_over_s3 = h * h / s;
955 let b_h2 = x2_over_s3 - s / 4.0;
956 let nu = (ln_beta - ln_b) * ln_b / ln_beta / bpob;
957 let lambda = 1.0 / ln_b;
958 let otl = 1.0 + 2.0 * lambda;
959 let h2 = b_h2 - bpob * otl;
960 let c = 3.0 * (x2_over_s3 / s);
961 let b_h3 = b_h2 * b_h2 - c - 0.25;
962 let mu = 6.0 * lambda * (1.0 + lambda);
963 let h3 = b_h3 + (bpob * bpob) * (2.0 + mu) - (b_h2 * bpob) * 3.0 * otl;
964
965 if theta_x < -190.0 {
966 let b_h4 = b_h2 * (b_h3 - 0.5) - (b_h2 - 2.0 / s) * 2.0 * c;
967 let bpob2 = bpob * bpob;
968 let bppob = b_h2 * bpob;
969 ds = nu * householder_factor_4(
970 nu,
971 h2,
972 h3,
973 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,
974 );
975 } else {
976 ds = nu * householder_factor_3(nu, h2, h3);
977 }
978 s += ds;
979 if s > s_right { s = s_right; } else if s < s_left { s = s_left; }
980 iterations += 1;
981 }
982 return s;
983 } else {
984 let inv_v_c = SQRT_TWO_PI / b_max;
985 let inv_v_l = inv_normalised_vega(theta_x, s_l);
986 let r_lm = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(
987 b_l, b_c, s_l, s_c, inv_v_l, inv_v_c, 0.0, false,
988 );
989 s = rational_cubic_interpolation(beta, b_l, b_c, s_l, s_c, inv_v_l, inv_v_c, r_lm);
990 s_left = s_l;
991 s_right = s_c;
992 }
993 } else {
994 let s_u = s_c + SQRT_PI_OVER_TWO * (2.0 - ome);
995 let b_u = b_u_over_b_max(s_c) * b_max;
996 if beta <= b_u {
997 let inv_v_c = SQRT_TWO_PI / b_max;
998 let inv_v_u = inv_normalised_vega(theta_x, s_u);
999 let r_um = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
1000 b_c, b_u, s_c, s_u, inv_v_c, inv_v_u, 0.0, false,
1001 );
1002 s = rational_cubic_interpolation(beta, b_c, b_u, s_c, s_u, inv_v_c, inv_v_u, r_um);
1003 s_left = s_c;
1004 s_right = s_u;
1005 } else {
1006 let (f_upper_map_h, d_f_upper_map_h_d_beta, d2_f_upper_map_h_d_beta2) =
1007 compute_f_upper_map_and_first_two_derivatives(theta_x, s_u);
1008 let r_uu = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(
1009 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,
1010 );
1011 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);
1012 if f <= 0.0 {
1013 let h = b_max - b_u;
1014 let t = (beta - b_u) / h;
1015 f = (f_upper_map_h * (1.0 - t) + 0.5 * h * t) * (1.0 - t);
1016 }
1017 s = inverse_f_upper_map(f);
1018 s_left = s_u;
1019 if beta > 0.5 * b_max {
1020 let beta_bar = b_max - beta;
1021 while iterations < n_iterations && ds.abs() > DBL_EPSILON * s {
1022 let h = theta_x / s;
1023 let t = s / 2.0;
1024 let gp = (2.0 / SQRT_TWO_PI) / (erfcx_cody((t + h) * (1.0 / SQRT_TWO)) + erfcx_cody((t - h) * (1.0 / SQRT_TWO)));
1025 let b_bar = normalised_vega(theta_x, s) / gp;
1026 if b_bar < beta_bar && s < s_right {
1027 s_right = s;
1028 } else if b_bar > beta_bar && s > s_left {
1029 s_left = s;
1030 }
1031
1032 let g = (beta_bar / b_bar).ln();
1033 let x2_over_s3 = (h * h) / s;
1034 let b_h2 = x2_over_s3 - s / 4.0;
1035 let c = 3.0 * (x2_over_s3 / s);
1036 let b_h3 = b_h2 * b_h2 - c - 0.25;
1037 let nu = -g / gp;
1038 let h2 = b_h2 + gp;
1039 let h3 = b_h3 + gp * (2.0 * gp + 3.0 * b_h2);
1040
1041 if theta_x < -580.0 {
1042 let b_h4 = b_h2 * (b_h3 - 0.5) - (b_h2 - 2.0 / s) * 2.0 * c;
1043 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));
1044 } else {
1045 ds = nu * householder_factor_3(nu, h2, h3);
1046 }
1047 s += ds;
1048 if s > s_right { s = s_right; } else if s < s_left { s = s_left; }
1049 iterations += 1;
1050 }
1051 return s;
1052 }
1053 }
1054 }
1055
1056 while iterations < n_iterations && ds.abs() > DBL_EPSILON * s {
1057 let b = normalised_black_internal(theta_x, s);
1058 let inv_bp = inv_normalised_vega(theta_x, s);
1059 let nu = (beta - b) * inv_bp;
1060 let h = theta_x / s;
1061 let x2_over_s3 = (h * h) / s;
1062 let h2 = x2_over_s3 - s * 0.25;
1063 let h3 = h2 * h2 - 3.0 * (x2_over_s3 / s) - 0.25;
1064 if b > beta && s < s_right {
1065 s_right = s;
1066 } else if b < beta && s > s_left {
1067 s_left = s;
1068 }
1069 ds = nu * householder_factor_3(nu, h2, h3);
1070 s += ds;
1071 if s > s_right { s = s_right; } else if s < s_left { s = s_left; }
1072 iterations += 1;
1073 }
1074 s
1075}
1076
1077pub fn normalised_black(x: f64, s: f64, theta: f64) -> f64 {
1078 if x == 0.0 {
1079 return erf_cody((0.5 / SQRT_TWO) * s);
1080 }
1081 let theta_x = if theta < 0.0 { -x } else { x };
1082 normalised_intrinsic(theta_x) + if s <= 0.0 { 0.0 } else { normalised_black_internal(-x.abs(), s) }
1083}
1084
1085pub fn black(f: f64, k: f64, sigma: f64, t: f64, theta: f64) -> f64 {
1086 let s = sigma * t.sqrt();
1087 if k == f {
1088 return f * erf_cody((0.5 / SQRT_TWO) * s);
1089 }
1090 let mu = if theta < 0.0 {
1091 (k - f).max(0.0)
1092 } else {
1093 (f - k).max(0.0)
1094 };
1095 mu + if s <= 0.0 {
1096 0.0
1097 } else {
1098 (f.sqrt() * k.sqrt()) * normalised_black_internal(-(f / k).ln().abs(), s)
1099 }
1100}
1101
1102pub fn implied_black_volatility(price: f64, f: f64, k: f64, t: f64, theta: f64) -> f64 {
1103 if price >= (if theta < 0.0 { k } else { f }) {
1104 return VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM;
1105 }
1106 let mu = if theta < 0.0 { k - f } else { f - k };
1107 let adjusted_price = if mu > 0.0 { price - mu } else { price };
1108 lets_be_rational(
1109 adjusted_price / (f.sqrt() * k.sqrt()),
1110 -(f / k).ln().abs(),
1111 2,
1112 ) / t.sqrt()
1113}
1114
1115pub fn normalised_implied_black_volatility(beta: f64, x: f64, theta: f64) -> f64 {
1116 let theta_x = if theta < 0.0 { -x } else { x };
1117 lets_be_rational(beta - normalised_intrinsic(theta_x), -x.abs(), 2)
1118}