Re: Program to calculate rational function with imbedded continued fraction
- To: mathgroup at smc.vnet.net
- Subject: [mg69957] Re: Program to calculate rational function with imbedded continued fraction
- From: Jean-Marc Gulliet <jeanmarc.gulliet at gmail.com>
- Date: Thu, 28 Sep 2006 06:15:32 -0400 (EDT)
- Organization: The Open University, Milton Keynes, UK
- References: <efdkq5$10v$1@smc.vnet.net>
Diana Mecum wrote:
> On 9/26/06, Andrzej Kozlowski <akoz at mimuw.edu.pl> wrote:
>> On 26 Sep 2006, at 18:21, Diana wrote:
>>
>>> Hello all,
>>>
>>> Would someone be able to help me write a program to calculate the
>>> rational function of (2x + 3)/(2x + 4), for example?
>> I have once already posted such a program and so did Daniel
>> Lichtblau. (Search the archive for the subject "How to do Continued
>> fraction of polynomials" and compare our posts. (I can't remember any
>> more what the difference between them was).
>>
>> Here is my version:
>>
>> F[f_, g_] :=
>> If[PolynomialReduce[f, g][[1, 1]] =!= 0, PolynomialReduce[f,
>> g][[
>> 1, 1]] + F[PolynomialReduce[f,
>> g][[2]], g], Module[{u = g, v = f, p, ls}, ls = Flatten[Last[
>> Reap[While[u =!= 0, p = PolynomialReduce[u, v]; u = v; v =
>> p[[2]]; Sow[
>> p[[1]]]]]]]; 1/Fold[Function[{x, y}, y + 1/x], Infinity,
>> Reverse[ls]]]]
>>
>> In the case of your example above you need to evaluate:
>>
>> In[50]:=
>> k=(2x+3)/(2x+4);
>>
>> In[51]:=
>> F[Numerator[k],Denominator[k]]
>>
>> Out[51]=
>> 1 + 1/(-2*x - 4)
[...snipped]
> I see your program. Thanks for forwarding it.
>
> To implement it, do I plug (e(1) + 1)(T^2e(1) in for k?
>
> When I run the routine, I get an error saying that the tags in In[50]
> and In[51] are protected.
Do not copy and paste the tags In[]/Out[] into Mathematica, just
whatever Mathematica expressions are below any In[] tag -- or between an
In[] and Out[] tags.
(* Start copy here *)
F[f_, g_] := If[PolynomialReduce[f, g][[1,1]] =!= 0, PolynomialReduce[f,
g][[1,1]] + F[PolynomialReduce[f, g][[2]], g], Module[{u = g, v = f, p,
ls}, ls = Flatten[Last[Reap[While[u =!= 0, p = PolynomialReduce[u, v]; u
= v; v = p[[2]]; Sow[p[[1]]]]]]]; 1/Fold[Function[{x, y}, y + 1/x],
Infinity, Reverse[ls]]]]
k = (2*x + 3)/(2*x + 4);
F[Numerator[k], Denominator[k]]
(* End copy here *)
will return
1
1 + --------
-4 - 2 x
Regards,
Jean-Marc