A135577 -id:A135577 - cp's OEIS frontend

A001576 a(n) = 1^n + 2^n + 4^n.

3, 7, 21, 73, 273, 1057, 4161, 16513, 65793, 262657, 1049601, 4196353, 16781313, 67117057, 268451841, 1073774593, 4295032833, 17180000257, 68719738881, 274878431233, 1099512676353, 4398048608257, 17592190238721, 70368752566273, 281474993487873, 1125899940397057

Offset: 0

N. J. A. Sloane

Equals A135576, except for the first term. - Omar E. Pol, Nov 18 2008
Conjecture: For n > 1, if a(n) = 1^n + 2^n + 4^n is a prime number then n is of the form 3^h. For example, for h=1, n=3, a(n) = 1^3 + 2^3 + 4^3 = 73 (prime); for h=2, n=9, a(n) = 1^9 + 2^9 + 4^9 = 262657 (prime); for h=3, n=27, a(n) is not prime. - Vincenzo Librandi, Aug 03 2010
The previous conjecture was proved by Golomb in 1978. See A051154. - T. D. Noe, Aug 15 2010
Another more elementary proof can be found in Liu link. - Bernard Schott, Mar 08 2019
Fills in one quarter section of the figurate form of the Sierpinski square curve. See illustration in links and A141725. - John Elias, Mar 29 2023

T. D. Noe, Table of n, a(n) for n = 0..200
John Elias, Illustration of Initial Terms: 1/4 Sierpinski Square Curve
Andy Liu, West German Mathematical Olympiad 1982 - Second round, Problem 4, Crux Mathematicorum, p. 105, Vol. 12, May. 86.
Index entries for linear recurrences with constant coefficients, signature (7,-14,8).

Subsequence of A002061.
Cf. A001550, A034513, A001579, A074501-A074580, A135576, A135577.
See also comments in A051154.
Cf. A006095, A022166.

Mathematica
```
Table[1^n + 2^n + 4^n, {n, 0, 24}]
```

a(n)=1+2^n+4^n \\ Charles R Greathouse IV, Jun 10 2011

[sigma(4,n)for n in range(0,23)] # Zerinvary Lajos, Jun 04 2009

a(n) = 6*a(n-1) - 8*a(n-2) + 3.
O.g.f.: -1/(-1+x) - 1/(-1+2*x) - 1/(-1+4*x) = ( -3+14*x-14*x^2 ) / ( (x-1)*(2*x-1)*(4*x-1) ). - R. J. Mathar, Feb 29 2008
E.g.f.: e^x + e^(2*x) + e^(4*x). - Mohammad K. Azarian, Dec 26 2008
a(n) = A024088(n)/A000225(n). - Reinhard Zumkeller, Feb 15 2009
Exp( Sum_{n >= 1} a(n)*x^n/n ) = 1 + 7*x + 35*x^2 + 155*x^3 + ... is the o.g.f. for the 2nd subdiagonal of triangle A022166, essentially A006095. - Peter Bala, Apr 07 2015

A333237 Numbers k such that 1/k contains at least one '9' in its decimal expansion.

Original entry on oeis.org

11, 13, 17, 19, 21, 23, 29, 31, 34, 38, 41, 42, 43, 46, 47, 49, 51, 52, 53, 57, 58, 59, 61, 62, 67, 68, 69, 71, 73, 76, 77, 81, 82, 83, 84, 85, 86, 87, 89, 91, 92, 94, 95, 97, 98, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 118

Offset: 1

Andrew Slattery, Mar 12 2020

Almost every prime appears in this sequence.
Among the first 10000 primes, only 2, 3, 5, 7, 37, 79, 239, 4649, and 62003 do not appear in the sequence. - Giovanni Resta, Mar 13 2020
The next primes not in the sequence are 538987, 35121409, and 265371653. - Robert Israel, Mar 18 2020

			5 is not in the sequence because 1/5 = 0.2 does not contain any 9s.

Robert Israel, Table of n, a(n) for n = 1..10000
Index entries for sequences related to decimal expansion of 1/n

Cf. A333236.
Subsequences (for terms > 1): A000533, A002275, A135577, A252491.
Cf. A216664 (a subsequence).
Cf. A187614.

f:= proc(n) local m,S,r;
   m:= 1; S:= {1};
   do
     r:= floor(m/n);
     if r = 9 then return true fi;
     m:= (m - r*n)*10;
     if member(m,S) then return false fi;
     S:= S union {m};
   od
end proc:
select(f, [$1..1000]); # Robert Israel, Mar 18 2020

Select[Range[120], MemberQ[ Flatten@ RealDigits[1/#][[1]], 9] &] (* Giovanni Resta, Mar 12 2020 *)

from itertools import count, islice
from sympy import n_order, multiplicity
def A333237_gen(startvalue=1): # generator of terms
    for m in count(max(startvalue,1)):
        m2, m5 = multiplicity(2,m), multiplicity(5,m)
        if max(str(10**(max(m2,m5)+n_order(10,m//2**m2//5**m5))//m)) == '9':
            yield m
A333237_list = list(islice(A333237_gen(), 10)) # Chai Wah Wu, Feb 07 2022

A333236(a(n)) = 9.

More terms from Giovanni Resta, Mar 12 2020

A138721 Concatenation of n digits 1, n digits 0 and n digits 1.

Original entry on oeis.org

101, 110011, 111000111, 111100001111, 111110000011111, 111111000000111111, 111111100000001111111, 111111110000000011111111, 111111111000000000111111111, 111111111100000000001111111111, 111111111110000000000011111111111, 111111111111000000000000111111111111

Offset: 1

Omar E. Pol, Mar 29 2008

a(n) is also A145641(n) written in base 2. - Omar E. Pol, Oct 15 2008

a(n) has 3n digits. - Omar E. Pol, Nov 12 2008

			From _Omar E. Pol_, Nov 12 2008: (Start)
n         Successive digits of a(n)
1                 ( 1 0 1 )
2              ( 1 1 0 0 1 1 )
3           ( 1 1 1 0 0 0 1 1 1 )
4        ( 1 1 1 1 0 0 0 0 1 1 1 1 )
5     ( 1 1 1 1 1 0 0 0 0 0 1 1 1 1 1 )
(End)

Paolo Xausa, Table of n, a(n) for n = 1..300
Index entries for linear recurrences with constant coefficients, signature (1111,-112110,1111000,-1000000).

Cf. A000533, A135577, A138120, A138144, A138145, A138146, A138826, A145641, A147757, A147759.

a:= n-> parse(cat(1$n,0$n,1$n)):
seq(a(n), n=1..14);  # Alois P. Heinz, Mar 03 2022

Table[(100^n + 1)*(10^n - 1)/9, {n, 15}] (* Paolo Xausa, Aug 02 2024 *)

Vec(x*(101000*x^2-2200*x+101)/((x-1)*(10*x-1)*(100*x-1)*(1000*x-1)) + O(x^100)) \\ Colin Barker, Sep 16 2013

G.f.: x*(101000*x^2 - 2200*x + 101) / ((x-1)*(10*x-1)*(100*x-1)*(1000*x-1)). - Colin Barker, Sep 16 2013

a(n) = (100^n+1)*(10^n-1)/9. - Paolo Xausa, Aug 02 2024

A333402 Numbers m such that the largest digit in the decimal expansion of 1/m is 1.

Original entry on oeis.org

1, 9, 10, 90, 99, 100, 900, 909, 990, 999, 1000, 9000, 9009, 9090, 9900, 9990, 9999, 10000, 90000, 90009, 90090, 90900, 90909, 99000, 99900, 99990, 99999, 100000, 900000, 900009, 900090, 900900, 909000, 909090, 990000, 990099, 999000, 999900, 999990, 999999, 1000000

Offset: 1

Bernard Schott, Mar 19 2020

If m is a term, 10*m is also a term.
If m is a term then m has only digits {1}, {9}, {1,0} or {9,0} in its decimal representation, but this is not sufficient to be a term (see examples).
Some subsequences below (not exhaustive, see crossrefs):
m = 10^k, k >= 0, hence m is in A011557 = {1, 10, 100, 1000, 10000, ...};
m = 9*10^k, k >= 0, hence m is in A052268 = {9, 90, 900, 9000, 90000, ...};
m = 10^k-1, k >= 1, hence m is in A002283 = {9, 99, 999, 9999, 99999, ...};
m = 9*(10^k+1), k >= 1, hence m is in 9*A000533 = {99, 909, 9009, 90009, ...};
m = 9+100*(100^k-1)/11, k >= 0, hence m is in 9*A094028 = {9, 909, 90909, 9090909, ...}.

			As 1/101 = 0.009900990099..., 101 is not a term.
As 1/909 = 0.001100110011..., 909 is a term.
As 1/9099 = 0.000109902187..., 9099 is not a term.
As 1/9999 = 0.000100010001..., 9999 is also a term.

Index entries for sequences related to decimal expansion of 1/n

Cf. A333236, A333237 (similar, with 9).
Subsequences: A002283, A011557, A052268.
Subsequences: 9*A000533, 9*A094028, 9*A135577, 9*A261544, 9*A330135.

Select[Range[10^4], Max @ RealDigits[1/#][[1]] == 1 &] (* Amiram Eldar, Mar 19 2020 *)

from itertools import count, islice
def A333402_gen(startvalue=1): # generator of terms >= startvalue
    for m in count(max(startvalue,1)):
        k = 1
        while k <= m:
            k *= 10
        rset = {0}
        while True:
            k, r = divmod(k, m)
            if max(str(k)) > '1':
                break
            else:
                if r in rset:
                    yield m
                    break
            rset.add(r)
            k = r
            while k <= m:
                k *= 10
A333402_list = list(islice(A333402_gen(),30)) # Chai Wah Wu, Feb 17 2022

A333236(a(n))= 1.

More terms from Jinyuan Wang, Mar 19 2020

A066138 a(n) = 10^(2*n) + 10^n + 1.

Original entry on oeis.org

3, 111, 10101, 1001001, 100010001, 10000100001, 1000001000001, 100000010000001, 10000000100000001, 1000000001000000001, 100000000010000000001, 10000000000100000000001, 1000000000001000000000001, 100000000000010000000000001, 10000000000000100000000000001, 1000000000000001000000000000001

Offset: 0

Henry Bottomley, Dec 07 2001

Palindromes whose digit sum is 3.
Essentially the same as A135577. - R. J. Mathar Apr 29 2008
From Peter Bala, Sep 25 2015: (Start)
For n >= 1, the simple continued fraction expansion of sqrt(a(n)) = [10^n; 1, 1, 2/3*(10^n - 1), 1, 1, 2*10^n, ...] has period 6. As n increases, the expansion has the large partial quotients 2/3*(10^n - 1) and 2*10^n.
For n >= 1, the continued fraction expansion of sqrt(a(2*n))/a(n) = [0; 1, 10^n - 1, 1, 1, 1/3*(10^n - 4), 1, 4, 1, 1/3*(10^n - 4), 1, 1, 10^n - 1, 2, 10^n - 1, ...] has pre-period of length 3 and period 12 beginning 1, 1, 1/3*(10^n - 4), ....  As n increases, the expansion has the large partial quotients 10^n - 1 and 1/3*(10^n - 4).
A theorem of Kuzmin in the measure theory of continued fractions says that large partial quotients are the exception in continued fraction expansions.
Empirically, we also see exceptionally large partial quotients in the continued fraction expansions of the m-th root of the numbers a(m*n), for m >= 3. For example, it appears that the continued fraction expansion of a(3*n)^(1/3), for n >= 2, begins [10^(2*n); 3*10^n - 1, 1, 0.5*10^(2*n) - 1, 1.44*10^n - 1, 1, ...]. Cf. A000533, A002283 and A168624. (End)

			From _Peter Bala_, Sep 25 2015: (Start)
Simple continued fraction expansions showing large partial quotients:
a(9)^(1/3) =[1000000; 2999, 1, 499999, 1439, 1, 2582643, 1, 1, 1, 2, 3, 3, ...].
a(20)^(1/4) = [10000000000; 39999999999, 1, 3999999999, 16949152542, 2, 1, 2, 6, 1, 4872106, 3, 9, 2, 3, ...].
a(25)^(1/5) = [10000000000; 4999999999999999, 1, 3333333332, 2, 1, 217391304347825, 2, 2, 1, 1, 1, 2, 1, 23980814, 1, 1, 1, 1, 1, 7, ...]. (End)

Harry J. Smith, Table of n, a(n) for n=0..100
Index entries for linear recurrences with constant coefficients, signature (111,-1110,1000).

Cf. A000533, A002283, A033934, A062397, A074992, A135577, A168624.

[10^(2*n) + 10^n + 1: n in [0..20]]; // Vincenzo Librandi, Sep 27 2015

Table[10^(2 n) + 10^n + 1, {n, 0, 15}] (* Michael De Vlieger, Sep 27 2015 *)
CoefficientList[Series[(3 - 222 x + 1110 x^2)/((1 - 100 x) (1 - 10 x) (1 - x)), {x, 0, 33}], x] (* Vincenzo Librandi, Sep 27 2015 *)

a(n) = { 10^(2*n) + 10^n + 1 } \\ Harry J. Smith, Feb 02 2010

Vec(-3*(370*x^2-74*x+1)/((x-1)*(10*x-1)*(100*x-1)) + O(x^20)) \\ Colin Barker, Sep 27 2015

A168624(n) = a(2*n)/a(n). - Peter Bala, Sep 24 2015
G.f.: (3 - 222*x + 1110*x^2)/((1 - 100*x)*(1 - 10*x)*(1 - x)). - Vincenzo Librandi, Sep 27 2015
From Colin Barker, Sep 27 2015: (Start)
a(n) = 111*a(n-1) - 1110*a(n-2) + 1000*a(n-3) for n > 2.
G.f.: -3*(370*x^2-74*x+1)/((x-1)*(10*x-1)*(100*x-1)). (End)
From Elmo R. Oliveira, Aug 27 2024: (Start)
E.g.f.: exp(x)*(exp(99*x) + exp(9*x) + 1).
a(n) = 3*A074992(n). (End)

Offset changed from 1 to 0 by Harry J. Smith, Feb 02 2010

More terms from Michael De Vlieger, Sep 27 2015

A135576 Numbers whose binary expansion has only the digit "1" as first, central and final digit.

Original entry on oeis.org

1, 7, 21, 73, 273, 1057, 4161, 16513, 65793, 262657, 1049601, 4196353, 16781313, 67117057, 268451841, 1073774593, 4295032833, 17180000257, 68719738881, 274878431233, 1099512676353, 4398048608257, 17592190238721

Offset: 1

Omar E. Pol, Feb 24 2008

This sequence is essentially identical to A001576.

a(n) is the number whose binary representation is A135577(n), (See example). - Omar E. Pol, Nov 18 2008

			--------------------------------------
n ........ a(n) ..... a(n) in base 2
--------------------------------------
1 .......... 1 ............ 1
2 .......... 7 ........... 111
3 ......... 21 .......... 10101
4 ......... 73 ......... 1001001
5 ........ 273 ........ 100010001
6 ....... 1057 ....... 10000100001
7 ....... 4161 ...... 1000001000001
8 ...... 16513 ..... 100000010000001
9 ...... 65793 .... 10000000100000001
10 .... 262657 ... 1000000001000000001

G. C. Greubel, Table of n, a(n) for n = 1..1000
Index entries for linear recurrences with constant coefficients, signature (7,-14,8).

Cf. A001576, A135577.

Subsequence of A006995.

nxt[n_]:=Module[{l=Floor[IntegerLength[n,2]/2]},FromDigits[Join[{1},Table[0,{l}],{1},Table[0,{l}],{1}],2]]
NestList[nxt,1,25] (* Harvey P. Dale, Dec 29 2010 *)
Join[{1},LinearRecurrence[{7,-14,8},{7,21,73},30]] (* Harvey P. Dale, Mar 22 2015 *)

a(n)=if(n--,4^n+2^n+1,1) \\ Charles R Greathouse IV, Dec 28 2012

a(1)=1. If n>1 then a(n) = A001576(n-1).

G.f.: -x*(16*x^3-14*x^2+1) / ((x-1)*(2*x-1)*(4*x-1)). - Colin Barker, Sep 16 2013

A147759 Palindromes formed from the reflected decimal expansion of the infinite concatenation of 1's and 0's.

Original entry on oeis.org

1, 11, 101, 1001, 10101, 101101, 1010101, 10100101, 101010101, 1010110101, 10101010101, 101010010101, 1010101010101, 10101011010101, 101010101010101, 1010101001010101, 10101010101010101, 101010101101010101

Offset: 1

Omar E. Pol, Nov 11 2008

a(k(n)) is divisible by 3 iff k(n) is defined by k(1) = 5 and k(n+1) - k(n) = A100285(n+2). - Altug Alkan, Dec 05 2015

			n .... Successive digits of a(n)
1 ............. ( 1 )
2 ............ ( 1 1 )
3 ........... ( 1 0 1 )
4 .......... ( 1 0 0 1 )
5 ......... ( 1 0 1 0 1 )
6 ........ ( 1 0 1 1 0 1 )
7 ....... ( 1 0 1 0 1 0 1 )
8 ...... ( 1 0 1 0 0 1 0 1 )
9 ..... ( 1 0 1 0 1 0 1 0 1 )
10 ... ( 1 0 1 0 1 1 0 1 0 1 )

Matthew Schulz, Table of n, a(n) for n = 1..1000
Index entries for linear recurrences with constant coefficients, signature (11,-20,110,-100).

Cf. A000533, A094028, A135577, A138120, A138144, A138145, A138146, A138721, A138826, A147757.

I:=[1,11,101,1001]; [n le 4 select I[n] else 11*Self(n-1)-20*Self(n-2)+110*Self(n-3)-100*Self(n-4): n in [1..30]]; // Vincenzo Librandi, Dec 05 2015

CoefficientList[Series[x/((1 - x) (1 - 10 x) (1 + 10 x^2)),{x, 0, 20}], x] (* Vincenzo Librandi, Dec 05 2015 *)
LinearRecurrence[{11,-20,110,-100},{1,11,101,1001},30] (* Harvey P. Dale, Apr 10 2022 *)

Vec(x/((1-x)*(1-10*x)*(1+10*x^2)) + O(x^30)) \\ Michel Marcus, Dec 05 2015

From R. J. Mathar, Feb 20 2009: (Start)
a(n) = 11*a(n-1)-20*a(n-2)+110*a(n-3)-100*a(n-4).
G.f.: x/((1-x)*(1-10*x)*(1+10*x^2)). (End)
E.g.f.: (exp(x)*(10*exp(9*x) - 1) - 9*cos(sqrt(10)*x))/99. - Stefano Spezia, Oct 12 2024

A147757 Palindromes formed from the reflected decimal expansion of the concatenation of 1, 0 and infinite digits 1.

Original entry on oeis.org

1, 11, 101, 1001, 10101, 101101, 1011101, 10111101, 101111101, 1011111101, 10111111101, 101111111101, 1011111111101, 10111111111101, 101111111111101, 1011111111111101, 10111111111111101, 101111111111111101

Offset: 1

Omar E. Pol, Nov 11 2008

a(n) is also A147758(n) written in base 2.

a(A016789(n)) is divisible by 3 for n > 0. - Altug Alkan, Dec 06 2015

			n .... Successive digits of a(n)
1 ............. ( 1 )
2 ............ ( 1 1 )
3 ........... ( 1 0 1 )
4 .......... ( 1 0 0 1 )
5 ......... ( 1 0 1 0 1 )
6 ........ ( 1 0 1 1 0 1 )
7 ....... ( 1 0 1 1 1 0 1 )
8 ...... ( 1 0 1 1 1 1 0 1 )
9 ..... ( 1 0 1 1 1 1 1 0 1 )
10 ... ( 1 0 1 1 1 1 1 1 0 1 )

Index entries for linear recurrences with constant coefficients, signature (11,-10).

Cf. A000533, A016789, A135577, A138120, A138144, A138145, A138146, A138721, A138826, A147758.

Cf. A144564 (a bisection).

f[n_] := Block[{w = {1, 0}}, Which[n == 1, w = {1}, n == 2, w = {1, 1}, n == 3, AppendTo[w, 1], n >= 4, w = Join[w, Table[1, {n - 4}], Reverse@ w]]; FromDigits@ w]; Array[f, 19] (* Michael De Vlieger, Dec 05 2015 *)
LinearRecurrence[{11,-10},{1,11,101,1001,10101},20] (* Harvey P. Dale, Aug 02 2017 *)

Vec( x+11*x^2+101*x^3 -91*x^4*(-11+10*x) / ( (10*x-1)*(x-1) ) + O(x^30)) \\ Michel Marcus, Dec 05 2015

G.f.: x+11*x^2+101*x^3-91*x^4*(-11+10*x) / ( (10*x-1)*(x-1) ). - R. J. Mathar, Aug 24 2011

a(n) = 11*a(n-1) - 10*a(n-2) for n>2. Wesley Ivan Hurt, Dec 06 2015

A152756 Bisection of A000533.

Original entry on oeis.org

1, 101, 10001, 1000001, 100000001, 10000000001, 1000000000001, 100000000000001, 10000000000000001, 1000000000000000001, 100000000000000000001, 10000000000000000000001

Offset: 1

Omar E. Pol, Dec 13 2008

a(1)=1, for n>1, a(n) is the concatenation of "1", 2(n-1)-1 digits "0" and "1". - Omar E. Pol, Dec 14 2008

			n ..... a(n)
1 ....... 1
2 ...... 101
3 ..... 10001
4 .... 1000001
5 ... 100000001

Index entries for linear recurrences with constant coefficients, signature (101,-100).

Cf. A000533, A135577, A138120, A138146.

[1] cat [10^(2*n)+1: n in [1..15]]; // Vincenzo Librandi, Jul 27 2014

A152756[n_] := If[n == 1, 1, 100^(n-1) + 1]; Array[A152756, 20] (* or *)
LinearRecurrence[{101, -100}, {1, 101, 10001}, 20] (* Paolo Xausa, Oct 05 2024 *)

Except the first term, a(n)=10^(2n-2)+1. - Robert G. Wilson v, Dec 14 2008

G.f.: x*(1-100*x^2)/(1-x)/(1-100*x). - Robert Israel, Jul 27 2014

A147816 Concatenation of n digits 1 and 2(n-1) digits 0.

Original entry on oeis.org

1, 1100, 1110000, 1111000000, 1111100000000, 1111110000000000, 1111111000000000000, 1111111100000000000000, 1111111110000000000000000, 1111111111000000000000000000, 1111111111100000000000000000000, 1111111111110000000000000000000000

Offset: 1

Omar E. Pol, Nov 13 2008

a(n) is also A016152(n) written in base 2.

			n ...... a(n)
1 ....... 1
2 ...... 1100
3 ..... 1110000
4 .... 1111000000
5 ... 1111100000000

Paolo Xausa, Table of n, a(n) for n = 1..300
Index entries for linear recurrences with constant coefficients, signature (1100,-100000).

Cf. A000533, A016152, A135577, A138119, A138120, A138144, A138145, A138146, A138721, A138826, A147757, A147759.

Array[(10^#-1)*10^(2*#-2)/9 &, 20] (* or *)
LinearRecurrence[{1100, -100000}, {1, 1100}, 20] (* Paolo Xausa, Feb 27 2024 *)

Vec(x/((100*x-1)*(1000*x-1))  + O(x^100)) \\ Colin Barker, Sep 16 2013

a(n) = A138119(n)/10.

a(n) = 1100*a(n-1)-100000*a(n-2). G.f.: x / ((100*x-1)*(1000*x-1)). - Colin Barker, Sep 16 2013

cp's OEIS Frontend

A001576 a(n) = 1^n + 2^n + 4^n.

Views

Author

Keywords

Comments

Links

Crossrefs

Programs

Mathematica

PARI

Sage

Formula

A333237 Numbers k such that 1/k contains at least one '9' in its decimal expansion.

Views

Author

Keywords

Comments

Examples

Links

Crossrefs

Programs

Maple

Mathematica

Python

Formula

Extensions

A138721 Concatenation of n digits 1, n digits 0 and n digits 1.

Views

Author

Keywords

Comments

Examples

Links

Crossrefs

Programs

Maple

Mathematica

PARI

Formula

A333402 Numbers m such that the largest digit in the decimal expansion of 1/m is 1.

Views

Author

Keywords

Comments

Examples

Links

Crossrefs

Programs

Mathematica

Python

Formula

Extensions

A066138 a(n) = 10^(2*n) + 10^n + 1.

Views

Author

Keywords

Comments

Examples

Links

Crossrefs

Programs

Magma

Mathematica

PARI

PARI

Formula

Extensions

A135576 Numbers whose binary expansion has only the digit "1" as first, central and final digit.

Views

Author

Keywords

Comments

Examples

Links

Crossrefs

Programs

Mathematica

PARI

Formula