y-cruncher - Language and Algorithms

By Alexander J. Yee

(Last updated: August 6, 2011)

Back To:

• Numberworld Home
• y-cruncher Home

Language:

• As of v0.5.5, y-cruncher is written in C with Intel SSE/AVX intrinsics.
  • The Windows versions use WinAPI for multithreading and I/O.
  • The Linux versions use OpenMP for multi-threading and a mix of standard C and Linux-specific libraries for I/O.
  • A small amount of inline assembly is used in the Linux versions to access the "cpuid" and "rdtsc" instructions. No assembly is used for anything else.
  • v0.6.x and later will use more inline assembly.
  • y-cruncher is compiled as C in Linux and C++ in Windows. (So it will compile as either C or C++.)

   

• y-cruncher is written completely from ground up using nothing but the most essential libraries.
  • Libraries:
    • Memory: malloc(), free(), memcpy(), memset()
    • Threads: WinAPI Threads + Mutex, OpenMP
    • Time: C "time.h", WinAPI Performance Counter, GNU "sys/time.h"
    • Math: C "math.h", rand()/srand()
    • I/O: C wprintf(), FILE, GNU file descriptors, C++ Streams (only for convenient debugging, there is no C++ in release builds)
      
      
  • Language Extensions:
    • Wide-Character/Unicode: wchar_t
    • C99 Extensions: "inline", intermingled declarations/code
    • 64-bit Integers: __int64, long long
    • Misc. Compiler Intrinsics
    • SSE and AVX intrinsics
      
      
• For performance reasons, y-cruncher does not use any 3rd-party libraries.
  There exist very well optimized libraries that perform the subroutines needed by y-cruncher. But its usually possible to do better with specialized code.
  
  • y-cruncher's custom FFT is up to 3x faster than FFTW 3.2.2 for the integer convolution operations that it needs.
    This is because FFTW does not support out-of-order FFTs which are faster than normal FFTs.
    
    
  • y-cruncher's custom arithmetic library can be more than 10x faster than GMP 5.0.1. (8-core Xeon server running Linux)
    This is because GMP focuses on small numbers - it was never meant for the gigantic sizes that y-cruncher needs. And as such, GMP lacks support for multi-threading. It also has very poor support for Windows.
    
    
    
• Virtually all arithmetic and low-level code is done using in-place functional programming. High-level code is done largely using lightweight objects.

• y-cruncher is compiled using different compilers.
  • Windows Versions: (y-cruncher v0.4.3 - v0.5.5)
    • x86 - Microsoft Visual Studio 2008
    • All versions with SSE - Intel C++ Compiler 11.0
      
      
  • Windows Versions: (y-cruncher v0.5.5 (fix 2) - current)
    • x86 - Microsoft Visual Studio 2010 SP1
    • All versions with SSE - Intel C++ Compiler 11.0
    • x64 AVX ~ Hina - Microsoft Visual Studio 2010 SP1
      
      
  • Linux Versions: (y-cruncher v0.5.4 - current)
    • x64 SSE3 - GNU C Compiler 4.4.3
      
      
• Note that I have recently received a lot of criticism for using the Intel Compiler because of its purported "unfairness" to non-Intel processors.
  Regardess of whether or not this is indeed the case, the Intel Compiler compiles the fastest code among all the compilers that I have tried - even on AMD processors. Therefore, I have no intention of changing compilers for now.

   

• y-cruncher will remain closed-source. For now at least... So please stop asking me.

• 99% of the source code is portable (standard C)
  • There is a special interface/library for any inherently non-portable function (such as multi-threading, I/O, and timers). Each of these interfaces contain multiple implementations specialized for whatever environment it will be compiled on. This allows the rest of y-cruncher to be fully portable.
    • Multi-threading (4 versions):
      • WinAPI threads- For Windows compilations.
      • Pthreads - For Linux and Unix compilations. (coming soon - not yet implemented)
      • OpenMP - For everything else. (The Linux version currently uses this.)
      • Multi-threading can also be completely disabled. In this case, the library will simply sequentialize everything. This will compile without OpenMP.
    • File I/O:
      • WinAPI I/O - For Windows compilations.
      • POSIX I/O - For Linux compilations.
      • Standard C - For everything else. However, there may be problems with support for large files (> 4 GB).
    • Timers:
      • WinAPI performance counter - For Windows compilations.
      • <sys/time.h> - For Linux compilations.
      • <time.h> - For everything else.
        
        
  • The only possible major wall for portability would be on systems that use big-endian instead of little-endian. y-cruncher does not use a fixed word-size. Some functions use 32-bit words, some 64-bit, and even 128-bit/256-bit SIMD. On little-endian machines, all word-sizes are interchangable as long as they are aligned properly. But this is not the case on big-endian machines - which will require the addition of swapping operations with additional overhead.
    
    
• Starting from v0.4.4, all new code that is vectorizable is written in a "vector-generic" manner.
  • Allows the code to be easily ported to any vector size.
  • No vectorization (via standard C/C++) is still supported as this is simply a vector of size 1.
  • Allows the code to be easily ported to any vector architecture that supports the basic functions. (packed addition, subtraction, and multiplication)

• y-cruncher's design is split into 3 major layers. There's an additional interface layer at the very top, but that's not really worth mentioning.
  
  
  • Math Layer: This is where all the constants and functions are implemented. This is done largely using very simple object-oriented programming structured in a mathematical way that most resembles the formulas/algorithms that are involved.
    Roughly 20% of the source code is in this layer.
    
    
  • Object Layer: This is where all the "number" objects are defined. (Large integer, large floating-point, etc...) This layer most resembles a C++/Java class library but without any fancy stuff like inheritance. (everything compiles in C)
    Roughly 2% of the source code is in this layer.
    
    
  • Modules Layer: Here's where it gets serious. More than 70% of the entire source code of y-cruncher resides here. The Modules layer implements all the core operations needed by the upper layers. This is where most of the performance critical code is. All that nasty stuff like FFTs, CPU-dispatching, loop-unrolling, etc... is all down here. Unlike the Math and Object layers, there is no structure in the Modules layer. It's exactly as the name implies: A collection of (mostly) independent modules that implement different low-level operations needed by y-cruncher.
    
    Together, these modules implement an interface of about 100 or so different functions. As a normal part of development, new modules are constantly being added and old ones phased out and retired. So at any one time, there are typically multiple modules that implement the same function(s). To prevent naming conflicts, name-mangling is used on all symbols and resolved by a compile-time version dispatcher.

Code-Size: Lines of Code in y-cruncher:

Although y-cruncher is completely close-sourced (and will likely remain that way for a long time), I'm fine with revealing some of the structure to satisfy the curious.

Below is the approximate code-size distribution of the y-cruncher project (July 2011).

Some of the numbers are more up-to-date than others. But for the most part, they are relatively accurate.

Notes:

• The line counts include inline commenting and external text-file documentation.
  More advanced documentation formats in the form of Excel spreadsheets and Mathematica ".nb" files are not included.

   

• There is a lot of code-replication among the different constants. All the algorithms that share the same Binary Splitting structure have the same near-identical code. However, slight differences in each algorithm prevent most of the code from being reused.
  
  In particular, the following algorithms all have the same 3-variable geometrically convergent Binary Splitting structure:
  • Pi - Chudnovsky
  • Pi - Ramanujan
  • Integer ArcCoth (used by Log(2) and Log(10))
  • Zeta(3) - Amdeberhan-Zeilberger Formula 1
  • Zeta(3) - Amdeberhan-Zeilberger Formula 2
  • Catalan - Lupas

The implementation for each of these is about 3000 lines - including multi-threading and swap. For each algorithm listed here, these 3000 lines are mostly copy/paste duplicates of each other with the necessary changes made for each algorithm.

• Much of the code-bloat in the performance-critical code come from a number of very dirty low-level optimizations:
  • Manual loop-unrolling.
  • Aggressive inlining of recursion end-points.
  • Multiple code-paths optimized for different instruction sets.

     

• I no longer keep track of how much deactivated code there is. I've probably thrown away close to a quarter million lines of code relating to y-cruncher since 2007.
  

Formulas:

y-cruncher has two algorithms for each major constant that it can compute - one for computation, and one for verification.

All complexities shown assume multiplication to be O(n log(n)). It is slightly higher than that, but for all practical purposes, O(n log(n)) is close enough.

Square Root of n and Golden Ratio

1st order Newton's Method - Runtime Complexity: O(n log(n))
Note that the final radix conversion from binary to decimal has a complexity of O(n log(n)²).

e

Taylor Series of exp(1): O(n log(n)²)

Taylor Series of exp(-1): O(n log(n)²)

 

Pi

Chudnovsky Formula: O(n log(n)³)

Ramanujan's Formula: O(n log(n)³)

 

Log(2)

Machin-Like Formulas: O(n log(n)³)

Despite the fact that the formulas are so similar, they are sufficient to qualify as computation and verification.
As long as two machin-like formula for the same constant have different coefficients for all shared terms, they are sufficent for computation and verification.

 

Log(10)

Machin-Like Formulas: O(n log(n)³)

 

Zeta(3) - Apery's Constant

Amdeberhan-Zeilberger Formula 2: O(n log(n)³)

Amdeberhan-Zeilberger Formula 1: O(n log(n)³)

Catalan's Constant

Lupas Formula: O(n log(n)³)

Ramanujan's Formula: O(n log(n)³)

Euler-Mascheroni Constant

Brent-McMillan with Refinement: O(n log(n)³)



                          

Brent-McMillan (alone): O(n log(n)³)



                 

Note that both of these formulas are essentially the same.
Therefore, in order for two computations to be independent enough to qualify as a verified record, they MUST be done using different n.

For simplicity of implementation, y-cruncher only uses n that are powers of two - which serves a dual purpose of allowing the use of (the easily computed) Log(2), as well as lending itself to shift optimizations.

Arithmetic Algorithms:

Multiplication:

y-cruncher uses a combination of grade-school, Karatsuba, Fast Fourier Transform, and Hybridized Number-Theoretic Transform multiplication algorithms.

This rather unique combination of integer and floating-point based algorithms has proven itself to be extremely useful. Not only is it fast and memory efficient, but it is also well suited for parallelism and vectorization.

Even without the use of assembly language optimizations, the current multiplication scheme used in y-cruncher is faster than that of GMP (as of the version used in Mathematica 6).

(GMP is an extremely fast arbitrary precison arithmetic library that is used by virtually every major math program including Mathematica, and Matlab. It is fast due in part to its use of highly optimized and specialized assembly code for different processors.)

The table below summarizes the approximate ranges where each algorithm is used: (for all versions prior to v0.4.3)

(v0.4.3 has a completely different tuning and is a lot more complicated than this.)

Algorithm	Complexity	Approximate Range (decimal digits)
		x86	x64
Basecase (grade-school multiplication)		1 - 130	1 - 170
Karatsuba		130 - 1000	170 - 560
3-Way Toom-Cook		Implemented - But Never Optimal to Use
Floating-Point Fast Fourier Transform (FFT)	w = # bits in mantissa	1000 - 4,800,000	560 - 120,000
Hybridized Number-Theoretic Transform (Hybrid NTT)	Better than FFT, but worse than SSA Pending Further Analysis...	> 4,800,000*	> 120,000
Schönhage-Strassen Algorithm (SSA)		Implemented - But Never Optimal to Use**

The current method of having a single threshold between each algorithm is not optimal as it does not take into account for multiple crossover points between adjacent algorithms.

*From my experience, it is never optimal to use Hybrid NTT on x86. FFT is faster for anything under the 32-bit limit.
But in the interest of reducing memory usage for large computations, I've forced y-cruncher to use Hybrid NTT for anything above 4.8 million digits - even if it means significantly slowing down a computation.

**The difference between the run-time complexities of Hybrid NTT and SSA is on the order of double-logarithms. However, SSA has a much larger Big-O constant - so large that Hybrid NTT is a lot faster for any attainable size despite having a slightly inferior complexity.
(Of course, this could change if there are any new ground-breaking optimizations that apply to SSA but not to Hybrid NTT.)

Other:

Nothing special here. These are all standard algorithms and approaches.

• Division: Division is done using first-order Newton's Method.
• Square Roots/Nth roots: All are done using first-order Newton's Method.
• Natural Logarithm: Logs that don't have a hard-coded Machin-like formula are done using Arithmetic-Geometric Mean (AGM).
• Radix Conversion: Prior to v0.5.3, radix conversions to base 10 were done using the standard Divide and Conquer algorithm.
  Starting from v0.5.3, they are done using Scaled Remainder Trees. This is the same approach that was used in the December 2009 world record for Pi.
• Series Summing: Series summing is done using standard Binary Splitting techniques with precision control to keep memory usage bounded.

As far as optimizations go, Division and Square Roots have been given very little attention and are therefore inefficient. (The current implementations aren't even fully paralleled.) Natural Log via AGM suffers as a result. Division and square roots account for such a tiny portion of the total run-time for the majority of constants that it has not been worth the effort to improve it. Most of the optimizations (as of version 0.5.3) have been concentrated in multiplication and Binary Splitting because they clearly dominate total run-time.