logicaffeine_base/
expansion.rs

1//! Expansions of the values of several infinitesimals in their smallest one.
2//!
3//! A value of `dx`, declared much smaller than `dt`, that is no rational function of a root of `dx`
4//! over the values of `dt` — `√(dt + dx)`, the state a loop stepped by `dx` leaves — is the sum of its
5//! expansion `Σ c_k·dx^(k/q)`, each coefficient a value of `dt` and of the infinitesimals larger than
6//! it. An [`Expansion`] holds such a value as the operations that built it and computes coefficients on
7//! demand, each once: a rational function of a root of `dx` by power-series division; sums and
8//! products term by term; a reciprocal and a real root by the recurrences [`crate::germ`] reads one
9//! infinitesimal's germs by. Every recurrence runs over any [`Coefficient`] ring, so the same
10//! expansion serves exact coefficients and coefficients that are themselves expansions.
11//!
12//! Orders are integers in units of `s^(1/q)`, `q` the expansion's ramification; two expansions meet on
13//! the least grid both lie on, each refined to it.
14
15use std::fmt;
16use std::sync::{Arc, Mutex};
17
18use crate::algebraic::Algebraic;
19use crate::multi::Multi;
20use crate::numeric::Rational;
21use crate::ordered_field::OrderedField;
22
23/// Why a coefficient or an operation of an expansion is not given.
24#[derive(Clone, Debug, PartialEq, Eq)]
25pub enum ExpansionError {
26    /// The reciprocal of zero.
27    DivisionByZero,
28    /// An even root of a negative value.
29    NegativeRoot(u32),
30    /// A coefficient's own refusal, in its own words.
31    Coefficient(String),
32}
33
34impl fmt::Display for ExpansionError {
35    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
36        match self {
37            ExpansionError::DivisionByZero => f.write_str("Division by zero"),
38            ExpansionError::NegativeRoot(n) => f.write_str(&crate::algebraic::negative_root(*n)),
39            ExpansionError::Coefficient(why) => f.write_str(why),
40        }
41    }
42}
43
44/// The values an expansion's coefficients are: a field with the operations the recurrences ask for,
45/// each refused in the coefficient's own words when it cannot be done.
46pub trait Coefficient: Clone + fmt::Debug + Send + Sync + 'static {
47    fn zero() -> Self;
48    fn one() -> Self;
49    fn of_rational(r: Rational) -> Self;
50    fn add(&self, o: &Self) -> Result<Self, ExpansionError>;
51    fn sub(&self, o: &Self) -> Result<Self, ExpansionError>;
52    fn neg(&self) -> Result<Self, ExpansionError>;
53    fn mul(&self, o: &Self) -> Result<Self, ExpansionError>;
54    /// Whether the coefficient is zero; refused when that cannot be told.
55    fn is_zero(&self) -> Result<bool, ExpansionError>;
56    fn recip(&self) -> Result<Self, ExpansionError>;
57    /// `−1`, `0` or `1`; refused when the sign cannot be told.
58    fn sign(&self) -> Result<i32, ExpansionError>;
59    /// The real `n`-th root, the nonnegative one for an even `n`.
60    fn nth_root(&self, n: u32) -> Result<Self, ExpansionError>;
61    /// `exp`, for an expansion's exponential.
62    fn exp(&self) -> Result<Self, ExpansionError>;
63    /// The natural logarithm of a positive coefficient, for an expansion's logarithm.
64    fn ln(&self) -> Result<Self, ExpansionError>;
65}
66
67/// A source of coefficients defined outside the operations here: a value of one infinitesimal read
68/// as a value of several, a flow's state.
69pub trait Source<C>: fmt::Debug + Send + Sync {
70    /// The coefficient of order `k`, `k` at least the floor the expansion was built with.
71    fn coefficient(&self, k: i64) -> Result<C, ExpansionError>;
72}
73
74/// A value of several infinitesimals is the coefficient of an expansion in a smaller one, exactly.
75impl Coefficient for Multi {
76    fn zero() -> Self {
77        <Multi as OrderedField>::zero()
78    }
79    fn one() -> Self {
80        <Multi as OrderedField>::one()
81    }
82    fn of_rational(r: Rational) -> Self {
83        Multi::Standard(Algebraic::from_rational(r))
84    }
85    fn add(&self, o: &Self) -> Result<Self, ExpansionError> {
86        Ok(OrderedField::add(self, o))
87    }
88    fn sub(&self, o: &Self) -> Result<Self, ExpansionError> {
89        Ok(OrderedField::sub(self, o))
90    }
91    fn neg(&self) -> Result<Self, ExpansionError> {
92        Ok(OrderedField::neg(self))
93    }
94    fn mul(&self, o: &Self) -> Result<Self, ExpansionError> {
95        Ok(OrderedField::mul(self, o))
96    }
97    fn is_zero(&self) -> Result<bool, ExpansionError> {
98        Ok(OrderedField::is_zero(self))
99    }
100    fn recip(&self) -> Result<Self, ExpansionError> {
101        OrderedField::recip(self).ok_or(ExpansionError::DivisionByZero)
102    }
103    fn sign(&self) -> Result<i32, ExpansionError> {
104        Ok(OrderedField::signum(self))
105    }
106    fn nth_root(&self, n: u32) -> Result<Self, ExpansionError> {
107        self.root(n).map_err(ExpansionError::Coefficient)
108    }
109    fn exp(&self) -> Result<Self, ExpansionError> {
110        Err(ExpansionError::Coefficient(
111            "the exponential of a value of several infinitesimals is a germ, not an exact value".to_string(),
112        ))
113    }
114    fn ln(&self) -> Result<Self, ExpansionError> {
115        Err(ExpansionError::Coefficient(
116            "the logarithm of a value of several infinitesimals is a germ, not an exact value".to_string(),
117        ))
118    }
119}
120
121#[derive(Debug)]
122enum Node<C: Coefficient> {
123    /// A value constant in `s`.
124    Constant(C),
125    /// `num/den`, polynomials in `t = s^(1/q)` lowest power first, with `den = t^d·u` and `u(0) ≠ 0`:
126    /// the power series of `num/u` shifted down by `d`.
127    Fraction { num: Vec<C>, unit: Vec<C>, d: i64 },
128    Add(Expansion<C>, Expansion<C>),
129    Sub(Expansion<C>, Expansion<C>),
130    Neg(Expansion<C>),
131    Mul(Expansion<C>, Expansion<C>),
132    /// `1/g` for `g` whose order and leading coefficient are known.
133    Recip(Expansion<C>, i64, C),
134    /// `g` on a grid `m` times finer.
135    Refine(Expansion<C>, u32),
136    /// The real `n`-th root of `g = c·t^v·(1 + …)`, `n` dividing `v`, with `ⁿ√c`.
137    Root { g: Expansion<C>, n: u32, v: i64, c: C, lead: C },
138    /// `exp(g)` for a limited `g`: `k·e_k = Σ_{j=1..k} j·g_j·e_{k−j}` from `e_0 = exp(g_0)`.
139    Exp(Expansion<C>),
140    /// `ln(g)` for `g` with no infinitesimal order and a positive leading coefficient `g_0`.
141    Ln(Expansion<C>, C),
142    Source(Arc<dyn Source<C>>),
143}
144
145/// An expansion in the smallest infinitesimal of a value. See the module documentation.
146#[derive(Clone, Debug)]
147pub struct Expansion<C: Coefficient> {
148    node: Arc<Node<C>>,
149    /// The lowest order that may be nonzero.
150    floor: i64,
151    /// Order `k` is the coefficient of `s^(k/q)`.
152    q: u32,
153    cache: Arc<Mutex<Vec<C>>>,
154}
155
156fn gcd(a: u64, b: u64) -> u64 {
157    if b == 0 {
158        a
159    } else {
160        gcd(b, a % b)
161    }
162}
163
164fn lcm(a: u32, b: u32) -> u32 {
165    let g = gcd(u64::from(a), u64::from(b));
166    u32::try_from(u64::from(a) / g * u64::from(b)).expect("an expansion's ramification fits")
167}
168
169impl<C: Coefficient> Expansion<C> {
170    fn of(node: Node<C>, floor: i64, q: u32) -> Expansion<C> {
171        Expansion {
172            node: Arc::new(node),
173            floor,
174            q,
175            cache: Arc::new(Mutex::new(Vec::new())),
176        }
177    }
178
179    /// The expansion of `num/den`, polynomials in `s^(1/q)` with `den` nonzero.
180    pub fn fraction(num: Vec<C>, den: Vec<C>, q: u32) -> Result<Expansion<C>, ExpansionError> {
181        let mut d = 0usize;
182        while d < den.len() && den[d].is_zero()? {
183            d += 1;
184        }
185        if d == den.len() {
186            return Err(ExpansionError::DivisionByZero);
187        }
188        let unit = den[d..].to_vec();
189        let mut n0 = 0usize;
190        while n0 < num.len() && num[n0].is_zero()? {
191            n0 += 1;
192        }
193        let floor = n0 as i64 - d as i64;
194        Ok(Expansion::of(Node::Fraction { num, unit, d: d as i64 }, floor, q))
195    }
196
197    /// A value constant in `s`: the coefficient itself at order `0`.
198    pub fn constant(c: C) -> Expansion<C> {
199        Expansion::of(Node::Constant(c), 0, 1)
200    }
201
202    /// The expansion a [`Source`] gives, lowest order `floor`, on the grid `q`.
203    pub fn source(source: Arc<dyn Source<C>>, floor: i64, q: u32) -> Expansion<C> {
204        Expansion::of(Node::Source(source), floor, q)
205    }
206
207    /// The ramification `q`: order `k` is the coefficient of `s^(k/q)`.
208    pub fn ramification(&self) -> u32 {
209        self.q
210    }
211
212    /// The lowest order that may be nonzero.
213    pub fn floor(&self) -> i64 {
214        self.floor
215    }
216
217    /// The same expansion on a grid `m` times finer.
218    pub fn refined(&self, m: u32) -> Expansion<C> {
219        if m == 1 {
220            return self.clone();
221        }
222        Expansion::of(Node::Refine(self.clone(), m), self.floor * i64::from(m), self.q * m)
223    }
224
225    /// Both expansions on the least grid both lie on.
226    fn met(&self, o: &Expansion<C>) -> (Expansion<C>, Expansion<C>, u32) {
227        let q = lcm(self.q, o.q);
228        (self.refined(q / self.q), o.refined(q / o.q), q)
229    }
230
231    /// The coefficient of `s^(k/q)`.
232    pub fn coefficient(&self, k: i64) -> Result<C, ExpansionError> {
233        if k < self.floor {
234            return Ok(C::zero());
235        }
236        let at = (k - self.floor) as usize;
237        if let Some(c) = self.cache.lock().expect("unpoisoned").get(at) {
238            return Ok(c.clone());
239        }
240        let c = match self.node.as_ref() {
241            Node::Constant(c) => {
242                if k == 0 {
243                    c.clone()
244                } else {
245                    C::zero()
246                }
247            }
248            Node::Fraction { num, unit, d } => {
249                let j = k + d;
250                if j < 0 {
251                    C::zero()
252                } else {
253                    let mut total = num.get(j as usize).cloned().unwrap_or_else(C::zero);
254                    for (i, u) in unit.iter().enumerate().skip(1) {
255                        let earlier = k - i as i64;
256                        if earlier < self.floor {
257                            break;
258                        }
259                        if !u.is_zero()? {
260                            total = total.sub(&u.mul(&self.coefficient(earlier)?)?)?;
261                        }
262                    }
263                    total.mul(&unit[0].recip()?)?
264                }
265            }
266            Node::Add(a, b) => a.coefficient(k)?.add(&b.coefficient(k)?)?,
267            Node::Sub(a, b) => a.coefficient(k)?.sub(&b.coefficient(k)?)?,
268            Node::Neg(a) => a.coefficient(k)?.neg()?,
269            Node::Mul(a, b) => {
270                let mut total = C::zero();
271                for i in a.floor..=(k - b.floor) {
272                    let x = a.coefficient(i)?;
273                    if x.is_zero()? {
274                        continue;
275                    }
276                    total = total.add(&x.mul(&b.coefficient(k - i)?)?)?;
277                }
278                total
279            }
280            Node::Recip(g, order, lead) => {
281                if k == -order {
282                    lead.recip()?
283                } else {
284                    let mut total = C::zero();
285                    for i in (-order)..k {
286                        let gi = g.coefficient(order + k - i)?;
287                        if gi.is_zero()? {
288                            continue;
289                        }
290                        total = total.add(&self.coefficient(i)?.mul(&gi)?)?;
291                    }
292                    total.neg()?.mul(&lead.recip()?)?
293                }
294            }
295            Node::Refine(g, m) => {
296                let m = i64::from(*m);
297                if k.rem_euclid(m) == 0 {
298                    g.coefficient(k.div_euclid(m))?
299                } else {
300                    C::zero()
301                }
302            }
303            Node::Root { g, n, v, c, lead } => {
304                let j = k - self.floor;
305                if j == 0 {
306                    lead.clone()
307                } else {
308                    let alpha = Rational::from_ratio_i64(1, i64::from(*n)).expect("n ≥ 1").add(&Rational::one());
309                    let mut total = C::zero();
310                    for i in 1..=j {
311                        let gi = g.coefficient(v + i)?;
312                        if gi.is_zero()? {
313                            continue;
314                        }
315                        let weight = alpha.mul(&Rational::from_i64(i)).sub(&Rational::from_i64(j));
316                        total = total.add(&C::of_rational(weight).mul(&gi)?.mul(&self.coefficient(k - i)?)?)?;
317                    }
318                    total.mul(&C::of_rational(Rational::from_i64(j)).mul(c)?.recip()?)?
319                }
320            }
321            Node::Exp(g) => {
322                if k == 0 {
323                    g.coefficient(0)?.exp()?
324                } else {
325                    let mut total = C::zero();
326                    for j in 1..=k {
327                        let gj = g.coefficient(j)?;
328                        if gj.is_zero()? {
329                            continue;
330                        }
331                        total = total.add(&C::of_rational(Rational::from_i64(j)).mul(&gj)?.mul(&self.coefficient(k - j)?)?)?;
332                    }
333                    total.mul(&C::of_rational(Rational::from_ratio_i64(1, k).expect("k ≥ 1")))?
334                }
335            }
336            Node::Ln(g, g0) => {
337                if k == 0 {
338                    g0.ln()?
339                } else {
340                    let mut inner = C::zero();
341                    for j in 1..k {
342                        let gk = g.coefficient(k - j)?;
343                        if gk.is_zero()? {
344                            continue;
345                        }
346                        inner = inner.add(&C::of_rational(Rational::from_i64(j)).mul(&self.coefficient(j)?)?.mul(&gk)?)?;
347                    }
348                    let inner = inner.mul(&C::of_rational(Rational::from_ratio_i64(1, k).expect("k ≥ 1")))?;
349                    g.coefficient(k)?.sub(&inner)?.mul(&g0.recip()?)?
350                }
351            }
352            Node::Source(source) => source.coefficient(k)?,
353        };
354        let mut cache = self.cache.lock().expect("unpoisoned");
355        if cache.len() == at {
356            cache.push(c.clone());
357        }
358        Ok(c)
359    }
360
361    pub fn add(&self, o: &Expansion<C>) -> Expansion<C> {
362        let (a, b, q) = self.met(o);
363        let floor = a.floor.min(b.floor);
364        Expansion::of(Node::Add(a, b), floor, q)
365    }
366
367    pub fn neg(&self) -> Expansion<C> {
368        Expansion::of(Node::Neg(self.clone()), self.floor, self.q)
369    }
370
371    /// The difference coefficient by coefficient, so a term both sides share cancels in the
372    /// coefficient's own arithmetic.
373    pub fn sub(&self, o: &Expansion<C>) -> Expansion<C> {
374        let (a, b, q) = self.met(o);
375        let floor = a.floor.min(b.floor);
376        Expansion::of(Node::Sub(a, b), floor, q)
377    }
378
379    pub fn mul(&self, o: &Expansion<C>) -> Expansion<C> {
380        let (a, b, q) = self.met(o);
381        let floor = a.floor + b.floor;
382        Expansion::of(Node::Mul(a, b), floor, q)
383    }
384
385    /// The first nonzero term at or below the order `through`: its order and coefficient; `None`
386    /// when every coefficient through it is zero.
387    pub fn leading_through(&self, through: i64) -> Result<Option<(i64, C)>, ExpansionError> {
388        for k in self.floor..=through {
389            let c = self.coefficient(k)?;
390            if !c.is_zero()? {
391                return Ok(Some((k, c)));
392            }
393        }
394        Ok(None)
395    }
396
397    /// `1/self` for a value whose leading term is `lead·s^(order/q)`.
398    pub fn recip(&self, order: i64, lead: C) -> Expansion<C> {
399        Expansion::of(Node::Recip(self.clone(), order, lead), -order, self.q)
400    }
401
402    /// `exp(self)` for a value with no negative order.
403    pub fn exp(&self) -> Expansion<C> {
404        Expansion::of(Node::Exp(self.clone()), 0, self.q)
405    }
406
407    /// `ln(self)` for a value with no infinitesimal order whose coefficient of order `0` is `g0`.
408    pub fn ln(&self, g0: C) -> Expansion<C> {
409        Expansion::of(Node::Ln(self.clone(), g0), 0, self.q)
410    }
411
412    /// The real `n`-th root (`n ≥ 2`) of a value whose leading term is `c·s^(v/q)`, the nonnegative
413    /// one for an even `n`; refused for a negative value and an even `n`. The grid is refined until `n`
414    /// divides the leading order.
415    pub fn root(&self, n: u32, v: i64, c: C) -> Result<Expansion<C>, ExpansionError> {
416        assert!(n >= 2, "a root of index 2 or more");
417        if n % 2 == 0 && c.sign()? < 0 {
418            return Err(ExpansionError::NegativeRoot(n));
419        }
420        let m = (u64::from(n) / gcd(u64::from(n), v.unsigned_abs())) as u32;
421        let g = self.refined(m);
422        let v = v * i64::from(m);
423        let lead = c.nth_root(n)?;
424        let floor = v / i64::from(n);
425        let q = g.q;
426        Ok(Expansion::of(Node::Root { g, n, v, c, lead }, floor, q))
427    }
428}
429
430#[cfg(test)]
431#[path = "../../logicaffeine_tests/unit/logicaffeine_base/expansion/tests.rs"]
432mod tests;
logicaffeine_base/expansion.rs

logicaffeine_base/
expansion.rs
