home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMBERTH.PAK / CONTINUE.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  3.5 KB  |  117 lines
  1. (* Copyright 1988, Wolfram Research, Inc. *)
  2.  
  3. (*:Version: Mathematica 2.0 *)
  4.  
  5. (*:Name: NumberTheory`ContinuedFractions` *)
  6.  
  7. (*:Author:
  8. *)
  9.  
  10. (*:History:
  11.     Modified Mar. 92 by John M. Novak:
  12.         allow rational numbers
  13.         fix behavior when partial fraction is finite, less than
  14.             order requested.
  15.     Modified May 92 by John M. Novak:
  16.         allow non-rational non-floating point numbers (guesses precision)
  17.         force machine numbers to bignums (machine precision) to reduce
  18.             or eliminate roundoff error (with suggestions from J. Keiper).
  19. *)
  20.  
  21. (*:Keywords:
  22.     number theory, continued fractions
  23. *)
  24.  
  25. (*:Requirements: none. *)
  26.  
  27. (*:Discussion:
  28.     ContinuedFraction[x,n] generates the first n quantities in the
  29.     continued fraction expansion for the number x.  The result is
  30.     expressed as the object ContinuedFractionForm[ list ].
  31.     All numbers are forced to bignums, minimizing roundoff error
  32.     problems.  Note that continued fraction may have order less
  33.     than n, if sequence terminates before n is reached.
  34.     If non-rational non-floating point number is given, the
  35.     routine tries to guess the precision needed to evaluate
  36.     the number to get a continued fraction of the correct
  37.     order.  From Jerry Keiper:
  38.        "It turns out that each term in a continued fraction expansion
  39.         requires an expected 0.429017 decimal digits (i.e.
  40.         Log[10, Kinchin], cf Ilan Vardi's book, p163).  Thus a
  41.         comfortable overestimate might be to start with one decimal
  42.         digit for every two continued fraction coefficients requested."
  43.     Normal[ ContinuedFractionForm[ list ] ] returns a number
  44.     very close to the original x given.
  45. *)
  46.  
  47. BeginPackage["NumberTheory`ContinuedFractions`"]
  48.  
  49. ContinuedFraction::usage = 
  50. "ContinuedFraction[x,n] generates the continued fraction 
  51. representation for the number x to order n.  Note that the
  52. order returned may be less than n if the continued fraction
  53. terminates before n steps."
  54.  
  55. ContinuedFractionForm::usage =
  56. "ContinuedFractionForm[{a0,a1,..}] is a continued-fraction;
  57. it is converted to an ordinary rational number using Normal."
  58.  
  59. Begin["`Private`"]
  60.  
  61. ContinuedFraction::incomp =
  62. "Warning: ContinuedFraction either terminated or argument was
  63. unable to be evaluated to sufficient precision in a reasonable
  64. number of attempts.";
  65.  
  66. ContinuedFraction[x_?MachineNumberQ,n_Integer?Positive] :=
  67.     ContinuedFraction[SetPrecision[x, $MachinePrecision],n]
  68.  
  69. ContinuedFraction[x:(_Real | _Rational), n_Integer?Positive] :=
  70.     ContinuedFractionForm[Floor[Prepend[cf[x,n - 1],x]]]
  71.  
  72. ContinuedFraction[x_?(Not[NumberQ[#]] && NumberQ[N[#]] &),
  73.         n_] :=
  74.     Module[{nx, prec = n/2 + 10, out = {}, i = 1},
  75.         While[Length[out] < n-1 && i++ < n,
  76.             nx = N[x, prec];
  77.             out = cf[nx, n-1];
  78.             prec += (n - Length[out])/2 + 10
  79.         ];
  80.         If[Length[out] < n-1,
  81.             Message[ContinuedFraction::incomp]
  82.         ];
  83.         ContinuedFractionForm[Floor[Prepend[out,nx]]]
  84.     ]
  85.  
  86. cf[num_,ord_] := {}/;((num - Floor[num]) == 0)
  87.  
  88. cf[num_,0] := {}
  89.  
  90. cf[num_, ord_] := 
  91.     Module[{tmp = 1/(num - Floor[num])},
  92.         Join[{tmp}, cf[tmp,ord - 1]]
  93.     ]
  94.  
  95. (*
  96. ContinuedFractionForm/: Normal[ContinuedFractionForm[a_List]] :=
  97.                 Last[Accumulate[ 1 / #1 + #2 &, Reverse[a]]]
  98. *)
  99.  
  100. ContinuedFractionForm/: Normal[ContinuedFractionForm[a_List]] :=
  101.                 Fold[ 1 / #1 + #2 &, Last[a], Rest[ Reverse[a]]]
  102.  
  103. End[]  (* NumberTheory`ContinuedFractions`Private` *)
  104.  
  105. Protect[ ContinuedFraction, ContinuedFractionForm ]
  106.  
  107. EndPackage[]  (* NumberTheory`ContinuedFractions` *)
  108.  
  109. (*:Limitations: none known. *)
  110.  
  111. (*:Examples:
  112.  
  113. ContinuedFraction[ N[Pi], 4 ]
  114.  
  115. Normal[ ContinuedFractionForm[{1, 1, 2, 2}] ]
  116. *)
  117.