Post Reply 
Some Continued Fraction Algorithms
10-23-2019, 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.
Find all posts by this user
Quote this message in a reply
10-31-2019, 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 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.)
Find all posts by this user
Quote this message in a reply
Post Reply 




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