From Newsgroup: comp.lang.forth
Up to 128 precision is obtainable if we approximate a real
number with a fraction n/d where n d is single precision 64 bit.
What is needed is an "oracle" that determines whether
n/d is too small or too large.
In this example the oracle say that p/q is larger than x,
but this oracle can be replaced to whatever you like.
`find-part finds the n-th convergent. If n is to large to accommodate
the largest ratio is returned.
is used to indicate overflow in adding positive number,
a technique acceptable in Forth.
\ -------------------- 8< --------------------------------------
WANT U>
VARIABLE x
\ For p1 p2 leave p2 . (pairs)
: PNIP >R >R 2DROP R> R> ;
\ Multiply double single : result double .
: DS*D >R R@ M* SWAP ROT R> UM* D+ ;
: D< DNEGATE D+ 0< NIP ;
\ For P Q :" The point is above the line" Q^2>RATIO*P^2
: oracle DUP 0< 1010 AND THROW
>R DUP M* x @ DS*D R> DUP M* D< ;
\ Leave START value for convergents.
: \\: 0 1 1 0 ;
\ For CONV1 CONV2 incorporate K (cf term), return CONV2 CONV3
: \\ >R 2SWAP 2OVER R@ * SWAP R> * SWAP D+ ;
\ For consecutive convergents C1 and C2 find next convergent.
\ Leave next C2 and C3.
VARIABLE CFI
: next-convergent 2SWAP 2OVER D+ 2DUP oracle >R 0 CFI !
BEGIN 1 CFI +! 2OVER D+ 2DUP oracle R@ <> UNTIL RDROP
2OVER DNEGATE D+ ;
: MAX-INT -1 1 RSHIFT ;
\ continued-fraction
\ For current values of 'a 'b 'c whatnot find
\ the nth convergent of as numerator denominator .
: find-part >R \\: 2SWAP
R> 0 DO
'next-convergent CATCH IF &* EMIT 2SWAP LEAVE THEN
DUP MAX-INT U> IF && EMIT 2SWAP LEAVE THEN
LOOP
PNIP ;
\ For A B print its continued fraction and the GCD.
: .CF SWAP BEGIN OVER /MOD . DUP WHILE SWAP REPEAT DROP &| EMIT . ;
\ ----------------------------------------------------------------
Example. Printing the continued fraction proves that the result
is correct.
\ -------------------- 8< --------------------------------------
INCLUDE cfsqrtd.frt
S[ ] OK 2 x !
S[ ] OK 3 find-part
S[ 5 7 ] OK .CF
0 1 2 2 |1
S[ ] OK 100 find-part
S[ 2015874949414289041 2850877693509864481 ] OK .CF
0 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 |1
29 x !
100 find-part .CF
*0 5 2 1 1 2 10 2 1 1 2 10 2 1 1 2 10 2 1 1 2 10 2 1 1 2 10 2 1 1
2 10 2 1 1 2 10 2 1 1 2 10 2 2 |1
\ -------------------- 8< --------------------------------------
Groetjes Albert
--
The Chinese government is satisfied with its military superiority over USA.
The next 5 year plan has as primary goal to advance life expectancy
over 80 years, like Western Europe.
--- Synchronet 3.21a-Linux NewsLink 1.2