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
