1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149

//! In addition to working out the error bounds algebraically, we can
//! also have the computer do this work for us as some computation is
//! being performed. This approach is known as *running error
//! analysis*.

// std
use std::ops::{Add, Div, Mul, Sub};
// pbrt
use crate::core::pbrt::MACHINE_EPSILON;
use crate::core::pbrt::{next_float_down, next_float_up};

// see efloat.h

/// Find solution(s) of the quadratic equation at<sup>2</sup> + bt + c = 0 using
/// *EFloat* instead of *Float* for error bounds.
pub fn quadratic_efloat(a: EFloat, b: EFloat, c: EFloat, t0: &mut EFloat, t1: &mut EFloat) -> bool {
    let discrim: f64 = b.v as f64 * b.v as f64 - 4.0f64 * a.v as f64 * c.v as f64;
    if discrim < 0.0 {
        false
    } else {
        let root_discrim: f64 = discrim.sqrt();
        let float_root_discrim: EFloat =
            EFloat::new(root_discrim as f32, MACHINE_EPSILON * root_discrim as f32);
        // compute quadratic _t_ values
        let q = if b.v < 0.0f32 {
            (b - float_root_discrim) * -0.5f32
        } else {
            (b + float_root_discrim) * -0.5f32
        };
        *t0 = q / a;
        *t1 = c / q;
        if t0.v > t1.v {
            std::mem::swap(&mut (*t0), &mut (*t1));
        }
        true
    }
}

/// **EFloat** keeps track of an interval that describes the
/// uncertainty of a value of interest.
#[derive(Debug, Default, Copy, Clone)]
pub struct EFloat {
    pub v: f32,
    pub low: f32,
    pub high: f32,
}

impl EFloat {
    pub fn new(v: f32, err: f32) -> Self {
        if err == 0.0 {
            EFloat { v, low: v, high: v }
        } else {
            EFloat {
                v,
                low: next_float_down(v - err),
                high: next_float_up(v + err),
            }
        }
    }
    pub fn lower_bound(&self) -> f32 {
        self.low
    }
    pub fn upper_bound(&self) -> f32 {
        self.high
    }
}

impl PartialEq for EFloat {
    fn eq(&self, rhs: &EFloat) -> bool {
        self.v == rhs.v
    }
}

impl Add for EFloat {
    type Output = EFloat;
    fn add(self, rhs: EFloat) -> EFloat {
        // TODO: r.Check();
        EFloat {
            v: self.v + rhs.v,
            low: next_float_down(self.lower_bound() + rhs.lower_bound()),
            high: next_float_up(self.upper_bound() + rhs.upper_bound()),
        }
    }
}

impl Sub for EFloat {
    type Output = EFloat;
    fn sub(self, rhs: EFloat) -> EFloat {
        // TODO: r.Check();
        EFloat {
            v: self.v - rhs.v,
            low: next_float_down(self.lower_bound() - rhs.upper_bound()),
            high: next_float_up(self.upper_bound() - rhs.lower_bound()),
        }
    }
}

impl Mul for EFloat {
    type Output = EFloat;
    fn mul(self, rhs: EFloat) -> EFloat {
        let prod: [f32; 4] = [
            self.lower_bound() * rhs.lower_bound(),
            self.upper_bound() * rhs.lower_bound(),
            self.lower_bound() * rhs.upper_bound(),
            self.upper_bound() * rhs.upper_bound(),
        ];
        // TODO: r.Check();
        EFloat {
            v: self.v * rhs.v,
            low: next_float_down(prod[0].min(prod[1]).min(prod[2].min(prod[3]))),
            high: next_float_up(prod[0].max(prod[1]).max(prod[2].max(prod[3]))),
        }
    }
}

impl Mul<f32> for EFloat {
    type Output = EFloat;
    fn mul(self, rhs: f32) -> EFloat {
        EFloat::new(rhs, 0.0) * self
    }
}

impl Div for EFloat {
    type Output = EFloat;
    fn div(self, rhs: EFloat) -> EFloat {
        let div: [f32; 4] = [
            self.lower_bound() / rhs.lower_bound(),
            self.upper_bound() / rhs.lower_bound(),
            self.lower_bound() / rhs.upper_bound(),
            self.upper_bound() / rhs.upper_bound(),
        ];
        // TODO: r.Check();
        if rhs.low < 0.0 && rhs.high > 0.0 {
            // the interval we're dividing by straddles zero, so just
            // return an interval of everything
            EFloat {
                v: self.v / rhs.v,
                low: -std::f32::INFINITY,
                high: std::f32::INFINITY,
            }
        } else {
            EFloat {
                v: self.v / rhs.v,
                low: next_float_down(div[0].min(div[1]).min(div[2].min(div[3]))),
                high: next_float_up(div[0].max(div[1]).max(div[2].max(div[3]))),
            }
        }
    }
}