Martin, one could pursue the continued square root approach by writing a
square root verb that takes and produces rational numbers. Let us
backtrack a little first, the continued fraction approach is simple but
converges slowly,
pp=. (j.@:[ ": ]) NB. Print precision
v0=. 1 + % NB. Continued fraction iteration
66 pp (v0^:_ ) 1 NB. Floating
1.618033988749908700000000000000000000000000000000000000000000000000
J decided that was good enough; but, we can force more iterations,
66 pp (v0^:99) 1 NB. Floating
1.618033988749894900000000000000000000000000000000000000000000000000
That is the floating number you also got with your approach (apparently the
last digit is wrong even after taking rounding into account).
I mentioned that converges slowly, this shows consecutive rational
approximations,
66 pp ,. (v0^:(i.11)) 1x NB. Rational
1.000000000000000000000000000000000000000000000000000000000000000000
2.000000000000000000000000000000000000000000000000000000000000000000
1.500000000000000000000000000000000000000000000000000000000000000000
1.666666666666666666666666666666666666666666666666666666666666666667
1.600000000000000000000000000000000000000000000000000000000000000000
1.625000000000000000000000000000000000000000000000000000000000000000
1.615384615384615384615384615384615384615384615384615384615384615385
1.619047619047619047619047619047619047619047619047619047619047619048
1.617647058823529411764705882352941176470588235294117647058823529412
1.618181818181818181818181818181818181818181818181818181818181818182
1.617977528089887640449438202247191011235955056179775280898876404494
...
66 pp ,. (v0^:(98 99)) 1x
1.618033988749894848204586834365638117720299848718891653939793515920
1.618033988749894848204586834365638117720312743963795685753591851088
How about the square root approach? As Pascal pointed out, %: converts
rationals to floating numbers,
datatype %: 2x
floating
One way to write a rational square root verb is using Newton's method,
u=. *: - 2:
( v1=. ] - u f. % (u f. d.1) ) NB. Newton's iteration for finding the
square root of 2
] - (*: - 2:) % +:
66 pp ,. (v1^:(i.11)) 1x NB. First 10 iterations...
1.000000000000000000000000000000000000000000000000000000000000000000
1.500000000000000000000000000000000000000000000000000000000000000000
1.416666666666666666666666666666666666666666666666666666666666666667
1.414215686274509803921568627450980392156862745098039215686274509804
1.414213562374689910626295578890134910116559622115744044584905019200
1.414213562373095048801689623502530243614981925776197428498289498623
1.414213562373095048801688724209698078569671875377234001561013133113
1.414213562373095048801688724209698078569671875376948073176679737991
1.414213562373095048801688724209698078569671875376948073176679737991
1.414213562373095048801688724209698078569671875376948073176679737991
1.414213562373095048801688724209698078569671875376948073176679737991
Which is not bad at all (comparing these digits to those in [0]); more
precision can be obtained by replacing 11 by a higher number. Moreover, a
square root verb follows somewhat easily,
root=. (] - (*:@:] - [) % +:@:])^:11&1
66 pp root 2 NB. Floating
1.414213562373095100000000000000000000000000000000000000000000000000
66 pp root 2x NB. Rational
1.414213562373095048801688724209698078569671875376948073176679737991
66 pp root 3x
1.732050807568877293527446341505872366942805253810380628055806979452
datatype %: 2x
floating
datatype root 2x
rational
So, you could go further with the continued square root approach by using
root (with a suitable number of iterations) instead of %: in your
iteration. Then again, if one is going to use Newton's method, one can
apply it directly to find the Golden Section,
u=. *: + - + _1:
( v2=. ] - u f. % (u f. d.1) ) NB. Newton's iteration
] - (*: + - + _1:) % _1 2x&p.
NB. Notice the x in _1 2x&p. :)
66 pp ,. (v2^:(i.11)) 1
1.000000000000000000000000000000000000000000000000000000000000000000
2.000000000000000000000000000000000000000000000000000000000000000000
1.666666666666666666666666666666666666666666666666666666666666666667
1.619047619047619047619047619047619047619047619047619047619047619048
1.618034447821681864235055724417426545086119554204660587639311043566
1.618033988749989097047296779290725053240839568674600343661069205517
1.618033988749894848204586838338166878717703891187710377694114386726
1.618033988749894848204586834365638117720309179805762869192919556392
1.618033988749894848204586834365638117720309179805762862135448622705
1.618033988749894848204586834365638117720309179805762862135448622705
1.618033988749894848204586834365638117720309179805762862135448622705
datatype ,. (v2^:(i.11)) 1 NB. Just, checking...
rational
As usual, when Newton's method converges, it converges relatively fast.
I hope it helps.
Reference
[0] The Square Root of Two to 1 Million Digits
http://apod.nasa.gov/htmltest/gifcity/sqrt2.1mil
On Sat, Apr 16, 2016 at 12:49 PM, Martin Kreuzer <[email protected]>
wrote:
> Thanks Pascal -
> for pointing that out: square root as (a C language double float) will
> have a maximum precision of about that number of significant digits.
> (I have become rusty - it's been a long time since.)
> -M
>
>
> At 2016-04-16 13:43, you wrote:
>
>> 0j47": (+`%)/ 1000$1x
>> 1.61803398874989484820458683436563811772030917981 0j35": %: 1r2
>> 0.70710678118654757000000000000000000 The division approach maintains a
>> rational throughout the expression and so can do unlimited precision when
>> converting to display. Square root of a rational is a floating point number
>> that cuts off precision. ----- Original Message ----- From: Martin Kreuzer <
>> [email protected]> To: [email protected] Sent: Saturday, April
>> 16, 2016 5:07 AM Subject: Re: [Jprogramming] Golden Ratio (sqrt approach)
>> ... Thanks Devon - I've experimented with extended precision in different
>> contexts and am able to reproduce your results here (which use the
>> "continued fraction" approach). Using the "continued root" approach, and
>> having set print precision to 20 digits (which seems to be the maximum, as
>> (9!:11 (21)) throws a limit error) I get (+&%:)/ 0 ,1000#1
>> 1.6180339887498949 which shows less than 20 digits; formatting this line
>> shows where it breaks off 0j35": 1-~ (+&%:)/ 1000#1x
>> 1.61803398874989490000000000000000000 I then tried (not knowing whether
>> this makes any sense as the implementation might be the same in the first
>> place) -- square root defined via the primitive: rtp=. %: -- square root
>> defined via fractional exponent rtf=. 1r2^~] with the same result:
>> 0j35": (+&rtp)/ 0 ,1000$1x 1.61803398874989490000000000000000000
>> 0j35": (+&rtf)/ 0 ,1000$1x 1.61803398874989490000000000000000000
>> Comparing 9!:10 '' 6 0j35": (+&rtf)/ 0 ,1000$1x
>> 1.61803398874989490000000000000000000 0j35": 1-~(+&rtf)/ 1000$1x
>> 1.61803398874989490000000000000000000 0j35": (+`%)/ 1000$1x
>> 1.61803398874989484820458683436563812 seems to indicate that the choice of
>> print precision doesn't effect it. That's where I got stuck ... -M At
>> 2016-04-16 03:23, you wrote: > (+`%)/10$1 1.625 (+`%)/100$1 1.61803
>> (+`%)/100$1x > 32951280099r20365011074 0j25":(+`%)/100$1x >
>> 1.6180339887498948482035085 0j25":(+`%)/1000$1x >
>> 1.6180339887498948482045868 0j55":(+`%)/1000$1x >
>> 1.6180339887498948482045868343656381177203091798057628621 >
>> 0j55":(+`%)/2000$1x >
>> 1.6180339887498948482045868343656381177203091798057628621 On Thu, > Apr 14,
>> 2016 at 7:02 PM, Martin Kreuzer <[email protected]> > wrote: > Thanks
>> Raul, Kip -- > > (1) Am now comfortable (again) > with foreign print
>> precision; just for the > record: > > def=. > 6 > ext=. 13 > 9!:10
>> '' NB. check > 6 > 9!:11 (ext) NB. > set > > 9!:10 '' NB. check >
>> 13 > 9!:11 ] def NB. > reset > > 9!:10 '' > 6 > > (2) The tacit
>> version looks pretty > straight forward (and works fine). > > (3) Will
>> follow up your > ideas/suggestions on "injection" and verb (cr) -- >
>> fodder for the > weekend perhaps (no chance to stay with it on a regular >
>> basis > yet, at least there's the habit of reading through the posts). > >
>> > (4) Meanwhile, looking at your parameters of the (a cr b) call, I > did
>> this > small modification, getting rid of the "cleanup" > (1-~): > > ]
>> v=. 0 ,13#1 > 0 1 1 1 1 1 1 1 1 1 1 1 1 > 1 > (+&%:)/ v > 1.618033473928
>> > > -M > > At 2016-04-14 21:43, > you wrote: > >> The following is not
>> simpler but invents a > "continued root" in which >> dyadic root %: plays
>> the role > dyadic % plays in a continued fraction. >> Here is the picture
>> of > a "continued root": a1 > b0 + %: >> b1 + a2
>> %: > b2 + a3 >> %:
>> b3 > + . >> . > . Verb cr >>
>> below produces "convergents" which stop with a > "diagonal" element. The
>> >> basic idea in verb cr belongs to Raul > Miller. cr =: {.@] , {.@] + [:
>> ([: >> %:`+/ ,)\ [: |: [ ,: }.@] NB. > Usage a cr b (13#2) cr 0,13#1 NB.
>> The 2 >> means square roots > are used 0 1 1.414213562 1.553773974
>> 1.598053182 >> 1.611847754 > 1.616121207 1.617442799 1.617851291
>> 1.617977531 1.618016542 >> > 1.618028597 1.618032323 1.618033474 --Kip
>> Murray On Thursday, April > 14, >> 2016, Raul Miller <
>> [email protected]> wrote: > One thing > you could >> do is get rid
>> of the intermediate names in gr: > > > gr=: monad define > >> 1-~ (+&%:)/
>> y$1 > ) > > And you might want > to make this tacit, for example: >> > >
>> 13 :'1-~ (+&%:)/ y$1' > > 1 -~ [: +&%:/ 1 $~ ] > > Or, depending on >> your
>> preferences, you > might want to use induction > rather than insertion: >>
>> > > > gri=: 1-~(1+%:@])^:([-1:)&1 > > I guess it's really a matter >
>> of >> what your idea of "elegance" is... > > Personally, when I am >
>> fiddling with >> small expressions, I like to set > up a line that >
>> evaluates and then tweak >> the expression and watch to > make sure > the
>> result does not change. For >> this example, I'd have lines > > like: > >
>> gr 10 > 1.61798 > >> 1-~(1+%:)^:9]1 > 1.61798 > 13 > :'((+&%:)/ y$1) -
>> 1' 10 > 1.61798 > >> (1-~(1+%:@])^:([-1:)&1) 10 > > 1.61798 > > (with
>> lots of other lines, >> including some errors, > mixed in) > > But the
>> precision issue you are seeing >> is really > the print precision > global
>> parameter. See > >> > http://www.jsoftware.com/help/dictionary/dx009.htm
>> for how to > change > >> that. > > I hope this helps. > > -- > Raul > > On
>> Thu, > Apr 14, 2016 at 3:13 >> PM, Martin Kreuzer <[email protected] >
>> > <javascript:;>> wrote: > > >> Moving from continued fraction to >
>> continued square root, I arrived at > >> this: > > > > NB. > modelling
>> gr=. rt(1+rt(1+rt(1+rt(1+...)))) > > > > >> gr=. monad > define > > ps=. +
>> > > rt=. %: > > v=. y $ 1 > > r=. 1-~ (ps&rt)/ >> > v > > ) > > gr 10 >
>> > 1.61798 > > gr 13 > > 1.61803 > > > > > Q1: > > >> What would be (more
>> elegant and/or concise) ways to do > this, especially > >> the > > line
>> with the return value (r)..? > > > Q2: > > What should I do to >> get
>> higher precision (more digits) > in the result (and > > still having a >>
>> floating point number); > does that need a "foreign"..? > > (I'm sure that
>> I >> have seen > this before, but can't remember where.) > > > > Thanks > >
>> > -M > > >> > > >
>> ---------------------------------------------------------------------- >
>> >> > > For information about J forums see >> >
>> http://www.jsoftware.com/forums.htm > >> >
>> ---------------------------------------------------------------------- >
>> > >> For information about J forums see >
>> http://www.jsoftware.com/forums.htm >> -- Sent from Gmail Mobile >> >
>> ---------------------------------------------------------------------- >
>> For >> information about J forums see >
>> http://www.jsoftware.com/forums.htm >> > > >
>> ---------------------------------------------------------------------- >
>> > For information about J forums see >
>> http://www.jsoftware.com/forums.htm > -- Devon McCormick, CFA >
>> Quantitative Consultant >
>> ---------------------------------------------------------------------- >
>> For information about J forums see http://www.jsoftware.com/forums.htm
>> ---------------------------------------------------------------------- For
>> information about J forums see http://www.jsoftware.com/forums.htm
>> ---------------------------------------------------------------------- For
>> information about J forums see http://www.jsoftware.com/forums.htm
>>
>
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
