123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140 |
- #!/usr/bin/perl
- func f(n) {
- #moebius(rad(n))
- (-1)**omega(n)
- }
- func g(n) {
- n.divisors.sum{|d|
- f(d)
- }
- }
- func g2(n) {
- n.divisors.sum {|d|
- moebius(rad(d))
- }
- }
- func g3(n) {
- #if (n.is_perfect_power) {
- #return (n.factor.lcm / n.factor_exp{.tail}.lcm)
- #return n.factor_exp.map{.tail}.max
- #assert_eq(n.perfect_root, n.perfect_root.rad)
- #}
- #n.is_perfect_power ? (moebius(n.rad) * ((n.perfect_power - 1) ** omega(n.rad))) : 0
- #moebius(n.perfect_root) *
- #(n.perfect_power-1)**omega(n) *
- (n.perfect_power - 1)**omega(n) * moebius(n.rad)
- }
- func gsum(n) {
- sum(1..n, {|k|
- (-1)**bigomega(k) * faulhaber(floor(n/k), 1)
- #k.divisors.sum {|d|
- # d**2 * liouville(d)
- #}
- })
- }
- func gsum2(n) {
- sum(1..n, {|k|
- moebius(k)**2 * faulhaber(floor(n/k), 1)
- #k.divisors.sum {|d|
- # d**2 * liouville(d)
- #}
- })
- }
- say g(2**3 * 3**3)
- say g3(2**3 * 3**3)
- for n in (2..100) {
- say "Testing: #{n}"
- #assert_eq(g3(n), g(n), "Failed for: n = #{n}")
- # next if n.is_squarefree
- #~ if (n.is_squarefree) {
- #~ assert_eq(g(n), 0)
- #~ }
- #~ next
- #n.is_prime_power || next
- #n.is_perfect_power || next
- #n.perfect_root == n.perfect_root.rad || next
- #n = n.pn_primorial
- #n.is_squarefree || next
- if (g3(n) != g(n)) {
- say "Not matching for #{n} -- #{g(n)} != #{g3(n)} -> #{n.factor_exp}"
- }
- }
- __END__
- #say 20.of(g)
- #say 20.of(g2)
- #var n = (43*13)**5
- say g(n)
- say g2(n)
- say g3(n)
- __END__
- # 282488053
- #say 20.of { gsum(_) }
- #say 20.of { gsum(_) }.map_cons(2, {|a,b| b-a })
- #say 20.of(g)
- for n in (1..1e4) {
- #var k = 1e20.irand
- #if (k.is_squarefree) {
- #say "Testing: #{k}"
- # assert_eq(g(k), euler_phi(k))
- #}
- if (g(n) != euler_phi(n)) {
- say n.factor
- }
- }
- #If n is squarefree (A005117), then a(n) = A000010(n) where A000010 is the Euler totient function.
- __END__
- var n = 10000
- say gsum(n)
- say gsum2(n)
- say ''
- say (n**3 * zeta(6)/(3*zeta(3)))
- say (n**3 * zeta(3)/(3*zeta(6)))
- say ''
- say (faulhaber(n, 2) * zeta(6)/(zeta(3)))
- say (faulhaber(n, 2) * zeta(3)/(zeta(6)))
- #Sum_{k=1..n} |a(k)| ~ n^3 * zeta(6)/(3*zeta(3)).
- #(n^3 * zeta(6)) / (3 * zeta(3))
- #say 60.of { gsum(_) } #.map_cons(2, {|a,b| b-a })
- #1, 4, 12, 25, 49, 73, 121, 172, 245, 317, 437, 541, 709, 853, 1045, 1250, 1538, 1757, 2117, 2429, 2813, 3173, 3701, 4109, 4710, 5214, 5870, 6494, 7334, 7910, 8870, 9689, 10649, 11513, 12665, 13614, 14982, 16062, 17406, 18630, 20310, 21462, 23310, 24870, 26622, 28206, 30414, 32054, 34407, 36210
- #1, 4, 12, 25, 49, 73, 121, 172, 245, 317, 437, 541, 709, 853, 1045, 1250, 1538, 1757, 2117, 2429, 2813, 3173, 3701, 4109, 4710, 5214, 5870, 6494, 7334, 7910, 8870, 9689, 10649, 11513, 12665, 13614, 14982, 16062, 17406, 18630, 20310, 21462, 23310, 24870, 26622, 28206, 30414, 32054, 34407, 36210, 38514, 40698, 43506, 45474, 48354, 50802, 53682, 56202, 59682
|