Write a short program on your favorite RPN calculator to display the Fibonacci sequence:

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, ...

No manual initialization allowed. On the HP-15C LE you might want to use R/S instead of PSE. Currently I have a couple of 7-step HP-12C programs and one 6-step WP 34S program, but this can always be improved (except by creating a new firmware just for this task - FIB can of course be used if needed :-).

Happy programming!

Two completely untest 34S solutions both six steps.

01   iC 0        Clx

02   iC 1        FIB

03   x[<->] Y    PSE 10

04   RCL+ Y      RCL L

05   PSE 10      INC X

06   BACK 003    BACK 004

The left version is more suitable to other RPN calculators of course. Both are untested and likely broken but the ideas are probably okay.

- Pauli

It works. It's similar to the one I came up with when using FIB. I forgot to mention the step count should include LBL when applicable, sorry!

001 LBL 99

002 0

003 FIB

004 PSE 10

005 INC L

006 RCL L

007 BACK 04

Gerson.

Edited: 7 Jan 2012, 9:20 p.m.

For 42S:

00 { 10-Byte Prgm }

01>LBL 01

02 CLST

03 SIGN

04>LBL 02

05 STO+ ST Y

06 X<>Y

07 PSE

08 GTO 02

09 .END.

Edited: 7 Jan 2012, 9:46 p.m.

10 bytes! One step shorter here... but 1 byte longer.

This is supposed to be just an exercise, not a challenge. Feel free to submit your solutions. Different approaches are always interesting, regardless of code size or register usage.

Here are my solutions:

WP-34S

001 LBL 88

002 1

003 LN

004 STOP

005 RCL+ L

006 BACK 02

HP-12C

01-  0

02-  n!

03-  LASTx

04-  PSE

05-  x<>y

06-  +

07-  GTO 03

01-  1

02-  LN

03-  PSE

04-  LASTx

05-  x<>y

06-  +

07-  GTO 03

Here are a few for the HP41C:

01 LBL 00

02 CLST

03 E

04 ENTER

05 LBL 01

06 X<>Z

07 ST+ Z

08 VIEW Z

09 PSE

10 GTO 01

11 END

16 Bytes

01 LBl 00

02 .

03 E

04 LBL 01

05 X<>Y

06 ST+ Y

07 VIEW Y

08 PSE

09 GTO 01

10 END

16 Bytes

01 LBL 00

02 E

03 .

04 LBL 01

05 +

06 PSE

07 LASTX

08 X<>Y

09 GTO 01

10 END

14 Bytes

I used local labels and synthetic instructions to make the programs shorter.

. = 0

E = 1

Cheers,

Fouad

Thanks for your participation!

I'd completely forgotten about this old thread. Not exactly the same problem, but some interesting ideas there, like Kioshi Ashima's as implemented on the WP 34S:

001 LBL 11

002 FIX 00

003 # Phi

004 FILL

005 2

006 +

007 /

008 PSE 10

009 *

010 BACK 02

This is based on the formula F(n) = round(phi^n/sqrt(5)) (See https://oeis.org/A000045). The first two Fibonacci numbers are not

properly rounded by FIX 0 though.

Gerson.

Edited: 8 Jan 2012, 3:21 p.m.

On an RPL machine, these two instructions

DUP2 +

So, a full program becomes

\<< 0 1

    2 100 START

        DUP2 +

    NEXT

\>>

(Though, on an RPL machine you'd want to collect the resulting numbers in a list, of course, and not just dump them on the stack.)

With this being too easy for an RPL machine, here's the advanced RPL programming exercise:

if F(10) is 55, what are the last 4 digits of F(10^14)?

(You're welcome to try this on your RPN machine, too, and I will be multo impressed if you can solve this...)

Quote:

---

if F(10) is 55, what are the last 4 digits of F(10^14)?

---

Just a guess: 6875.

lol. For a guess, that's pretty good! Did a calculator help you with your guess?

Yes. I used your program to check it up to F(10^5) on my HP 50g:

%%HP: T(3)A(R)F(,);

\<< 0 1 2 4 ROLL

  START DUP2 + ROT DROP

  NEXT

\>>

  F(10^2)    354224848179261915075

  F(10^3)    final digits:  228875

  F(10^4)    final digits:  366875

  F(10^5)    final digits:  746875

These appear to follow a pattern.

P.S.: The latter has 20,899 digits (~10^5*log₁₀(phi)). It took about one hour or so on the real 50g (I haven't timed it).

Edited: 8 Jan 2012, 7:06 p.m.

Brilliant, Gerson!

Ok.

So, for the tenacious ones: how about the last 14 decimal digits of F(10^14)?

EDIT: you could have used the MOD instruction to keep the number of digits down. That would compute F(10^5) *a lot* faster. But, still, 10^14 is out of reach with this approach.

Edited: 8 Jan 2012, 7:21 p.m.

The last four digits of F(10^4) coincide with the last four digits of F(10000) because 10^14 mod 15000 equals 10000.So it suffices to compute F(10000), easily done on the real 50g. See Pisano periods here. There is a calculator for the cycle lengths (scroll down halfway through the page). For modulus 10000 the result is

Fibs mod 10000: cycle length is 15000

It will not compute the cycle length for modulus 10^14, but it can be extrapolated from lower power of 10. The length appears to be 1.5*10^14, which means this approach won't be of help in this case.

Cheers,

Gerson.

P.S.: The fourth and higher powers of 10 modulus 15000 equals 10000 - so my guess was indeed right :-)

Edited: 8 Jan 2012, 10:07 p.m.

My approach was slightly different: It's clear that the last digit of F_n has to repeat since there are only 10x10 = 100 possibilities for the two previous values. I used the following program with initial values ST X: 1 and ST Y: 0:

X<>Y

RCL+ ST Y

10

MOD

It turns out the length of the period is 60. So I added a loop to find the period for the last two digits:

00 { 27-Byte Prgm }

01 LBL "Fib"

02 X<>Y

03 60

04 X<> ST Z

05 LBL 00

06 X<>Y

07 RCL+ ST Y

08 100

09 MOD

10 DSE ST Z

11 GTO 00

12 END

This gave me the following sequence:

y:  20   40   60   80   0

x:  61   21   81   41   1

Thus the period is: 60 x 5 = 300

Next step was to change line 03 to 300 and line 08 to 1E3 to find another sequence of 5:

y:  600   200   800   400   0

x:  801   601   401   201   1

After changing line 03 to 1500 and line 08 to 1E4 I got a sequence of 10:

y:  8000   6000   4000   2000   0      8000   6000   4000   2000   0

x:  9001   8001   7001   6001   5001   4001   3001   2001   1001   1

Thus I ended up with a period of 60 x 5 x 5 x 10 = 15,000.
From here I did the same as Gerson: 10¹⁴ MOD 15,000 = 10,000.
Changed line 03 to 10000 and line 08 to 1E5 to find:

y: 30626

x: 66875

Cheers

Thomas

Edited: 9 Jan 2012, 3:07 a.m.

Wow, Thomas, very nice! I *am* impressed.

Yes, your guess was right. I volunteer your result is even right to 5 digits.

(Which still doesn't help much toward 14 digits...)

If you look up the Fibonacci page on wikipedia, you'll see various methods for direct computation of F(n).

One is matrix exponentiation: F(n) = [[0 1][1 1]] ^ n

This program for the 50g will tell you what F(100) mod 10^14 is, for example:

\<< 1.E14 R\->I MODSTO

    [[0 1][1 1]] 1.E2 R\->I POWMOD

    {1 2} GET

\>>

Now, if you change 1.E2 to 1.E3 it will produce a negative result. No idea why.

If you try 1.E14 it will tell you: "POWMOD Error: Integer too large"

I conclude that on the 50g POWMOD first does the POW, then the MOD, which, uhm, renders this instruction pretty much useless for matrices. (You want to MOD while you're POW'ing, or else you don't avoid the giant intermediate numbers.)

Well, it works on ND1, and the result is: 88299560546875

Modular matrix exponentiation can be made very fast with the binary method, where you're, essentially, moving in powers of 2 with every iteration. There're only ~50 (~log2(10^14)) matrix multiplications necessary to arrive at 10^14.

This could also be made to work on the 50g with ~10 lines of code to write a proper version of POWMOD. I'm pretty sure this could be made to run in under a minute, for F(10^14).

There's a Project Euler problem (#304) where finding this number (almost this number, there's a different modulus) is the first part of the simple but brutal task of summing up the Fibonacci numbers for each of 100,000 primes following 10^14.

On the 50g, 10 of those primes

\<< 1.E14 R\->I 1 10 START NEXTPRIME NEXT \>>

(I'm not sure what's done on the 50g, but this speed suggests to me that they're not using Miller-Rabin for the probabilistic primality check that kicks in for large numbers, used by NEXTPRIME.)

Once you have the F(NEXTPRIME(10^14)) and F(PREVPRIME(10^14)), you can use normal Fibonacci iteration (with MOD) to compute the spans between primes and eliminate the need to use matrix exponentiation to find each F(n). With an average of 33 numbers from prime to prime in this range, you're looking at having to compute 3.3 million Fibonacci numbers, though.

If you're interested, here's RPL+ code to do all that:

\<< 1234567891011 =m

   1e14 =bias

   1e5 =count

   \<< [[0 1][1 1]] \->big  SWAP

      m modpow

      DUP {1 1} GET

      SWAP {2 1} GET \>> =F2

   0 =sum

   bias =:prev =a

   a F2 EVAL

   1 count START

      1

      a NEXTPRIME =:a prev -

      START

         DUP ROT +

         m MOD

      NEXT

      a =prev

      DUP sum + m MOD =sum 

   NEXT

   DROP DROP

   sum

\>>

Note that the matrix exponentiation method actually computes F(n-1), F(n), F(n+1) at once. So, to continue with Fib iteration you can pick two of them from the result matrix, and you're ready to go.

Quote:

---

So, for the tenacious ones: how about the last 14 decimal digits of F(10^14)?

---

We can combine the Fibonacci- and the Lucas-Sequence to find a fast algorithm to calculate both sequences using the following identities:

The idea is similar to binary exponentiation (a.k.a square-and-multiply algorithm). Here's a Python-program:

#!/usr/bin/python

K = 10**14

B = 2**47

N = 0

M = 10**14

F = 0

L = 2

while B > 1:

    (F, L) = (F*L % M, (5*F*F+L*L)/2 % M)

    B /= 2

    N *=2

    print "%20d %20d %20d %20d %20d" % (N, B, K, F, L)

    if K >= B:

        (F, L) = ((F+L)/2 % M,(5*F+L)/2 % M)

        K -= B

        N += 1

        print "%20d %20d %20d %20d %20d" % (N, B, K, F, L)

And that's the output:

                   N                    B                    K                    F                    L

                   0       70368744177664      100000000000000                    0                    2

                   1       70368744177664       29631255822336                    1                    1

                   2       35184372088832       29631255822336                    1                    3

                   4       17592186044416       29631255822336                    3                    7

                   5       17592186044416       12039069777920                    5                   11

                  10        8796093022208       12039069777920                   55                  123

                  11        8796093022208        3242976755712                   89                  199

                  22        4398046511104        3242976755712                17711                39603

                  44        2199023255552        3242976755712            701408733           1568397607

                  45        2199023255552        1043953500160           1134903170           2537720636

                  90        1099511627776        1043953500160       67194370816120       26026380244498

                 180         549755813888        1043953500160       38521399707760       95774259272002

                 181         549755813888         494197686272       67147829489881       44190628905401

                 362         274877906944         494197686272       73003235747281       62383406970803

                 363         274877906944         219319779328       67693321359042       13699792853604

                 726         137438953472         219319779328       22604627687368       59229375788818

                 727         137438953472          81880825856       40917001738093       86126257112829

                1454          68719476736          81880825856       68797008295097       27914836383243

                1455          68719476736          13161349120       48355922339170       85949938929364

                2910          34359738368          13161349120       29880280387880        1222581444498

                5820          17179869184          13161349120       19590131884240       30816254472002

               11640           8589934592          13161349120       19866585048480       83263801888002

               11641           8589934592           4571414528       51565193468241       91298363565201

               23282           4294967296           4571414528       94702926281441       21451378170403

               23283           4294967296            276447232       58077152225922       47483004788804

               46566           2147483648            276447232       53214104177288       83596643750418

               93132           1073741824            276447232       75272696106384         856675174722

              186264            536870912            276447232       32522299625248       84969227777282

              372528            268435456            276447232       10340608015936       31006195307522

              372529            268435456              8011776       20673401661729       41354617693601

              745058            134217728              8011776       58418769896129       96892716347203

             1490116             67108864              8011776       48860609677187       39467245923207

             2980232             33554432              8011776       67118061778709       99561741164847

             5960464             16777216              8011776       35902403842523       89664428533407

            11920928              8388608              8011776       20398272665861       51376915027647

            23841856              4194304              8011776       15043208059067       60138774356607

            23841857              4194304              3817472       37590991207837       67677407325971

            47683714              2097152              3817472       79233668834727       93180651092843

            47683715              2097152              1720320       86207159963785       44674497633239

            95367430              1048576              1720320       68859452249615       77012557631123

            95367431              1048576               671744       72936004940369       10654909439599

           190734862               524288               671744       76957202272031       55876229280803

           190734863               524288               147456       66416715776417       20331120320479

           381469726               262144               147456       89852350343743       94115666789443

           762939452               131072               147456       10390253505149       64051296250247

           762939453               131072                16384       37220774877698       58001281887996

          1525878906                65536                16384        1539434313208       72834288896018

          3051757812                32768                16384       34530356005744       10933216256322

          6103515624                16384                16384       86304808313568        3648804967682

          6103515625                16384                    0       44976806640625       17586423267761

         12207031250                 8192                    0       74957275390625       45489501953123

         24414062500                 4096                    0       14520263671875       81231689453127

         48828125000                 2048                    0       71563720703125       64288330078127

         97656250000                 1024                    0       90838623046875        8721923828127

        195312500000                  512                    0       85528564453125       49542236328127

        390625000000                  256                    0       60955810546875       86651611328127

        781250000000                  128                    0       57012939453125       87432861328127

       1562500000000                   64                    0       78924560546875       85870361328127

       3125000000000                   32                    0       92950439453125       88995361328127

       6250000000000                   16                    0       50799560546875       82745361328127

      12500000000000                    8                    0       36700439453125       95245361328127

      25000000000000                    4                    0       38299560546875       70245361328127

      50000000000000                    2                    0       11700439453125       20245361328127

     100000000000000                    1                    0       88299560546875       20245361328127

My guess is that this algorithm should be fast on a HP-50g as well. Due to the lack of arbitrary precision the HP-42S doesn't qualify for this.

Kind regards

Thomas

How about another:

01 iC 0

02 iC 1

03 [cmplx]ENTER

04 +

05 PSE 10

06 BACK 003

Six steps, seven with an initial label.

- Pauli

Nice use of ^CENTER as DUP2 in RPL!

%%HP: T(3)A(R)F(,);

\<< \-> n

  \<< 0 1 1 n 2 -

    START DUP2 +

    NEXT n \->LIST

  \>>

\>>

Gerson.

And another:

01 iC 0

02 iC 1

03 +

04 PSE 10

05 y[<->] L

06 BACK 003

Same length though.

- Pauli

There's more than a way to skin a cat. On the HP-15C:

001- f LBL A

002- g CLx

003- f MATRIX 1

004- f LBL 0

005- f PSE

006-   RCL + 0

007- f x<> 0

008-   GTO 0

Edited: 9 Jan 2012, 11:20 p.m.

Nice work, again, Thomas.

Quote:

---

Due to the lack of arbitrary precision the HP-42S doesn't qualify for this.

---

That's what I used to think too, until, almost exactly one year ago, an intrepid bunch (some active posters here) showed the way: Making Fib bits visible (on Silicium.org)

(I think this is my all-time favorite calculator-related discussion. It starts slow but quickly shows real ingenuity. I still can't get over those HP-28C print-outs.)

Now, I'd be willing to bet that F(10^14) mod 10^14 can even be solved on a 15C (LE).

Edited: 10 Jan 2012, 9:08 a.m.

Thanks for link: it's an interesting post! I didn't mean to say it's impossible to do it using the HP-42S. Or we could use Peter's MCODE Multi-Precision Library for HP-41. I might give it a try but this is definitely not my first choice.

Thanks for the challenge as well!

Thomas

I'm very tempted to add some more complex constants to the 34s "[cmplx] 1" e.g.

- Pauli

That would be useful, not just for this task, so go ahead if there's room enough. By the way, is it possible to enter i (1 0) in one step only? Thanks!

Gerson.

I can't think of a way to get i on the stack in one step. Can do it in two though.

We do have a complex version of the constants catalogue which puts zero in Y already. This doesn't include any small integers though. "CPX h-shift CONST"

I've been thinking about possible functions to help with the commands in xrom stuff. If any more get added, some will be exposed to the user but don't expect them to necessarily be pleasant to use.

- Pauli

Quote:

---

We do have a complex version of the constants catalogue which puts zero in Y already. This doesn't include any small integers though. "CPX h-shift CONST"

---

Yes, I had checked that but then I realized it wouldn't make sense to include small integers there. Changing the behavior of the CPX key before numeric keys and [.] might be a option, for instance CPX 123 would put 0 in Y and 123 in X. [CPX] [R/S] might be a non-intuitive shortcut for entering i. I don't know the impact of these changes in the executable size though. Also, this may not be useful enough to be worth the trouble.

Gerson.

A function that does pretty much nothing requires of the order of 100 bytes. This is function table overheads, catalogue overheads and the actual function prelude/epilogue that seems almost unavoidable.

Making small changes to the keyboard handler isn't so bad -- much of this is table based now so adding a command to an unused location is usually cheap.

Adding a complex prefix to digit entry is quite a different change. No idea how large this would be.

The current work is moving as much functionality into internal keystroke programs as we can. Very slow but saving lots of space.

- Pauli

XLBL"i"

iC 1

iC 0

RTN

This should do the trick with minimum overhead. We have to add a catalogue entry and a command table entry but that's not too expensive.

Quote:

---

This should do the trick with minimum overhead. We have to add a catalogue entry and a command table entry but that's not too expensive.

---

Not quite -- the stack dynamics would be incorrect :-(

XLBL"i"

xIN NILADIC_1_COMPLEX

iC 1

iC 0

xOUT xOUT_NORMAL

will work however. The code is 8 bytes, the command table is 10 - 12 bytes (I forget which) and the catalogue entry is 10 bits but can't live in the complex constants catalogue which is compiled automatically -- the complex functions catalogue seems like an okay place to put it.

Pauli

Edited: 11 Jan 2012, 5:21 p.m.

We've got a shorter solution on the 34s now:

01 [cmplx]i

02 RCL+ Y

03 PSE 10

04 x[<->] Y

05 BACK 003

- Pauli

That's what [cmplx]i was for! If the steps 02 and 03 are swapped the sequence will begin with F(0). I'd searched for an i constant in the HP-42S catalog but there was none. I think this first appeared in the HP-28C. If now WP 34S has it, then it was worth posting this simple RPN exercise :-)

BTW, how [cmplx]i is entered?

Gerson.

Quote:

---

BTW, how [cmplx]i is entered?

---

It's in the complex catalogue (CPX 3) but I just mapped it to CPX dot as well. I'm sure Walter will find a better place on the keyboard.

Remember that WP 34S instructions are generally 16 bits long. This allows to cramp so much functionality into a single instruction but costs more bytes than other designs. We are down to 10 bytes here, still pretty decent...

We couldn't put it in the complex constants catalogue since that is automatically produced alongside the real version. Still, the complex catalogue is a decent place and Walter will likely find somewhere decent on the keyboard to locate it. "CPX 1" might be okay which would allow "CPX 0" as well :-)

- Pauli

We have way more than 256 different instructions kind of required this and by doing so the argument can be squeezed into the instruction as well.

We currently have 571 instructions and that is not counting integer, real & complex as different functions (there are 660 if you do that). We've almost got as many different functions as the 15C has total op-codes.

- Pauli

Edited: 12 Jan 2012, 5:38 p.m.

HP 35s

F001 LBL F

F002 0

F003 1

F004 PSE

F005 eqn REGX+REGY

F006 GTO F004

Great! The first 6-step program for an unmodified HP calculator!

Since it does show F(0) and F(1) before looping (thanks to the 2-line display) that's okay. Also, this allows us to match Allen's 10-byte step solution for the HP-42S above:

00 {10-Byte Prgm }

01 LBL 00

02 CLST

03 N!

04 LBL 01

05 PSE

06 RCL+ ST L

07 GTO 01

08 .END.

Gerson.

The 42S not, no, but Free42 decimal (as on the iPhone) can determine the last 12 digits this way.

Cheers, Werner

I rather doubt the 35s solution is less than 10 bytes.

- Pauli

Sorry! I think I didn't make myself clear. I meant showing F(0) and F(1) at the same time like Morten did would allow me to achieve the 10-byte count in the HP-42S. My post was really ambiguous, I realize now.

Gerson.

Challenge accepted:

1E12

STO 00

1E14

XEQ "Fib"

gives:

y: 245361328127

x: 299560546875

Cheers

Thomas

00 { 90-Byte Prgm }

01 LBL "Fib"

02 STO 01

03 LOG

04 2

05 LOG

06 /

07 IP

08 1

09 +

10 2

11 X<>Y

12 Y^X

13 STO 02

14 CLST

15 2

16 X<>Y

17 LBL 00

18 1

19 RCL 02

20 X<=Y?

21 GTO 02

22 RDN

23 RDN

24 RCL ST Y

25 X^2

26 X<>Y

27 STO* ST Z

28 X^2

29 STO+ ST Y

30 2

31 STO/ 02

32 STO/ ST Z

33 *

34 +

35 RCL 00

36 MOD

37 X<>Y

38 RCL 00

39 MOD

40 RCL 01

41 RCL 02

42 X<=Y?

43 GTO 01

44 RDN

45 RDN

46 GTO 00

47 LBL 01

48 -

49 STO 01

50 RDN

51 STO+ ST Y

52 2

53 STO/ ST Z

54 *

55 RCL+ ST Y

56 RCL 00

57 MOD

58 X<>Y

59 RCL 00

60 MOD

61 GTO 00

62 LBL 02

63 RDN

64 RDN

65 END

Edited: 14 Jan 2012, 3:11 a.m.