[VA] Short & Sweet Math Challenge #25 "San Valentin's Special: Weird Math"

[Valentin Albillo]

      
Hi all,

Thank you very much for your overwhelming participation in my S&SMC#25, much, much appreciated. I hope you've enjoyed its two main tasks: sleuthing your way through the challenges, using your brains and your math intuition, then writing efficient code for your favorite HP calcs to let them take away the drudgery. Seeing the many posted solutions and insights, it's quite clear you did well !

Now for my original solutions. Giving the solutions to all six Concoctions at once would make for a tremendously long post which would be most unwieldy to write and to read, so (as I did in some previous S&SMC) I'll give the solutions one Concoction at a time, beginning with the first one, immediately below: first my code and Sleuthing process, then my Results, and finally my extensive Comments. Last but not least, The Hall of Fame !

Note: My HP-71B code may use keywords from the JPC ROM, MATH ROM, HP-IL ROM and STRINGLX LEX file, executed on go71b, while RPN code is for the HP-42S, executed on a DM42.


My original solution for "Concoction the First: Weird Limit"

The Sleuthing

This simple 2-line, 96-byte HP-71B program asks for the sum to exceed, then simulates 10, 100, ..., 100,000 tests and outputs the average counts :

      1  DESTROY ALL @ INPUT N @ FOR L=1 TO 5 @ RANDOMIZE 1 @ T=10^L @ M=0 @ FOR I=1 TO T
      2  S=0 @ WHILE S<N @ S=S+RND @ M=M+1 @ END WHILE @ NEXT I @ DISP T;M/T @ NEXT L

      >RUN
      ? 1
            10      2.7            { error:  0.01828, see below }
            100     2.79           { error: -0.07172 }
            1000    2.729          { error: -0.01072 }
            10000   2.7084         { error:  0.00988 }
            100000  2.71959        { error: -0.00131 }

This simpler 24-step, 46-byte RPN program accepts the sum and the number of tests and produces the average count to exceed it:

      LBL "SUMRND"  SEED      X<Y?
      "Wait..."     LBL 01    GTO 00
      AVIEW         CLX       DSE ST Z
      STO 01        LBL 00    GTO 01
      X<>Y          ISG 00    RCL 00
       0            LBL 00    RCL/ 01
      STO 00        RAN       CLD
      E^X            +        END

         1 [ENTER] 1E6 [XEQ] "SUMRND" -> 2.717405     { error: 0.0008768 }
         1 [ENTER] 1E7 [XEQ] "SUMRND" -> 2.7180730    { error: 0.0002088 }

Running one or the other (whichever is faster), we get the following results:

      N    S = 1      S = 2      S = 2.021  S = 3      S = Pi     S = 4
      ---------------------------------------------------------------------
      10   2.7        4.8        4.8        6.6        7.1        9.1
      100  2.79       4.72       4.75       6.81       7.13       8.85
      1E3  2.729      4.701      4.737      6.701      6.978      8.708
      1E4  2.7084     4.6751     4.7132     6.6757     6.9571     8.6718
      1E5  2.71959    4.67827    4.71806    6.66808    6.95027    8.66601
      1E6  2.717405   4.671500   4.712124   6.666287   6.949634   8.664481
      1E7  2.7180730  4.6707911  4.7118212     -       6.9500464     -

      N    S = 5      S = 5+1/6  S = 10     S = 15     S = 20      
      ------------------------------------------------------------
      10   11.2       11.5       20.9       30.6       41.1
      100  10.77      11.22      20.66      30.72      40.82  
      1E3  10.719     11.07      20.623     30.604     40.68
      1E4  10.673     11.0006    20.6876    30.7079    40.7027
      1E5  10.66641   10.99947   20.65914   30.66700   40.66927
      1E6     -          -       20.664967     -          -


The Results

Considering the data obtained above, these are my answers:

• a. What do you think is the average count of generated random numbers for their sum to exceed 1 ? Can you recognize what the exact value would be ?
  
        The limit seems to be ~ 2.7181, which I recognize as the constant e = 2.7182+
   
• b. What would the average count be for the sum to exceed 2 ? To exceed 2.021 ? To exceed Pi ?
  
        Count(2) = ~ 4.67079, Count(2.021) = ~ 4.71182, Count(Pi) = ~ 6.95004
   
• c. What do you think is the asymptotic expression for the average count needed to exceed large values of the sum ?
  
        We have Count(5) = ~ 10.666, Count(10) = ~ 20.665, Count(15) = ~ 30.667, Count(20) = ~ 40.669
  
        so the asymptotic expression for large sums S seems to be: Count(S) ~ 2*(S + 1/3)
   
• d. Can you explain the constant component of said expression ?
  
  The constant in the asymptotic expression seems to be 1/3, and there's a formal explanation but I'll give here an easier, informal one: when the last random number added up results in a sum exceeding the given threshold, the overshoot is also uniformly distributed in [0,1), and the expected value of this overshoot is the expected value of the minimum of two uniform random variables (the random number and the overshoot), which theoretically is 1/3 (i.e.: average of MIN(RND,RND) = 1/3).
  
  If you don't know that piece of theory,this modification of the above HP-71B program will produce it for a given sum by running 10-10,000 tests. It just simulates the process looking at the average overshoot (instead of the average of the count of random numbers generated), and computing the average of the fractional parts (overshoot) of the sum (for the case of integer thresholds), like this:
  
  1  DESTROY ALL @ INPUT N @ FOR L=1 TO 4 @ RANDOMIZE 1 @ T=10^L @ M=0 @ FOR I=1 TO T
  2  S=0 @ WHILE S<N @ S=S+RND @ END WHILE @ M=M+FP(S) @ NEXT I @ DISP T;M/T @ NEXT L
  
  >RUN
  ?           S = 5     S = 10    S = 15       
      ------------------------------------
      10      0.356...  0.262...  0.289...
      100     0.305...  0.321...  0.276...
      1000    0.323...  0.335...  0.333...
      10000   0.333...  0.332...  0.333...
  

The Comments

The limit average count for the sum of a series of [0,1) uniformly distributed random numbers to exceed 1 is exactly e = 2.71828182845904523536+, the base of the natural logarithms, which is pretty "weird" and can be considered an analog of Buffon's Needle experiment to estimate the value of Pi. Here we don't throw needles on a grid but merrily add up random numbers keeping count and we get e instead.

These are the exact theoretical values for the sum to exceed S:

      S  Average Count                               20-decimal value
      ----------------------------------------------------------------------
      1  e                                           2.71828182845904523536
      2  e² - e                                      4.67077427047160499187
      3  (2 e³ - 4 e² + e)/2                         6.66656563955588990415
      4  (6 e⁴ - 18 e³ + 12 e² - e)/6                8.66660449003269543723
      5  (24 e⁵ - 96 e⁴ + 108 e³ - 32 e² + e)/24    10.66666206862241185802 

This is the general formula to numerically compute the theoretically exact value ...



... and this is my simple 1-line, 53-byte HP-71B program to instantly compute them given the sum to exceed:

      1  DESTROY ALL @ INPUT X @ S=0 @ FOR K=0 TO IP(X) @ S=S+(K-X)^K/FACT(K)*EXP(X-K) @ NEXT K @ DISP S

      >RUN
      ?        1      2.71828182846         2.021      4.71182750642
               2      4.67077427047            Pi      6.94950388760
               3      6.66656563953         5+1/6      11.0000029914
               4      8.66660448999
               5      10.6666620697
              10      20.6666664745

As the successive terms have alternating signs and go on increasing, the precision for large sums (say > 10) degrades very quickly and we can see that for Sum > 100 the result is but garbage:

             100     -2.69702821806E43     { correct value: 200.666666666... }

Implementing the above exact formula in RPN (24 steps, 37 bytes) and running it in the 34-digit DM42 we get the following:

      LBL "FSUMRN"   E^X         x
      STO 00         LASTX       +
      IP             +/-        DSE 01
      STO 01         RCL 01     GTO 00
       0             Y^X        RCL 00
      LBL 00         LASTX      E^X
      RCL 00          N!         +
      RCL- 01         /         END

A sample run would be:

      1 [XEQ] "FSUMRN"   ->    2.71828182846

and assorted results truncated to 20 decimal digits (use [SHOW] to see them in full) would be:

     Sum   Average count               Sum     Average count
    -------------------------------------------------------------------
      1    2.71828182845 904523536    2.021    4.71182750642 763255399
      2    4.67077427047 160499187        e    6.10400234136 375415166
      3    6.66656563955 588990414       Pi    6.94950388752 954473480
      4    8.66660449003 269543722    5+1/6   11.00000299090 420020529
      5   10.66666206862 241185801   20+1/6   40.99999999999 999999992
     10   20.66666666647 631880061    20.21   41.08666666666 666666663
     15   30.66666666666 666034379   21+1/6   43.00000000000 000000000
     20   40.66666666666 666666648

which are fully correct to the digits shown. Though the precision attained using the 34-digit DM42 is much greater than using the 12-digit HP-71B program, for large sums it will still quickly degrade. For instance:

     50   100.666667031        { only ~ 6 correct decimals, about 9 digits in all }
    100   1.12982914443E21     { utter garbage, should be 200.66666666666666666...}


The Hall of Fame

Some of you did bravely tackle this Concoction the First: Weird Limit, namely these four experts:

• robve posted code for the HP Prime, generic BASIC and Python, as well as correct results for questions a and b, gave c a try and left d unanswered. A fine achievement which would be finer had he followed the rules re ony HP models and HP languages, as everyone else did.
   
• ramon_ea1gth posted RPL code and correctly guessed e as the limit.
   
• Werner posted RPN code with results identical to robve's post, but correctly guessed the asymptotic expression.
   
• Nihotte(lma) posted RPN code for the HP-35s and RPL code for the HP48G, correctly guessed the limit to be e and stated that his results match those of robve. He also used a novel (and bold!) thinking-out-of-the-box attempt to guess the asymptotic expression using the HP10BII+ to fit the data to a linear regression, getting correctly the 2*S part, and later trying out yet another approach.

That's all for now, thanks a lot to those who contributed, I really hope you enjoyed it. I'll post my original solutions for "Concoction the Second: Weird Sum" in a couple of days. Stay tuned !  

Best regards.
V.