√m≈√x+√y

[waxtep]

Andrey A в сообщении #1569847 писал(а):

Вам не кажется, что тут телега впереди лошади немножко заехала? Четные откинуть не проблема, а нечетное кратное квадрату $>1$ — как раз хороший тест для итогового алгоритма, если он отыщется.

Этот мой кусок кода - не алгоритм-кандидат, а фиксация уже известных знаний о природе Чтобы просто иметь возможность увидеть корректные результаты для отдельных $m$ или их диапазонов. У меня был исходно реализован только банальный перебор, но там я столкнулся с забавной и раздражающей проблемой погрешности округления; а таблички уважаемого juna мне показалось мало, к тому же, для нескольких $m$ оттуда, если я правильно понимаю, есть лучшие нижние приближения.

Andrey A в сообщении #1569847 писал(а):

P.S. И проверить заодно $\gcd (x,y)$, которое, кстати, вовсе не обязано быть свободным от квадратов.

Не обязано, просто произвольно выбираю один из вариантов, в котором $x=\gcd(x,y)$ минимально и равно свободной от квадратов части $m$, т.е. $\sqrt{75}=\sqrt3+\sqrt{48}$. Может быть, в поисках закономерностей, интереснее выводить вариант с максимальным $x$, или какой-нибудь еще промежуточный, не знаю.
Пока наблюдений не густо:

• если идти по всем $m$ подряд, конечно, встретятся все возможные значения $y$, без пропусков;
• при этом чуть меньше половины ($24$ из $55$ рассмотренных) нечетных значений $y$ впервые встречаются только для $m=4y-2$, вот начало их последовательности: $1,3,15,17,21,29,35,37,43,47,53,59,65,67,69,73,77,83,85,91,93,97,107,109\ldots$; они свободны от квадратов, но это не все свободные от квадратов нечетные в диапазоне;
• четные $y$, кроме $y=2$, впервые встречаются для $m$ меньших, чем $4y-2$.

[juna]

waxtep в сообщении #1569795 писал(а):

Тут, видимо, некоторое количество неточностей

Вы правы.

Оказывается, искать

$$d_y=m+x-\sqrt{4mx}-\lfloor m+x-\sqrt{4mx}\rfloor\rightarrow\min$$

недостаточно. Отсюда и неточности.

Пусть $y_d=\lfloor m+x-\sqrt{4mx}\rfloor$

$$d_x=m+y_d-\sqrt{4my_d}-\lfloor m+y_d-\sqrt{4my_d}\rfloor$$

Будем искать:

$$d_x+d_y\rightarrow\min$$

при $x=1\ldots \lfloor\frac{m}{2}\rfloor-1$
Вот соответствующий код в maxima.

gen(m):=block(
n:floor(m/2)-1,
min:m, 
for x:1 thru n do 
  block([y,dy,yd,dx], 
        y:m+x-2*sqrt(m*x),yd:floor(m+x-2*sqrt(m*x)),
        dy:y-yd,dx:m+yd-2*sqrt(m*yd)-floor(m+yd-2*sqrt(m*yd)),
       /* print(float(dx+dy), float(dx), float(dy), yd, x), */
        if min>(dx+dy) then 
           block(min:(dx+dy), xn:x,yn:y)),
return([xn,floor(yn)]))$

Выполнял сравнение с простым двухпараметрическим перебором. Отличий в диапазоне от 4 до 200 не было.

(Оффтоп)

Прикрепил бы сюда расчетную таблицу, но не знаю, как

[waxtep]

juna в сообщении #1569860 писал(а):

Прикрепил бы сюда расчетную таблицу, но не знаю, как

Можно наверное просто текстом с разделителями:

(осторожно, 397 строк)

m;x;y
4;1;1
5;1;1
6;1;2
7;1;2
8;2;2
9;1;4
10;2;3
11;1;5
12;3;3
13;3;3
14;3;4
15;2;6
16;1;9
17;2;7
18;2;8
19;1;11
20;5;5
21;2;10
22;5;6
23;1;14
24;6;6
25;1;16
26;6;7
27;3;12
28;7;7
29;1;19
30;7;8
31;5;11
32;2;18
33;4;14
34;8;9
35;6;12
36;1;25
37;6;13
38;9;10
39;4;18
40;10;10
41;7;14
42;10;11
43;3;23
44;11;11
45;5;20
46;11;12
47;3;26
48;3;27
49;1;36
50;2;32
51;6;22
52;13;13
53;6;23
54;6;24
55;2;36
56;14;14
57;6;26
58;14;15
59;11;19
60;15;15
61;5;31
62;15;16
63;7;28
64;1;49
65;10;24
66;16;17
67;4;38
68;17;17
69;8;30
70;17;18
71;1;55
72;2;50
73;14;23
74;18;19
75;3;48
76;19;19
77;6;40
78;19;20
79;11;31
80;5;45
81;1;64
82;20;21
83;16;26
84;21;21
85;14;30
86;21;22
87;10;38
88;22;22
89;7;46
90;10;40
91;2;66
92;23;23
93;10;42
94;23;24
95;4;60
96;6;54
97;2;71
98;2;72
99;11;44
100;1;81
101;5;61
102;25;26
103;14;41
104;26;26
105;12;46
106;26;27
107;3;74
108;3;75
109;22;33
110;27;28
111;12;50
112;7;63
113;23;34
114;28;29
115;18;42
116;29;29
117;13;52
118;29;30
119;10;60
120;30;30
121;1;100
122;30;31
123;14;54
124;31;31
125;5;80
126;14;56
127;17;51
128;2;98
129;14;58
130;32;33
131;27;39
132;33;33
133;24;44
134;33;34
135;15;60
136;34;34
137;7;82
138;34;35
139;29;41
140;35;35
141;16;62
142;35;36
143;30;42
144;1;121
145;6;92
146;36;37
147;3;108
148;37;37
149;11;79
150;6;96
151;5;101
152;38;38
153;17;68
154;38;39
155;6;100
156;39;39
157;33;46
158;39;40
159;18;70
160;10;90
161;30;52
162;2;128
163;15;79
164;41;41
165;26;60
166;41;42
167;23;66
168;42;42
169;1;144
170;42;43
171;19;76
172;43;43
173;9;103
174;43;44
175;7;112
176;11;99
177;20;78
178;44;45
179;31;61
180;5;125
181;17;87
182;45;46
183;20;82
184;46;46
185;30;66
186;46;47
187;6;126
188;47;47
189;21;84
190;47;48
191;41;55
192;3;147
193;14;103
194;48;49
195;42;56
196;1;169
197;42;57
198;22;88
199;19;95
200;2;162
201;22;90
202;50;51
203;4;150
204;51;51
205;8;132
206;51;52
207;23;92
208;13;117
209;28;84
210;52;53
211;5;151
212;53;53
213;24;94
214;53;54
215;34;78
216;6;150
217;18;110
218;54;55
219;24;98
220;55;55
221;12;130
222;55;56
223;32;86
224;14;126
225;1;196
226;56;57
227;11;138
228;57;57
229;29;95
230;57;58
231;42;76
232;58;58
233;31;94
234;26;104
235;38;84
236;59;59
237;26;106
238;59;60
239;41;82
240;15;135
241;11;149
242;2;200
243;3;192
244;61;61
245;5;180
246;61;62
247;6;176
248;62;62
249;28;110
250;10;160
251;55;71
252;7;175
253;2;210
254;63;64
255;56;72
256;1;225
257;14;151
258;64;65
259;48;84
260;65;65
261;29;116
262;65;66
263;46;89
264;66;66
265;42;96
266;66;67
267;30;118
268;67;67
269;19;145
270;30;120
271;34;113
272;17;153
273;40;104
274;68;69
275;11;176
276;69;69
277;61;78
278;69;70
279;31;124
280;70;70
281;62;79
282;70;71
283;9;191
284;71;71
285;46;102
286;71;72
287;6;210
288;2;242
289;1;256
290;72;73
291;32;130
292;73;73
293;59;89
294;6;216
295;12;188
296;74;74
297;33;132
298;74;75
299;28;144
300;3;243
301;6;222
302;75;76
303;34;134
304;19;171
305;12;196
306;34;136
307;68;86
308;77;77
309;34;138
310;77;78
311;11;205
312;78;78
313;46;119
314;78;79
315;35;140
316;79;79
317;15;194
318;79;80
319;24;168
320;5;245
321;36;142
322;80;81
323;72;90
324;1;289
325;13;208
326;81;82
327;36;146
328;82;82
329;60;108
330;82;83
331;41;139
332;83;83
333;37;148
334;83;84
335;54;120
336;21;189
337;69;101
338;2;288
339;38;150
340;85;85
341;70;102
342;38;152
343;7;252
344;86;86
345;38;154
346;86;87
347;61;117
348;87;87
349;59;121
350;14;224
351;39;156
352;22;198
353;79;98
354;88;89
355;14;228
356;89;89
357;66;116
358;89;90
359;31;179
360;10;250
361;1;324
362;90;91
363;3;300
364;91;91
365;58;132
366;91;92
367;17;226
368;23;207
369;41;164
370;92;93
371;30;190
372;93;93
373;46;157
374;93;94
375;15;240
376;94;94
377;36;180
378;42;168
379;85;105
380;95;95
381;42;170
382;95;96
383;49;158
384;6;294
385;80;114
386;96;97
387;43;172
388;97;97
389;59;145
390;97;98
391;12;266
392;2;338
393;44;174
394;98;99
395;16;252
396;11;275
397;49;167
398;99;100
399;90;110
400;1;361

[juna]

Нашел, как добавить.

[waxtep]

juna в сообщении #1569863 писал(а):

Нашел, как добавить.

Ага, совпадаем

[Andrey A]

Спасибо за списки и за статистику. Всё это любопытно в сочетании с ожидаемо всплывшей задачей описания последовательностей радикалов. Нечто похожее уже было, причем с решением — задача восстановления целого радикала по одному, двум и более начальных знаков дроби. Там образуются прогрессии с "почти целой разностью", которая стремится к целому. Теперь что-то и не найти. А хотя, кое-что есть здесь (перед заголовком "Радикалы степени k"). Возвращаясь к началу, замечу только, что для составного нечетного задача в некотором смысле проще, ведь получив любое квадратное $R$, уже имеем пару множителей $m$ и некоторое решение с $R=1$ или $R=0$. Решить же ее для простого $m$, исходя из величины, означало бы решить полностью.

[waxtep]

juna в сообщении #1569860 писал(а):

Будем искать:

$$d_x+d_y\rightarrow\min$$

Тут до меня с некоторым запозданием дошло, в каком смысле $R$ должно быть мало: для $4mx=u^2-R, 4my=v^2-R, u+v=2m$ можно записать

$$2\sqrt{m}\left(\sqrt{m}-\sqrt{x}-\sqrt{y}\right)=u-\sqrt{u^2-R}+v-\sqrt{v^2-R}\approx\frac{R}2\left(\frac1{u^2}+\frac1{v^2}\right)+\frac{R^2}8\left(\frac1{u^4}+\frac1{v^4}\right)+\ldots$$

т.е. оптимум соответствует случаю, когда среднее арифметическое относительного отставания $4mx$ и $4my$ от квадрата целого минимально. Например, для $\sqrt{683}\approx\sqrt{121}+\sqrt{229}$, само по себе $R=53$ не сказать, чтобы "мало", но $R\left(\frac1x+\frac1y\right)$ действительно минимально для всех возможных $x$.
Заодно небольшая модификация кода основной функции, чтобы не связываться с нехваткой стека + возвращает и $R$:

fra_base1(m)={mf=factor(m); np=matsize(mf)[1]; p=vector(np,i,mf[i,1]); ed=vector(np,i,floor(mf[i,2]/2)); msq=factorback(p,ed); msqf=m/msq^2;
if(msq>1, return([0,msqf,msqf*(msq-1)^2,0]));
sm=floor(m/4);
if(mf[1,1]==2, return([fra_d(m,sm,sm+1),sm,sm+1,1])); 
dxy=vector(3,i,0); x=1; y=fra_y(m,x); d=fra_d(m,x,y); dxy=[d,x,y]; 
for(x=2,sm,if(fra_d(m,x,fra_y(m,x))<dxy[1],dxy=[fra_d(m,x,fra_y(m,x)),x,fra_y(m,x)])); 
d=dxy[1]; x=dxy[2]; y=dxy[3]; 
return([d,x,y,ceil(sqrt(4*x*y))^2-4*x*y])};

[Andrey A]

waxtep в сообщении #1569926 писал(а):

... но $R\left(\frac1x+\frac1y\right)$ действительно минимально для всех возможных $x$.

Жаль только, что судить о качестве приближения можем только перебрав все остальные варианты.

Andrey A в сообщении #1569101 писал(а):

Короче, тут лишь способ получить относительно хорошее приближение для начала некоторого (загадочного) процесса.

Как раз об этом. Но закономерности копятся не в пустую, глядишь и просветлеет. Вы не упускаете из виду, что $v$ и $R$ — числа разной четности?

[waxtep]

Andrey A в сообщении #1569949 писал(а):

Вы не упускаете из виду, что $v$ и $R$ — числа разной четности?

Для нечетных $m$, да. Я в предыдущем сообщении неудачно ввел обозначения, то что у Вас здесь $m\pm v$, у меня $u$ и $v$

[waxtep]

Andrey A в сообщении #1569081 писал(а):

Такая схема действительно выражает общее решение для составного $m$ и хорошо описывается в цепных дробях. Любая пара множителей $m$ взаимно проста, поскольку $m$ свободно от квадратов.
Возьмем $m=77=11 \cdot 7$, разложим $\dfrac{11}{7}=1,1,1,3$ и выпишем подходящие дроби, уменьшив последний знак на единицу: $1,1,1,2=\dfrac{1}{1},\dfrac{2}{1},\dfrac{3}{2},\dfrac{8}{5}.$ Элементы двух последних дробей возвращают нужные $a,b,c,d$ из описанной выше схемы: $\sqrt{3 \cdot 2}+\sqrt{8 \cdot 5}=\sqrt{(8+3)(5+2)}.$

А у Вас есть понимание, какое именно разбиение на множители составного свободного от квадратов нечетного дает оптимум? Почему, например, для $105$ это именно $3\cdot35$, а для $1001$ - $11\cdot91$?

[Andrey A]

waxtep в сообщении #1569968 писал(а):

Почему, например, для $105$ это именно $3\cdot35$, а для $1001$ - $11\cdot91$?

О! Наконец пошли вопросы, а то я думал уже не выберусь из-под ответов ) Для составного $m$ должно быть $R=1$ и, значит, $v$ четно. Уже говорилось, что при наличии вариантов оно должно быть наименьшим. Если оно добывается перебором, то первое по возрастанию. В противном случае, боюсь, должны получить все возможные. При разложении в цепную дробь, думаю, следует выбрать с наименьшим последним знаком, но тогда уже не очень важно — работа сделана. Вот и проверим. Для $m=1001$ имеем $77/13=5,1,12;\ 91/11=8,3,1,2;\ 143/7=20,2,3.$

$20,2,2=\dfrac{20}{1},\dfrac{41}{2},\dfrac{102}{5};\ x=82,y=510;\ v=102 \cdot 5 - 41 \cdot 2=428.$

$8,3,1,1=\dfrac{8}{1},\dfrac{25}{3},\dfrac{33}{4},\dfrac{58}{7};\ x=132,y=406;\ v=58 \cdot 7 - 33 \cdot 4=274.$

$5,1,11=\dfrac{5}{1},\dfrac{6}{1},\dfrac{71}{12};\ x=6,y=852;\ v=71 \cdot 12 - 6 \cdot 1=846.$

Действительно решение берется из второй дроби. Но ведь последние знаки могут совпадать в различных разложениях. Тогда, возможно, по предпоследнему.

[Andrey A]

P.S.

Andrey A в сообщении #1569974 писал(а):

Если оно добывается перебором, то первое по возрастанию.

Не факт. Первый по возрастанию квадрат сравнимый с единицей может оказаться нечетным. А перебирать все равносильно факторизации. Неинтересно. Да, и по поводу многомерных чисел: $m$ по предположению свободно от квадратов, в каноническом разложении только первые степени. Положим, оно имеет $N$ простых делителей. Тогда общее число делителей $2^N$, общее число пар множителей — $2^{N-1}$, без учета тривиальной пары $m \cdot 1$ их $2^{N-1}-1.$ Восстанавливать $31$-у дробь действительно скучновато. В связи с этим
$1$) Хорошо бы доказать, что сортировка дробей в обратном движении равносильна сортировке по качеству приближений. Или опровергнуть.
$2$) Предположу, что перебор по убыванию до первого $x$ такого, что $(x-y)^2 \equiv 1 \mod m$ окажется в подобных обстоятельствах продуктивней факторизации. А он уж точно соответствует лучшему приближению. Верхней границей, кстати, нужно брать не $m/4$, а $(\sqrt{m}-1)^2/4$ и учитывать только четные $x,y$. Для нашего опыта это будет $(\sqrt{1001}-1)^2/4 \approx 234$, а искомый $x=274.$ Всего $20$ итераций — уже сравнимая работа, а это только для трех простых $7 \cdot 11 \cdot 13.$
$3$) Похоже, в этой процедуре заложена некая модернизация метода Ферма. Об эффективности судить пока не берусь, но плюсы на поверхности: квадраты выбираются не любые, а соответствующие малым вычетам, причем целый класс заведомо неквадратных вычетов $2,3 \mod 4$ отсеиваются автоматически. Оставим на будущее.

[waxtep]

Здесь кажется коварнее: иногда решение соответствует разложению на два составных множителя, например для $23023=7\cdot11\cdot13\cdot23$ наилучшее приближение $\sqrt{23023}\approx\sqrt{4674}+\sqrt{6950}$ соответствует разложению $\dfrac{11\cdot23}{7\cdot13}=2, 1, 3, 1, 1, 4, 2$. Подручный инструмент для борьбы с цепными дробями:

fra_chfr_main(p)={sp=#p; numpmax=floor(sp/2); nmax=2^(sp-numpmax)*(2^numpmax-1); for(cnt=1,nmax,if(vecsum(binary(cnt))<=numpmax,print(fra_chfr(p,binary(cnt)))))};

fra_chfr(p,np0)={m=factorback(p); sp=#p; snp0=#np0; np=vector(sp,i,if(i>=sp-snp0+1,np0[i-sp+snp0],0)); m1=factorback(p,np); m2=m/m1; if(m1>m2,[m1,m2]=[m2,m1]);
a=contfrac(m2/m1); sa=#a; a[sa]=a[sa]-1; pq=contfracpnqn(a); x=pq[1,2]*pq[2,2]; y=pq[1,1]*pq[2,1]; d=fra_d(m,x,y); return([m1,m2,x,y,d]);};

Пример использования (на входе - массив простых делителей $m$, на выходе - печать всех возможных разложений на множители (иногда с дубликатами, тут я пока поленился сделать аккуратнее) и соответствующие им $x,y,d$):

(23:49) gp > fra_chfr_main([7,11,13,23])
[23, 1001, 174, 19194, 4.5078730747833521573646101671686627188 E-7]
[13, 1771, 2180, 11034, 1.6797100425581026541502997663988349877 E-7]
[77, 299, 1122, 13980, 2.0800752218743054735222972986637040989 E-7]
[11, 2093, 3044, 9324, 1.5463421680150141156305087900261879303 E-7]
[91, 253, 4674, 6950, 1.4454141411589071551796186041444902486 E-7]
[143, 161, 72, 20520, 6.7775715710877483065151211079561307366 E-7]
[7, 3289, 470, 16914, 2.9218435717189735088011575289986968511 E-7]
[143, 161, 72, 20520, 6.7775715710877483065151211079561307366 E-7]
[91, 253, 4674, 6950, 1.4454141411589071551796186041444902486 E-7]
[77, 299, 1122, 13980, 2.0800752218743054735222972986637040989 E-7]

[waxtep]

Так поаккуратнее, хотя такой кусок кода на выставку тоже не отнесешь Без дубликатов, отсортировано по качеству приближения

fra_chfr_main(p)={sp=#p; numpmax=floor(sp/2); nmax=2^(sp-numpmax)*(2^numpmax-1); chfr=vector(5,i,-1);
for(cnt=1,nmax,if(vecsum(binary(cnt))<=numpmax,chfr=matconcat([chfr;fra_chfr(p,binary(cnt))])));
chfr=vecsort(chfr~,,8)~; chfr=vecextract(matrix(matsize(chfr)[1]-1,5,m,n,chfr[m+1,n]),[2,3,4,5,1]); return(chfr)};

fra_chfr(p,np0)={m=factorback(p); sp=#p; snp0=#np0; np=vector(sp,i,if(i>=sp-snp0+1,np0[i-sp+snp0],0)); m1=factorback(p,np); m2=m/m1; if(m1>m2,if(sp==2*vecsum(np0),return(vector(5,i,-1)),[m1,m2]=[m2,m1]));
a=contfrac(m2/m1); sa=#a; a[sa]=a[sa]-1; pq=contfracpnqn(a); x=pq[1,2]*pq[2,2]; y=pq[1,1]*pq[2,1]; d=fra_d(m,x,y); return([d,m1,m2,x,y]);};

(02:03) gp > fra_chfr_main([7,11,13,23])
%287 =
[ 91  253 4674  6950 1.4454141411589071551796186041444902486 E-7]

[ 11 2093 3044  9324 1.5463421680150141156305087900261879303 E-7]

[ 13 1771 2180 11034 1.6797100425581026541502997663988349877 E-7]

[ 77  299 1122 13980 2.0800752218743054735222972986637040989 E-7]

[  7 3289  470 16914 2.9218435717189735088011575289986968511 E-7]

[ 23 1001  174 19194 4.5078730747833521573646101671686627188 E-7]

[143  161   72 20520 6.7775715710877483065151211079561307366 E-7]

-- 15.11.2022, 02:17 --

(02:14) gp > fra_chfr_main([7,11,13,23,41,47])
%290 =
[2093   21197 8665400 13816272 1.7151355717169643923691576460627093444 E-12]...

Здесь решение соответствует разложению на множители о трех простых делителях каждый; в общем, похоже, всякое бывает

[Andrey A]

waxtep в сообщении #1570038 писал(а):

Здесь кажется коварнее:

Никакого колдунства, их ровно $2^3-1=7,$ как и положено для $4$-х простых делителей. Теперь выписываю последний и предпоследний знаки дробей в указанном Вами порядке:
$2,4$
$2,1$
$3,4$
$4,1$
$6,1$
$11,1$
$17,1$
Как видим, последний знак действительно не убывает, свидетельствуя о качестве приближения. Для предпоследнего мало данных.
А в чем смысл борьбы с цепными дробями, не поделитесь? Или тут что-то личное )