Some Continued Fraction Algorithms

10232019, 04:06 AM
Post: #1




Some Continued Fraction Algorithms
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. 

10312019, 02:51 AM
Post: #2




RE: Some Continued Fraction Algorithms
Another number theory program using HP50g's big integer capability.
Given a prime, M, and a nonzero 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.) 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)