moebius_transform_fast.sf 887 B

123456789101112131415161718192021222324252627282930313233343536373839404142
  1. #!/usr/bin/ruby
  2. # Decently efficient implementations of the Möbius inversion formula.
  3. # Formula definition:
  4. # g(n) = Sum_{d|n} mu(d) * f(n/d)
  5. # f(n) = Sum_{d|n} g(d)
  6. # See also:
  7. # https://en.wikipedia.org/wiki/M%C3%B6bius_inversion_formula
  8. func moebius_transform(n, f) {
  9. var a = n.factor_prod {|p,e| ipow(p, e-1) }
  10. var b = idiv(n, a)
  11. b.divisors.sum {|d|
  12. moebius(idiv(b, d)) * f(a * d)
  13. }
  14. }
  15. func moebius_transform_alt(n, f) {
  16. n.squarefree_divisors.sum {|d|
  17. moebius(d) * f(idiv(n, d))
  18. }
  19. }
  20. var f = { .phi }
  21. say 30.of {|n| moebius_transform(n, f) } # OEIS: A007431
  22. say 30.of {|n| moebius_transform_alt(n, f) } # OEIS: A007431
  23. assert_eq(
  24. 30.of {|n| f(n) },
  25. 30.of {|n| n.divisors.sum{|d| moebius_transform(d,f) } }
  26. )
  27. assert_eq(
  28. 30.of {|n| f(n) },
  29. 30.of {|n| n.divisors.sum{|d| moebius_transform_alt(d,f) } }
  30. )