pub trait FormalPowerSeriesCoefficient:
Sized
+ Clone
+ Zero
+ PartialEq
+ One
+ From<usize>
+ From<isize>
+ Add<Output = Self>
+ Sub<Output = Self>
+ Mul<Output = Self>
+ Div<Output = Self>
+ for<'r> Add<&'r Self, Output = Self>
+ for<'r> Sub<&'r Self, Output = Self>
+ for<'r> Mul<&'r Self, Output = Self>
+ for<'r> Div<&'r Self, Output = Self>
+ AddAssign<Self>
+ SubAssign<Self>
+ MulAssign<Self>
+ DivAssign<Self>
+ for<'r> AddAssign<&'r Self>
+ for<'r> SubAssign<&'r Self>
+ for<'r> MulAssign<&'r Self>
+ for<'r> DivAssign<&'r Self>
+ Neg<Output = Self> {
type Base: MIntConvert<usize> + MIntConvert<isize>;
// Required methods
fn pow(self, exp: usize) -> Self;
fn memorized_fact(mf: &MemorizedFactorial<Self::Base>) -> &[Self];
fn memorized_inv_fact(mf: &MemorizedFactorial<Self::Base>) -> &[Self];
fn memorized_inv(mf: &MemorizedFactorial<Self::Base>, n: usize) -> Self;
// Provided methods
fn signed_pow(self, exp: isize) -> Self { ... }
fn memorized_factorial(n: usize) -> MemorizedFactorial<Self::Base> { ... }
}Required Associated Types§
type Base: MIntConvert<usize> + MIntConvert<isize>
Required Methods§
fn pow(self, exp: usize) -> Self
fn memorized_fact(mf: &MemorizedFactorial<Self::Base>) -> &[Self]
fn memorized_inv_fact(mf: &MemorizedFactorial<Self::Base>) -> &[Self]
fn memorized_inv(mf: &MemorizedFactorial<Self::Base>, n: usize) -> Self
Provided Methods§
Sourcefn signed_pow(self, exp: isize) -> Self
fn signed_pow(self, exp: isize) -> Self
Examples found in repository?
crates/competitive/src/math/formal_power_series/formal_power_series_impls.rs (line 516)
462 pub fn mul_of_pow_sparse(&self, q: &Self, exp_p: isize, exp_q: isize, deg: usize) -> Self {
463 if deg == 0 {
464 return Self::zero();
465 }
466 if exp_p == 0 && exp_q == 0 {
467 return Self::from_vec(
468 once(T::one())
469 .chain(repeat_with(T::zero))
470 .take(deg)
471 .collect(),
472 );
473 }
474 if exp_p != 0 && self.iter().all(|x| x.is_zero()) {
475 assert!(exp_p > 0);
476 return Self::zeros(deg);
477 }
478 if exp_q != 0 && q.iter().all(|x| x.is_zero()) {
479 assert!(exp_q > 0);
480 return Self::zeros(deg);
481 }
482
483 let normalize = |f: &Self, exp: isize| {
484 if exp == 0 {
485 return (0usize, T::one(), Self::from_vec(vec![T::one()]));
486 }
487 let k = f.iter().position(|value| !value.is_zero()).unwrap();
488 assert!(
489 exp >= 0 || k == 0,
490 "Negative exponent with zero constant term"
491 );
492 let c = f[k].clone();
493 let f = (f.clone() >> k) / &c;
494 (k, c, f)
495 };
496 let (sp, cp, mut p) = normalize(self, exp_p);
497 let (sq, cq, mut q) = normalize(q, exp_q);
498
499 let shift = exp_p
500 .saturating_mul(sp as _)
501 .saturating_add(exp_q.saturating_mul(sq as _)) as usize;
502 if shift >= deg {
503 return Self::zeros(deg);
504 }
505 p.truncate(deg - shift);
506 q.truncate(deg - shift);
507
508 let mut f = Self::solve_sparse_differential2(
509 &p,
510 &q,
511 &p,
512 T::from(exp_p),
513 T::from(exp_q),
514 deg - shift,
515 );
516 f *= cp.signed_pow(exp_p) * cq.signed_pow(exp_q);
517 if shift > 0 {
518 f <<= shift;
519 }
520 f.prefix(deg)
521 }Sourcefn memorized_factorial(n: usize) -> MemorizedFactorial<Self::Base>
fn memorized_factorial(n: usize) -> MemorizedFactorial<Self::Base>
Examples found in repository?
crates/competitive/src/math/formal_power_series/formal_power_series_impls.rs (line 305)
293 pub fn exp(&self, deg: usize) -> Self {
294 debug_assert!(self[0].is_zero());
295 if self.data.iter().filter(|x| !x.is_zero()).count()
296 <= deg.next_power_of_two().trailing_zeros() as usize * 16
297 {
298 let diff = self.clone().diff();
299 let pos: Vec<_> = diff
300 .data
301 .iter()
302 .enumerate()
303 .filter_map(|(i, x)| if x.is_zero() { None } else { Some(i) })
304 .collect();
305 let mf = T::memorized_factorial(deg);
306 let mut f = Self::zeros(deg);
307 f[0] = T::one();
308 for i in 1..deg {
309 let mut tot = T::zero();
310 for &j in &pos {
311 if j > i - 1 {
312 break;
313 }
314 tot += f[i - 1 - j].clone() * &diff[j];
315 }
316 f[i] = tot * T::memorized_inv(&mf, i);
317 }
318 return f;
319 }
320 let mut f = Self::one();
321 let mut i = 1;
322 while i < deg {
323 let mut g = -f.log(i * 2);
324 g[0] += T::one();
325 for (g, x) in g.iter_mut().zip(self.iter().take(i * 2)) {
326 *g += x.clone();
327 }
328 f = (f * g).prefix(i * 2);
329 i *= 2;
330 }
331 f.prefix(deg)
332 }
333 pub fn log(&self, deg: usize) -> Self {
334 (self.inv(deg) * self.clone().diff()).integral().prefix(deg)
335 }
336 pub fn pow(&self, rhs: usize, deg: usize) -> Self {
337 if rhs == 0 {
338 return Self::from_vec(
339 once(T::one())
340 .chain(repeat_with(T::zero))
341 .take(deg)
342 .collect(),
343 );
344 }
345 if let Some(k) = self.iter().position(|x| !x.is_zero()) {
346 if k >= deg.div_ceil(rhs) {
347 Self::zeros(deg)
348 } else {
349 let deg = deg - k * rhs;
350 let x0 = self[k].clone();
351 let mut f = (self >> k) / &x0;
352 if f.data.iter().filter(|x| !x.is_zero()).count()
353 <= deg.next_power_of_two().trailing_zeros() as usize * 12
354 {
355 f = f.pow_sparse1(T::from(rhs), deg);
356 } else {
357 f = (f.log(deg) * &T::from(rhs)).exp(deg);
358 }
359 f *= x0.pow(rhs);
360 f <<= k * rhs;
361 f
362 }
363 } else {
364 Self::zeros(deg)
365 }
366 }
367 fn pow_sparse1(&self, rhs: T, deg: usize) -> Self {
368 debug_assert!(!self[0].is_zero());
369 let pos: Vec<_> = self
370 .data
371 .iter()
372 .enumerate()
373 .skip(1)
374 .filter_map(|(i, x)| if x.is_zero() { None } else { Some(i) })
375 .collect();
376 let mf = T::memorized_factorial(deg);
377 let mut f = Self::zeros(deg);
378 f[0] = T::one();
379 for i in 1..deg {
380 let mut tot = T::zero();
381 for &j in &pos {
382 if j > i {
383 break;
384 }
385 tot += (T::from(j) * &rhs - T::from(i - j)) * &self[j] * &f[i - j];
386 }
387 f[i] = tot * T::memorized_inv(&mf, i);
388 }
389 f
390 }
391
392 fn sparse_fold(&self, sparse: impl IntoIterator<Item = (usize, T)>, deg: usize) -> T {
393 sparse
394 .into_iter()
395 .take_while(|&(i, _)| i <= deg)
396 .fold(T::zero(), |sum, (i, x)| sum + x * self.coeff(deg - i))
397 }
398
399 /// solve: $X(QF)'=\alpha P'(QF)+\beta P(Q'F)$ in $O(deg * max(nz(P), nz(Q), nz(X)))$
400 pub fn solve_sparse_differential2(
401 p: &Self,
402 q: &Self,
403 x: &Self,
404 alpha: T,
405 beta: T,
406 deg: usize,
407 ) -> Self {
408 if deg == 0 {
409 return Self::zero();
410 }
411 let collect_sparse = |p: &Self| -> Vec<(usize, T)> {
412 p.iter()
413 .enumerate()
414 .filter(|&(_, x)| !x.is_zero())
415 .map(|(i, x)| (i, x.clone()))
416 .collect()
417 };
418 assert!(q.coeff(0).is_one());
419 assert!(x.coeff(0).is_one());
420 let p = collect_sparse(p);
421 let q = collect_sparse(q);
422 let x = collect_sparse(x);
423 let diff = |p: &[(usize, T)]| -> Vec<(usize, T)> {
424 p.iter()
425 .filter(|&&(i, _)| i > 0)
426 .map(|&(i, ref x)| (i - 1, x.clone() * T::from(i)))
427 .collect()
428 };
429 let dp = diff(&p);
430 let dq = diff(&q);
431
432 let mf = T::memorized_factorial(deg);
433 let mut f = Self::zeros(deg);
434 let mut qf = Self::zeros(deg);
435 let mut dq_f = Self::zeros(deg);
436 let mut d_qf = Self::zeros(deg);
437 f[0] = T::one();
438 for i in 0..deg - 1 {
439 qf[i] = f.sparse_fold(q.iter().cloned(), i);
440 dq_f[i] = f.sparse_fold(dq.iter().cloned(), i);
441 let dp_qf_i = qf.sparse_fold(dp.iter().cloned(), i);
442 let p_dq_f_i = dq_f.sparse_fold(p.iter().cloned(), i);
443 let x_d_qf_i = d_qf.sparse_fold(
444 x.iter()
445 .map(|&(i, ref x)| (i, x.clone() - T::from((i == 0) as usize))),
446 i,
447 );
448 d_qf[i] = alpha.clone() * dp_qf_i + beta.clone() * p_dq_f_i - x_d_qf_i;
449
450 let mut f_ip1 = d_qf[i].clone();
451 for &(j, ref q) in q.iter().take_while(|&&(j, _)| j <= i) {
452 if j > 0 {
453 f_ip1 -= q.clone() * &f[i - (j - 1)] * T::from(i - (j - 1));
454 }
455 }
456 f[i + 1] = f_ip1 * T::memorized_inv(&mf, i + 1);
457 }
458 f
459 }
460
461 /// P^exp_p * Q^exp_q
462 pub fn mul_of_pow_sparse(&self, q: &Self, exp_p: isize, exp_q: isize, deg: usize) -> Self {
463 if deg == 0 {
464 return Self::zero();
465 }
466 if exp_p == 0 && exp_q == 0 {
467 return Self::from_vec(
468 once(T::one())
469 .chain(repeat_with(T::zero))
470 .take(deg)
471 .collect(),
472 );
473 }
474 if exp_p != 0 && self.iter().all(|x| x.is_zero()) {
475 assert!(exp_p > 0);
476 return Self::zeros(deg);
477 }
478 if exp_q != 0 && q.iter().all(|x| x.is_zero()) {
479 assert!(exp_q > 0);
480 return Self::zeros(deg);
481 }
482
483 let normalize = |f: &Self, exp: isize| {
484 if exp == 0 {
485 return (0usize, T::one(), Self::from_vec(vec![T::one()]));
486 }
487 let k = f.iter().position(|value| !value.is_zero()).unwrap();
488 assert!(
489 exp >= 0 || k == 0,
490 "Negative exponent with zero constant term"
491 );
492 let c = f[k].clone();
493 let f = (f.clone() >> k) / &c;
494 (k, c, f)
495 };
496 let (sp, cp, mut p) = normalize(self, exp_p);
497 let (sq, cq, mut q) = normalize(q, exp_q);
498
499 let shift = exp_p
500 .saturating_mul(sp as _)
501 .saturating_add(exp_q.saturating_mul(sq as _)) as usize;
502 if shift >= deg {
503 return Self::zeros(deg);
504 }
505 p.truncate(deg - shift);
506 q.truncate(deg - shift);
507
508 let mut f = Self::solve_sparse_differential2(
509 &p,
510 &q,
511 &p,
512 T::from(exp_p),
513 T::from(exp_q),
514 deg - shift,
515 );
516 f *= cp.signed_pow(exp_p) * cq.signed_pow(exp_q);
517 if shift > 0 {
518 f <<= shift;
519 }
520 f.prefix(deg)
521 }
522
523 /// exp(P/Q)
524 pub fn exp_of_div_sparse(&self, q: &Self, deg: usize) -> Self {
525 if deg == 0 {
526 return Self::zero();
527 }
528 let shift_q = q
529 .iter()
530 .position(|value| !value.is_zero())
531 .expect("Zero denominator");
532 let shift_p = self.iter().position(|value| !value.is_zero()).unwrap_or(!0);
533 assert!(shift_p > shift_q);
534
535 let mut p = self >> shift_q;
536 let mut q = q >> shift_q;
537 assert!(!q.coeff(0).is_zero());
538
539 let c = q[0].clone();
540 p /= c.clone();
541 q /= c;
542
543 Self::solve_sparse_differential2(&p, &q, &q, T::one(), -T::one(), deg)
544 }
545}
546
547impl<T, C> FormalPowerSeries<T, C>
548where
549 T: FormalPowerSeriesCoefficientSqrt,
550 C: ConvolveSteps<T = Vec<T>>,
551{
552 pub fn sqrt(&self, deg: usize) -> Option<Self> {
553 if self[0].is_zero() {
554 if let Some(k) = self.iter().position(|x| !x.is_zero()) {
555 if k % 2 != 0 {
556 return None;
557 } else if deg > k / 2 {
558 return Some((self >> k).sqrt(deg - k / 2)? << (k / 2));
559 }
560 }
561 } else {
562 let s = self[0].sqrt_coefficient()?;
563 if self.data.iter().filter(|x| !x.is_zero()).count()
564 <= deg.next_power_of_two().trailing_zeros() as usize * 4
565 {
566 let t = self[0].clone();
567 let mut f = self / t;
568 f = f.pow_sparse1(T::one() / T::from(2usize), deg);
569 f *= s;
570 return Some(f);
571 }
572
573 let mut f = Self::from(s);
574 let inv2 = T::one() / (T::one() + T::one());
575 let mut i = 1;
576 while i < deg {
577 f = (&f + &(self.prefix_ref(i * 2) * f.inv(i * 2))).prefix(i * 2) * &inv2;
578 i *= 2;
579 }
580 f.truncate(deg);
581 return Some(f);
582 }
583 Some(Self::zeros(deg))
584 }
585}
586
587impl<T, C> FormalPowerSeries<T, C>
588where
589 T: FormalPowerSeriesCoefficient,
590 C: ConvolveSteps<T = Vec<T>>,
591{
592 pub fn count_subset_sum<F>(&self, deg: usize, mut inverse: F) -> Self
593 where
594 F: FnMut(usize) -> T,
595 {
596 let n = self.length();
597 let mut f = Self::zeros(n);
598 for i in 1..n {
599 if !self[i].is_zero() {
600 for (j, d) in (0..n).step_by(i).enumerate().skip(1) {
601 if j & 1 != 0 {
602 f[d] += self[i].clone() * &inverse(j);
603 } else {
604 f[d] -= self[i].clone() * &inverse(j);
605 }
606 }
607 }
608 }
609 f.exp(deg)
610 }
611 pub fn count_multiset_sum<F>(&self, deg: usize, mut inverse: F) -> Self
612 where
613 F: FnMut(usize) -> T,
614 {
615 let n = self.length();
616 let mut f = Self::zeros(n);
617 for i in 1..n {
618 if !self[i].is_zero() {
619 for (j, d) in (0..n).step_by(i).enumerate().skip(1) {
620 f[d] += self[i].clone() * &inverse(j);
621 }
622 }
623 }
624 f.exp(deg)
625 }
626 /// [x^n] P(x) / Q(x)
627 pub fn bostan_mori(mut self, mut rhs: Self, mut n: usize) -> T
628 where
629 C: NttReuse<T = Vec<T>>,
630 {
631 let mut res = T::zero();
632 rhs.trim_tail_zeros();
633 if self.length() >= rhs.length() {
634 let r = &self / &rhs;
635 if n < r.length() {
636 res = r[n].clone();
637 }
638 self -= r * &rhs;
639 self.trim_tail_zeros();
640 }
641 let k = rhs.length().next_power_of_two();
642 let mut p = C::transform(self.data, k * 2);
643 let mut q = C::transform(rhs.data, k * 2);
644 while n > 0 {
645 let t = C::even_mul_normal_neg(&q, &q);
646 p = if n.is_multiple_of(2) {
647 C::even_mul_normal_neg(&p, &q)
648 } else {
649 C::odd_mul_normal_neg(&p, &q)
650 };
651 q = t;
652 n /= 2;
653 if n != 0 {
654 if C::MULTIPLE {
655 p = C::transform(C::inverse_transform(p, k), k * 2);
656 q = C::transform(C::inverse_transform(q, k), k * 2);
657 } else {
658 p = C::ntt_doubling(p);
659 q = C::ntt_doubling(q);
660 }
661 }
662 }
663 let p = C::inverse_transform(p, k);
664 let q = C::inverse_transform(q, k);
665 res + p[0].clone() / q[0].clone()
666 }
667 /// return F(x) where [x^n] P(x) / Q(x) = [x^d-1] P(x) F(x)
668 pub fn bostan_mori_msb(self, n: usize) -> Self {
669 let d = self.length() - 1;
670 if n == 0 {
671 return (Self::one() << (d - 1)) / self[0].clone();
672 }
673 let q = self;
674 let mq = q.clone().parity_inversion();
675 let w = (q * &mq).even().bostan_mori_msb(n / 2);
676 let mut s = Self::zeros(w.length() * 2 - (n % 2));
677 for (i, x) in w.iter().enumerate() {
678 s[i * 2 + (1 - n % 2)] = x.clone();
679 }
680 let len = 2 * d + 1;
681 let ts = C::transform(s.prefix(len).data, len);
682 mq.reversed().middle_product(&ts, len).prefix(d + 1)
683 }
684 /// x^n mod self
685 pub fn pow_mod(self, n: usize) -> Self {
686 let d = self.length() - 1;
687 let q = self.reversed();
688 let u = q.clone().bostan_mori_msb(n);
689 let mut f = (u * q).prefix(d).reversed();
690 f.trim_tail_zeros();
691 f
692 }
693 fn middle_product(self, other: &C::F, deg: usize) -> Self {
694 let n = self.length();
695 let mut s = C::transform(self.reversed().data, deg);
696 C::multiply(&mut s, other);
697 Self::from_vec((C::inverse_transform(s, deg))[n - 1..].to_vec())
698 }
699 pub fn multipoint_evaluation(self, points: &[T]) -> Vec<T> {
700 let n = points.len();
701 if n <= 32 {
702 return points.iter().map(|p| self.eval(p.clone())).collect();
703 }
704 let mut subproduct_tree = Vec::with_capacity(n * 2);
705 subproduct_tree.resize_with(n, Zero::zero);
706 for x in points {
707 subproduct_tree.push(Self::from_vec(vec![-x.clone(), T::one()]));
708 }
709 for i in (1..n).rev() {
710 subproduct_tree[i] = &subproduct_tree[i * 2] * &subproduct_tree[i * 2 + 1];
711 }
712 let mut uptree_t = Vec::with_capacity(n * 2);
713 uptree_t.resize_with(1, Zero::zero);
714 subproduct_tree.reverse();
715 subproduct_tree.pop();
716 let m = self.length();
717 let v = subproduct_tree.pop().unwrap().reversed().resized(m);
718 let s = C::transform(self.data, m * 2);
719 uptree_t.push(v.inv(m).middle_product(&s, m * 2).resized(n).reversed());
720 for i in 1..n {
721 let subl = subproduct_tree.pop().unwrap();
722 let subr = subproduct_tree.pop().unwrap();
723 let (dl, dr) = (subl.length(), subr.length());
724 let len = dl.max(dr) + uptree_t[i].length();
725 let s = C::transform(uptree_t[i].data.to_vec(), len);
726 uptree_t.push(subr.middle_product(&s, len).prefix(dl));
727 uptree_t.push(subl.middle_product(&s, len).prefix(dr));
728 }
729 uptree_t[n..]
730 .iter()
731 .map(|u| u.data.first().cloned().unwrap_or_else(Zero::zero))
732 .collect()
733 }
734 pub fn product_all<I>(iter: I, deg: usize) -> Self
735 where
736 I: IntoIterator<Item = Self>,
737 {
738 let mut heap: BinaryHeap<_> = iter
739 .into_iter()
740 .map(|f| PartialIgnoredOrd(Reverse(f.length()), f))
741 .collect();
742 while let Some(PartialIgnoredOrd(_, x)) = heap.pop() {
743 if let Some(PartialIgnoredOrd(_, y)) = heap.pop() {
744 let z = (x * y).prefix(deg);
745 heap.push(PartialIgnoredOrd(Reverse(z.length()), z));
746 } else {
747 return x;
748 }
749 }
750 Self::one()
751 }
752 pub fn sum_all_rational<I>(iter: I, deg: usize) -> (Self, Self)
753 where
754 I: IntoIterator<Item = (Self, Self)>,
755 {
756 let mut heap: BinaryHeap<_> = iter
757 .into_iter()
758 .map(|(f, g)| PartialIgnoredOrd(Reverse(f.length().max(g.length())), (f, g)))
759 .collect();
760 while let Some(PartialIgnoredOrd(_, (xa, xb))) = heap.pop() {
761 if let Some(PartialIgnoredOrd(_, (ya, yb))) = heap.pop() {
762 let zb = (&xb * &yb).prefix(deg);
763 let za = (xa * yb + ya * xb).prefix(deg);
764 heap.push(PartialIgnoredOrd(
765 Reverse(za.length().max(zb.length())),
766 (za, zb),
767 ));
768 } else {
769 return (xa, xb);
770 }
771 }
772 (Self::zero(), Self::one())
773 }
774 pub fn kth_term_of_linearly_recurrence(self, a: Vec<T>, k: usize) -> T
775 where
776 C: NttReuse<T = Vec<T>>,
777 {
778 if let Some(x) = a.get(k) {
779 return x.clone();
780 }
781 let p = (Self::from_vec(a).prefix(self.length() - 1) * &self).prefix(self.length() - 1);
782 p.bostan_mori(self, k)
783 }
784 pub fn kth_term(a: Vec<T>, k: usize) -> T
785 where
786 C: NttReuse<T = Vec<T>>,
787 {
788 if let Some(x) = a.get(k) {
789 return x.clone();
790 }
791 Self::from_vec(berlekamp_massey(&a)).kth_term_of_linearly_recurrence(a, k)
792 }
793 /// sum_i a_i exp(b_i x)
794 pub fn linear_sum_of_exp<I, F>(iter: I, deg: usize, mut inv_fact: F) -> Self
795 where
796 I: IntoIterator<Item = (T, T)>,
797 F: FnMut(usize) -> T,
798 {
799 let (p, q) = Self::sum_all_rational(
800 iter.into_iter()
801 .map(|(a, b)| (Self::from_vec(vec![a]), Self::from_vec(vec![T::one(), -b]))),
802 deg,
803 );
804 let mut f = (p * q.inv(deg)).prefix(deg);
805 for i in 0..f.length() {
806 f[i] *= inv_fact(i);
807 }
808 f
809 }
810 /// sum_i (a_i x)^j
811 pub fn sum_of_powers<I>(iter: I, deg: usize) -> Self
812 where
813 I: IntoIterator<Item = T>,
814 {
815 let mut n = T::zero();
816 let prod = Self::product_all(
817 iter.into_iter().map(|a| {
818 n += T::one();
819 Self::from_vec(vec![T::one(), -a])
820 }),
821 deg,
822 );
823 (-prod.log(deg).diff() << 1) + Self::from_vec(vec![n])
824 }
825
826 pub fn power_projection(&self, w: &[T], m: usize) -> Self
827 where
828 C: NttReuse<T = Vec<T>>,
829 {
830 if w.is_empty() {
831 return Self::zeros(m);
832 }
833 if m <= 1 {
834 return Self::from_vec(vec![w[0].clone(); m]);
835 }
836
837 let n0 = w.len();
838 let mut n = n0.next_power_of_two();
839 let mut f = self.prefix_ref(n);
840 f.resize(n);
841
842 let base = n * 2;
843 let mut p_flat = vec![T::zero(); base];
844 for (i, wi) in w.iter().enumerate() {
845 p_flat[n - 1 - i] = wi.clone();
846 }
847 let mut q_flat = vec![T::zero(); base * 2];
848 q_flat[0] = T::one();
849 let q_offset = base;
850 for (i, fi) in f.iter().enumerate() {
851 q_flat[q_offset + i] = -fi.clone();
852 }
853 let mut py = 1usize;
854 let mut qy = 2usize;
855
856 let y_limit = m;
857 while n > 1 {
858 let base = n * 2;
859 let len_p = base * py;
860 let len_q = base * qy;
861 let len = (len_p + len_q - 1).max(len_q + len_q - 1);
862 let len_pot = len.max(1).next_power_of_two();
863 let half = len_pot / 2;
864
865 let p_fft = C::transform(p_flat, len);
866 let q_fft = C::transform(q_flat, len);
867 let pr_fft = C::odd_mul_normal_neg(&p_fft, &q_fft);
868 let qr_fft = C::even_mul_normal_neg(&q_fft, &q_fft);
869 let mut pr = C::inverse_transform(pr_fft, half);
870 let mut qr = C::inverse_transform(qr_fft, half);
871
872 let new_py = (py + qy - 1).min(y_limit);
873 let new_qy = (qy + qy - 1).min(y_limit);
874 let new_base = n;
875 let need_p = new_base * new_py;
876 if pr.len() < need_p {
877 pr.resize_with(need_p, T::zero);
878 } else if pr.len() > need_p {
879 pr.truncate(need_p);
880 }
881 let need_q = new_base * new_qy;
882 if qr.len() < need_q {
883 qr.resize_with(need_q, T::zero);
884 } else if qr.len() > need_q {
885 qr.truncate(need_q);
886 }
887
888 let n2 = n / 2;
889 if n2 < new_base {
890 for y in 0..new_py {
891 let start = y * new_base + n2;
892 let end = y * new_base + new_base;
893 for v in pr[start..end].iter_mut() {
894 *v = T::zero();
895 }
896 }
897 for y in 0..new_qy {
898 let start = y * new_base + n2;
899 let end = y * new_base + new_base;
900 for v in qr[start..end].iter_mut() {
901 *v = T::zero();
902 }
903 }
904 }
905
906 p_flat = pr;
907 q_flat = qr;
908 py = new_py;
909 qy = new_qy;
910 n = n2;
911 }
912
913 let base = 2;
914 let mut p_y = Vec::with_capacity(py);
915 for y in 0..py {
916 p_y.push(p_flat[base * y].clone());
917 }
918 let mut q_y = Vec::with_capacity(qy);
919 for y in 0..qy {
920 q_y.push(q_flat[base * y].clone());
921 }
922 (Self::from_vec(p_y) * Self::from_vec(q_y).inv(m)).prefix(m)
923 }
924
925 pub fn compositional_inverse(&self, deg: usize) -> Self
926 where
927 C: NttReuse<T = Vec<T>>,
928 {
929 if deg == 0 {
930 return Self::zero();
931 }
932 if deg == 1 {
933 return Self::from_vec(vec![T::zero()]);
934 }
935 debug_assert!(self[0].is_zero());
936 debug_assert!(!self[1].is_zero());
937
938 let mut f = self.prefix_ref(deg);
939 f.resize(deg);
940 let c = f[1].clone();
941 f /= c.clone();
942
943 let mut w = vec![T::zero(); deg];
944 w[deg - 1] = T::one();
945 let s = f.power_projection(&w, deg);
946
947 let n = deg - 1;
948 let n_t = T::from(n);
949 let mut h = vec![T::zero(); n];
950 for i in 1..=n {
951 h[n - i] = s[i].clone() * &n_t / T::from(i);
952 }
953
954 let h_fps = Self::from_vec(h);
955 let inv_n = T::one() / n_t;
956 let mut t = h_fps.log(n);
957 t *= -inv_n;
958 let g_over_x = t.exp(n);
959 let mut g = (g_over_x << 1).prefix(deg);
960
961 let inv_c = T::one() / c;
962 let mut pow = T::one();
963 for coef in g.iter_mut() {
964 *coef *= pow.clone();
965 pow *= inv_c.clone();
966 }
967 g
968 }
969 /// f(x) <- f(x + a)
970 pub fn taylor_shift(mut self, a: T) -> Self {
971 let f = T::memorized_factorial(self.length());
972 let n = self.length();
973 for i in 0..n {
974 self.data[i] *= T::memorized_fact(&f)[i].clone();
975 }
976 self.data.reverse();
977 let mut b = a.clone();
978 let mut g = Self::from_vec(T::memorized_inv_fact(&f)[..n].to_vec());
979 for i in 1..n {
980 g[i] *= b.clone();
981 b *= a.clone();
982 }
983 self *= g;
984 self.truncate(n);
985 self.data.reverse();
986 for i in 0..n {
987 self.data[i] *= T::memorized_inv_fact(&f)[i].clone();
988 }
989 self
990 }Dyn Compatibility§
This trait is not dyn compatible.
In older versions of Rust, dyn compatibility was called "object safety", so this trait is not object safe.