The Museum of HP Calculators

HP Forum Archive 14

 geometric progressionMessage #1 Posted by Jim Kimes on 15 June 2004, 3:09 p.m. I'm trying to establish a geometric progression formula for the following series but I'm not sure it's a geometric progression. I need to be able to solve the nth term of the following series: (364/365) * (363/365) * (362/365)... I want to use 'a' for the changing numerator and 'b' for unchanging denominator and r for the ratio. I know the formula for a geometric progression but it's not that simple. What is the answer here? I need a formula that will fit all similar conditions that I can program on an HP 41. Right now I'm getting the answer by looping but I think there's a better way. Thanks.

 Re: geometric progressionMessage #2 Posted by hugh steers on 15 June 2004, 6:49 p.m.,in response to message #1 by Jim Kimes how about (aPr)/b^r for r terms. where aPr means permutations of a taking r. ? Edited: 15 June 2004, 6:51 p.m.

 Simplified formMessage #3 Posted by Tizedes Csaba [Hungary] on 16 June 2004, 3:48 a.m.,in response to message #1 by Jim Kimes Hello, try this: ```FACTORIAL(n-1)/n^(n-1) = GAMMA(n)*n^(1-n) in your case n=365 ``` Csaba

 Simplified form - I think I was crazy...Message #4 Posted by Tizedes Csaba [Hungary] on 16 June 2004, 3:56 a.m.,in response to message #3 by Tizedes Csaba [Hungary] `FACTORIAL(n-1)/n^(n-1) = GAMMA(n)*n^(1-n)` I'm so sorry, I think my mind is not workin today... And I can't to delete my message... Sorry again! Cs.

 Re: Simplified form - I think I was crazy...Message #5 Posted by Jim Kimes on 16 June 2004, 2:33 p.m.,in response to message #4 by Tizedes Csaba [Hungary] The problem with factorials is that calculators and computers can't handle super large factoials. They overflow. Besides I'm looking for a partial or a truncated factorial. Logs of factorials may be able to handle that. With your proposed solution, when a = 364, b = 365 and n = 23, do you get .493? If not the real problem isn't solved. Try the problem manually, doesn't take that long, and verify the .493 answer and then formulate it for any value of a & b & n, and then you will feel comfortable with working on a solution.

 Re: geometric progressionMessage #6 Posted by bill platt on 16 June 2004, 3:25 p.m.,in response to message #1 by Jim Kimes Quote: (364/365) * (363/365) * (362/365)... Won't this series go on for 364 cycles, decreasing all the time, and at an increasingly rapid rate, until you get to (0/365) at whcih point it goes disjointly to zero? You could very easily write a loop, using i as a counter variable, and loop 364 times.....what am I missing here? Regards,Bill

 Re: geometric progressionMessage #7 Posted by Jim Kimes on 16 June 2004, 8:16 p.m.,in response to message #6 by bill platt "You could very easily write a loop, using i as a counter variable, and loop 364 times.....what am I missing here?" You're under the impression that I haven't solved my problem. Oh, I've solved it. I have a loop routine that works very well on a 41 but I can't loop with a 19BII and so I want to use some kind of progression formula that will find partial factorials. At n = 23 the ratio becomes .493. I'm beginning to think there is no formula. I haven't heard back from Dr. Math yet.

 A solution to Jim's problem [LONG]Message #8 Posted by Valentin Albillo on 17 June 2004, 4:58 a.m.,in response to message #7 by Jim Kimes Hi, Jim: Jim posted: "I can't loop with a 19BII and so I want to use some kind of progression formula that will find partial factorials. At n = 23 the ratio becomes .493. I'm beginning to think there is no formula." Yes, there is. A "partial" factorial is technically called a Pochhammer Symbol, which itself is nothing more than a Gamma function divided over another gamma function, like this: ``` (x)n = x*(x+1)*...*(x+n-1) = Gamma(x+n)/Gamma(x) ``` This is the clue for solving your problem non-iteratively. In HP-71B BASIC (i.e.: nearly plain-vanilla English), your problem would be solved like this: ``` 10 B=365 @ N=23 20 DISP FACT(B)/FACT(B-N)/B^N ``` where FACT is the factorial function. However, upon running this results in an overflow, as FACT(365) exceeds the range for real numbers in the HP-71B and in most any other calculator or system. What can we do, then ? Easy. Just use the asymptotic series for the factorial function, aka Stirling's series, to compute the natural logarithm of x! instead of x!. This computation won't overflow for your arguments, and once we have the logarithm, we just use the exponential function EXP (aka e^x) to deliver the final result. Stirling series is: ``` x! = (x/e)^x*Sqrt(2*Pi*x)*(1+1/(12*x)+1/(288*x^2)-139/(51840*x^3)+ ...) ``` so its logarithm will be: ``` z = ln(x!) = x*ln(x)-x+ln(2*Pi*x)/2+ln(1+1/(12*x)+1/(288*x^2)-139/(51840*x^3))+ ... ``` and of course your final result is simply e^z. The resulting, non-iterative HP-71B BASIC program is: ``` 10 B=365 @ N=23 20 DISP EXP(FNF(B)-FNF(B-N)-N*LN(B)) 30 END 90 DEF FNF(X)=X*LN(X)-X+LN(2*PI*X)/2+LN(1+1/12/X+1/288/X^2-139/51840/X^3) ``` where DISP is the statement used to show a result on the display (aka PRINT, VIEW, etc), EXP is the exponential function e^x, LN is the natural logarithm function, PI is 3.14159265359, and FNF is a user-defined function, which can be substituted for a subroutine in other languages, passing it the appropriate parameter. Upon running, it promptly returns: ```>RUN .492702772644 ``` which agrees well with your .493 and is accurate to nearly 8 decimal places (actually, the exact result is 0.49270276567601459277458277+). If still more accuracy is desired, you can use some additional terms to the Stirling's series, have a look at this link: and its twin for the denominators. You can get easily 12+ decimals by using more terms. On the other hand, if you don't need that much accuracy, you can save RAM and execution time by deleting terms from the Stirling's series above, like this: ``` 90 DEF FNF(X)=X*LN(X)-X+LN(2*PI*X)/2+LN(1+1/12/X) ``` though you should test if the accuracy you then get is adequate for your intended purposes, particularly for large N. For your specific numeric example, you still get 8 decimals ! Hope that helps. Best regards from V. Edited: 17 June 2004, 10:26 a.m.

 Re: A solution to Jim's problem [LONG]Message #9 Posted by Jim Kimes on 17 June 2004, 5:36 p.m.,in response to message #8 by Valentin Albillo This is rich stuff! We're definitely on the right track. I've got a monthly newsletter to get out after which I'll try to digest this. Thank you, Valentin. You just may have solved my problem. I never in my life have heard of a Pochhammer Symbol. It sounds like a novel by Frederick Forsyth.

 The "Birthday Problem"Message #10 Posted by Karl Schneider on 17 June 2004, 12:32 a.m.,in response to message #6 by bill platt This problem has an application in the well-known statistical exercise called the "birthday problem". "Among 25 randomly-selected people in a room (with no twins), what is the probability that at least two people in the room have the same birthday?" To answer this question, determine the probability of no one sharing a birthday. (Ignore leap year): ```prob = (364/365)*(363/365)*...*(341/365) = Perm(364,24) / 365^24 = 0.43130 ``` So, the chance that two of 25 random people have the same birthday is (1-0.43130) * 100% = 56.87% Rather counter-intuitive, no? Of course, for larger values in the quotient, a term-by-term looped approach must be programmed. -- Karl S.

 Re: The "Birthday Problem"Message #11 Posted by John Mosand (Norway) on 17 June 2004, 11:40 a.m.,in response to message #10 by Karl Schneider When I saw the original posting, it immediately reminded me of the 'birthday problem'. About 20 years ago, I tried it out at our office. We were 20 persons at that time. And yes, there were two with the same birthday. Happened to be a girl who had recently started to work with us and myself. The counter-intuitive side of it appears because one may tend to think of the chances that two particular persons share birthdays - which is something quite different!

 Re: The "Birthday Problem" (drifting OT)Message #12 Posted by Cameron on 17 June 2004, 1:26 p.m.,in response to message #11 by John Mosand (Norway) My first sweetheart's dad pulled that trick at a dinner party when I was in my mid teens (he was a statistician who worked for ICI). The occasion is etched on my mind to this day. I was convinced that professionals that played with numbers were either demigods or magicians--a feeling reinforced when I read some of the postings on this forum. Suffice it to say that I never could cut the mustard in the numbers game whic is probably why I remain in awe of those that do. ;-) Cameron xyzzy

 Re: The "Birthday Problem"Message #13 Posted by Jim Kimes on 17 June 2004, 5:29 p.m.,in response to message #11 by John Mosand (Norway) It is indeed. I've written a very provocative manuscript I'm trying market. I can't tell you all much about it now for the sake of confidentiality but when it's published, everyone is welcome to it.

 Re: The "Birthday Problem"Message #14 Posted by Jim Kimes on 25 June 2004, 8:45 a.m.,in response to message #10 by Karl Schneider Karl, this is the best approach for solving this problem. Not only is it easy to program for a calculator with a SOLVER, it is also easy to explain to groups as to what is going on mathematically. I'm grateful to you and everyone for your interest in this topic. It has fascinated me for years. I've had one paper published on it and I have a newspaper interested in an update. I think this is really a great forum. Extremely helpful!

 Re: geometric progressionMessage #15 Posted by bill platt on 17 June 2004, 10:12 a.m.,in response to message #6 by bill platt Hi, I couldn't resist, and so I wrote a loop lat night, both on the 32s and the 11c. The 11c is much slower. The 32s prompts for the number of terms (i) and the starting numerator (A). The 11c must be hand-laoded to I and 0 respectively before GSB A to start. 32s starts with xeq M. Here is the 11c: hERE IS 32S: ``` LBL B LBL B 365 365 / / RTN RTN LBL E LBL E RCL 2 RCL B RTN RTN LBL A LBL M RCL 0 INPUT A GSB B INPUT i STO 2 RCL A 1 XEQ B STO-0 STO B RCL I 1 X<>Y STO-A - STO-i STO I GTO C GTO C LBL C LBL C RCL I RCL i X=0? X=0? GTO E GTO E RCL 0 RCL A GSB B XEQ B STO * 2 STO * B 1 1 STO - 0 STO - A RCL I STO -i X<>Y GTO C - STO I GTO C ``` regards, Bill Edited: 17 June 2004, 10:29 a.m.

 A better 11/15/34 routine for "geometric progression"Message #16 Posted by Karl Schneider on 19 June 2004, 3:29 a.m.,in response to message #15 by bill platt Bill -- I just couldn't resist writing a tighter program for the HP-11C/15C/34C for the "birthday problem" that is shorter and simpler in flow, requires no initialization of memory, uses fewer labels, and utilizes the built-in loop-counter function "DSE": ```LBL "A" STO 0 store number of people X<->Y INT STO 1 store integer number of days LST x FRAC STO 2 store fractional number of days RCL 1 RCL 0 - EEX 3 / RCL 1 + STO I loop counter variable is set RCL 2 STO + 1 replace days in year with full floating value 1 STO 3 initialize answer variable LBL "0" start of processing loop DSE I GTO "1" RCL 3 RTN display final result and stop LBL "1" calculate term RCL I INT RCL 2 + RCL 1 / STO x 3 GTO "0" Usage: (number of days in year -- may be non-integer) ENTER (number of people) GSB "A" : : (read probability of all unique birthdays) ``` This program can accommodate leap years by representing 365.25 days (or a more precise value) in a year. This is somewhat useful, because you can't do ```Perm(365.25,25) / (365.25^25) ``` Edited: 22 June 2004, 2:08 a.m. after one or more responses were posted

 Re: A better 11/15/34 routine for "geometric progression"Message #17 Posted by bill platt on 21 June 2004, 5:35 p.m.,in response to message #16 by Karl Schneider Hi Karl, Thank you very much for your post. I was really hoping to get a critique / response! I figured I would get a response regarding the possibility of using DSE. The thing was, I can never remember the format for DSE without checking the manual--so I made my own counter on the spot. I discovered that you can run my 32s version on the 20s---with no modification except the Loop being " / 365 =" instead of "365 /" (There is no DSE/ISG on the 20s). The other conundrum (my weakest part with loops) was getting the count right---you know--the "when do you count--before or after the action" and so when I tried later to convert to DSE and use a single loop, I ended up with a program that ran fine, except that the counter was off by one. So thank you for the chance to see how to make it work correctly. Logic and problem solving are interesting. My solution I posted was very much the first attack---and I programmed it in that classically useful way of doing it manually and recording what you do. It is for that reason that I ended up with an initial solve, store to a result variable, then go into a loop and add to the result variable. I knew it should be possible to do it in one fell swoop---but I couldn't see it right then and there. But programming is a skill, and you only get better with practice. and to that end, this was a very good exercise. The problem for me is that I am not a programmer---yet every year, or sometimes after 2 or three years, I find myself needing to write programs---in a very specialized languege----and so while rusty, I do need to know what I am doing. (The language is written for developing stability caclulations for ships---and only infrequently do I need to write a custom program to attack an unusual problem).) So, thanks again! Best regards, Bill

 Thanks, Bill...Message #18 Posted by Karl Schneider on 22 June 2004, 2:13 a.m.,in response to message #17 by bill platt ... Check my previous post (to which you commented) for an edited version that can calculate correct results with non-integer days in a year.

 Re: geometric progressionMessage #19 Posted by Jim Kimes on 25 June 2004, 8:39 a.m.,in response to message #15 by bill platt OK, my newsletter is mailed so I can get back to the ‘birthday paradox.’ Below is my contribution to the recent spate of programs to solve the birthday paradox. My program appears to be beefy at first glance, but it does a lot of things automatically, including prompts. First, in this problem, what you’re really solving for is the following: 1 ~ [(364/365) * (363/365) * 362/365)…] So that leaves us with three variables: 365, or what I call ‘max’; number of iterations, which I call ‘group’: and the probability at a certain group count. You would never want to solve for ‘max’ but you have to input it. This program can possibly be done with fewer steps by taking ‘max’ to the group number power. You can try that; it would be interesting, but I don’t think you’re going to cut down very many steps. I’m not a proficient BASIC programmer but I’m sure a shorter program can be written in either in Virtual or old BASIC. There is no pride of authorship here so have at it. We can all learn. Here’s the way the program looks on the HP 19BII, using the solver: 1 – (Perm ( X:Y) / (365^Y)) = P where X = universe (365), Y= no. in group (23, and P = target or resultant probability of a matching birthday. Step Action Explanation for action 01 LBL BDP Label BirthDay Paradox 02 CLRG Clear all registers 03 CF01=GRP Reminder of which variable to solve 04 PROMPT 05 1 Explained in text above 06 STO 3 07 MAX? What is size of universerse? 08 PROMPT 09 STO 2 10 PROB? What is target probability, answer will be in group size 11 PROMPT 12 STO 04 13 GRP? What is group size, answer will be in probability 14 PROMPT 15 STO 05 16 LBL 01 LBL 01 begins the routine to solve for group size 17 1 18 STO +01 19 RCL 02 20 RCL 01 21 --- 22 RCL 02 23 / 24 RCL 03 25 * 26 STO 03 27 1 28 --- 29 CHS 30 STO 06 31 FS? 01 Sets flag to branch to probability routine, otherwise continue 32 GTO 02 33 RCL 04 Recalls target probability to compare with value in Step 30 34 X>Y? 35 GOTO 01 If probability doesn't satisfy, return to LBL 01 & iterate 36 STOP 37 VIEW 01 View the group size since the routine has stopped 38 STOP 39 VIEW 06 View exact probability since it won't exactly = target 40 LBL 02 Starts new routine solvling for Probability 41 RCL 05 Steps 39 & 40 compare the group counter and 42 RCL 01 compares in Step 41 43 X=Y? Here can use X =Y because number in group will be an integer 44 GTO 03 If they equal then branch to Step 44. Read probability 45 GTO 01 46 LBL 03 47 VIEW 06 48 END

Go back to the main exhibit hall