# Calculating (3 + sqrt(5))^n exactly

Today your goal is to find integers a and b given non-negative integer n such that:

You should write a program or a function that takes parameter n and outputs a and b in a format of your choice.

Standard loopholes apply. Additionally, it's intended that you implement the above problem using basic arithmetic yourself. So you may not use built-in exact algebra functionality, rationals, or functions implementing non-trivial mathematical constructs (for example the Lucas sequence).

Shortest code in bytes wins.

Example input/output:

0 → 1, 0
1 → 3, 1
2 → 14, 6
3 → 72, 32
4 → 376, 168
5 → 1968, 880
6 → 10304, 4608
7 → 53952, 24128
8 → 282496, 126336
9 → 1479168, 661504

# Octave, 26 bytes

[3 5;1 3]**input('')*[1;0]


Because (a+b*sqrt(5)) * (3+sqrt(5)) = (3a+5b) + (a+3b) * sqrt(5),

multiplying input vector

| 1 |    /* a = 1 */
| 0 |    /* b = 0 */


which stands for 1 = (3+sqrt(5))^0 by matrix

| 3 5 |
| 1 3 |


seems natural. Instead of looping n times, we rather raise the matrix to the power of n and then multiply it by input vector.

• You're selling yourself short, [3 5;1 3]**input('')*[1;0] is 26 bytes, not 41.
– orlp
Apr 12, 2015 at 21:49
• @(n)[3 5;1 3]^n*[1;0] (function handle) would save you five characters, mut nice idea! Apr 12, 2015 at 22:36

# Python 2, 50

a=1;b=0
exec"a,b=3*a+5*b,3*b+a;"*input()
print a,b


Multiplies by 3+sqrt(5) repeatedly by its action on the pair (a,b) representing a+b*sqrt(5). Equivalent to starting with the column vector [1,0] and left-multiplying n times by the matrix [[3,5],[1,3]].

# Julia, 22 20 bytes

n->[3 5;1 3]^n*[1;0]


This creates a lambda function which takes a single integer as input and returns a 2-element vector of integers corresponding to the solution [a, b]. To call it, give it a name, e.g. f=n->....

Start by multiplying

We can then translate the right hand side of this equation into a 2-column matrix, where the first corresponds to the coefficient of a and the second to the coefficient of b:

Multiply this matrix by itself n times, then right multiply by the column vector (1, 0), and POOF! Out pops the solution vector.

Examples:

julia> println(f(0))
[1,0]

julia> println(f(5))
[1968,880]


# J, 20 bytes

+/@:*(3 5,.1 3&)&1 0


Multiplicate the vector [1 0] with the matrix [[3 5] [1 3]] n times.

2 bytes saved thanks to @algorithmshark.

Usage and test:

   (+/@:*(3 5,.1 3&)&1 0) 5
1968 880

(+/@:*(3 5,.1 3&)&1 0) every i.6
1   0
3   1
14   6
72  32
376 168
1968 880

• You can get down to 20 by exploiting tacit adverb parsing: +/ .*(3 5,:1 3&)&1 0. Apr 13, 2015 at 0:11
• @algorithmshark Thanks, although why (+/@:*&(3 5,.1 3)&1 0) works and (+/@:*&1 0&(3 5,.1 3)) not? Shouldn't the second one bond correctly and the first one swapped? Apr 13, 2015 at 0:42
• Got it, they bond as I expected but the outer & makes the powering/looping so you modify the left side input during powering (opposite to the normal right-side modification). Apr 13, 2015 at 0:50

# Pyth, 20 bytes

u,+*3sGyeG+sGyeGQ,1Z


u which is reduce in general, is used here as an apply repeatedly loop. The updating function is G -> ,+*3sGyeG+sGyeG, where G is a 2 tuple. That function translates to 3*sum(G) + 2*G[1], sum(G) + 2*G[1]. s is sum, y is *2.

• I chose @randomra's answer over yours because his/hers was posted 16 minutes earlier, sorry.
– orlp
May 7, 2015 at 0:41

# APL (22)

{⍵+.×⍨2 2⍴3 5 1}⍣⎕⍨2↑1


Explanation:

• {...}⍣⎕⍨2↑1: read a number, and run the following function that many times, using [1,0] as the initial input.
• 2 2⍴3 5 1: the matrix [[3,5],[1,3]]
• ⍵+.×⍨: multiply the first number in ⍵ by 3, the second by 5, and sum them, this is the new first number; then multiply the first number in ⍵ by 1, the second by 3, and sum those, that is the new second number.
• Awww yiss, APL. Apr 13, 2015 at 7:38

# Jelly, 13 bytes

5W×U++Ḥ
2BÇ³¡


Try it online!

### How it works

5W×U++Ḥ    Helper link. Argument: [a, b]

5W         Yield [5].
×U       Multiply it by the reverse of [a, b]. This yields [5b, a].
+      Hook; add the argument to the result. This yields [a + 5b, a + b].
+Ḥ    Fork; add the doubled argument ([2a, 2b]) to the result.
This yields [3a + 5b, a + 3b].

2B         Convert 2 to binary, yielding [1, 0].
¡      Repeat:
³           n times.

• No, I'm pretty sure Jelly was around a long time before the creation of the internet :P Jan 23, 2016 at 4:04
• @Doᴡɴɢᴏᴀᴛ For non-competing answers, I prefer keeping the byte count on the second line. This keeps the answer from rising to the top in leaderboards and userscripts, which seems unfair. Jan 23, 2016 at 4:33

## Dyalog APL, 18 bytes

((3∘×+5 1×⌽)⍣⎕)1 0


This is a program that takes input through ⎕.

 (         )         Monadic train:
3∘×                3 times argument
+               Plus
5 1×⌽          (5 1) times the reverse
(           ⍣⎕)      Apply that function (input) times
1 0   starting with (1 0)


The features used here were implemented well before April 2015, making this answer valid.

Try it here. Note that tryapl.org is a limited subset of Dyalog and does not support ⎕.

# Mathematica, 31

Nest[{{3,5},{1,3}}.#&,{1,0},#]&


# CJam, 21 bytes

0X{_2$3*+@5*@3*+}li*p  Try it online. ### How it works 0X " Stack: [ 0 1 ] "; li{ " Do int(input()) times: "; _2$    " Stack: [ a b ] -> [ a b b a ]                 ";
3*+    " Stack: [ a b b a ] -> [ a b (b+3a) ]          ";
@5*@3* " Stack: [ a b (b+3a) ] -> [ (b+3a) 5a 3b ]     ";
+      " Stack: [ (b+3a) 5a 3b ] -> [ (b+3a) (5a+3b) ] ";
}*       "                                               ";
p        " Print topmost stack item plus linefeed.       ";
" Print remaining stack item (implicit).        ";


# Javascript, 63 61 bytes

I am using a recursive evaluation of the binomial: (x+y)^n = (x+y)(x+y)^{n-1}

New(thanks to @edc65)

F=n=>{for(i=y=0,x=1;i++<n;)[x,y]=[3*x+5*y,x+3*y];return[x,y]}


Old

F=n=>{for(i=y=0,x=1;i<n;i++)[x,y]=[3*x+5*y,x+3*y];return [x,y]}

• Might want to consider editing your formula. We don't have MathJax anymore. Apr 12, 2015 at 21:35
• I thought it was just introduced a few days ago? Apr 12, 2015 at 22:39
• Yeah, but it messed up the stack snippets, so it had to be disabled. Apr 12, 2015 at 23:30
• I count 63 as is, and can be shortened to 61 F=n=>{for(i=y=0,x=1;i++<n;)[x,y]=[3*x+5*y,x+3*y];return[x,y]} Apr 13, 2015 at 6:45
• n=>[...Array(n)].map(_=>[x,y]=[3*x+5*y,x+3*y],y=0,x=1)[n-1] same length
– l4m2
Dec 27, 2017 at 0:36

# C, 114 bytes

g(n){int i,a[2]={1,0},b[2];for(i=0;i<n;i++)*b=*a*3+5*a[1],b[1]=*a+3*b[1],*a=*b,a[1]=b[1];printf("%d,%d",*a,a[1]);}


This implements matrix multiplication the boring way. For a more fun (quote: "awesomely horrific") 238 byte solution, look no further!

f(n){int p[2][n+3],i,j,k=0,a[2]={0};for(j=0;j<n+3;j++)p[0][j]=0;*p[1]=0;(*p)[1]=1;for(j=0;j<n;j++,k=!k)for(i=1;i<n+3;i++)p[!k][i]=p[k][i-1]+p[k][i];for(i=1;i<n+2;i++)a[!(i%2)]+=p[k][i]*pow(3,n+1-i)*pow(5,(i-1)/2);printf("%d,%d",*a,a[1]);}


Unraveled:

g(n){
int i,a[2]={1,0},b[2];
for(i=0;i<n;i++)
*b=3**a+5*a[1],b[1]=*a+3*b[1],*a=*b,a[1]=b[1];
printf("%d,%d",*a,a[1]);
}


This could probably be shortened a bit. Try a test program online!

• This is uses a rather overcomplicated algorithm :P
– orlp
Apr 12, 2015 at 22:10
• @orlp I couldn't think of a shorter algorithm for this language. I thought this one would work out, but it kind of got out of hand, haha. Implementing matrix multiplication by hand could very well be shorter. Apr 12, 2015 at 22:16
• Upvote because this is awesomely horrific. Apr 12, 2015 at 22:43
• 82 bytes Sep 22, 2020 at 1:26

# k2 - 22 char

Function taking one argument.

_mul[(3 5;1 3)]/[;1 0]


_mul is matrix multiplication so we curry it with the matrix (3 5;1 3) and then hit it with the functional power adverb: f/[n;x] applies f to x, n times. Again we curry it, this time with the starting vector 1 0.

  _mul[2 2#3 5 1]/[;1 0] 5
1968 880
f:_mul[2 2#3 5 1]/[;1 0]
f'!8  /each result from 0 to 7 inclusive
(1 0
3 1
14 6
72 32
376 168
1968 880
10304 4608
53952 24128)


This will not work in Kona, because for some reason f/[n;x] isn't implemented correctly. Only the n f/x syntax works, so the shortest fix is {x _mul[(3 5;1 3)]/1 0} at 23 char.

• Wow. This use of currying is so smart I feel like my K answer is stupid. Anyway, I raised the issue you found in Kona on their bug tracker. Apr 14, 2015 at 18:02
• Apr 22, 2015 at 19:28

# ised, 25 bytes (20 characters)

({:{2,4}·x±Σx:}$1)∘1  I hoped for better, but there are just too many braces needed in ised to make it competent, the operator precedence is not optimal for golfing. It expects the input to be in$1 memory slot, so this works:

ised '@1{9};' '({:{2,4}·x±Σx:}\$1)∘1'


For n=0, the zero is skipped (outputs 1, instead of 1 0). If that's an issue, replace the final 1 with ~[2].

## Seriously, 32 bytes, non-competing

,╗43/12╜";)@4*≈(6*-"£n.X4ì±0n


Hex Dump:

2cbb34332f313260bd223b2940342af728362a2d229c6e2e58348df130606e7f


Try It Onlline

Obviously not a contender for shortest, but at least the method is original. (Noting that such a problem necessarily indicates a Lucas sequence, as mentioned in the description, this program generates successive terms of the sequences using the recurrence relation

a_n = 6*a_{n-1} - 4*a_{n-2}.)

(iterate(\(a,b)->(3*a+5*b,a+3*b))(1,0)!!)


Usage example: (iterate(\(a,b)->(3*a+5*b,a+3*b))(1,0)!!) 8 -> (282496,126336).

C/C++ 89 bytes

void g(int n,long&a,long&b){if(n){long j,k;g(n-1,j,k);a=3*j+5*k;b=j+3*k;}else{a=1;b=0;}}


Formatted:

    void g(int n, long&a, long&b) {
if (n) {
long j, k;
g(n - 1, j, k);
a = 3 * j + 5 * k;
b = j + 3 * k;
} else {
a = 1;
b = 0;
}}


Same concept:

void get(int n, long &a, long& b) {
if (n == 0) {
a = 1;
b = 0;
return;
}
long j, k;
get(n - 1, j, k);
a = 3 * j + 5 * k;
b = j + 3 * k;
}


The test bench:

#include <iostream>
using namespace std;
int main() {
long a, b;
for (int i = 0; i < 55; i++) {
g(i, a, b);
cout << i << "-> " << a << ' ' << b << endl;
}
return 0;
}


The output:

0-> 1 0
1-> 3 1
2-> 14 6
3-> 72 32
4-> 376 168
5-> 1968 880
6-> 10304 4608
7-> 53952 24128
8-> 282496 126336
9-> 1479168 661504
10-> 7745024 3463680
11-> 40553472 18136064
12-> 212340736 94961664
13-> 1111830528 497225728
14-> 5821620224 2603507712
15-> 30482399232 13632143360
16-> 159607914496 71378829312
17-> 835717890048 373744402432
18-> 4375875682304 1956951097344
19-> 22912382533632 10246728974336
20-> 119970792472576 53652569456640
21-> 628175224700928 280928500842496
22-> 3289168178315264 1470960727228416
23-> 17222308171087872 7702050360000512
24-> 90177176313266176 40328459251089408
25-> 472173825195245568 211162554066534400
26-> 2472334245918408704 1105661487394848768

• Welcome to the site, and nice first answer! Apr 28, 2016 at 20:33

# Jelly, 10 bytes

31,53Dæ*³Ḣ


Try it online!

Uses matrix. Computes ([[3,1],[5,3]]**input())[0].

# K, 37 bytes

f:{:[x;*(1;0)*_mul/x#,2 2#3 1 5;1 0]}


or

f:{:[x;*(1;0)*_mul/x#,(3 1;5 3);1 0]}


They're both the same thing.

# Python 3, 49 bytes

w=5**0.5;a=(3+w)**int(input())//2+1;print(a,a//w)


although on my machine, this only gives the correct answer for inputs in the range 0 <= n <= 18.

This implements the closed form formula

w = 5 ** 0.5
u = 3 + w
v = 3 - w
a = (u ** n + v ** n) / 2
b = (u ** n - v ** n) / (2 * w)


and takes advantage of the fact that the v ** n part is small, and can be computed by rounding rather than direct calculation.

• This is not a valid solution (you must support any n), but since you're nowhere near being the shortest I don't see a reason to downvote. It's a cool solution.
– orlp
Apr 12, 2015 at 23:37

# Scheme, 97 bytes

(define(r n)(let s([n n][a 1][b 0])(if(= 0 n)(cons a b)(s(- n 1)(+(* a 3)(* b 5))(+ a(* b 3))))))


C 71 bytes (60 with pre-initialised variables)

Scope for golfing yet but just to prove that C does not have to be "awesomely horrific".

f(int n,int*a){for(*a=1,a[1]=0;n--;a[1]=*a+3*a[1],*a=(5*a[1]+4**a)/3);}


If the values in a are initialised to {1,0}, we do better.

f(int n,int*a){for(;n--;a[1]=*a+3*a[1],*a=(5*a[1]+4**a)/3);}


This is iteratively using the mappings a->3a+5b, b->a+3b but avoiding a temporary variable by calculating a from the new value of b instead.

• Your solution overflows integers for large inputs :)
– orlp
Apr 13, 2015 at 16:10
• @orlp - That's C for you. Granted this solution fails earlier than others because of the interim calculation in brackets but it would only manage a couple of extra steps anyway unless I change the datatype. Is it worth explicitly changing the question to give the range you expect to support? Probably too late now though. Apr 14, 2015 at 8:20
• There is no range to support, a proper solution should work for any input. In C that means you'll have to implement arbitrary width integers, sorry =/
– orlp
Apr 14, 2015 at 9:08
• Suggest a[*a=1]=0 instead of *a=1,a[1]=0 Feb 23, 2019 at 3:37