델파이에서 파이(Pi) 더 많은 숫자값 확인하는 방법

[RAD]

• 원문제목: Getting Big with Pi in Delphi
• 원문: https://blogs.embarcadero.com/getting-big-with-pi-in-delphi/
• 작성자: Jim McKeeth (2023.4)

 

When I was very young and first learned about Pi they told me 3.14 was a good approximation, but it was an irrational number that went on forever, and people with computers were able to calculate many digits. I had a computer (Commodore Vic 20) and wanted to see how many digits of Pi I could calculate. This is when I discovered floating point number precision limitations, which foiled my first attempt. The idea of calculating more digits of Pi has always itched in the back of my mind, but CPU based floating point precision was never enough. 아주 어릴 적 파이(Pi)를 처음 배웠을 때, 선생님들은 나에게 3.14가 근사치라고 가르쳐줬지만 사실은 무한하게 이어지는 비이성적인 수였다. 그리고 컴퓨터를 이용해 사람들은 많은 자릿수까지 계산할 수가 있었다. 당시 나는 컴퓨터 (Commodore Vic 20)가 있었고, 파이를 몇 자리까지 계산할 수 있는지 확인해보고 싶었다. 이 때 부동 소수점 숫자 정밀도의 한계를 발견했고, 첫 번째 시도는 실패로 돌아갔다. 더 많은 자릿수까지 파이를 계산하고 싶다는 생각은 늘 마음 한구석에 있었지만, CPU 기반의 부동 소수점 정밀도로는 충분하지가 않았다.

In my 2023 Pi day post I experimented with Rudy Velthuis’ arbitrary precision big numbers library and still only got 20 digits, but I wanted more. As you may know, Rudy Velthuis passed away a few years ago, so his library has remained unmaintained for four years. There were a few forks, and a pull request, but no new canonical home of his library. First thing I did was merge all the pull requests and forks into a new fork in the Turbo Pack organization on GitHub, then I collected up Rudy’s articles on his libraries, and added them to the wiki, and did some additional general clean up and updating. I worked with Roman Kassebaum, the maintainer of Turbo Pack, to make sure it built and installed on more than just my machine. 2023년 파이의 날을 기념한 글에서 루디 벨투이스(Rudy Velthuis)의 임의 정밀도 라이브러리를 사용해 파이 값을 계산해보았는데, 지금도 20자리까지밖에 확인할 수 없었고 난 더 많은 자릿수를 확인해보고 싶었다. 알다시피 루디 벨투이스는 몇 년 전 별세하신 이후, 그가 만든 라이브러리는 4년간 업데이트되지 않고 있다. 몇 번의 fork와 pull request가 있었지만, 그 라이브러리 관련 새로운 공식 홈은 없는 상황이었다. 내가 가장 먼저 했던 일은 모든 pull request와 fork를 깃허브(GitHub)의 Turbo Pack 조직에 있는 새로운 fork에 병합하는 것이었다. 그리고 이를 위키에 추가한 뒤, 추가로 정리 및 업데이트를 진행했다. 이 작업을 로만(Roman Kassebaum)과 같이 진행했는데, 그는 Turbo Pack의 유지보수 담당자로, 내 컴퓨터 뿐 아니라 다른 환경에도 해당 라이브러리를 빌드하고 설치할 수 있도록 작업했다.

[루디의 빅 넘버 라이브러리를 위한 새 공식 홈.
github.com/TurboPack/RudysBigNumbers ]

 

Once I was sure Rudy’s Big Numbers libraries were working correctly, I went about rewriting my Pi calculation implementations. I used ChatGPT and Bing at first, but they never got as much accuracy as I wanted. Instead I ended up searching for more details and examples in other languages, and attempted to translate them to Delphi. I ended up with two results, both calculate a verified 100k digits of Pi, and with a respectable speed.

루디의 빅 넘버 라이브러리가 제대로 작동한다는 것을 확인을 한 뒤에는, 파이 계산 구현을 다시 작성했다. 처음에는 ChatGPT와 Bing을 사용했지만, 내가 원하는 만큼의 정확도가 나오지 않았다. 대신 다른 언어로 된 더 자세한 내용과 예제를 찾아서 델파이로 변경해보기로 했다. 결론적으로 두 가지 결과를 얻었는데, 두 결과 모두 10만 자리까지 계산됐고 속도 또한 굉장히 빨랐다.

 

 

Chudnovsky와 BigDecimals

The Chudnovsky algorithm was published by the Chudnovsky brothers in 1988. It was used in the world record calculation of 100 trillion digits of Pi on March 21, 2022.

추드노브스키(Chudnovsky) 알고리즘은 1988년 추드노브스키 형제가 발표한 것으로, 2022년 3월 21일에는 파이의 100조 자리까지 계산해내는 세계 기록 신기록을 세우는데 사용되기도 했다. 

 

You specify the number of places you want it to calculate too, and then it starts with an estimate, and continues improving it until there are no changes out to that many decimal places. Internally is uses an additional six digits of precision, but then rounds it down to the specified precision once it settles on a number. The 32-bit version ran out of memory, but with 64-bit Windows version I reached 100K digits.

계산하려는 자릿수를 지정하면, 프로그램은 이에 따른 추정치를 시작해 해당 자릿수까지 계산하면서 동일한 결과가 나올 때까지 개선해나아간다. 내부적으로는 추가 6자리의 정밀도를 확인하며, 최종 결과에서는 지정된 정밀도에 맞추기 위해 반올림한다. 32-bit 버전은 메모리가 부족해 사용할 수 없었지만, 64-bit 윈도우 버전에서는 10만 자릿수까지 결과를 얻어낼 수 있었다.

function Chudnovsky(Places: Integer):BigDecimal;
begin
  // Use +6 internally for calculations
  var internalPrecision := MaxInt;
  if Places <= MaxInt - 6 then
    internalPrecision := Places + 6;
 
  var lastSum: BigDecimal;
  var t := BigDecimal.Create(3);
  var sum := BigDecimal.Create(3);
  sum.DefaultPrecision := internalPrecision;
  sum.DefaultRoundingMode := rmFloor;
  var n := BigInteger.One;
  var d := BigInteger.Zero;
  var na: UInt64 := 0;
  var da: UInt64 := 24;
  while sum <> lastSum do
  begin
    lastSum := sum;
    n := n + na;
    na := na + 8;
    d := d + da;
    da := da + 32;
    t := ((t * n)/d);
    sum := (sum + t).RoundToPrecision(internalPrecision);
  end;
  Result := sum.RoundToPrecision(Places);
end;

 

Rudy’s implementation of BigDecimals points out it isn’t very efficient, so I was curious if I could calculate faster by switching to BigIntegers.

루디가 구현한 BigDecimals는 그다지 효율적이지 않다는 평가가 있었으므로, BigIntegers로 전환하면 더 빠르게 계산할 수 있는지 궁금했다.

 

베일리-보베인-플루프(Bailey-Borwein_Plouffe) 그리고 BigIngeters

Discovered in 1995, the BBP algorithm works with integers to generate successive digits of pi. So it doesn’t calculate Pi as a whole, but just each individual digit. My implementation adds those digits to an array.

1995년 선보인 BBP 알고리즘은 정수를 사용해 파이의 연속된 자릿수를 생성하는 알고리즘이다. 따라서 해당 알고리즘은 파이 전체를 계산하는 것이 아니라 각 개별 숫자만 계산하는 것이다. 내 구현에서는 이와 같은 자릿수를 배열에 추가한 것다.

function BBPpi(Places: UInt64): TDigits;
// Bailey-Borwein-Plouffe
begin
  SetLength(Result, Places);
 
  var idx: Uint64 := 0;
  var q := BigInteger.One;
  var r := BigInteger.Zero;
  var t := BigInteger.One;
  var k := BigInteger.One;
  var n := BigInteger.Create(3);
  var l := BigInteger.Create(3);
 
  while true do
  begin
    if 4*q+r-t < n*t then
    begin
      result[idx] := n.AsInt64; // It is just a byte
      inc(idx);
      //write(n.ToString[1]);
      if idx >= places then break;
      var newR := 10 * (r - n * t);
      n := ((10 * (3 * q + r)) div t) - 10 * n;
      q := q * 10;
      r := newR;
    end
    else
    begin
      var newR := (2 * q + r) * l;
      var newN := (q * (7 * k)+2+(r * l)) div (t * l);
      q := q * k;
      t := t * l;
      l := l + 2;
      k := k + 1;
      n := newN;
      r := newR;
    end;
  end;
end;

You can easily modify it to just output each digit as it goes, instead of storing it in the array. If you do then you will see the rate of generation slows down the deeper you go. I thought about adding a call back for each digit, but more likely I will modify it to generate digits in a range, so you don’t need to start at the top each time. But I’ll leave that for another day.

배열에 저장하지 않고 각 자릿수를 그대로 출력하도록 변경하는 것은 쉽다. 그렇게 하면, 더 깊이 들어갈수록 생성 속도는 느려지는 것을 확인할 수 있다. 각 자릿수에 대해 콜백을 추가하는 것도 고려해봤지만, 범위 내에서 자릿수를 생성하도록 수정하는 것이 더 적합하다고 보았다. 하지만 이 부분은 일단 다음 기회로 미루기로 했다.

 

연산 시간

Performance will vary based on your computer hardware, but I wanted to see how it scaled and which algorithm was faster. I ran each on Win32, Win64, and Linux64, using the same hardware (Linux was under Ubuntu with WSL2).

성능은 컴퓨터 하드웨어에 따라 다르겠지만, 확장 방식과 어떤 알고리즘이 더 빠른지를 확인해보고 싶었다. 동일한 하드웨어에서 Win32, Win64, Linux64에서 각각 실행해 비교해보았다 (리눅스는 WSL2가 설치된 우분투에서 사용해보았다).

Win64 was the fastest, with BBP implementation using BigIntegers being the leader.

Win64가 가장 빨랐는데, 그 중에서도 BigIntegers를 사용한 BBP구현이 선두를 이끌었다.

 

모바일에서의 파이 연산

Sometimes you need to calculate a few thousand digits of Pi while you are away from your computer. So to make sure this all works under Android, I created a FMX version. It was actually faster on Android than I expected, but I didn’t bother benchmarking it or seeing how many digits it could generate. In theory it is only limited by your available memory.

가끔은 컴퓨터와 떨어져서 몇 천 자리의 파이 값을 계산해야 할 때가 있다. 그래서 안드로이드에서도 이 모든 작업이 제대로 작동하는지 확인해보기 위해 FMX 버전을 만들었다. 실제로 안드로이드에서의 연산은 생각했던 것보다 더 빨랐지만, 벤치마킹 여부나 몇 자리까지 생성할 수 있는지는 확인하지 않았다. 이론적으로는 사용 가능한 메모리에 따라서 달라지게 될 것이다.

          

 

 

코드 확인하기

I added all of my code to the samples folder on GitHub. There are also unit tests that verify both algorithms work as expected. It uses a series of hashes to test it to various numbers of digits, as well as a a saved 100K digits of Pi. If you don’t want to install from GitHub then stay tuned to GetIt where we will have an installer soon.

작성한 모든 코드는 깃허브 샘플 폴더에 올려놓았다. 두 알고리즘이 예상대로 작동하는지를 확인할 수 있는 단위 테스트도 깃허브에서 확인할 수 있다. 해당 코드에서는 여러 자릿수 테스트를 위한 일련의 해시(hash)를 사용하고, 100,000  자리의 파이를 저장한다. 깃허브에서 인스톨을 하고 싶지 않다면, 조만간 깃허브에서 인스톨러를 제공할 예정이니 조금만 기다려주기를 바란다.

 

 

파이 외 그 이상

Once I figured out Rudy’s libraries I wanted to try other irrational numbers besides Pi. They were all pretty easy to generate, but Euler’s number required a little work.

루디의 라이브러리를 이해하고 나니, 파이 외에 다른 무리수도 시도해보고 싶었다. 다른 수들은 그래도 생성하기 꽤 쉬웠지만, 오일러(Euler) 수는 작업을 조금 해야했다.

Precision: 200
√2, Square root of 2, Pythagoras constant:
1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727
√3, Square root of 3, Theodorus' constant:
1.7320508075688772935274463415058723669428052538103806280558069794519330169088000370811461867572485756
√5, Square root of 5:
2.2360679774997896964091736687312762354406183596115257242708972454105209256378048994144144083787822749
φ, Phi, Golden ratio, (1 + √5)/2:
1.61803398874989484820458683436563811772030917980576286213544862270526046281890244970720720418939113745
Common logarithm of 2:
0.301029995663981
Natural logarithm of 2:
0.693147180559945
------
Euler's Number to 500 digits
2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952605956307381323286279434907632338298807531952510190115738341879307021540891499348841675092447614606680822648001684774118537423454424371075390777449920695517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416928368190255151086574637721112523897844250569536967707854499699679468644549059879316368892300987931

 

After a few different algorithms and multiple iterations I ended up with the following for Euler’s number

몇 가지 여러 알고리즘과 여러번의 반복 끝에 오일러 수에 대해 다음과 같은 결과를 얻을 수 있었다.

 

function Euler(const Digits: Integer): BigDecimal;
begin
  Result := 1;
  var InternalPrecision := MaxInt;
  if Digits <= MaxInt - 2 then
    InternalPrecision := Digits + 2;
 
  Result.DefaultPrecision :=InternalPrecision;
  var lastFactorial := BigInteger.One;
  var lastResult := BigDecimal.Zero;
  var iteration: UInt64 := 1;
  while true do
  begin
    lastResult := Result;
    lastFactorial := lastFactorial * iteration;
    Result := (Result + BigDecimal.One / lastFactorial).RoundToPrecision(InternalPrecision);
    if lastResult = Result then break;
    inc(Iteration);
  end;
  //Writeln(Format('%d digits took %d iterations',[Digits, iteration]));
  Result.DefaultRoundingMode := rmFloor;
  Result := Result.RoundToPrecision(Digits);
end;

 

The process is similar to my Chudnovsky implementation in that it continues to improve it’s estimate until there are no changes within the specified number of digits. The formula does use factorials, but does them in sequence, which is much faster. I verified it to 10K digits as well.

이 프로세스는 지정된 자릿수 내에서 더 변경 사항이 없을 때까지 추정치를 계속해서 개선해 나아가는 방식이라는 점에서 Chudnovsky 구현과 유사하다. 이 공식은 계승(factorials)을 사용하지만 순차적으로 수행하기 때문에 더 빠르다. 10,000자리까지는 직접 검증을 완료했다.