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

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 .

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 '15 at 21:49
• @(n)[3 5;1 3]^n*[1;0] (function handle) would save you five characters, mut nice idea! – flawr Apr 12 '15 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. – algorithmshark Apr 13 '15 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? – randomra Apr 13 '15 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). – randomra Apr 13 '15 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, sum(G) + 2*G. 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 '15 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. – Nit Apr 13 '15 at 7:38

Jelly, 13 bytes

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

Try it online!

How it works

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

5W         Yield .
×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Ç³¡      Main link. Argument: n

2B         Convert 2 to binary, yielding [1, 0].
¡      Repeat:
Ç            Apply the helper link...
³           n times.
• No, I'm pretty sure Jelly was around a long time before the creation of the internet :P – Conor O'Brien Jan 23 '16 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. – Dennis Jan 23 '16 at 4:33

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. – Alex A. Apr 12 '15 at 21:35
• I thought it was just introduced a few days ago? – flawr Apr 12 '15 at 22:39
• Yeah, but it messed up the stack snippets, so it had to be disabled. – Alex A. Apr 12 '15 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]} – edc65 Apr 13 '15 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 '17 at 0:36

C, 114 bytes

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

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

f(n){int p[n+3],i,j,k=0,a={0};for(j=0;j<n+3;j++)p[j]=0;*p=0;(*p)=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);}

Unraveled:

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

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

• This is uses a rather overcomplicated algorithm :P – orlp Apr 12 '15 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. – BrainSteel Apr 12 '15 at 22:16
• Upvote because this is awesomely horrific. – kirbyfan64sos Apr 12 '15 at 22:43

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. – kirbyfan64sos Apr 14 '15 at 18:02
• – kirbyfan64sos Apr 22 '15 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 ~.

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! – DJMcMayhem Apr 28 '16 at 20:33

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 '15 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=0;n--;a=*a+3*a,*a=(5*a+4**a)/3);}

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

f(int n,int*a){for(;n--;a=*a+3*a,*a=(5*a+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 '15 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. – Alchymist Apr 14 '15 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 '15 at 9:08
• Suggest a[*a=1]=0 instead of *a=1,a=0 – ceilingcat Feb 23 at 3:37

(non-competing) Jelly, 10 bytes

31,53Dæ*³Ḣ

Try it online!

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