Well, that was kind of embarrassing. Apparently, we did not have any tests on converting bulky rational numbers to double.
Fix was simple, though. (And, I have added this example as a test.) Thanks, -- Raul On Mon, Feb 20, 2023 at 8:14 AM Raul Miller <[email protected]> wrote: > > Today I stumbled over a problem. The short form is that I goofed in > my implementation of jtDfromQ in k.c. I might have this fixed in an > hour or so, or it might take me longer. (I'm going to get breakfast, > first.) > > Long form: > > phi=: (+%)/1e3#1x > 0.0+phi > 1.61803 > 0.0+2*phi > 3.23607 > 0.0+_1+2*phi > 2.23607 > 0.0+*~_1+2*phi > 4.25353e38 > > So I built a program in C to mimic this calculation: > > #include <gmp.h> > #include <stdio.h> > > int main() { > // golden ratio as continued fraction > mpq_t one; > mpq_init(one); > mpq_set_ui(one, 1, 1); > mpq_t phi; > mpq_init(phi); > mpq_set_ui(phi, 1, 1); > for (int j= 0; j<999; j++) { > mpq_div(phi, one, phi); > mpq_add(phi, one, phi); > } > printf("phi: %g -- %s\n", mpq_get_d(phi), mpq_get_str(0, 10, phi)); > mpq_t two; > mpq_init(two); > mpq_set_ui(two, 2, 1); > mpq_t phi2; > mpq_init(phi2); > mpq_mul(phi2, two, phi); > printf("2*phi: %g -- %s\n", mpq_get_d(phi2), mpq_get_str(0, 10, phi2)); > mpq_t sqr5; > mpq_init(sqr5); > mpq_sub(sqr5, phi2, one); > printf("square root of 5: %g -- %s\n", mpq_get_d(sqr5), > mpq_get_str(0, 10, sqr5)); > mpq_t five; > mpq_init(five); > mpq_mul(five, sqr5, sqr5); > printf("5: %g -- %s\n", mpq_get_d(five), mpq_get_str(0, 10, five)); > } > > and built it with cc -g sq5.c -o sq5 -lgmp > > It displayed: > phi: 1.61803 -- > 70330367711422815821835254877183549770181269836358732742604905087154537118196933579742249494562611733487750449241765991088186363265450223647106012053374121273867339111198139373125598767690091902245245323403501/43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875 > 2*phi: 3.23607 -- > 140660735422845631643670509754367099540362539672717465485209810174309074236393867159484498989125223466975500898483531982176372726530900447294212024106748242547734678222396278746251197535380183804490490646807002/43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875 > square root of 5: 2.23607 -- > 97194177735908175207981982079326473737797879155345685082728081084772518818444815269080619149045968297679578305403209347401163036907660573971740862463751801641201490284097309096322681531675707666695323797578127/43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875 > 5: 5 -- > 9446708185759308415384067495999677431530963218480368032804826598281856324445977322684945038267086094364761366000137291348836189673785457326607903364013465483957273836804336595888397782139002535468799414419546535346394066447256463745311310661259359973909189379826722425332112242554370313063917929424669185186291673823764654829513873821477637371237697744102254002802127905427315493403711022179894479121632130910668828129/1889341637151861683076813499199935486306192643696073606560965319656371264889195464536989007653417218872952273200027458269767237934757091465321580672802693096791454767360867319177679556427800507093759882883909307069278813289451292749062262132251871994781837875965344485066422448510874062612783585884933837037258334764752930965902774764295527474247539548820450800560425581085463098680742204435978895824326426182133765625 > > Comparing with my J session: > > phi > 70330367711422815821835254877183549770181269836358732742604905087154537118196933579742249494562611733487750449241765991088186363265450223647106012053374121273867339111198139373125598767690091902245245323403501r4346655768693745643568852767504062580256466051... > 2*phi > 140660735422845631643670509754367099540362539672717465485209810174309074236393867159484498989125223466975500898483531982176372726530900447294212024106748242547734678222396278746251197535380183804490490646807002r434665576869374564356885276750406258025646605... > (2*phi)-1 > 97194177735908175207981982079326473737797879155345685082728081084772518818444815269080619149045968297679578305403209347401163036907660573971740862463751801641201490284097309096322681531675707666695323797578127r4346655768693745643568852767504062580256466051... > *~(2*phi)-1 > 9446708185759308415384067495999677431530963218480368032804826598281856324445977322684945038267086094364761366000137291348836189673785457326607903364013465483957273836804336595888397782139002535468799414419546535346394066447256463745311310661259359973909189... > > (And if I take the C value for five, replace '/' with 'r' and compare > it to ":*~(2*phi)-1 the values are identical.) > > So, anyways... I need to figure out where I went wrong in jtDfromQ. > > Thanks, > > -- > Raul ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
