A161942 -id:A161942 - cp's OEIS frontend

A336937 The 2-adic valuation of sigma(n), the sum of divisors of n.

0, 0, 2, 0, 1, 2, 3, 0, 0, 1, 2, 2, 1, 3, 3, 0, 1, 0, 2, 1, 5, 2, 3, 2, 0, 1, 3, 3, 1, 3, 5, 0, 4, 1, 4, 0, 1, 2, 3, 1, 1, 5, 2, 2, 1, 3, 4, 2, 0, 0, 3, 1, 1, 3, 3, 3, 4, 1, 2, 3, 1, 5, 3, 0, 2, 4, 2, 1, 5, 4, 3, 0, 1, 1, 2, 2, 5, 3, 4, 1, 0, 1, 2, 5, 2, 2, 3, 2, 1, 1, 4, 3, 7, 4, 3, 2, 1, 0, 2, 0, 1, 3, 3, 1, 6

Offset: 1

Antti Karttunen, Aug 12 2020

nonn

Equally, the 2-adic valuation of A000593(n), the sum of odd divisors of n.

Proof for the given additive formula: It's easy to see that for all powers of 2 and all even powers of odd primes the result is zero. Thus assuming p is an odd prime, factorize sigma(p^(2e-1)) = (1 + p + p^2 + ... + p^(2e-1)) as (1+p)*(1 + u + u^2 + u^3 + ... + u^(e-1)), where u=p^2. Note that u [and its powers] are always of the form 4k+1, thus the 2-adic valuation of that sum is A007814(e) [see my Aug 15 2020 comment there] which when added to the 2-adic valuation of 1+p then gives the 2-adic valuation for whole sigma(p^(2e-1)).

Antti Karttunen, Table of n, a(n) for n = 1..65537
Index entries for sequences related to sigma(n)

One less than A286357.

Cf. A000035, A000203, A000265, A000593, A007814, A028982 (positions of zeros), A161942, A295664, A336932.

a[n_] := IntegerExponent[DivisorSigma[1, n], 2]; Array[a, 100] (* Amiram Eldar, Jul 04 2022 *)

PARI
```
A336937(n) = valuation(sigma(n),2);
```

A007814(n) = valuation(n,2);
A336937(n) = { my(f=factor(n)); sum(i=1,#f~,(f[i,1]%2) * (f[i,2]%2) * (A007814(1+f[i,1])+A007814(1+f[i,2])-1)); };

from sympy import divisor_sigma
def A336937(n): return (~(m:=int(divisor_sigma(n))) & m-1).bit_length() # Chai Wah Wu, Jul 01 2022

a(n) = A007814(A000203(n)) = A007814(A000593(n)).
Additive with a(2^e) = 0, a(p^2e) = 0, a(p^(2e-1)) = A007814(1+p) + A007814(e).
For n with factorization p_1^e_1 * ... * p_k^e_k, a(n) = Sum_{i=1..k} A000035(p_i) * A000035(e_i) * (A007814(1+p_i)+A007814(1+e_i)-1).
a(n) = A286357(n)-1.
a(n) = a(2*n) = a(A000265(n)).
For any odd n, a(n) >= A295664(n).
a(A003961(n)) = A336932(n).
A161942(n) = A000203(n) / (2^a(n)).

A350073 a(n) = A064989(sigma(n)), where A064989 is multiplicative with a(2^e) = 1 and a(p^e) = prevprime(p)^e for odd primes p.

Original entry on oeis.org

1, 2, 1, 5, 2, 2, 1, 6, 11, 4, 2, 5, 5, 2, 2, 29, 4, 22, 3, 10, 1, 4, 2, 6, 29, 10, 3, 5, 6, 4, 1, 20, 2, 8, 2, 55, 17, 6, 5, 12, 10, 2, 7, 10, 22, 4, 2, 29, 34, 58, 4, 25, 8, 6, 4, 6, 3, 12, 6, 10, 29, 2, 11, 113, 10, 4, 13, 20, 2, 4, 4, 66, 31, 34, 29, 15, 2, 10, 3, 58, 49, 20, 10, 5, 8, 14, 6, 12, 12, 44, 5, 10

Offset: 1

Antti Karttunen, Dec 12 2021

Cf. A000203, A064989, A161942, A337343 (a(n) >= n), A348748, A348749.

Cf. also A326042, A350072.

f[2, e_] := 1; f[p_, e_] := NextPrime[p, -1]^e; a[1] = 1; a[n_] := Times @@ f @@@ FactorInteger[DivisorSigma[1, n]]; Array[a, 100] (* Amiram Eldar, Dec 12 2021 *)

A064989(n) = { my(f = factor(n)); for (i=1, #f~, f[i,1] = if(2==f[i, 1],1,precprime(f[i, 1]-1))); factorback(f); };
A350073(n) = A064989(sigma(n));

Multiplicative with a(p^e) = A064989(1 + p + p^2 + ... + p^e).

a(n) = A064989(A000203(n)) = A064989(A161942(n)).

A347240 a(n) is the largest prime factor (A006530) of all terms encountered when iterating the map x -> A000593(x), when starting from x = n, but excluding the n itself. If n is a power of 2, then a(n) = 1. If 1 is never reached, then a(n) = -1.

Original entry on oeis.org

1, 1, 2, 1, 3, 2, 2, 1, 13, 3, 3, 2, 7, 2, 3, 1, 13, 13, 5, 3, 2, 3, 3, 2, 31, 7, 5, 2, 5, 3, 2, 1, 3, 13, 3, 13, 19, 5, 7, 3, 7, 2, 11, 3, 13, 3, 3, 2, 19, 31, 13, 7, 5, 5, 13, 2, 5, 5, 5, 3, 31, 2, 13, 1, 7, 3, 17, 13, 3, 3, 13, 13, 37, 19, 31, 5, 3, 7, 5, 3, 19, 7, 7, 2, 5, 11, 5, 3, 13, 13, 7, 3, 2, 3, 5, 2, 19

Offset: 1

Antti Karttunen, Aug 28 2021

nonn

			For n = 17, the iteration proceeds as follows 17 -> 18 (= 2*3*3), 18 -> 13 (13 is a prime), 13 -> 14 (= 2*7), 14 -> 8 (= 2*2*2), 8 -> 1. The largest prime factor present after the initial step is 13, thus a(17) = 13.

Antti Karttunen, Table of n, a(n) for n = 1..20000

Cf. A000593, A006530, A336361, A347241, A347242, A347243, A347244.

Cf. also A161942.

A006530(n) = if(1==n, n, my(f=factor(n)); f[#f~, 1]);
A000265(n) = (n >> valuation(n, 2));
A000593(n) = sigma(A000265(n));
A347240(n) = { my(m=1); while(n>1, n = A000593(n); m = max(m, A006530(n))); (m); };

a(n) = A347241(A000593(n)). - Antti Karttunen, Feb 10 2022

A348992 a(n) = A000265(sigma(n)) / gcd(sigma(n), A003961(n)), where A003961(n) is fully multiplicative with a(prime(k)) = prime(k+1), and sigma is the sum of divisors function.

Original entry on oeis.org

1, 1, 1, 7, 3, 1, 1, 5, 13, 3, 3, 7, 7, 1, 3, 31, 9, 13, 5, 1, 1, 3, 3, 1, 31, 7, 1, 7, 15, 3, 1, 7, 3, 9, 3, 91, 19, 5, 7, 5, 21, 1, 11, 7, 39, 3, 3, 31, 57, 31, 9, 49, 27, 1, 9, 5, 1, 15, 15, 1, 31, 1, 13, 127, 3, 3, 17, 7, 3, 3, 9, 13, 37, 19, 31, 35, 3, 7, 5, 31, 121, 21, 21, 7, 27, 11, 3, 5, 45, 39, 7, 7, 1, 3

Offset: 1

Antti Karttunen, Nov 10 2021

Denominator of ratio A003961(n) / A161942(n).

Odd part of A349162.
Cf. A000203, A000265, A003961, A161942, A342671, A348993, A349169 (where equal to A348990).
Cf. A349161 (numerators).

Array[#1/(2^IntegerExponent[#1, 2]*GCD[##]) & @@ {DivisorSigma[1, #], Times @@ Map[NextPrime[#1]^#2 & @@ # &, FactorInteger[#]]} &, 94] (* Michael De Vlieger, Nov 11 2021 *)

A000265(n) = (n >> valuation(n, 2));
A003961(n) = { my(f = factor(n)); for (i=1, #f~, f[i, 1] = nextprime(f[i, 1]+1)); factorback(f); };
A348992(n) = { my(s=sigma(n)); (A000265(s)/gcd(s,A003961(n))); };

a(n) = A161942(n) / A342671(n) = A000265(A349162(n)).

a(n) = A003961(A348993(n)).

A336701 Numbers k for which A000265(1+A000265(sigma(k))) is equal to A000265(1+k).

Original entry on oeis.org

1, 3, 7, 15, 31, 127, 1023, 8191, 34335, 57855, 131071, 524287, 2147483647

Offset: 1

Antti Karttunen, Aug 02 2020

Numbers k such that A336698(k) [= A000265(1+A161942(k))] is equal to A000265(1+k).
Numbers k such that A337194(k) = 2^e * A000265(1+k), for some e >= 1, where that e = A337195(k).
Any odd perfect number would trivially satisfy this condition.
Also, all hypothetical quasiperfect numbers, numbers k that satisfy sigma(k) = 2k+1, would be members.
Question: Is A066175 a subsequence of this sequence?
From Antti Karttunen, Aug 23 2020: (Start)
Numbers k such that (1+k) = 2^e * A336698(k), for some e >= 0.
Thus numbers k such that for some e >= 0, (1+k) = 2^(e-A337195(k)) * A337194(k), or equally, that A337194(k) = 2^(A337195(k)-e) * (1+k).
Conjecture: There are no even terms. This is equivalent to claim that there are no k such that A336698(k) = 1+k: If we assume that k is even, then in above equations we set e=0, and the requirement will then become that A337194(k) = 2^A337195(k)*(1+k), thus 1+k = A336698(k) = A000265(1+A000265(sigma(k))).
(End)

Paolo Cattaneo, Sui numeri quasiperfetti, Bollettino dell’Unione Matematica Italiana, Serie 3, Vol.6(1951), n.1, p. 59-62.
P. Hagis and G. L. Cohen, Some Results Concerning Quasiperfect Numbers, J. Austral. Math. Soc. Ser. A 33, 275-286, 1982.
V. Siva Rama Prasad and C. Sunitha, On quasiperfect numbers, Notes on Number Theory and Discrete Mathematics, Vol. 23, 2017, No. 3, 73-78.
Eric Weisstein's World of Mathematics, Quasiperfect Number
Index entries for sequences where odd perfect numbers must occur, if they exist at all
Index entries for sequences related to sigma(n)

Subsequence of A336700.
Cf. A000203, A000265, A066175, A161942, A336698, A336699, A336937, A337194, A337195.
Cf. A000668 (a subsequence).
See also comments in A326042, A332223.

Block[{f}, f[n_] := n/2^IntegerExponent[n, 2]; Select[Range[2^20], f[1 + f[DivisorSigma[1, #]]] == f[1 + #] &] ] (* Michael De Vlieger, Aug 22 2020 *)

A000265(n)  = (n>>valuation(n,2));
isA336701(n) = (A000265(1+A000265(sigma(n))) == A000265(1+n));

A082903 Highest power of two that divides the sum of divisors of n.

Original entry on oeis.org

1, 1, 4, 1, 2, 4, 8, 1, 1, 2, 4, 4, 2, 8, 8, 1, 2, 1, 4, 2, 32, 4, 8, 4, 1, 2, 8, 8, 2, 8, 32, 1, 16, 2, 16, 1, 2, 4, 8, 2, 2, 32, 4, 4, 2, 8, 16, 4, 1, 1, 8, 2, 2, 8, 8, 8, 16, 2, 4, 8, 2, 32, 8, 1, 4, 16, 4, 2, 32, 16, 8, 1, 2, 2, 4, 4, 32, 8, 16, 2, 1, 2, 4, 32, 4, 4, 8, 4, 2, 2, 16, 8, 128, 16, 8, 4, 2

Offset: 1

Labos Elemer, Apr 22 2003

a(n) = gcd(2^n, sigma_1(n)) = gcd(A000079(n), A000203(n)) also a(n) = gcd(2^n, sigma_3(n)) = gcd(A000079(n), A001158(n)). (The original, equivalent definition for this sequence).
a(n) = gcd(2^n, sigma_k(n)) when k is an odd positive integer. Proof: It suffices to show that v_2(sigma_k(n)) does not depend on k, where v_2(n) is the 2-adic valuation of n. Since v_2(ab) = v_2(a)+v_2(b) and sigma_k(n) is an arithmetic function, we need only prove it for n=p^e with p prime. If p is 2 or e is even, sigma_k(p^e) is odd, so we can disregard those cases. Otherwise, we sum the geometric series to obtain v_2(sigma_k(p^e)) = v_2(p^(k(e+1))-1)-v_2(p-1). Applying the well-known LTE Lemma (see Hossein link) Case 4 arrives at v_2(p^(k(e+1))-1)-v_2(p-1) = v_2(p+1)+v_2(k(e+1))-1. But v_2(k(e+1)) = v_2(k)+v_2(e+1), and k is odd, so we conclude that v_2(sigma_k(p^e)) = v_2(p+1)+v_2(e+1)-1, a result independent of k. - Rafay A. Ashary, Oct 15 2016
Also the highest power of two that divides the sum of odd divisors of n. - Antti Karttunen, Mar 27 2022

Robert Israel, Table of n, a(n) for n = 1..10000
Amir Hossein, Lifting The Exponent Lemma (Containing PDF file)
Jon Maiga, Computer-generated formulas for A082903, Sequence Machine.

Cf. A000203, A000265, A000593, A001157, A001158, A006519, A007814, A082902, A161942, A286357.

Cf. A028982 (positions of 1's).

seq(2^min(n, padic:-ordp(numtheory:-sigma(n),2)), n=1..100); # Robert Israel, Oct 23 2016

Array[2^IntegerExponent[DivisorSigma[1, #], 2] &, 97] (* Michael De Vlieger, Apr 03 2022 *)

a(n) = gcd(2^n, sigma(n)); \\ Michel Marcus, Oct 15 2016

A082903(n) = (2^valuation(sigma(n), 2)); \\ Antti Karttunen, Mar 27 2022

from sympy import divisor_sigma
def A082903(n): return 1<<(~(m:=int(divisor_sigma(n))) & m-1).bit_length() # Chai Wah Wu, Jul 02 2022

From Antti Karttunen, Mar 27 2022: (Start)
(Some of these formulas were found by Sequence Machine.)
a(n) = a(A000265(n)) = a(2*n).
a(n) = A006519(A000593(n)) = A006519(A000203(n)) = A000203(n) / A161942(n).
a(n) = 2^(A286357(n)-1).
(End)
a(n) = 2^A336937(n). - Chai Wah Wu, Jul 02 2022

Name replaced with a simpler one and the original definition moved to the Comments section by Antti Karttunen, Apr 03 2022

A286460 Compound filter (2-adic valuation & sum of the divisors): a(n) = P(A001511(n), A000203(n)), where P(n,k) is sequence A000027 used as a pairing function.

Original entry on oeis.org

1, 8, 7, 39, 16, 80, 29, 157, 79, 173, 67, 438, 92, 302, 277, 600, 154, 782, 191, 949, 497, 668, 277, 1957, 466, 905, 781, 1656, 436, 2630, 497, 2284, 1129, 1487, 1129, 4281, 704, 1832, 1541, 4282, 862, 4658, 947, 3658, 3004, 2630, 1129, 8133, 1597, 4373, 2557, 4953, 1432, 7262, 2557, 7507, 3161, 4097, 1771, 14368, 1892, 4658, 5357, 8785, 3487, 10442, 2279

Offset: 1

Antti Karttunen, May 10 2017

nonn

Antti Karttunen, Table of n, a(n) for n = 1..10000
Eric Weisstein's World of Mathematics, Pairing Function
Index entries for sequences related to sigma(n)

Cf. A000027, A286360.

Cf. A000593, A146076 (sequences matching to this filter), also A000203, A161942, A286260, A286357.

A000203(n) = sigma(n);
A001511(n) = (1+valuation(n,2));
A286460(n) = (1/2)*(2 + ((A001511(n)+A000203(n))^2) - A001511(n) - 3*A000203(n));
for(n=1, 10000, write("b286460.txt", n, " ", A286460(n)));

from sympy import divisor_sigma as D
def a001511(n): return bin(n)[2:][::-1].index("1") + 1
def T(n, m): return ((n + m)**2 - n - 3*m + 2)//2
def a(n): return T(a001511(n), D(n))
print([a(n) for n in range(1, 101)]) # Indranil Ghosh, May 12 2017

(define (A286460 n) (* (/ 1 2) (+ (expt (+ (A001511 n) (A000203 n)) 2) (- (A001511 n)) (- (* 3 (A000203 n))) 2)))

a(n) = (1/2)*(2 + ((A001511(n)+A000203(n))^2) - A001511(n) - 3*A000203(n)).

A336699 a(n) = A000265(1+A000265(sigma(A000265(n)))), where A000265(k) gives the odd part of k, and sigma is the sum of divisors function.

Original entry on oeis.org

1, 1, 1, 1, 1, 1, 1, 1, 7, 1, 1, 1, 1, 1, 1, 1, 5, 7, 3, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 1, 1, 1, 5, 1, 7, 5, 3, 1, 1, 11, 1, 3, 1, 5, 1, 1, 1, 29, 1, 5, 1, 7, 3, 5, 1, 3, 1, 1, 1, 1, 1, 7, 1, 11, 1, 9, 5, 1, 1, 5, 7, 19, 5, 1, 3, 1, 1, 3, 1, 61, 11, 11, 1, 7, 3, 1, 1, 23, 5, 1, 1, 1, 1, 1, 1, 25, 29, 5, 1, 13, 5, 7, 1, 1

Offset: 1

Antti Karttunen, Aug 02 2020

See the "lacunae" in the scatter plot. - Antti Karttunen, Mar 27 2022

Antti Karttunen, Table of n, a(n) for n = 1..16384
Antti Karttunen, Data supplement: n, a(n) computed for n = 1..65537
Jon Maiga, Computer-generated formulas for A336699, Sequence Machine.
Index entries for sequences related to sigma(n)

Cf. A000203, A000265, A000593, A002131, A082903, A161942, A171435, A322250, A336698, A336700, A336701, A336844 [= a(A003961(n))], A351565.

A000265(n) = (n>>valuation(n,2));
A336699(n) = A000265(1+A000265(sigma(A000265(n))));

a(n) = A000265(1+A000265(A000593(n))) = A000265(1+A161942(A000265(n))).
a(n) = A336698(A000265(n)).
From Antti Karttunen, Mar 27 2022: (Start)
a(n) = A351565(A000593(n)).
[The following formulas were discovered by Sequence Machine]:
a(n) = A351565(A002131(n)) = A000265(1+A000265(A002131(n))).
a(n) = A336698(1+A322250(n)).
a(n) = A171435(A000593(n)+A082903(n)).
(End)

A347242 Numbers k such that when iterating the map x -> A000593(x), at some point before 1 is reached (after starting from x=k), a term is encountered whose largest prime factor is at least as large as A006530(k).

Original entry on oeis.org

9, 18, 25, 27, 36, 45, 49, 50, 54, 55, 63, 72, 75, 81, 90, 98, 99, 100, 108, 110, 117, 121, 125, 126, 135, 144, 147, 150, 162, 165, 169, 175, 180, 196, 198, 200, 216, 220, 225, 234, 242, 243, 245, 250, 252, 270, 275, 288, 289, 294, 300, 315, 324, 325, 330, 338, 343, 350, 360, 361, 363, 375, 385, 392, 396, 400

Offset: 1

Antti Karttunen, Aug 28 2021

nonn

Provided there do not exist any odd perfect numbers, these are numbers k for which A347240(k) >= A006530(k), as for any odd perfect number x, A347240(x) = -1 by its escape clause.
If k is included as a term, then 2*k is also present.
Not all odd squares of primes are present. For example, 67^2 and 79^2 are not included. See also A091490, which seems to be a subsequence of those exceptions.
Conjecture: There are no primes in this sequence. Checked up to the 2^20-th prime, 16290047.

			For n = 55 = 5*11, on the first iteration we get A000593(55) = 72 = 2^3 * 3^2, but both 2 and 3 are less than 11; therefore we iterate a second time to get A000593(72) = 13, which is the first value whose largest prime factor is larger than that of 55 (13 > 11), thus 55 is included in the sequence.

Index entries for sequences where odd perfect numbers must occur, if they exist at all

Cf. A000593, A006530, A091490, A161942, A336353, A336361, A347240, A347241, A347243 (complement), A347244 (characteristic function).

Positions of nonzero terms in A347245.

A006530(n) = if(1==n, n, my(f=factor(n)); f[#f~, 1]);
A000265(n) = (n >> valuation(n, 2));
A000593(n) = sigma(A000265(n));
A347244(n) = { my(gpf=A006530(n)); while(n>1, n = A000593(n); if(A006530(n)>=gpf,return(1))); (0); };
isA347242(n) = A347244(n);

A348739 Numbers k for which A326042(k) > k.

Original entry on oeis.org

4, 9, 12, 16, 18, 25, 32, 36, 44, 48, 49, 64, 72, 81, 96, 99, 100, 108, 124, 144, 147, 162, 169, 176, 180, 192, 196, 225, 236, 243, 252, 256, 279, 284, 288, 300, 320, 324, 361, 372, 396, 400, 405, 432, 441, 448, 450, 468, 484, 486, 496, 507, 512, 529, 531, 567, 576, 588, 604, 612, 625, 639, 648, 675, 676, 700, 704

Offset: 1

Antti Karttunen, Nov 02 2021

nonn

Terms that occur also in A337386 are: 180, 300, 720, 900, 960, 1008, 1200, 1440, 1620, 1800, 2016, 2400, ...

Positions of negative terms in A348736.
Cf. A326182 (subsequence after its initial 1), A348738.
Cf. A000203, A003961, A064989, A161942, A191218, A326042, A337386, A348742, A348749 (corresponding odd numbers), A348942.

f1[2, e_] := 1; f1[p_, e_] := NextPrime[p, -1]^e; s1[1] = 1; s1[n_] := Times @@ f1 @@@ FactorInteger[n]; f2[p_, e_] := NextPrime[p]^e; s2[1] = 1; s2[n_] := Times @@ f2 @@@ FactorInteger[n]; Select[Range[700], s1[DivisorSigma[1, s2[#]]] > # &] (* Amiram Eldar, Nov 04 2021 *)

A003961(n) = my(f = factor(n)); for (i=1, #f~, f[i, 1] = nextprime(f[i, 1]+1)); factorback(f); \\ From A003961
A064989(n) = {my(f); f = factor(n); if((n>1 && f[1,1]==2), f[1,2] = 0); for (i=1, #f~, f[i,1] = precprime(f[i,1]-1)); factorback(f)};
A326042(n) = A064989(sigma(A003961(n)));
isA348739(n) = (A326042(n)>n);

cp's OEIS Frontend

A336937 The 2-adic valuation of sigma(n), the sum of divisors of n.

Views

Author

Keywords

Comments

Links

Crossrefs

Programs

Mathematica

PARI

PARI

Python

Formula

A350073 a(n) = A064989(sigma(n)), where A064989 is multiplicative with a(2^e) = 1 and a(p^e) = prevprime(p)^e for odd primes p.

Views

Author

Keywords

Links

Crossrefs

Programs

Mathematica

PARI

Formula

A347240 a(n) is the largest prime factor (A006530) of all terms encountered when iterating the map x -> A000593(x), when starting from x = n, but excluding the n itself. If n is a power of 2, then a(n) = 1. If 1 is never reached, then a(n) = -1.

Views

Author

Keywords

Examples

Links

Crossrefs

Programs

PARI

Formula

A348992 a(n) = A000265(sigma(n)) / gcd(sigma(n), A003961(n)), where A003961(n) is fully multiplicative with a(prime(k)) = prime(k+1), and sigma is the sum of divisors function.

Views

Author

Keywords

Comments

Links

Crossrefs

Programs

Mathematica

PARI

Formula

A336701 Numbers k for which A000265(1+A000265(sigma(k))) is equal to A000265(1+k).

Views

Author

Keywords

Comments

Links

Crossrefs

Programs

Mathematica

PARI

A082903 Highest power of two that divides the sum of divisors of n.

Views

Author

Keywords

Comments

Links

Crossrefs

Programs

Maple

Mathematica

PARI

PARI

Python

Formula

Extensions

A286460 Compound filter (2-adic valuation & sum of the divisors): a(n) = P(A001511(n), A000203(n)), where P(n,k) is sequence A000027 used as a pairing function.

Views

Author

Keywords

Links

Crossrefs

Programs

PARI

Python

Scheme