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

23

2

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

(3 + sqrt(5))^n = a + b * sqrt(5)

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

orlp

Posted 2015-04-12T20:35:38.987

Reputation: 37 067

Answers

3

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 .

lirtosiast

Posted 2015-04-12T20:35:38.987

Reputation: 20 331

16

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.

pawel.boczarski

Posted 2015-04-12T20:35:38.987

Reputation: 1 243

You're selling yourself short, [3 5;1 3]**input('')*[1;0] is 26 bytes, not 41. – orlp – 2015-04-12T21:49:25.873

3@(n)[3 5;1 3]^n*[1;0] (function handle) would save you five characters, mut nice idea! – flawr – 2015-04-12T22:36:22.087

14

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]].

xnor

Posted 2015-04-12T20:35:38.987

Reputation: 115 687

12

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

Initial expand

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:

Matrix

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]

Alex A.

Posted 2015-04-12T20:35:38.987

Reputation: 23 761

8

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

randomra

Posted 2015-04-12T20:35:38.987

Reputation: 19 909

You can get down to 20 by exploiting tacit adverb parsing: +/ .*(3 5,:1 3&)&1 0. – algorithmshark – 2015-04-13T00:11:14.783

@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 – 2015-04-13T00:42:50.487

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 – 2015-04-13T00:50:21.417

7

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.

isaacg

Posted 2015-04-12T20:35:38.987

Reputation: 39 268

I chose @randomra's answer over yours because his/hers was posted 16 minutes earlier, sorry. – orlp – 2015-05-07T00:41:17.117

5

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.

marinus

Posted 2015-04-12T20:35:38.987

Reputation: 30 224

1Awww yiss, APL. – Nit – 2015-04-13T07:38:09.540

5

Jelly, 13 bytes

5W×U++Ḥ
2Bdz¡

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].

2Bdz¡      Main link. Argument: n

2B         Convert 2 to binary, yielding [1, 0].
    ¡      Repeat:
  Ç            Apply the helper link...
   ³           n times.

Dennis

Posted 2015-04-12T20:35:38.987

Reputation: 196 637

No, I'm pretty sure Jelly was around a long time before the creation of the internet :P – Conor O'Brien – 2016-01-23T04:04:13.060

1@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 – 2016-01-23T04:33:32.973

3

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]}

flawr

Posted 2015-04-12T20:35:38.987

Reputation: 40 560

1Might want to consider editing your formula. We don't have MathJax anymore. – Alex A. – 2015-04-12T21:35:01.733

I thought it was just introduced a few days ago? – flawr – 2015-04-12T22:39:17.110

Yeah, but it messed up the stack snippets, so it had to be disabled. – Alex A. – 2015-04-12T23:30:38.367

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 – 2015-04-13T06:45:13.057

n=>[...Array(n)].map(_=>[x,y]=[3*x+5*y,x+3*y],y=0,x=1)[n-1] same length – l4m2 – 2017-12-27T00:36:41.717

F=(n,x=1,y=0)=>n?F(n-1,3*x+5*y,x+3*y):[x,y] lot shorter – l4m2 – 2017-12-27T00:40:57.293

3

Mathematica, 31

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

alephalpha

Posted 2015-04-12T20:35:38.987

Reputation: 23 988

3

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).        ";

Dennis

Posted 2015-04-12T20:35:38.987

Reputation: 196 637

2

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!

BrainSteel

Posted 2015-04-12T20:35:38.987

Reputation: 5 132

1This is uses a rather overcomplicated algorithm :P – orlp – 2015-04-12T22:10:52.167

@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 – 2015-04-12T22:16:26.373

1Upvote because this is awesomely horrific. – kirbyfan64sos – 2015-04-12T22:43:14.430

2

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.

algorithmshark

Posted 2015-04-12T20:35:38.987

Reputation: 8 144

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 – 2015-04-14T18:02:15.820

FYI, this was recently fixed in Kona.

– kirbyfan64sos – 2015-04-22T19:28:33.397

2

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].

orion

Posted 2015-04-12T20:35:38.987

Reputation: 3 095

2

Seriously, 32 bytes, non-competing

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

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}.)

quintopia

Posted 2015-04-12T20:35:38.987

Reputation: 3 899

1

Haskell, 41 bytes

(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).

nimi

Posted 2015-04-12T20:35:38.987

Reputation: 34 639

1

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

philn

Posted 2015-04-12T20:35:38.987

Reputation: 111

Welcome to the site, and nice first answer! – James – 2016-04-28T20:33:25.777

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.

kirbyfan64sos

Posted 2015-04-12T20:35:38.987

Reputation: 8 730

0

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.

user16488

Posted 2015-04-12T20:35:38.987

Reputation:

1This 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 – 2015-04-12T23:37:05.487

0

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))))))

Penguino

Posted 2015-04-12T20:35:38.987

Reputation: 231

0

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.

Alchymist

Posted 2015-04-12T20:35:38.987

Reputation: 544

Your solution overflows integers for large inputs :) – orlp – 2015-04-13T16:10:32.757

@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 – 2015-04-14T08:20:16.523

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 – 2015-04-14T09:08:48.713

Suggest a[*a=1]=0 instead of *a=1,a[1]=0 – ceilingcat – 2019-02-23T03:37:47.870

0

(non-competing) Jelly, 10 bytes

31,53Dæ*³Ḣ

Try it online!

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

Leaky Nun

Posted 2015-04-12T20:35:38.987

Reputation: 45 011