2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему На страницу Пред.  1, 2
 
 Re: найти корешок по двум элементам цепной дроби
Сообщение24.10.2020, 15:49 
Заслуженный участник
Аватара пользователя


21/11/12
1870
Санкт-Петербург
P.S. Во второй версии оплошность. Ограничив процесс пятью итерациями, урезал заодно и конечную проверку, а этого делать не следовало.
На всякий случай: если на заключительном этапе (наименьшее $m=...$) в столбце "J" под заголовком $a_2$ все пять строк заполнены числами, лучше проверить результат в первой версии. Бывает редко, не хочется переписывать таблицу.

 Профиль  
                  
 
 Re: найти корешок по двум элементам цепной дроби
Сообщение25.10.2020, 10:27 
Заслуженный участник
Аватара пользователя


21/11/12
1870
Санкт-Петербург
Andrey A в сообщении #1488834 писал(а):
Бывает редко, не хочется переписывать таблицу.

Вроде бы и не редко. Переписал. Ссылка прежняя https://yadi.sk/i/D_CzGy4qt8reiw

В случае одного знака $ задача решается устно:
Для четного $u\ \ m=(u/2)^2+1.$
Для нечетного $u\ \ m=u^2+2.$

Для $3$-х знаков постараюсь всё же написать, хотя тема не для Excel.

 Профиль  
                  
 
 Re: найти корешок по двум элементам цепной дроби
Сообщение26.10.2020, 20:41 
Заслуженный участник
Аватара пользователя


21/11/12
1870
Санкт-Петербург
Andrey A в сообщении #1488990 писал(а):
Для $3$-х знаков постараюсь всё же написать

Готово. Ссылка та же.

 Профиль  
                  
 
 Re: найти корешок по двум элементам цепной дроби
Сообщение27.10.2020, 03:46 
Заслуженный участник
Аватара пользователя


21/11/12
1870
Санкт-Петербург
P.S. Исправлена ошибка 27.10.2020. :oops:

 Профиль  
                  
 
 Re: найти корешок по двум элементам цепной дроби
Сообщение27.10.2020, 06:50 
Заслуженный участник


20/08/14
11067
Россия, Москва

(Тупой перебор)

Сделал втупую перебор $a_0$ на PARI/GP в соответствии с наблюдением mihaild и Null хотя бы просто для проверки:
Код:
? sq(b,c)=my(a,m,x,y); for(a=0,oo, x=(a+c/(b*c+1))^2; y=(a+(c+1)/(b*c+b+1))^2; m=ceil(min(x,y)); if(m<max(x,y), printf("sqrt(%d)=[%d,%d,%d,...]\n", m,a,b,c); break) );
? sq(5,4)
sqrt(174)=[13,5,4,...]
? sq(7,11)
sqrt(2129)=[46,7,11,...]
? sq(997,999)
sqrt(249000003002)=[498999,997,999,...]
time = 1,295 ms.
? sq(999,997)
sqrt(994014984005)=[997003,999,997,...]
time = 2,606 ms.
Не быстро конечно. Вроде линейно по $a_0$, которое очевидно квадратично по $a_1, a_2$ (точнее скорее линейно по их произведению $a_1 a_2$), так что давать $a_1,a_2 \approx  10^9$ не стоит.

-- 27.10.2020, 07:24 --

Более умный перебор (для совсем малых аргументов может ошибаться UPD: проверил, не ошибается) на PARI/GP:
Код:
? sq(b,c)=my(a,k,m,x,y,bc1); bc1=b*c+1; for(k=1,oo, a=floor((k*bc1/c-c/bc1)/2); x=(a+c/(b*c+1))^2; y=(a+(c+1)/(b*c+b+1))^2; m=ceil(min(x,y)); if(m<max(x,y), printf("sqrt(%d)=[%d,%d,%d,...]\n", m,a,b,c); break) );
? sq(999,997)
sqrt(994014984005)=[997003,999,997,...]
? sq(9999,9997)
sqrt(9994001499840005)=[99970003,9999,9997,...]
time = 62 ms.
? sq(99999,99997)
sqrt(99994000149998400005)=[9999700003,99999,99997,...]
time = 765 ms.
? sq(999999,999997)
sqrt(999994000014999984000005)=[999997000003,999999,999997,...]
time = 8,019 ms.
Кажется уже линейно по коэффициентам ... С коэффициентами под $10^9$ (каждый) должно справиться за пару-тройку часов ...
В принципе можно даже и к операциям в целых перейти, тогда вопрос точности вычислений и возможной ошибки отпадёт, но анализировать что-то лень ... С другой стороны PARI похоже и сам всё в целых (точнее конечно рациональных) считает, поленился проверить типы данных, но раз не ошибается и время не зависит от установленной точности вычислений с плавающей точкой, то очень похоже.

-- 27.10.2020, 08:06 --

Не знаю надо ли пояснять идею, она проста: перебирать только те $a_0$, которые после возведения в квадрат дают максимально близкие снизу к целым числам, для чего оценивается слагаемое $k=\frac{2 a_0 a_2}{a_1 a_2 + 1}+(\frac{a_2}{a_1 a_2 + 1})^2$ из раскрытого $(a_0+\frac{a_2}{a_1 a_2 +1})^2$.

 Профиль  
                  
 
 Re: найти корешок по двум элементам цепной дроби
Сообщение28.10.2020, 12:11 


02/04/18
240
Все не получается выделить достаточно времени, поэтому весьма коряво и не для слишком больших значений, зато для произвольной (тоже ограниченной, впрочем, если нужна точность) длины цепи.
Код тут:
http://cpp.sh/85mur
Экзешник тут:
https://yadi.sk/d/pX5RloDX8ja6aw

 Профиль  
                  
 
 Re: найти корешок по двум элементам цепной дроби
Сообщение08.11.2020, 20:39 
Заслуженный участник


20/08/14
11067
Россия, Москва
По подсказкам Andrey A и небольшой своей доработкой реализовал его алгоритм на PARI/GP:
Код:
? {sq3(b,c)=my(x,pq,r,m,p,q2,v,i,pf);
   pf=default(realprecision);default(realprecision,300);\\Для больших чисел может понадобиться и большая точность
   v=contfracpnqn([b,c,b]); if((v[2,2]*v[2,1])%2==0, v[2,2]=-v[2,2], v=contfracpnqn([b,c,c,b]));
   m=round(((v[2,2]*v[2,1]/2)%v[1,1]+v[2,1]/v[1,1])^2);\\Инициализация m каким-то большим точно подходящим числом
   while(1,\\Цикл максимального уменьшения найденного числа сохраняя его "подходящесть"
      p=0; pq=contfracpnqn(contfrac(sqrt(4*m),12),10+2);\\Разложение sqrt(4m) в 10 (+2 запас) подходящих дробей, пока 10-ти хватало
      forstep(i=10,1,-1,\\Перебираем подходящие дроби в обратном порядке, от больших
         r=m-pq[1,i]+pq[2,i]^2;\\Это будет будущее m=m-P+Q^2
         if(r<2 || r>=m, next);\\Такие не интересуют
         v=contfrac(sqrt(r),5); if(#v>3 && v[2..3]==[b,c], m=r; p=pq[1,i]; q2=pq[2,i]^2; break);\\Если новое подходит, то его и используем
      );
      if(p==0, break);\\Нашли ли новое m
      r=p%(2*q2);\\Попытаемся уменьшить p до минимума
      while(r<p,\\Пока оно меньше уже найденного проверяем подходит ли такое p
         x=round(r^2/q2/4); if(x>=m, break);\\Это новое m, если больше уже найденного, то не интересно
         v=contfrac(sqrt(x),5); if(#v>3 && v[2..3]==[b,c], m=x; break);\\Проверяем новое m и если подходит, то используем и выходим
         r+=2*q2;\\Продолжаем перебор в сторону увеличения до нахождения наименьшего подходящего, пока хватало 3-4 циклов
      );
   );
   default(realprecision,pf);\\Вернём точность вычислений обратно
   m
};
? m=sq33(7,11);printf("sqrt(%d)=%d\n", m,contfrac(sqrt(m),5))
sqrt(2129)=[46,7,11,3]
? m=sq33(22,45);printf("sqrt(%d)=%d\n", m,contfrac(sqrt(m),5))
sqrt(1048669)=[1024,22,45,2,7]
? m=sq33(222,11);printf("sqrt(%d)=%d\n", m,contfrac(sqrt(m),5))
sqrt(6522939)=[2554,222,11,2,9]
Для всех тестовых векторов оказалось достаточно 10-ти элементов в цепной дроби разложения корня, 12 элементов берётся исключительно из-за возможного усечения результата contfrac если последнее значение окажется 1.
Опять же для всех тестовых векторов достаточно 5-ти итераций цикла "спуска" по $m$, для больших чисел обычно хватает 2-3 итерации.
Использование встроенной функции contfrac потребовало увеличения точности операций с плавающей точкой, поставил с запасом, на скорость это влияет мало.
Отдельно с помощью Andrey A было проверено что возможна реализация только в длинных рациональных числах, уже не зависящая от точности вычислений с плавающей точкой. Но она медленнее (на не слишком больших числах типа как вон ниже где-то вдвое) и код вдвое объёмнее, потому не привожу, на алгоритм оно не влияет (за одним исключением: перебирать подходящие дроби можно в прямом порядке, причём тогда и разложение в цепную дробь тоже становится итерационным и можно оборвать процесс ранее 10-ти членов, но на скорость это влияет слабо).

Отличия этой реализации от алгоритма Andrey A:
подходящие дроби перебираются от конца к началу;
они формируются сразу все 10 (+2 для обхода оптимизации PARI/GP), а не итерационно;
вместо предложенного им двоичного поиска (деления пополам) наименьшего подходящего $P$ берётся сразу минимальное положительное $P$ и перебирается подряд (с шагом арифметической прогрессии) в сторону увеличения, обычно хватает 1-2 итераций, больше 4 на тестовых векторах не понадобилось ни разу.

Замеры скорости.
Сравнивал со своей реализацией полного перебора выше (оптимизированной примерно вдвое).
98901 комбинацию $[?,1,1]-[?,999,99]$ мой перебирает за 11с, этот за 10.5с.
98901 комбинацию $[?,1,1]-[?,99,999]$ мой перебирает за 105с, этот за 9.7с.
9801 комбинацию $[?,9901,9901]-[?,9999,9999]$ мой перебирает за 224с, этот за 1.0с.
9801 комбинацию $[?,999999901,999999901]-[?,999999999,999999999]$ этот перебирает за 1.2с (тут не хватило точности в 200 значащих цифр, поставил 300, PARI/GP позволяет хоть миллион значащих цифр использовать).
Т.е. зависимость от величины чисел практически отсутствует (в пределах десятка процентов) и в среднем решение занимает порядка 100мкс (3.7ГГц CPU).

 Профиль  
                  
 
 Re: найти корешок по двум элементам цепной дроби
Сообщение08.11.2020, 21:30 
Заслуженный участник
Аватара пользователя


21/11/12
1870
Санкт-Петербург
Dmitriy40 Спасибо! Очень интересно и поучительно.
Фразу ... зависимость от величины чисел практически отсутствует я для себя перевожу так: программа написана толково.

 Профиль  
                  
 
 Re: найти корешок по двум элементам цепной дроби
Сообщение08.11.2020, 22:51 
Заслуженный участник


20/08/14
11067
Россия, Москва

(Оффтоп)

Andrey A в сообщении #1491242 писал(а):
Фразу ... зависимость от величины чисел практически отсутствует я для себя перевожу так: программа написана толково.
Нет, правильным переводом будет "Алгоритм выполняет много сложных операций (корень, цепная и подходящие дроби) и мало простых (арифметика), потому накладные расходы на операции с длинными числами теряются на общем фоне.". Это достоинство не реализации программы, а наличия встроенных функций в PARI/GP. И разумеется Вашего алгоритма, что он не $O(a^n)$, не $O(n^2)$, не $O(n \log n)$, и даже не $O(n)$, а всего лишь $O(\log n)$ или менее (навскидку, просто столько занимают операции с длинными числами, а тут они необходимы в том или ином виде). Т.е. зависимость времени работы от длины чисел присутствует лишь в реализации чисел, но не в алгоритме, он сам по себе занимает близкое к константе время, хотя это чисто эмпирический факт, без доказательства (независимости количества итераций спуска по $m$ от величины чисел). Хорошей аналогией будет операция взятия остатка (дробной части числа), её можно вычислить как последовательность вычитаний с линейной зависимостью от величины чисел, а можно одной операцией деления (с округлением и вычитанием), вообще не зависящей от величины чисел (в частном случае их реализации).
Я программу почти не оптимизировал, больше боролся за простоту кода.

Да, хотелось бы обратить внимание всех, я детально в алгоритме не разбирался, лишь по минимуму для написания работающего кода, потому не понял есть ли нормальное доказательство что он выдаёт именно минимальное значение. По-моему это лишь эмпирический факт, но возможно это гарантируется формулой спуска $m=m-P+Q^2$, не знаю (в остальных местах гарантия есть).

 Профиль  
                  
 
 Re: найти корешок по двум элементам цепной дроби
Сообщение08.11.2020, 23:41 
Заслуженный участник
Аватара пользователя


21/11/12
1870
Санкт-Петербург
Dmitriy40 в сообщении #1491260 писал(а):
... есть ли нормальное доказательство

Доказательства нет и контрпримера нет. Тонкое место.

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 25 ]  На страницу Пред.  1, 2

Модераторы: Модераторы Математики, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group