Approximate the plastic number

Question

The challenge

The plastic number is a number related to the golden ratio, with many interesting mathematical properties. As such, there are many approaches that can be used to calculate the number.

In order to precisely specify the number for the purposes of this challenge, we'll use the following definition (although there are plenty of equivalent definitions, and you can use any definition you wish as long as it comes to the same number):

The plastic number is a real number ρ such that ρ³=ρ+1.

Your challenge is to write a program or function which takes an integer x as input (with x > 1), and produces an approximation to ρ as output, such that the larger the value of x gets, the closer the output gets to ρ (with at most finitely many exceptions; staying at the same value counts as "closer" for this purpose), and for any positive number δ, there's some input x to your program that produces an output that's within δ of ρ.

Clarifications

• If you're outputting via a method that inherently outputs strings (e.g. the standard output stream), you can format output either in decimal (e.g. 1.3247179572), or as a ratio of two integers with a / character between them.
• If you're outputting as a value within your programming language (e.g. returning from a function), it must be of a fixed-point, floating-point, or rational type. (In particular, you can't use data types that store numbers symbolically, unless they're used only to hold the ratio of two integers. So if you're using Mathematica or a similar language, you'll need to include the extra code to actually generate the digits of the output.)
• Your answer must work in a hypothetical variant of your language in which integers can be arbitrarily large, and memory (including stack) is unlimited. You may not assume that floating-point arithmetic in your language is arbitrarily accurate, but must instead use its actual accuracy (meaning that outputting a floating-point number is only going to be possible in languages where the accuracy of floating-point numbers can be controlled at runtime).
• x can have any meaning you want (so long as increasing it gives more accurate outputs). I imagine that most submissions will have it control the number of digits of output to produce, or the number of iterations of the algorithm used by your program to converge on the plastic number, but other meanings are acceptable.

Testcase

Here are the first few digits of the plastic number:

1.32471795724474602596090885

More digits are available on OEIS.

Victory condition

As usual for code-golf, shorter is better, measured in bytes. However, feel free to post answers even if they don't win, so long as they add something (e.g. a different language, or a different algorithm) to the existing answers.

Answer 1

Python 2, 49 bytes

n=x=input()
while n**3/x/x<n+x:n+=1
print n,'/',x

Try it online!

The idea is to express the ρ with ρ³=ρ+1 as a fraction n/x whose denominator x is the input accuracy parameter. We take (n/x)³=n/x+1 and clear denominators to get n³=x²(x+n).

Since the LHS increases in n faster than the RHS, we can approximate the equality point n as the smallest with n³≥x²(x+n). The code counts up n until this is the case, starting at x which is smaller.

A small byte save is to divide both sides by x² to write n³/x²≥x+n (negated in the while condition). This is floor division in the code, but the fractional part lost is negligible.

A same-length alternative instead puts x as the numerator:

Python 2, 49 bytes

n=x=input()
while x**3/n/n<n+x:n-=1
print x,'/',n

Try it online!

Answer 2

Mathematica, 20 bytes

#^3-#-1&~Root~1~N~#&

Mathematica's builtin Root function gives the solutions to a polynomial equation f[x] == 0.

Explanation

#^3-#-1&~Root~1~N~#&
                   &  (* Function *)
#^3-#-1&              (* A pure-function polynomial, x^3-x-1 *)
        ~Root~1       (* Find the first root *)
               ~N~#   (* approximate to (input) digits *)

Sample I/O

In[1]:= f=#^3-#-1&~Root~1~N~#&;
        f[1]

Out[1]= 1.

In[2]:= f[9]

Out[2]= 1.32471796

In[3]:= f[100]

Out[3]= 1.324717957244746025960908854478097340734404056901733364534015050302827851245547594054699347981787280

Answer 3

Mathematica, 27 bytes

x/.Solve[x^3==x+1>2,x]~N~#&

-1 byte from Martin
-2 bytes from ovs

input

[27]

output

{1.32471795724474602596090885}

Answer 4

sed, 67 60 (59+1) bytes

s,^,1/1/1 ,
:;s,(1*/(1*)/(1*).*)1$,\2\3/\1,
t
s,(/1*).*,\1,

Try it online!

+1 for the -E flag (ERE instead of BRE). Input and output are both unary: input 11111 for x=5 e.g. Output is a fraction of two unary numbers: the aforementioned 11111 input yields output 11111/1111 (5/4 in decimal).

Approximates the plastic number as a fraction between to consecutive elements of the Padovan sequence.

Answer 5

Mathematica, 27 bytes

Nest[(1+#)^(1/3)&,1,#]~N~#&

Uses a truncated approximation of the nested cubic radical form ³√(1+³√(1+³√(1+...))). While the output will always have x-1 decimal places, the result is actually less accurate than that, because the expression converges more slowly than one digit per iteration (x is also used as the number of nested radicals that are computed). For example x = 100 gives

_________________________________________________________________________
1.324717957244746025960908854478097340734404056901733364534015050302827850993693624204577670741656151

where the overlined part is correct.

Answer 6

Octave, 50 bytes

@(n)char(digits(n)*0+vpasolve(sym('r^3-r-1'))(1));

Try it online!

Defines an anonymous function, withn the desired number of digits of output.

This answer abuses that digits returns the current setting for the number of digits in variable precision arithmetic. This means we can just use it in an anonymous function without errors about 'Too many output arguments'.

Other than that, it's really straightforward: vpasolve is short for Variable-Precision Arithmetic Solve, with the precision set by the last call of digits. Since vpa is a Symbolic data type in Octave, which is banned per the spec, we just wrap the whole function in char(...) to get string output. Note that in solve and vpasolve, the f==0 is implied, so r^3==r+1 has been replaced by r^3-r-1 (==0)

Answer 7

MATL (27 28 bytes)

7BG:"t@)y@Q)+h]tG3+)yG2+)/

My first solution (27 bytes)

Try it online!

It's certainly not optimal, I'm still getting used to MATL.

Explanation:

I create a Padovan sequence up to input+3 then find the ratio of the last two numbers.

7B     % Turn 7 into binary to give 1 1 1 
G:"    % For k=1:input do...
t@)    % Existing sequence member k
y@1+)  % Existing sequence member k+1
+h     % Add them together and concatenate to the sequence array
]      % End loop
tG3+)  % Final sequence member
yG2+)  % Second last sequence member
/      % Divide to approximate ρ

Proper fraction output (35 bytes) (28 bytes, @Sanchises):

However, the first solution doesn't fulfill the need for arbitrary precision being the floating point limit of default MATL settings. So rather than adding several bytes to extend this precision, it's simpler to take the proper fraction route and write a fraction of the final two integers in the (N-1)^th and N^th elements of the truncated Padovan sequence.

e.g "114/86"

7BG:"t@)y@1+)+h]tG3+)V'/'YcyG2+)VYc

7BG:"t@tQh)sh]tJ)V47hyJq)Vh&

Courtesy of user @Sanchises. :)

Try it online!

Non-iterative evaluation:

Notably, my shortest code for the 'exact' version is (23 bytes):

1-1h69X^12**108+1I/^6/s

Try it online!

...but doesn't give arbitrary precision. I wonder if anyone can adjust this to fulfill the rules (use the input etc) and still add less than 5 bytes? :P

Answer 8

M, 15 14 bytes

²×3’
*3Ḥ‘÷Ç
Ç¡

Try it online!

Algorithm

This uses rationals and Newton's method. Specifically, for input x, the first x iterations with starting value x are applied.

We're trying to find a specific root of the polynomial p(t) = t³ - t - 1. Newton's method achieves this by taking a starting value t₀ – sufficiently close to ρ – and recursively defining a sequence by
t_n+1 = t_n - p(t_n) / p'(t_n).

Since p'(t) = 3t² -1, we get
t_n+1 = t_n - (t_n³ - t_n - 1)/(3t_n² - 1) = (3t_n³ - t_n - t_n³ + t_n + 1) / (3t_n² - 1) = (2t_n³ + 1) / (3t_n² - 1).

Note that the initial approximation x gets progressively worse as x increases. While the output for x = 3 is slightly less precise than the output for x = 2, since Newton's method converges quadratically to ρ, this shouldn't be an issue for large values of x.

How it works

Ç¡    Main link. Argument: x

Ç¡    Call the second helper link x times, which initial argument x.


*3Ḥ‘÷Ç  Second helper link. Argument: t

*3      Compute t³.
  Ḥ     Unhalve; yield 2t³.
   ‘    Increment; yield 2t³+1.
     Ç  Call the first helper link with argument t.
    ÷   Divide the left result by the right one.


²×3’    First helper link. Argument: t

²       Compute t².
 ×3     Compute 3t².
   ’    Decrement; yield 3t²-1.

Answer 9

Julia 0.5,  44  40 bytes

|(x,r=1)=x<1?r:~-x|big(2r^3+1)//(3r^2-1)

Uses rationals and Newton's method.

Try it online!

Answer 10

APL (Dyalog Extended), 19 bytes

{{(⍺+⍵*÷3)*÷3}/⍵/1}

Try it online!

Uses the cube root approximation method, which very, very slowly creeps towards the value of the plastic number.

Explanation

{{(⍺+⍵*÷3)*÷3}/⍵/1} ⍵ → input
               ⍵/1} replicate 1 ⍵ times
 {           }/     reduce by inner fn:
          *÷3       cube root of
     ⍵*÷3           cube root of right arg
   ⍺+               plus left arg

Answer 11

Rust, 66 bytes

|n|{let mut i=(1,1,1);for _ in 0..n{i=(i.1,i.2,i.0+i.1)}(i.1,i.2)}

Try it online!

Uses the padovan sequence to approximate the plastic number. Returns a tupled (num,denom). A readable version that works with arbitrary precision integers can be found here.

Answer 12

AWK, 34 bytes

{for(;$1^3>$1+1;)$1=($1+1)^(1/3)}1

Step by step:

{                  # $1 is the input (a positive integer)
for(;$1^3>$1+1;)   # as long as $1 cubed is bigger than the input plus one
  $1=($1+1)^(1/3)  # assigns $1 this new value: cube root of $1 plus one
}
1                  # prints $1

Try it online!

Answer 13

05AB1E, 23 bytes

U[X3m¹/¹/¹X+›#X>U]X'/¹J

Try it online!

---

Direct port of https://codegolf.stackexchange.com/a/126822/59376 by xnor.

Answer 14

Charcoal, 28 bytes

ＡＩθθＡθνＷ‹∕∕Ｘν³θθ⁺νθＡ⁺ν¹νＩ∕νθ

Try it online! Link to verbose mode. Also I apparently messed up Divide and IntDivide :|
Uses the same method as the Python and JavaScript answers.

Answer 15

NewStack, 14 bytes

¹Fᵢ{E2x³⁺÷3x²⁻

Break down:

¹                Add arbitrary number 1 to the stack.
 Fᵢ{             Define for loop with a user's input amount of itterations.
    E            Define new edit for element 0 (element 0 being the 1 added. earlier).
     2x³⁺÷3x²⁻   update x to equal (2x^3+1)/(3x^2-1). (x = element 0).

How it works:

The formula (2x³+1)/(3x²-1) comes from the simplification of Newton's method for the equasion x³=x+1. You can find it here. Repeating this process an infinite amoune of times converges to the plastic number. It's rate of convergence is rather quick at around 2.6 decimals per iteration.

INPUT ITERATION >> VALUE
0 >> 1
1 >> 1.5
2 >> 1.3478260869565217
3 >> 1.325200398950907
4 >> 1.3247181739990537
5 >> 1.3247179572447898
6 >> 1.324717957244746    <- 16 decimal precision in 6 iterations!
...
100 >> 1.324717957244746

Padovan sequence alternative, 27 25 17 bytes

¹Fᵢ{[ƨ2+ƨ3]ℲƤƨ/ƨ2

Break down:

¹                  Append first element of Padovan sequence.
 Fᵢ{       Ⅎ       Define for loop of user's input amount of iterations.
    [ƨ2+ƨ3]        Append sum second and third to last elements.
            Ƥƨ/ƨ2  Print ratio of last two elements.

-2 bytes by choosing better print strategy

-8 bytes by choosing better way to index stack

How it works:

As the Padovan sequence continues, the ratio of the last two elements converge to the plastic number.

INPUT ITERATION >> VALUE
0 >> 1
1 >> 2
...
10 >> 1.3157894736842106
...
89 >> 1.324717957244746    <- 16 decimal precision in 89 iterations
...
100> > 1.324717957244746

Answer 16

Ruby, 35 bytes

->x{n=1;n+=0.1**x while n**3-n<1;n}

Try it online!

Answer 17

Clojure, 46 bytes

#(nth(iterate(fn[i](Math/pow(inc i)(/ 3)))1)%)

Uses the iterated cube-root formula. This is a bit more interesting but longer:

(def f #(apply comp(repeat %(fn[i](Math/pow(inc i)(/ 3))))))

((f 10)1)
1.3247179361449652

Answer 18

Javascript, 36 bytes

f=(x,n=x)=>n**3/x/x<n+x?f(x,++n):n/x

Works the same as the top python answer. No console.log was included because if you run f(x) in console it will be logged automatically.

f=(x,n=x)=>n**3/x/x<n+x?f(x,++n):n/x
console.log(f(300))

Answer 19

><>, 38+3 = 41 bytes

11\n;
?!\}2,:01{::::}**-+0(?$-{+{1-:}

Expects the input to be present on the stack at program start, so +3 bytes for the -v flag.

Try it online!

Effectively performs a binary search to narrow in on the output value. Increasing x increases the number of iterations to perform.

Edit: refactored calculation slightly to save 1 byte, previous version:

11\n;
?!\}2,:0{::::}**$-1-0)?$-{+{1-:}

Answer 20

k, 27 bytes

{"/"/$1_|x({y,z,x+y}.)/3#1}

Try it online! This assumes infinite integers (which, alas, is not true). It uses the Padovan sequence.

Answer 21

TI-BASIC, 21 bytes

:Prompt X //Prompt for input, 3 bytes
:While X  //While X, 3 bytes
:³√(1+Y→Y //Calculate cube root of 1+Y and store to Y, 7 bytes
:DS<(X,0  //Decrement X and skip next command (which doesn't do anything), 5 bytes
:End      //End While loop, 2 bytes
:Y        //Display Y, 1 byte

Uses this recursive formula.

Interestingly, hard-coding the number and rounding it gives the same byte-count:

TI-BASIC, 21 bytes

:Prompt X    //Prompt for input, 3 bytes
:.5√(3       //Store √(3)/2 to Ans, 5 bytes
:Ansֿ¹cosh(3ֿ¹coshֿ¹(3Ans //Store the plastic number to Ans, 9 bytes
:round(Ans,X //Round the plastic number X decimal digits, 4 bytes

Uses this trigonometric formula.

Answer 22

C#, 317 bytes

using m=System.Math;a=x=>{if(a==0)return "1/1";var d=a(x-1).Split('/');var b=int.Parse(d[0]);var c=int.Parse(d[1]);return string.Format("{0}/{1}",(2*m.Pow(b,3)+m.Pow(c,3)).ToString(new string('#',int.MaxValue.ToString().Length)),(3*m.Pow(b,2)*c-m.Pow(c,3)).ToString(new string('#',int.MaxValue.ToString().Length)));};

It returns the result as a fraction.

Explanation

It uses the Newton's method with x iterations for finding the root of the polynomial p^3-p-1=0. The formula is x_n=1-(f(x_(n-1)))/(f'(x_(n-1))), and x_0 is a starting point.

The polynomials derivative is 3p^2-1, and let's say x_(n-1)=b/c. Then, by using the above formula we get, that x_n=(2 b^3+c^3)/(3 b^2 c-c^3). Let's also say, that we start from 1, this will happen, when x=2, because x>1, and is an integer. Idented, and commented code:

using System;
string PlasticNumber(int x)
{
    if (x == 2) 
        return "1/1";                 

//If x=2, we return our starting value, but we need to return it as a fraction

    var d = PlasticNumber(x - 1).Split('/');
    var b = System.Convert.ToInt32(d[0]);
    var c = int.Parse(d[1]);

//We parse the previous value of the fraction, and put it into two variables

    return string.Format("{0}/{1}", 
        (2 * Math.Pow(b, 3) + Math.Pow(c, 3))
        .ToString(new string('#', int.MaxValue.ToString().Length)),
        (3 * Math.Pow(b, 2) * c - Math.Pow(c, 3))
        .ToString(new string('#', int.MaxValue.ToString().Length)));

//We return the result as a fraction, but it's important not to return it in
  scientific notation, because that will cause issues in the parsing process 

}

Answer 23

PHP, 86 bytes

for($p=[1,1,1];$i++<$argn;)$p[]=bcadd($p[$i],$p[$i-1]);echo bcdiv($p[$i+1],$p[$i],$i);

PHP Sandbox Online

Creates the Padovan Spiral and print the ratio of the last two numbers.

Answer 24

Axiom, 96 bytes

h(n:NNI):Float==(n>1.E5=>-1;n:=n+1;j:=digits(n::PI);r:=solve(x^3-x=1,10.^-n);digits(j);rhs(r.1))

results

(31) -> [h(i) for i in 0..10]
   (31)
   [1.0, 1.3, 1.33, 1.325, 1.3247, 1.32472, 1.324718, 1.324718, 1.32471796,
    1.324717957, 1.3247179572]
                                                         Type: List Float

how you can see h(2) should be 1.32 and not 1.33 so there is some error in last digits

Then there would be this one of 110 bytes

g(n:NNI):Float==(n>1.E5=>-1;n:=n+1;j:=digits(n::PI);x:=sqrt(23./108);r:=(.5+x)^(1/3)+(.5-x)^(1/3);digits(j);r)

It use the formula for resolve equation of III grade of type x^3-3*p*x-2*q=0 in the case q^2-p^3>=0 that is m=sqrt(q^2-p^3) and x=(q+m)^(1/3)+(q-m)^(1/3)

In our case r^3-r-1=0 this can be written as r^3-3*(1/3)r-2*(1/2)=0 so p=1/3 q=1/2 m=1/4-1/27=23/108 x=(0.5+m)^(1/3)+(0.5-m)^(1/3)

this one that use Newton iteration with start point r=1

f(n:NNI):Float==(n>1.E5=>-1;n:=n+1;j:=digits(n::PI);e:=10^-n;r:=1.;repeat(v:=(r^3-r-1)/(3*r^2-1);abs(v)<e=>break;r:=r-v);digits(j);r)

it change in the function, digits value for obtain one obj of n+1 digits afther the float point. At end the digits() value is assigned back to preciding value.