Реализация рекуррентного соотношения Мюллера

@XZentus · Регистрация: 06.10.2013

Студворк — интернет-сервис помощи студентам

Вопрос про реализацию рекуррентного соотношения Мюллера.
Полная статья здесь.

Вкратце: если считать эту штуку через арифметику с плавающей точкой, довольно быстро получите результат, что это соотношение стремится к 100 или 3. Хотя реально оно стремится к 5.
$https://www.cyberforum.ru/cgi-bin/latex.cgi?f\left(y,z\right) = 108 - \frac{815 - \frac{1500}{z}}{y} {x}_{0} = 4 {x}_{1} = 4.25 {x}_{i} = f({x}_{i-1}, {x}_{i-2})$

Собственно возникла идея попробовать посчитать это дело через Haskell - модуль Data.Ratio.

Программа:

Haskell
1
2
3
4
5
6
7
import Data.Ratio
 
f :: Rational -> Rational -> Rational
f y z = 108 - (815 - 1500 / z) / y
 
xlist :: [Rational]
xlist = 4 : (17 % 4) : [f y x | x <- xlist, y <- tail xlist]

Собственно, то ли я что не так написал, то ли Data.Ratio что-то не так считает, но у меня выдает следующие результаты:

Haskell
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
λ:> mapM putStrLn $ map show $ zip [1..] $ take 30 xlist
(1,4 % 1)
(2,17 % 4)
(3,76 % 17)
(4,182 % 19)
(5,5648 % 91)
(6,71243 % 706)
(7,7383604 % 71243)
(8,191520578 % 1845901)
(9,9936012992 % 95760289)
(10,128869359497 % 1242001624)
(11,13371410111116 % 128869359497)
(12,346852443455462 % 3342852527779)
(13,17994604390483568 % 173426221727731)
(14,233388717076502963 % 2249325548810446)
(15,24216278202785723764 % 233388717076502963)
(16,628166752596799215698 % 6054069550696430941)
(17,32589109339073942840672 % 314083376298399607849)
(18,422678390381086249917377 % 4073638667384242855084)
(19,43856865147508248134839756 % 422678390381086249917377)
(20,1137640736040803212149761942 % 10964216286877062033709939)
(21,59020472163090419808670958288 % 568820368020401606074880971)
(22,765491253960598579082939483483 % 7377559020386302476083869786)
(23,79426929458774673451480561510324 % 765491253960598579082939483483)
(24,2060323057451250339490851817595618 % 19856732364693668362870140377581)
(25,106888963982134911292674567267095552 % 1030161528725625169745425908797809)
(26,1386342129678911918115108233121910457 % 13361120497766863911584320908386944)
(27,143846056986305067035334587977476073996 % 1386342129678911918115108233121910457)
(28,3731345904365556498961371969748443847622 % 35961514246576266758833646994369018499)
(29,193581145701493272256970684027654783701808 % 1865672952182778249480685984874221923811)
(30,2510733454600106371747666505205257374164803 % 24197643212686659032121335503456847962726)

Т.е. вычисления скатились к 103.

Аналогичных проблем на С++ (использовал boost) не замечено.
Прошу просмотреть код и указать на косяки.

@korvin_ · 30.05.2015, 16:06

Сообщение от XZentus

Аналогичных проблем на С++ (использовал boost) не замечено.

В C++ ты тоже использовал ленивые списки и list comprehension?

http://ideone.com/QMcRZj

Haskell
1
2
3
4
5
6
7
8
9
10
11
12
import Data.Ratio
import Control.Monad
 
f :: Rational -> Rational -> Rational
f y z = 108 - (815 - 1500 / z) / y
 
x :: Int -> Rational
x 0 = 4
x 1 = 17 % 4
x i = f (x (i-1)) (x (i-2))
 
main = forM_ [0..29] $ \i -> print (i, x i)

Bash
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
(0,4 % 1)
(1,17 % 4)
(2,76 % 17)
(3,353 % 76)
(4,1684 % 353)
(5,8177 % 1684)
(6,40156 % 8177)
(7,198593 % 40156)
(8,986404 % 198593)
(9,4912337 % 986404)
(10,24502636 % 4912337)
(11,122336033 % 24502636)
(12,611148724 % 122336033)
(13,3054149297 % 611148724)
(14,15265963516 % 3054149297)
(15,76315468673 % 15265963516)
(16,381534296644 % 76315468673)
(17,1907542343057 % 381534296644)
(18,9537324294796 % 1907542343057)
(19,47685459212513 % 9537324294796)
(20,238423809278164 % 47685459212513)
(21,1192108586037617 % 238423809278164)
(22,5960511549128476 % 1192108586037617)
(23,29802463602463553 % 5960511549128476)
(24,149012035582781284 % 29802463602463553)
(25,745059330625296977 % 149012035582781284)
(26,3725294111260656556 % 745059330625296977)
(27,18626462930705797793 % 3725294111260656556)
(28,93132291776736534004 % 18626462930705797793)
(29,465661390253305305137 % 93132291776736534004)

@XZentus · 30.05.2015, 16:18 **[ТС]**

Сообщение от korvin_

В C++ ты тоже использовал ленивые списки и list comprehension?

Нет, там решение "в лоб".

Видимо, я что-то недопонимаю в list comprehensions...
Где именно проблема в моём варианте?

@korvin_ · 30.05.2015, 16:34

Сообщение от XZentus

я что не так написал

Да, почитай про синтаксис list comprehensions. http://ideone.com/kF7R9v

Haskell
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
{-# LANGUAGE ParallelListComp #-}
 
import Data.Ratio
import Control.Monad
 
f :: Rational -> Rational -> Rational
f y z = 108 - (815 - 1500 / z) / y
 
x :: Int -> Rational
x 0 = 4
x 1 = 17 % 4
x i = f (x (i-1)) (x (i-2))
 
xlist :: [Rational]
xlist = 4 : 17 % 4 : zipWith f (tail xlist) xlist
 
xlist2 :: [Rational]
xlist2 = 4 : 17 % 4 : [ f y z | y <- tail xlist2 | z <- xlist2 ]
 
main = forM_ (zip3 [0..29] xlist xlist2) $ \(i, xitem, xitem2) -> do
    print (i, x i)
    print (i, xitem)
    print (i, xitem2)
    putStrLn ""

Bash
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
(0,4 % 1)
(0,4 % 1)
(0,4 % 1)
 
(1,17 % 4)
(1,17 % 4)
(1,17 % 4)
 
(2,76 % 17)
(2,76 % 17)
(2,76 % 17)
 
(3,353 % 76)
(3,353 % 76)
(3,353 % 76)
 
(4,1684 % 353)
(4,1684 % 353)
(4,1684 % 353)
 
(5,8177 % 1684)
(5,8177 % 1684)
(5,8177 % 1684)
 
(6,40156 % 8177)
(6,40156 % 8177)
(6,40156 % 8177)
 
(7,198593 % 40156)
(7,198593 % 40156)
(7,198593 % 40156)
 
(8,986404 % 198593)
(8,986404 % 198593)
(8,986404 % 198593)
 
(9,4912337 % 986404)
(9,4912337 % 986404)
(9,4912337 % 986404)
 
(10,24502636 % 4912337)
(10,24502636 % 4912337)
(10,24502636 % 4912337)
 
(11,122336033 % 24502636)
(11,122336033 % 24502636)
(11,122336033 % 24502636)
 
(12,611148724 % 122336033)
(12,611148724 % 122336033)
(12,611148724 % 122336033)
 
(13,3054149297 % 611148724)
(13,3054149297 % 611148724)
(13,3054149297 % 611148724)
 
(14,15265963516 % 3054149297)
(14,15265963516 % 3054149297)
(14,15265963516 % 3054149297)
 
(15,76315468673 % 15265963516)
(15,76315468673 % 15265963516)
(15,76315468673 % 15265963516)
 
(16,381534296644 % 76315468673)
(16,381534296644 % 76315468673)
(16,381534296644 % 76315468673)
 
(17,1907542343057 % 381534296644)
(17,1907542343057 % 381534296644)
(17,1907542343057 % 381534296644)
 
(18,9537324294796 % 1907542343057)
(18,9537324294796 % 1907542343057)
(18,9537324294796 % 1907542343057)
 
(19,47685459212513 % 9537324294796)
(19,47685459212513 % 9537324294796)
(19,47685459212513 % 9537324294796)
 
(20,238423809278164 % 47685459212513)
(20,238423809278164 % 47685459212513)
(20,238423809278164 % 47685459212513)
 
(21,1192108586037617 % 238423809278164)
(21,1192108586037617 % 238423809278164)
(21,1192108586037617 % 238423809278164)
 
(22,5960511549128476 % 1192108586037617)
(22,5960511549128476 % 1192108586037617)
(22,5960511549128476 % 1192108586037617)
 
(23,29802463602463553 % 5960511549128476)
(23,29802463602463553 % 5960511549128476)
(23,29802463602463553 % 5960511549128476)
 
(24,149012035582781284 % 29802463602463553)
(24,149012035582781284 % 29802463602463553)
(24,149012035582781284 % 29802463602463553)
 
(25,745059330625296977 % 149012035582781284)
(25,745059330625296977 % 149012035582781284)
(25,745059330625296977 % 149012035582781284)
 
(26,3725294111260656556 % 745059330625296977)
(26,3725294111260656556 % 745059330625296977)
(26,3725294111260656556 % 745059330625296977)
 
(27,18626462930705797793 % 3725294111260656556)
(27,18626462930705797793 % 3725294111260656556)
(27,18626462930705797793 % 3725294111260656556)
 
(28,93132291776736534004 % 18626462930705797793)
(28,93132291776736534004 % 18626462930705797793)
(28,93132291776736534004 % 18626462930705797793)
 
(29,465661390253305305137 % 93132291776736534004)
(29,465661390253305305137 % 93132291776736534004)
(29,465661390253305305137 % 93132291776736534004)

@Catstail · 30.05.2015, 17:02

Решил проверить это код на Лиспе:

Lisp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
(defun muller (n &optional (c 17/4) (p 4))
  (printline (rat2flo c)) 
  (if (zerop n) c (muller (- n 1) (- 108 (/ (- 815 (/ 1500 p)) c)) c)))
 
==> muller
(muller 30)
4.25
4.47058823529412
4.64473684210526
4.77053824362606
4.85570071258907
4.91084749908279
4.94553740412392
4.9669625817627
4.98004570135563
4.98797944847839
4.99277028806207
4.99565589150663
4.99739126838134
4.99843394394482
4.99906007197089
4.99943593714684
4.99966152410377
4.99979690071342
4.99987813547793
4.9999268795046
4.99995612706116
4.99997367600571
4.99998420552027
4.99999052328223
4.99999431395856
4.99999658837126
4.99999795302136
4.99999877181231
4.99999926308721
4.99999955785226
4.99999973471133
 
==> 11641533109203575871233/2328306745375394431036

@Curry · 30.05.2015, 18:31

Я просто добавил наглядности.

Haskell
1
2
3
4
5
6
7
8
9
10
11
import Data.Ratio
import Control.Monad
 
f :: Rational -> Rational -> Rational
f y z = 108 - (815 - 1500 / z) / y
 
xlist :: [Rational]
xlist = 4 : 17 % 4 : zipWith f (tail xlist) xlist
 
main = forM_ (zip [0..29] xlist) $ \(i, xitem) -> 
    print (i, xitem, fromRational xitem :: Double)

@XZentus · 30.05.2015, 19:06 **[ТС]**

Разобрался, всем спасибо.

Добавлено через 1 минуту
Catstail, Очевидно, при наличии поддержки Rational, проблем быть по определению не должно, раз таз там округления нет.

Добавлено через 13 минут
Да, кстати, любопытная статистика:

Haskell:

Кликните здесь для просмотра всего текста

Haskell
1
2
3
4
5
6
7
8
9
10
11
12
13
14
import Data.Ratio
import System.Environment
 
f y z = 108 - (815 - 1500 / z) / y
 
xlist = 4 : (17 % 4) : zipWith f (tail xlist) xlist
 
main = do
         args <- getArgs
         case args of
                    [arg1] -> do
                                let x = read arg1
                                putStrLn $ show $ xlist !! x
                    _      -> error "ERR: Main.Args"

Bash
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
$ ghc -O2 -Wall main.hs
[ворнинги]
$ time -l ./main 10000
[вывод программы]
        6,82 real         6,82 user         0,00 sys
     28748  maximum resident set size
       832  average shared memory size
       100  average unshared data size
       128  average unshared stack size
      6466  page reclaims
         0  page faults
         0  swaps
         0  block input operations
         0  block output operations
         0  messages sent
         0  messages received
       682  signals received
         4  voluntary context switches
        73  involuntary context switches

C++

Кликните здесь для просмотра всего текста

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
#include <iostream>
#include <limits>
#include <cmath>
#include <string>
#include <cstdlib>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/rational.hpp>
 
using namespace boost;
 
typedef boost::multiprecision::cpp_int LL;
 
const char MSG_ERR[] = "Error\n";
const size_t MSG_ERR_LEN = std::strlen(MSG_ERR);
 
rational<LL> *x;
 
rational<LL> fun(const rational<LL> & y, const rational<LL> & z) {
    rational<LL> ret = 108 - (815 - (1500 / z))/y;
    return ret;
}
 
int main(int argn, char *argv[]) {
    if(argn != 2) {
        std::cout.write(MSG_ERR, MSG_ERR_LEN);
        return 1;
    }
    int n = std::atoi(argv[1]);
    if(n < 1) {
        std::cout.write(MSG_ERR, MSG_ERR_LEN);
        return 1;
    }
    x = new rational<LL>[n];
    x[0] = 4;
    x[1] = rational<LL>(17,4);
    for(int i = 2; i < n; ++i)
        x[i] = fun(x[i-1], x[i-2]);
    std::cout << x[n-1] << std::endl;
    return 0;
}

Bash
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
$ clang++36 -Wall -Wextra -O3 -I/usr/local/include/ main.cpp
[ворнинги из-за буста]
$ time -l ./a.out 10000
[вывод программы]
       33,64 real        33,63 user         0,00 sys
     43764  maximum resident set size
        63  average shared memory size
         3  average unshared data size
       127  average unshared stack size
     10397  page reclaims
         0  page faults
         0  swaps
         0  block input operations
         0  block output operations
         0  messages sent
         0  messages received
         0  signals received
         4  voluntary context switches
       395  involuntary context switches

Итого без ручных оптимизаций Хаскель быстрее почти в 5 раз.

@Curry · 30.05.2015, 20:57

Сообщение от XZentus

Итого без ручных оптимизаций Хаскель быстрее почти в 5 раз.

Что то, как то ... хм. А за счёт чего?

@XZentus · 30.05.2015, 21:20 **[ТС]**

KolodeznyDiver, возможно, из-за сочетания rational + cpp_int.

Добавлено через 13 минут
В boost есть ещё boost::multiprecision::cpp_rational. В общем то логичнее использовать его, вместо вложенной конструкции rational<cpp_int>. Но почему-то это ещё более тормозной вариант. Те же 10000 у меня он считал больше 2 минут, до того как я его выключил.

Если что, вот код:

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#include <iostream>
#include <cstdlib>
#include <boost/multiprecision/cpp_int.hpp>
//#include <boost/rational.hpp>
 
typedef boost::multiprecision::cpp_rational LL;
 
const char MSG_ERR[] = "Error\n";
const size_t MSG_ERR_LEN = std::strlen(MSG_ERR);
 
LL *x;
 
LL fun(const LL & y, const LL & z) {
    return 108 - (815 - (1500 / z))/y;
}
 
int main(int argn, char *argv[]) {
    if(argn != 2) {
        std::cout.write(MSG_ERR, MSG_ERR_LEN);
        return 1;
    }
    int n = std::atoi(argv[1]);
    n += 1; //досадная ошибка, из-за которой выводился предыдущий элемент. Но на общее время для кол-ва элементов порядка 10000 влияет в пределах погрешности.
    if(n < 1) {
        std::cout.write(MSG_ERR, MSG_ERR_LEN);
        return 1;
    }
    x = new LL[n];
    x[0] = 4;
    x[1] = LL(17,4);
    for(int i = 2; i < n; ++i)
        x[i] = fun(x[i-1], x[i-2]);
    std::cout << x[n-1] << std::endl;
    return 0;
}

Новые блоги и статьи Все статьи Все блоги /
SDL3 для Desktop (MinGW): Создаём пустое окно с нуля для 2D-графики на SDL3, Си и C++ 8Observer8 10.03.2026 Содержание блога Финальные проекты на Си и на C++: hello-sdl3-c. zip hello-sdl3-cpp. zip Результат:	Установка CMake и MinGW 13.1 для сборки С и C++ приложений из консоли и из Qt Creator в EXE 8Observer8 10.03.2026 Содержание блога MinGW - это коллекция инструментов для сборки приложений в EXE. CMake - это система сборки приложений. Здесь описаны базовые шаги для старта программирования с помощью CMake и. . .	Как дизайн сайта влияет на конверсию: 7 решений, которые реально повышают заявки Neotwalker 08.03.2026 Многие до сих пор воспринимают дизайн сайта как “красивую оболочку”. На практике всё иначе: дизайн напрямую влияет на то, оставит человек заявку или уйдёт через несколько секунд. Даже если у вас. . .	Модульная разработка через nuget packages DevAlt 07.03.2026 Сложившийся в . Net-среде способ разработки чаще всего предполагает монорепозиторий в котором находятся все исходники. При создании нового решения, мы просто добавляем нужные проекты и имеем. . .
Модульный подход на примере F# DevAlt 06.03.2026 В блоге дяди Боба наткнулся на такое определение: В этой книге («Подход, основанный на вариантах использования») Ивар утверждает, что архитектура программного обеспечения — это структуры,. . .	Управление камерой с помощью скрипта OrbitControls.js на Three.js: Вращение, зум и панорамирование 8Observer8 05.03.2026 Содержание блога Финальная демка в браузере работает на Desktop и мобильных браузерах. Итоговый код: orbit-controls-threejs-js. zip. Сканируйте QR-код на мобильном. Вращайте камеру одним пальцем,. . .	SDL3 для Web (WebAssembly): Синхронизация спрайтов SDL3 и тел Box2D 8Observer8 04.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-sync-physics-sprites-sdl3-c. zip На первой гифке отладочные линии отключены, а на второй включены:. . .	SDL3 для Web (WebAssembly): Идентификация объектов на Box2D v3 - использование userData и событий коллизий 8Observer8 02.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-collision-events-sdl3-c. zip Сканируйте QR-код на мобильном и вы увидите, что появится джойстик для управления главным героем. . . .