Some Continued Fraction Algorithms - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html) +--- Forum: General Forum (/forum-4.html) +--- Thread: Some Continued Fraction Algorithms (/thread-13840.html) |
Some Continued Fraction Algorithms - ttw - 10-23-2019 04:06 AM Four algorithms I've been using. Convert partial quotients to fractional and vice versa. Also convert surd to partial quotients and vice versa. These do work with big integers. I tested the fraction and partial quotients with numbers bigger than 2^64. ->Frac Input: list of integers representing the partial quotients of a continued fraction. HP50g Output: a list of the continued fractions exhibited for each partial quotient. (All fractional, no improper fractions.) Input is echoed. « DUP SIZE LASTARG { 0 1} {1 0} -> N A P Q « 1 N FOR I A I GET DUP P HEAD * P 2 GET + P + 'P' STO Q HEAD * Q 2 GET + Q + 'Q' STO NEXT P Q / 1 N SUB » » This just uses the usual recursion. Making P and Q a single list with two elements runs slower (even though one can do the arithmetic using list operations.) Frac-> Input: fraction. (Like '30/101'). Input is echoed. Trivial to change to pair of numbers or even to sort the pair. Output: list of partial quotients; first entry is the integer part. « DUP FXND -> A B « {} WHILE B REPEAT A B IQUOT + A B MOD B 'A' STO 'B' STO END » » Repeated Euclid's Algorith. IDIV2 is very slow. IREMAINDER slower than MOD. Also / and FLOOR slower than IQUOT. Perhaps different size numbers would change the results. Lagrange's therorems show that every periodic continued fraction is a quadratic irrational (old name was "surd"). Also every quadratic irrationsl (A+Sqrt(D))/C yields a continued fraction. The next to programs compute these results. ->Surd Input is a list of partial quotients making up the periodic part of the quadratic irrational. (Doing the aperiodic stuff was too hard and I didn't need it.) Output is a number of the form: (A+Sqrt(B))/C. « DUPDUP SIZE IF 1 == THEN DUP + END ->Frac 1 2 SUB FXND EVAL -> P2 Q2 P1 Q1 « P1 Q1 - DUP SQ Q1 P2 * 4 * + SQRT + 2 Q1 * / EXPAND » » Surd-> Input: quadratic irrational ( A + Sqrt(D) / B ) Output: list of partial quotients. First entry is integer part. Then, list of aperiodic part (no list if none). Followed by periodic part. « DUP 0 {} 0 0 -> X A S N M « EXPAND DUP FLOOR DUP 1 ->LIST 'A' STO - EXPAND DUP DO 'S' STO+ INV DUP FLOOR DUP 'A' STO+ - EXPAND DUPDUP S SWAP POS DUP 'N' STO UNTIL END DROP2 X A REVLIST HEAD A SIZE DUP 'M' STO N - IF 1 > THEN A N 1 + M 1 - SUB REVLIST 1 -> LIST + N END A 1 N SUB REVLIST + » » A bit slow as the program checks for a cycle completion at each new generated partial quotient. Algorithm is to iterate the Gauss Map: X = FRAC( 1 / X ). There are some faster stuff but they just get partial quotients ignoring the aperiodic stuff. RE: Some Continued Fraction Algorithms - ttw - 10-31-2019 02:51 AM Another number theory program using HP50g's big integer capability. Given a prime, M, and a non-zero number, A, "pRoot" determines if A is a quadratic residue for M, and if not, if A is a primitive root of M. A being a primitive root means A can be used a multiplier for a Linear Congruential Generator, with M as modulus (with no additive term.) While these generators are not cryptgraphcally secure, the lower order bits have much less structure than the low order bits of M being a power of two. « MIN LASTARG MAX DUP 1 - 0 0 -> A M M1 D N « 1 SF M MODSTO M1 DIVIS REVLIST DUP 'D' STO SIZE 'N' STO A D 2 GET POWMOD IF 1 == THEN 1 CF M A "Residue" END IF 1 FS? THEN 3 N 1 - FOR I A D I GET POWMOD IF 1 == THEN N 3 + 'I' STO 1 CF END NEXT M A IF 1 FS? THEN "Primitive Root" ELSE "Non Root" END »» If A is a quadratic residue and M = 1 Mod 4, then these can be used to construct a quadratic random number generator that has nice properties (though I haven't really tried to prove anything.) |