Some Continued Fraction Algorithms

[ttw]

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.

[ttw]

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.)