Search code examples
mathwolfram-mathematicadiscrete-mathematicsmathematica-8

How do I define a function in mathemathica that uses a complicated recurrence relation?


I am trying to write a little script to compute an integer sequence. The function I'm trying to write in code is the one int the blackboard, a(n). The problem is that I was expecting the function h(n) I defined in the script to give a number as a result, but it is giving something else: for h(2) it gives ArgMax[{p, HarmonicNumber[p] <= 1}, p, Integers] How can I correct that? (You have to understand that I am by no means a programmer, nor know much about mathematica. Thanks in advnace. The script I wrote is this one:

    h[n_] := (ArgMax[{p, 
      Sum[1/s, {s, 1 + Sum[h[k], {k, 1, (n - 1)}], p}] <= 1}, p, 
     Integers]) - Sum[h[k], {k, 1, (n - 1)}]; h[1] = 1;

Image of the definition written by hand a(n)=(maximum p such that the sum from s equals r to p is less or equal than one)-r+1, where r=1+the sum from k=1 to (n-1) of a(k), and a(1)=1

PD: Those letters that look like v's are r's. Sorry.


Solution

  • a[1] = 1;
    
    a[n_] := Module[{sum = 0},
      r = 1 + Sum[a[k], {k, n - 1}];
      x = r;
      While[sum <= 1, sum += 1/x++];
      p = x - 2;
      p - r + 1]
    
    Table[a[n], {n, 6}]
    

    {1, 2, 6, 16, 43, 117}

    The result for a[4] is 16 not 14.

    To illustrate, when n = 4

    r = 1 + Sum[a[k], {k, 4 - 1}]
    
      = 1 + a[1] + a[2] + a[3]       (* refer to established results for a[n] *)
    
      = 1 +  1  +  2  +  6  =  10
    
    sum = 0;
    x = r;
    While[sum <= 1, sum += 1/x++];
    p = x - 2;
    p - r + 1
    

    16

    Or in another form

    Total[Table[1/s, {s, 10, 25}]] <= 1   (* True *)
    
    p - r + 1 = 25 - 10 + 1 = 16
    

    Using memoisation, as mentioned by ogerard

    Clear[a]
    
    a[1] = 1;
    
    a[n_] := a[n] = Module[{sum = 0},
       r = 1 + Sum[a[k], {k, n - 1}];
       x = r;
       While[sum <= 1, sum += 1/x++];
       p = x - 2;
       p - r + 1]
    

    only reduces the following run time by 9 seconds

    Timing[Table[a[n], {n, 14}]]
    

    {40.8906, {1, 2, 6, 16, 43, 117, 318, 865, 2351, 6391, 17372, 47222, 128363, 348927}}