'********************************************************************************
'Subject: Toy version of R. W. Gosper's continued fraction arithmetic.
'Ref. : HAKMEM, A.I. Lab Memo #239, M.I.T., 1972
'Quote : 'Sad to say, the integer variables in these algorithms do not
' usually shrink on outputs as much as they grow on inputs.'
' If they did, we would indeed have 'unlimited significance arithmetic
' without multiprecision multiplication or division.' Utilizing 64-bit
' floating point, significance more than doubles to roughly 40 decimals,
' but transcendental functions of nice rational arguments do give a much
' longer range of correct terms.
'Author : Sjoerd.J.Schaper - vspickelen
'Date : 10-25-2009
'Code : PDS, VBdos and FreeBasic
'********************************************************************************
'This program is copyright (c) 2009 by the author. It is made available as is,
'and no warranty - about the program, its performance, or its conformity to any
'specification - is given or implied. It may be used, modified, and distributed
'freely, as long as the original author is credited as part of the final work.
'********************************************************************************
REM $LANG: "qb"
DEFINT A-Z
CONST Verb = 0 ' verbosity 0,1,2
CONST Eps = 1D-40 ' set Eps = 0 for full (noisy) Cf
CONST Cfx = 102 ' fixed Cf length
TYPE cfa ' Cf arithmetic structure:
u(3) AS DOUBLE ' GL(2,Z) matrices
v(3) AS DOUBLE
x(Cfx) AS INTEGER ' continued fractions
y(Cfx) AS INTEGER
flag AS INTEGER ' break signal
END TYPE
DECLARE FUNCTION nextcfd (n AS cfa)
'compute next Cf digit z(x, y)
DECLARE FUNCTION digin (n AS cfa, p)
'update tensor for input x(p=0) or y(p=1)
DECLARE SUB Cget (q, c())
'read q from Cf_c
DECLARE SUB Digout (n AS cfa, q)
'update tensor for output digit q
DECLARE SUB Cset (c(), q)
'write q in Cf_c
DECLARE SUB Cdo (c(), n AS cfa)
'safe Cf arithmetic loop
DECLARE SUB Cdup (a(), b())
'copy Cf_a to Cf_b
DECLARE FUNCTION comp (t, a(), b())
'compare Cf's; mismatch pointer t
DECLARE SUB ReadCf (c(), w AS STRING)
'input continued fraction "q0.q1.q2..."
DECLARE SUB OutCf (s AS STRING, c())
'print Cf_c and convergent
DECLARE SUB Trans (n AS cfa)
'transpose matrices
DECLARE SUB Turn (n AS cfa)
'turn coefficient cube
DECLARE SUB Tmul (tr() AS DOUBLE, n AS cfa)
'premultiply trafo * tensor
DECLARE SUB Tset (n AS cfa, a, b, c, d, e, f, g, h)
'assign (a b c d) / (e f g h) to tensor,
'reset read / write pointers and flag
DECLARE SUB Lnorm (n AS cfa)
'normalize (c d) / (g h) registers
'when computing generalized Cf's
DECLARE FUNCTION homocfd (n AS cfa)
'compute homographic digit z(y)
DECLARE SUB OutFrm (n AS cfa)
'print bilinear fractional form
'powers and roots
DECLARE SUB Cpwr (a(), k, sw)
'Cf_a:= Cf_a ^ k, integer k
'set switch for tangent sum
DECLARE SUB Csqrt (a())
'Cf_a:= Cf_a ^ (1 / 2)
DECLARE SUB Cubrt (a())
'Cf_a:= Cf_a ^ (1 / 3)
'reciprocation and negation
DECLARE SUB Cinv (a())
'Cf_a:= 1 / Cf_a
DECLARE SUB Cneg (a())
'Cf_a:= -Cf_a
'streaming functions with rational argument,
'non-reentrant, to be initialized with t = 0
DECLARE FUNCTION cquad (p, q, r, t)
'solve p*x^2 + q*x - r = 0; p <> 0
DECLARE FUNCTION clog (x, y, t)
'x / y > 0
DECLARE FUNCTION cexp (x, y, t)
DECLARE FUNCTION csin (x, y, t)
DECLARE FUNCTION ccos (x, y, t)
DECLARE FUNCTION ctan (x, y, t)
DECLARE FUNCTION ctanh (x, y, t)
DECLARE FUNCTION catan (x, y, t)
DECLARE FUNCTION catanh (x, y, t)
'Abs(x / y) < 1
DECLARE FUNCTION cpi (t)
'special purpose adaptations w/ integer argument:
DECLARE FUNCTION csin2 (k, t)
'sin(big argument)
DECLARE FUNCTION ctanh2 (k, t)
'tanh(square root(k))
'
DECLARE FUNCTION gencf (n AS cfa, x2 AS DOUBLE, y, t, sw)
'generalized Cf for transcendental functions,
'evaluates tanh(sw=0) or atan(sw=1)
DECLARE FUNCTION gcd (x, y)
'highest common divisor
DECLARE SUB Halve (x, y)
'x / y:= x / 2y
'some demos
DECLARE SUB Allbasic ()
'try basic functions for fi and e
DECLARE SUB Alltrans (x, y)
'run all transcendental functions
DECLARE SUB Appxtest ()
'Knuth's exercise and Gosper's Appendix 2
DECLARE SUB Gospernum ()
'Gosper's numbers in items 101A and B
CONST Inf = &H7FFF ' infinity
CONST Pend = NOT Inf, Done = -Inf ' markers
STACK 4096 ' for DOS-Basic only, Rem out in FB
Tim# = TIMER: CLS
Allbasic ' main
Alltrans -1, 4
Appxtest
Gospernum
eind:
PRINT "Timer:"; CSNG(TIMER - Tim#); "s"
SYSTEM
SUB Allbasic
DIM n AS cfa
DIM a(Cfx), b(Cfx), c(Cfx)
a(1) = 2: t = 0
b(1) = 2: t2 = 0
DO
Cset a(), cquad(1, -1, 1, t) ' fi
Cset b(), cexp(1, 1, t2) ' e
LOOP UNTIL b(b(1)) = Done
OutCf "x", a(): OutCf "y", b()
Cdup a(), n.x(): Cdup b(), n.y()
Tset n, 0, 1, 1, 0, 0, 0, 0, 1 ' sum
Cdo c(), n: OutCf "x+y", c()
Cdup a(), n.x(): Cdup b(), n.y()
Tset n, 0, 1, -1, 0, 0, 0, 0, 1 ' difference
Cdo c(), n: OutCf "x-y", c()
Cdup a(), n.x(): Cdup b(), n.y()
Tset n, 1, 0, 0, 0, 0, 0, 0, 1 ' product
Cdo c(), n: OutCf "x*y", c()
Cdup a(), n.x(): Cdup b(), n.y()
Tset n, 0, 1, 0, 0, 0, 0, 1, 0 ' quotient
Cdo c(), n: OutCf "x/y", c()
Cdup a(), c(): Cpwr c(), 10, 0
OutCf "x^10", c()
Cdup b(), c(): Cpwr c(), 10, 0 ' power
OutCf "y^10", c()
Cdup a(), c(): Csqrt c()
OutCf "x^1/2", c()
Cdup b(), c(): Csqrt c() ' roots
OutCf "y^1/2", c()
Cdup a(), c(): Cubrt c()
OutCf "x^1/3", c()
Cdup b(), c(): Cubrt c()
OutCf "y^1/3", c()
c(1) = 2: t = 0
DO: Cset c(), cquad(349, -500, 0, t) ' rational to Cf
LOOP UNTIL c(c(1)) = Done
OutCf "349x = 500:", c()
c(1) = 2: t = 0
DO: Cset c(), cquad(5, 18, 15, t) ' solve quadratic
LOOP UNTIL c(c(1)) = Done
OutCf "5xx + 18x - 15 = 0:", c()
Tset n, 1, 0, 0, 0, 0, 0, 0, 1
c(1) = 2: t = 0
DO
d = cquad(1, 0, 2, t)
Cset n.x(), d: Cset n.y(), d ' accuracy test
IF t = 1 THEN n.x(2) = 0: n.y(2) = 2
Cset c(), nextcfd(n)
LOOP UNTIL c(c(1)) = Done
OutCf "(sqrt(2) - 1)(sqrt(2) + 1):", c()
ReadCf a(), "-2.1.7.9.1.6171.1.15.2.1.3.10000"
OutCf "x ", a()
Cpwr a(), 81, 0: OutCf "y = x^81", a()
FOR t = 1 TO 4: Cubrt a(): NEXT ' Newton stress-test
OutCf "y^1/81", a()
END SUB
SUB Alltrans (x, y)
DIM c(Cfx), q AS STRING
FOR tl = 1 TO 8
PRINT
q = "(" + LTRIM$(STR$(x)) + "/"
q = q + LTRIM$(STR$(y)) + ")"
c(1) = 2: t = 0
DO: Cset c(), clog(x, y, t)
LOOP UNTIL c(c(1)) = Done
OutCf "log" + q, c()
c(1) = 2: t = 0
DO: Cset c(), cexp(x, y, t)
LOOP UNTIL c(c(1)) = Done
OutCf "exp" + q, c()
c(1) = 2: t = 0
DO: Cset c(), csin(x, y, t)
LOOP UNTIL c(c(1)) = Done
OutCf "sin" + q, c()
c(1) = 2: t = 0
DO: Cset c(), ccos(x, y, t)
LOOP UNTIL c(c(1)) = Done
OutCf "cos" + q, c()
c(1) = 2: t = 0
DO: Cset c(), ctan(x, y, t)
LOOP UNTIL c(c(1)) = Done
OutCf "tan" + q, c()
c(1) = 2: t = 0
DO: Cset c(), ctanh(x, y, t)
LOOP UNTIL c(c(1)) = Done
OutCf "tanh" + q, c()
c(1) = 2: t = 0
DO: Cset c(), catan(x, y, t)
LOOP UNTIL c(c(1)) = Done
OutCf "atan" + q, c()
c(1) = 2: t = 0
DO: Cset c(), catanh(x, y, t)
LOOP UNTIL c(c(1)) = Done
OutCf "atanh" + q, c()
x = x * 4 + y: y = y * 4 ' step 1 / 4
d = gcd(x, y)
IF d > 1 THEN x = x \ d: y = y \ d
NEXT tl
c(1) = 2: t = 0
DO: Cset c(), cpi(t)
LOOP UNTIL c(c(1)) = Done
OutCf "pi", c()
END SUB
SUB Appxtest
DIM m AS cfa, n AS cfa, c(Cfx)
PRINT : n.x(2) = Done
ReadCf n.y(), "-1.5.1.1.1.2.1.2" ' Knuth's exercise 4.5.3.15
Tset n, 97, 39, 0, 0, -62, -25, 0, 0
Cdo c(), n: OutCf "Knuth:", c()
PRINT : n.x(2) = Done ' From Gosper's Appendix 2:
n.y(2) = 0: n.y(3) = Done
Tset n, 70, 29, 0, 0, 12, 5, 0, 0 ' drain out Cf-matrix
Cdo c(), n: OutCf "terms:", c()
Tset n, 0, 2, 0, 0, -1, 3, 0, 0
n.x(2) = Done: t1 = 0
DO
Cset n.y(), cquad(1, 0, 2, t1) ' sqrt(2)
Cset n.x(), nextcfd(n)
LOOP UNTIL n.x(n.x(1)) = Done
Cinv n.x(): Cinv n.x() ' strip zeroes
OutCf "2 / (3 - sqrt(2)):", n.x()
Tset m, 1, -1, 0, 0, 1, 1, 0, 0
Tset n, -4, 4, 0, 0, 1, 1, 0, 0
m.x(2) = Done: n.x(2) = Done
t1 = 0
DO
Cset m.y(), cexp(1, 1, t1)
Cset n.y(), nextcfd(m) ' y:= (e - 1) / (e + 1)
Cset n.x(), nextcfd(n) ' 4 * (1 - y) / (1 + y)
LOOP UNTIL n.x(n.x(1)) = Done
OutCf "tanh(1/2):", n.y()
OutCf "4 / e:", n.x()
Tset m, 1, 0, 0, 1, 1, 0, 0, -1 ' (xy + 1) / (xy - 1)
Tset n, 2, 1, 0, 0, 1, 0, 1, 0
c(1) = 2: t1 = 0: t2 = 0
DO
d = cexp(1, 1, t1)
Cset m.x(), d: Cset m.y(), d
Cset n.x(), nextcfd(m) ' x = coth(1)
Cset n.y(), cquad(1, 0, 6, t2) ' y = sqrt(6)
Cset c(), nextcfd(n)
LOOP UNTIL c(c(1)) = Done
OutCf "x", n.x(): OutCf "y", n.y()
OutCf "(2xy + x) / (xy + y):", c()
c(1) = 2: t1 = 0
DO: Cset c(), cquad(10, 0, 17, t1) ' rational sqrt warmup
LOOP UNTIL c(c(1)) = Done
OutCf "sqrt(17/10)", c()
Cset n.x(), Done
OutCf "coth(1)", n.x() ' compute coth(1/2) from coth(1)
n.y(2) = 2: n.y(3) = Done ' initial guess Y0 = 2
t1 = 3
DO
Tset n, 1, 0, 0, -1, 0, -1, 1, 0 ' operations
Tset m, 0, 1, 1, 0, 0, 0, 0, 2
Cdup n.y(), m.y(): c(1) = 2
DO ' Newton iteration
Cset m.x(), nextcfd(n) ' mx:= (xy - 1) / (y - x)
Cset c(), nextcfd(m): d = c(1) 'average (mx + y) / 2
LOOP UNTIL c(d) = Done OR d = t1
Cset c(), Done: Cdup c(), n.y() ' feed back
PRINT "t ="; t1 - 2
OutCf "Y", c() ' approximate coth(1/2)
t1 = (t1 - 1) * 2
LOOP UNTIL t1 > Cfx
ReadCf c(), "5.1": Csqrt c() ' sqrt with Cf-argument
OutCf "sqrt(6)", c()
END SUB
FUNCTION catan (x, y, t)
STATIC n AS cfa, x2 AS DOUBLE
catan = Done
IF t = 0 THEN
IF y = 0 THEN EXIT FUNCTION
Tset n, 0, 0, x, 0, 0, 0, y, 1
x2 = CDBL(x) * x
END IF
catan = gencf(n, x2, y, t, 1)
END FUNCTION
FUNCTION catanh (x, y, t)
STATIC n AS cfa, x2 AS DOUBLE
catanh = Done
IF t = 0 THEN
IF y = 0 THEN EXIT FUNCTION
IF 1 - ABS(x / y) < 1E-14 THEN ' undefined
EXIT FUNCTION
END IF
Tset n, 0, 0, x, 0, 0, 0, y, 1
x2 = -CDBL(x) * x
END IF
catanh = gencf(n, x2, y, t, 1)
END FUNCTION
FUNCTION ccos (x, y, t)
STATIC m AS cfa, n AS cfa, a2 AS DOUBLE, b
ccos = Done
IF t = 0 THEN
IF y = 0 THEN EXIT FUNCTION
a = x: b = y: Halve a, b ' A&S 4.3.25
Tset m, 0, 0, a, 0, 0, 0, b, 1 ' x = y = tan a / 2b
Tset n, -1, 0, 0, 1, 1, 0, 0, 1 ' (1 - xy) / (1 + xy)
a2 = -CDBL(a) * a
END IF
DO
d = gencf(m, a2, b, t, 0)
Cset n.x(), d: Cset n.y(), d
d = nextcfd(n)
LOOP UNTIL d > Pend
ccos = d
END FUNCTION
SUB Cdo (c(), n AS cfa)
c(1) = 2 ' recipient
DO
Cset c(), nextcfd(n)
LOOP UNTIL c(c(1)) = Done OR n.flag < 0
IF n.flag < 0 THEN
Cset c(), Done: PRINT " break" ' handle .999 problem
END IF
END SUB
SUB Cdup (a(), b())
FOR r = 0 TO Cfx: b(r) = a(r): NEXT
END SUB
FUNCTION cexp (x, y, t)
STATIC n AS cfa, a2 AS DOUBLE, b
cexp = Done
IF t = 0 THEN
IF y = 0 THEN EXIT FUNCTION ' 2x x2 x2 x2
a = x: b = y ' 1 + ------ --- --- --- ...
Halve a, b: a2 = CDBL(a) * a ' (y-x)+ 3y+ 5y+ 7y+
Tset n, 0, 0, b + a, 1, 0, 0, b - a, 1
END IF
cexp = gencf(n, a2, b, t, 0)
END FUNCTION
SUB Cget (q, c())
r = c(0): q = c(r) ' c(Cfx)==Done
IF q > Done THEN c(0) = r + 1 ' advance read ptr
END SUB
SUB Cinv (a())
IF a(2) <= Done THEN EXIT SUB
IF a(2) = 0 THEN
FOR t = 2 TO Cfx - 1 ' remove initial 0
a(t) = a(t + 1)
NEXT t
ELSE
sw = a(2) < 0
IF sw THEN Cneg a()
FOR t = Cfx - 1 TO 3 STEP -1
a(t) = a(t - 1)
NEXT t: a(2) = 0 ' prefix 0
a(Cfx) = Done
IF sw THEN Cneg a()
END IF
END SUB
FUNCTION clog (x, y, t)
STATIC m AS cfa, n AS cfa, s AS cfa
STATIC k, a2 AS DOUBLE, b, t2
clog = Done
IF t = 0 THEN
IF y = 0 THEN EXIT FUNCTION
IF x / y < 1E-14 THEN EXIT FUNCTION ' undefined
'
a = x: b = y: k = 0
WHILE CDBL(a) / b > 2
Halve a, b: k = k + 1 ' scaling factor 2^k
WEND
WHILE CDBL(b) / a > 2
Halve b, a: k = k - 1
WEND
t2 = a - b: b = b + a: a = t2
IF Verb > 1 THEN
PRINT "arg"; a; "/"; b; " k ="; k
END IF
Tset m, 0, 0, a, 0, 0, 0, b, 1 ' 2*atanh((x-y) / (x+y))
Tset n, 0, 0, k, 0, 0, 0, 3, 1 ' + k*log(2)
Tset s, 0, 2, 2, 0, 0, 0, 0, 1
IF k = 0 THEN m.u(2) = a * 2
a2 = -CDBL(a) * a: t2 = 0
END IF
DO
IF k THEN
Cset s.x(), gencf(m, a2, b, t, 1)
Cset s.y(), gencf(n, -1, 3, t2, 1)
d = nextcfd(s)
IF d > Pend THEN EXIT DO
ELSE
d = gencf(m, a2, b, t, 1): EXIT DO
END IF
LOOP
clog = d
END FUNCTION
SUB Cneg (a())
DIM n AS cfa
Cdup a(), n.y()
Tset n, -1, 0, -1, 0, 0, 1, 0, 1
Cdo a(), n
END SUB
FUNCTION comp (r, a(), b())
FOR t = 2 TO Cfx - 1 ' first discrepant terms
IF CLNG(a(t)) - b(t) THEN EXIT FOR
IF a(t) <= Done THEN EXIT FOR ' both finite
NEXT t
r = a(t): IF r <= Done THEN r = Inf
s = b(t): IF s <= Done THEN s = Inf
IF t AND 1 THEN SWAP r, s ' invert decision
comp = SGN(CLNG(r) - s): r = t
END FUNCTION
FUNCTION cpi (t)
STATIC m AS cfa, n AS cfa, s AS cfa, t2
IF t = 0 THEN
Tset m, 0, 0, 2, 0, 0, 0, 7, 1 ' Charles Hutton's 1776 formula
Tset n, 0, 0, 4, 0, 0, 0, 3, 1 ' 4*atan(1/7) + 8*atan(1/3)
Tset s, 0, 2, 2, 0, 0, 0, 0, 1 ' 6E-46
t2 = 0
END IF
DO
Cset s.x(), gencf(m, 1, 7, t, 1)
Cset s.y(), gencf(n, 1, 3, t2, 1)
d = nextcfd(s)
LOOP UNTIL d > Pend
cpi = d
END FUNCTION
SUB Cpwr (a(), k, sw)
DIM b(Cfx), n AS cfa
IF k = 0 THEN
a(2) = 1: IF sw THEN a(2) = 0
a(3) = Done: EXIT SUB
END IF
'
Cdup a(), b() ' b:= initial a
c = ABS(k): c2 = 1 + Inf \ 2
DO UNTIL c AND c2
c2 = c2 \ 2: LOOP ' highest set bit(k)
c2 = c2 \ 2
WHILE c2 ' L=>R binary building
Cdup a(), n.x(): Cdup a(), n.y() ' x = y = a
IF sw THEN
Tset n, 0, 1, 1, 0, -1, 0, 0, 1 ' tan(2a) 4.3.26
ELSE
Tset n, 1, 0, 0, 0, 0, 0, 0, 1 ' square
END IF
Cdo a(), n
IF c AND c2 THEN
Cdup a(), n.x(): Cdup b(), n.y() ' x = a, y = b
IF sw THEN
Tset n, 0, 1, 1, 0, -1, 0, 0, 1 'tan(a + b) 4.3.18
ELSE
Tset n, 1, 0, 0, 0, 0, 0, 0, 1 ' mult
END IF
Cdo a(), n
END IF
c2 = c2 \ 2
WEND
END SUB
FUNCTION cquad (a, b, c, t)
STATIC p() AS DOUBLE, q() AS DOUBLE, rD, i, j
DIM disc AS DOUBLE
cquad = Done
IF t THEN
SWAP i, j: t = t + 1
IF q(j) = 0 THEN EXIT FUNCTION
d = INT((rD + p(j)) / q(j)) ' next digit
p(i) = d * q(j) - p(j)
q(i) = q(i) + d * (p(j) - p(i)) ' recurrence
ELSE
disc = CDBL(b) * b + CDBL(a) * c * 4
IF disc < 0 OR a = 0 THEN EXIT FUNCTION
'
REDIM p(1) AS DOUBLE, q(1) AS DOUBLE
i = 0: j = 1: t = 1
rD = INT(SQR(disc))
p(j) = b: q(j) = a * 2 ' initialize
d = INT((rD - p(j)) / q(j)) ' first digit
p(i) = d * q(j) + p(j)
q(i) = (disc - p(i) * p(i)) / q(j)
END IF
cquad = d
END FUNCTION
SUB Cset (c(), q)
r = c(1)
IF r < Cfx THEN
IF q > Done THEN c(1) = r + 1 ' advance write ptr
c(c(1)) = Pend: c(r) = q
ELSE
c(Cfx) = Done
END IF
END SUB
FUNCTION csin (x, y, t)
STATIC m AS cfa, n AS cfa, a2 AS DOUBLE, b
csin = Done
IF t = 0 THEN
IF y = 0 THEN EXIT FUNCTION
a = x: b = y: Halve a, b ' 4.3.24
Tset m, 0, 0, a, 0, 0, 0, b, 1 ' x = y = tan a / 2b
Tset n, 0, 1, 1, 0, 1, 0, 0, 1 ' (x + y) / (xy + 1)
a2 = -CDBL(a) * a
END IF
DO
d = gencf(m, a2, b, t, 0)
Cset n.x(), d: Cset n.y(), d
d = nextcfd(n)
LOOP UNTIL d > Pend
csin = d
END FUNCTION
FUNCTION csin2 (k, t)
STATIC n AS cfa
IF t = 0 THEN
Tset n, 0, 0, 1, 0, 0, 0, 2, 1 ' tan 1 / 2
n.x(2) = Done
DO
Cset n.x(), gencf(n, -1, 2, t, 0)
LOOP UNTIL n.x(n.x(1)) = Done
Cpwr n.x(), k, -1 ' peasant's mult tan k / 2
Cdup n.x(), n.y(): t = 0
Tset n, 0, 1, 1, 0, 1, 0, 0, 1 ' sin
END IF
t = t + 1
csin2 = nextcfd(n)
END FUNCTION
SUB Csqrt (a())
DIM m AS cfa, n AS cfa
IF a(2) < 0 THEN
PRINT " illegal argument Csqrt"
a(2) = Done: EXIT SUB
END IF
sw = a(2) = 0
IF sw THEN Cinv a()
'
Cdup a(), m.x() ' radicand R > 1
a(2) = INT(SQR(a(2) - (a(3) = 1))) ' initial guess X
a(3) = Done
IF Verb THEN
OutCf "R", m.x()
OutCf "X0", a()
END IF
'
t = 3
DO
Cdup a(), m.y(): Cdup a(), n.x() ' feed back
Tset m, 0, 1, 0, 0, 0, 0, 1, 0
Tset n, 0, 1, 1, 0, 0, 0, 0, 2
a(1) = 2 ' Newton iteration
DO
Cset n.y(), nextcfd(m) ' y:= R / X
Cset a(), nextcfd(n): d = a(1) 'average (X + y) / 2
LOOP UNTIL a(d) = Done OR d = t
Cset a(), Done
IF Verb THEN
PRINT "t ="; t - 2
OutCf " q", n.y()
OutCf " X", a()
END IF
d = 1 + t \ 2
t = (t - 1) * 2
LOOP WHILE d < Cfx
IF sw THEN Cinv a()
END SUB
FUNCTION ctan (x, y, t)
STATIC n AS cfa, x2 AS DOUBLE
ctan = Done
IF t = 0 THEN
IF y = 0 THEN EXIT FUNCTION
Tset n, 0, 0, x, 0, 0, 0, y, 1
x2 = -CDBL(x) * x
END IF
ctan = gencf(n, x2, y, t, 0)
END FUNCTION
FUNCTION ctanh (x, y, t)
STATIC n AS cfa, x2 AS DOUBLE
ctanh = Done
IF t = 0 THEN
IF y = 0 THEN EXIT FUNCTION
Tset n, 0, 0, x, 0, 0, 0, y, 1
x2 = CDBL(x) * x
END IF
ctanh = gencf(n, x2, y, t, 0)
END FUNCTION
FUNCTION ctanh2 (k, t)
STATIC m AS cfa, n AS cfa, t2
ctanh2 = Done
IF k = 0 THEN
IF t = 0 THEN t = 1: ctanh2 = 0
EXIT FUNCTION
END IF
'
IF t = 0 THEN
Tset m, 0, 0, 1, 1, 0, 0, 1, 0 ' sqrt(k) k k k
Tset n, 0, 1, 0, 0, 0, 0, 1, 0 ' ------- -- -- -- ...
t2 = 0 ' 1+ 3+ 5+ 7+
END IF
DO
Cset n.x(), cquad(1, 0, k, t2) ' sqrt(k)
Cset n.y(), gencf(m, CDBL(k), 1, t, 0) ' / sqrt(k) * coth(sqrt(k))
d = nextcfd(n)
LOOP UNTIL d > Pend
ctanh2 = d
END FUNCTION
SUB Cubrt (a())
DIM m AS cfa, n AS cfa, s AS cfa
IF a(2) <= Done THEN
PRINT " illegal argument Cubrt"
a(2) = Done: EXIT SUB
END IF
fl = a(2) < 0
IF fl THEN Cneg a()
sw = a(2) = 0
IF sw THEN Cinv a()
'
Cdup a(), m.x() ' radicand R > 1
a(2) = INT((a(2) - (a(3) = 1)) ^ (1 / 3))
a(3) = Done ' initial guess X
IF Verb THEN
OutCf "R", m.x()
OutCf "X0", a()
END IF
'
t = 3
DO
Cdup a(), n.x() ' feed back
Cdup a(), s.x(): Cdup a(), s.y()
Tset s, 2, 0, 0, 0, 0, 0, 0, 1
Tset m, 0, 1, 0, 0, 0, 0, 1, 0
Tset n, 0, 2, 2, 0, 0, 0, 0, 3
a(1) = 2 ' Newton iteration
DO
Cset m.y(), nextcfd(s)
Cset n.y(), nextcfd(m) ' y:= R / 2X^2
Cset a(), nextcfd(n): d = a(1) 'average 2(X + y) / 3
LOOP UNTIL a(d) = Done OR d = t
Cset a(), Done
IF Verb THEN
PRINT "t ="; t - 2
OutCf " q", n.y()
OutCf " X", a()
END IF
d = 1 + t \ 2
t = (t - 1) * 2
LOOP WHILE d < Cfx
IF sw THEN Cinv a()
IF fl THEN Cneg a()
END SUB
FUNCTION digin (n AS cfa, p)
DIM tr(3) AS DOUBLE ' trafo
digin = 0
IF p THEN
Cget q, n.y() ' Cf digit
IF Verb > 1 THEN
PRINT " <-y" + LTRIM$(STR$(n.y(0) - 2)); q
END IF
ELSE
Cget q, n.x()
IF Verb > 1 THEN
PRINT " <-x" + LTRIM$(STR$(n.x(0) - 2)); q
END IF
END IF
'
IF q > Done THEN
tr(0) = q: tr(1) = 1
tr(2) = 1: tr(3) = 0
ELSE
IF q = Pend THEN EXIT FUNCTION ' wait for next digit
tr(0) = 1: tr(1) = 0
tr(2) = 1: tr(3) = 0
END IF
IF p THEN Trans n ' transpose for y
Tmul tr(), n ' premultiply trafo
IF p THEN Trans n ' restore
digin = -1
END FUNCTION
SUB Digout (n AS cfa, q)
DIM tr(3) AS DOUBLE ' trafo
n.flag = n.flag + (q = Inf) ' Inf count
'
tr(0) = 0: tr(1) = 1 ' inverse
tr(2) = 1: tr(3) = -q
Turn n: Tmul tr(), n: Turn n ' tilted cube
END SUB
FUNCTION gcd (x, y)
a = x: b = y
WHILE b
c = a MOD b
a = b: b = c
WEND
gcd = ABS(a)
END FUNCTION
FUNCTION gencf (n AS cfa, x2 AS DOUBLE, y, t, sw)
DIM tr(3) AS DOUBLE
DO
IF t THEN
d = homocfd(n)
IF d > Pend THEN EXIT DO
END IF
t = t + 1 ' x x2 4x2 9x2 16x2
IF sw THEN ' -- --- --- --- ----... atan 4.4.43
tr(1) = CDBL(t) * t * x2 ' y+ 3y+ 5y+ 7y+ 9y+
ELSE
tr(1) = x2 ' x x2 x2 x2 x2
END IF ' -- --- --- --- ---... tanh 4.5.70
tr(0) = CDBL(2 * t + 1) * y ' y+ 3y+ 5y+ 7y+ 9y+
tr(2) = 1: tr(3) = 0
IF Verb > 1 THEN
PRINT " <-p,q_" + LTRIM$(STR$(t)); tr(0); tr(1)
END IF
Trans n: Tmul tr(), n: Trans n
LOOP
gencf = d
END FUNCTION
SUB Gospernum
DIM c(Cfx), m AS cfa, n AS cfa
DIM p AS cfa, q AS cfa, s AS cfa
PRINT
c(1) = 2: t1 = 0 ' 3 full-precision numbers
DO: Cset c(), cquad(1, 0, 105, t1)
LOOP UNTIL c(c(1)) = Done
OutCf "sqrt(105)", c()
c(1) = 2: t1 = 0
DO: Cset c(), ctanh(1, 69, t1)
LOOP UNTIL c(c(1)) = Done
Cinv c(): OutCf "coth(1/69)", c()
Tset n, 4, -2, 0, 0, 1, -1, 0, 0
n.x(2) = Done: t1 = 0
DO
Cset n.y(), cexp(2, 3, t1)
Cset n.x(), nextcfd(n)
LOOP UNTIL n.x(n.x(1)) = Done
OutCf "exp(2/3)", n.y()
OutCf "(4*exp(2/3) - 2) / (exp(2/3) - 1):", n.x()
PRINT
Tset m, 0, 0, 0, 3, 1, 0, 0, 0
Tset n, 0, 1, 1, 0, 0, 0, 0, 1
Tset s, 0, 1, -1, 0, 0, 0, 0, 1
Tset p, 1, 0, 0, 0, 0, 0, 0, 1
Tset q, 0, 1, 0, 0, 0, 0, 1, 0
t1 = 0: t2 = 0: t3 = 0: t4 = 0
c(1) = 2
DO: d = cpi(t1)
Cset m.x(), d: Cset m.y(), d
Cset n.x(), nextcfd(m)
Cset n.y(), cexp(1, 1, t2)
Cset q.x(), nextcfd(n) ' x = 3 / pi^2 + e
Cset s.x(), ctanh2(5, t3)
Cset s.y(), csin2(69, t4)
d = nextcfd(s) ' y = tanh(sqrt(5)) - sin(69)
Cset p.x(), d: Cset p.y(), d
Cset q.y(), nextcfd(p) ' y^2
Cset c(), nextcfd(q) ' quotient
LOOP UNTIL c(c(1)) = Done
OutCf "3/pi^2", n.x() ' 6E-44
OutCf "exp(1)", n.y() ' 2E-115
OutCf "sum:", q.x()
OutCf "tanh(sqrt(5))", s.x() ' 2E-44
OutCf "sin(69)", s.y() ' 3E-40
OutCf "diff^2:", q.y()
Csqrt c() ' root 4E-40
OutCf "sqrt(3/pi^2 + e) / (tanh(sqrt(5)) - sin(69))", c()
END SUB
SUB Halve (x, y)
IF (x AND 1) = 0 THEN
x = x \ 2
ELSE
y = y * 2
END IF
END SUB
FUNCTION homocfd (n AS cfa)
DIM d(3), q AS DOUBLE
homocfd = Pend
FOR r = 2 TO 3
IF n.v(r) THEN
q = INT(n.u(r) / n.v(r)) ' both z-values
IF ABS(q) > Inf THEN q = SGN(q) * Inf
d(r) = q
ELSE
d(r) = Done
END IF
NEXT r
fl = d(2) = d(3)
'
IF Verb > 1 THEN
OutFrm n: PRINT " ";
FOR r = 2 TO 3: PRINT d(r); : NEXT
IF fl THEN PRINT " ->out"
END IF
'
IF fl THEN ' flat z-range
IF d(2) > Done THEN
Digout n, d(2): Lnorm n
END IF
homocfd = d(2)
END IF
END FUNCTION
SUB Lnorm (n AS cfa)
DIM d AS DOUBLE
DO
d = n.u(2) * n.v(3) - n.v(2) * n.u(3)
IF ABS(d) < 4 THEN EXIT DO
FOR r = 2 TO 3
n.u(r) = n.u(r) * .5#
n.v(r) = n.v(r) * .5#
NEXT r
LOOP
IF ABS(d) < 1D-19 THEN EXIT SUB
DO
IF ABS(d) > .25 THEN EXIT DO
FOR r = 2 TO 3
n.u(r) = n.u(r) * 2
n.v(r) = n.v(r) * 2
NEXT r
d = n.u(2) * n.v(3) - n.v(2) * n.u(3)
LOOP
END SUB
FUNCTION nextcfd (n AS cfa)
DIM d(3), q AS DOUBLE
DIM df(2) AS LONG
nextcfd = Pend
DO
fl = -1 ' fl =-1 for flat z-range
FOR r = 0 TO 3
IF n.v(r) THEN
q = INT(n.u(r) / n.v(r)) ' four corner values z
IF q > Inf THEN q = Inf ' piecewise transmission
d(r) = q
ELSE
d(r) = Done
END IF
IF r THEN fl = fl AND (d(r - 1) = d(r))
NEXT r
'
IF Verb > 1 THEN
OutFrm n: PRINT " ";
FOR r = 0 TO 3: PRINT d(r); : NEXT
IF fl THEN PRINT " ->out"
END IF
'
IF fl THEN
IF d(0) > Done THEN Digout n, d(0)
nextcfd = d(0): EXIT DO ' happens eventually
END IF
'
FOR r = 1 TO 2 ' choose Cf x or y
df(r) = ABS(CLNG(d(0)) - d(r))
df(0) = ABS(CLNG(d(3 - r)) - d(3))
IF df(0) > df(r) THEN df(r) = df(0)
NEXT r
p = 0: IF df(1) > df(2) THEN p = 1
'
LOOP WHILE digin(n, p) ' digit stream <-in
END FUNCTION
SUB OutCf (s AS STRING, c())
DIM p(1) AS DOUBLE, q(1) AS DOUBLE
DIM n AS cfa, t AS DOUBLE
i = 0: j = 1
p(j) = 1: q(i) = 1
Cdup c(), n.y()
'
PRINT s; " [";
FOR r = 2 TO Cfx - 1
IF c(r) > Done THEN
PRINT c(r);
p(i) = p(i) + c(r) * p(j) ' Cf recurrence
q(i) = q(i) + c(r) * q(j)
SWAP i, j
IF q(j) THEN
t = 1 / (q(j) * q(j))
IF t < Eps OR t = 0 THEN
r = r + 1: n.y(r) = Done
EXIT FOR
END IF ' fp-noise obliterates all
IF t < 1D-48 THEN n.y(r) = Done
END IF
ELSE
EXIT FOR
END IF
NEXT r
PRINT "] "; r - 2
'
PRINT p(j); "/"; q(j) ' convergent
IF q(j) THEN
r = c(2) < 0
IF r THEN
PRINT "-";
ELSE
r = 1: PRINT " ";
END IF
Tset n, r, 0, r, 0, 0, 1, 0, 1
PRINT LTRIM$(STR$(nextcfd(n))) + ".";
DO
FOR r = 0 TO 3 ' decimal expansion
t = n.u(r): n.u(r) = n.v(r) * 10: n.v(r) = t
NEXT r
r = n.y(n.y(0)) > Done
PRINT LTRIM$(STR$(nextcfd(n)));
LOOP WHILE r
PRINT
ELSE
PRINT " Inf"
END IF
END SUB
SUB OutFrm (n AS cfa)
DIM Fs AS STRING, Gs AS STRING
FOR r = 0 TO 3
Fs = Fs + STR$(n.u(r)) + " "
Gs = Gs + " " + STR$(n.v(r))
NEXT r
PRINT Fs + "/"; Gs
END SUB
SUB ReadCf (c(), w AS STRING)
DIM s AS STRING
s = LTRIM$(RTRIM$(w)) + "."
k = 2: i = 1: r = LEN(s)
DO
j = INSTR(i, s, ".")
IF k < Cfx THEN ' read partial quotient
c(k) = VAL(MID$(s, i, j - i))
ELSE
PRINT "Cf truncated"
EXIT DO
END IF
k = k + 1: i = j + 1
LOOP UNTIL i > r
c(k) = Done
END SUB
SUB Tmul (tr() AS DOUBLE, n AS cfa)
DIM tmp AS DOUBLE
FOR r = 0 TO 1 ' unrolled matrix mult
s = r + 2
tmp = tr(0) * n.u(r) + tr(1) * n.u(s)
n.u(s) = tr(2) * n.u(r) + tr(3) * n.u(s)
n.u(r) = tmp
tmp = tr(0) * n.v(r) + tr(1) * n.v(s)
n.v(s) = tr(2) * n.v(r) + tr(3) * n.v(s)
n.v(r) = tmp
NEXT r
END SUB
SUB Trans (n AS cfa)
SWAP n.u(2), n.u(1)
SWAP n.v(2), n.v(1)
END SUB
SUB Tset (n AS cfa, a, b, c, d, e, f, g, h)
n.u(0) = a: n.u(1) = b: n.u(2) = c: n.u(3) = d
n.v(0) = e: n.v(1) = f: n.v(2) = g: n.v(3) = h
n.x(0) = 2: n.x(1) = 2: n.y(0) = 2: n.y(1) = 2
n.flag = 25
END SUB
SUB Turn (n AS cfa)
SWAP n.u(2), n.v(0)
SWAP n.u(3), n.v(1)
END SUB