7 ms·

Outperforming LAPACK with C++ metaprogramming

157 points by okaleniuk 9y ago

vog 9y ago
I wonder if the similar tricks work in other languages with similarily optimizing compilers, but "nicer" metaprogramming, such as Common LISP, Rust and perhaps OCaml? (with camlp4) Even more interesting would be whether this leads to the same huge performance gains in those languages as well.
- fithisux 9y ago
  Or even Dlang .
- ufo 9y ago
  These kinds of tricks certainly are more pleasant to do in languages with real metaprogramming support. For numerical applications, one interesting project that comes to mind is the Terra programming language. It uses a flexible high level language (Lua) to metaprogram a fast and predictable C-like language (Terra). The advantage compared to the examples you mentioned is that each language can focus on its strength: the low level language being similar to C allows the programmer to better reason about performance and having a high level language for metaprogramming keep things simple. https://terralang.org https://terralang.org There are also systems out there specifically designed for numeric applications that let you specify an algorithm in a high level DSL and generate highly performing C code, some of which might even tesr various versions of the code to find one that runs better on your machine (for example, the best way to iterate the data might be affected by cache sizes). I don't remember the names of these off the top of my head but the Terra PhD thesis has lots of references to them.
- ska 9y ago
  It was playing around with exactly these sorts of tricks in c++ (with late 90's compilers) that drove me to learn common lisp. It's a much nicer environment to do this sort of thing in, but difficult to compete head to head with a good optimizing c++ compiler on code generation in general. Fwiw, in common lisp particularly having access to to compiler macros allows for interesting approaches. I was doing primarily research code at the time though, and the overall improvement in productivity more than made up for the fractional difference in run times for me. Of course, ymmv.
- ChrisRackauckas 9y ago
  It does lead to performance gains. This kind of stuff is done all the time in Julia's numerical codes. Flexible macros definitely make this easier, and multiple dispatch + automatic function specialization on types + inlining + interprocedural optimizations makes a lot of what's discussed here automatic. If it's not automatic, then libraries like StaticArrays.jl make functions that use small arrays do these kinds of tricks automatically (it overloads and unrolls operations specifically to match the size of the array which is compile time information). I am sure other languages can do this too, but having it semi-automatic is pretty cool if what you do all day is write numerical codes.
  - uadjet 9y ago
    Julia's flexibility in the numerical space is why it's my daily driver language for prototyping and implementation nowadays. Truly a fantastic language.
- sidlls 9y ago
  Just be wary of performance gains claims. Per other comments this implementation is highly optimized for a specific use and uses a numerically unstable algorithm.
barrkel 9y ago
Please check with array initialized from file or other dynamic source. With maximal inlining, there's a risk of constant propagation resulting in a lot of computation happening at compile time rather than run time.
- dpwm 9y ago
  As an anecdote of just how powerful this can be, I once wrote a prototype of a simulation (written in C) that used a lot of integer divisions and needed to be run many times with different parameters. As a quick and dirty hack that took about half hour to implement, I decided to just have the program as a template string and the parameters be substituted in by a python script that called GCC. To the amazement of a very sceptical friend, the result was a 5x speedup. A lot of this speedup was consistent with replacing the integer divisions with multiplicative inverses, but also merging constants in a way that would have made things quite unreadable. I was quite impressed with how far GCC had gone. Considering each instance of this simulation ran for tens of minutes, it was well worth compilation time of a few seconds. It kept the code really simple, too.
  - contravariant 9y ago
    What made you choose a python script rather than C++ templates or the C preprocessor?
    - mmlb 9y ago
      The way I read GP, the python script can very well do the equivalent of gcc -DPARAM=42 -o binary-param=42 && ./binary-param-42. edit: oh nm, reread GP...
      - dpwm 9y ago
        Pretty much, except with more parameters taken from a file. The substituting in of parameters was done in Python just so that I could manually check the output for a few cases. I think if I had to do something similar again today, any script would output a shell script that used command line macro definitions. That didn't really occur to me at the time.
sdeer 9y ago
It is somewhat unclear from his writing but did he compare this metaprogramming solution against OpenBlas or Atlas or just the reference BLAS/LAPACK implementation? If I was trying to do any significant linear algebra computations, I would definitely first try an optimized BLAS/LAPACK solution first.
- pletnes 9y ago
  Also, it seems the author solved a triangular system. LAPACK has special routines for that. Were they used? LAPACK is typically optimized for larger systems, not 5 unknowns. Also, 5 is not a great number for vectorized operations - it might even be beneficial to zero/one pad the matrix. Optimized LAPACK is often 5-10x faster than «basic LAPACK». BLAS/LAPACK was originally written in the 70s/80s and sometimes unroll loops to the tune of 5 or 7 - it made sense then. Not so much these days. I didn’t read the article in detail, but there appear to be a lot of holes in it.
  - sevensor 9y ago
    Personally, my biggest reason to prefer LAPACK in general is that its authors have already put a great deal of effort into correctness and numerical stability, so I don't have to. Even basic LAPACK is pretty fast, let alone the optimized libraries. Hand-optimizing my own special case is an absolute last resort.
    - pletnes 9y ago
      Yes, Cramer’s rule (referred to in the article) does not give numerically stable results AFAIK - which is yet another can of worms. Good catch.
      - thanatropism 9y ago
        I don't think I have heard Cramer's rule mentioned without the caveat "but don't use this for real problems".
virgilp 9y ago
[edit] I had a bug in my code; it does seem that the compiler is not picking up constants at compiletime. That is to say, it's not precomputing the final value at compiletime based on the value from those matrices. Here's how I modified the code to make sure: https://pastebin.com/iG21pTiP https://pastebin.com/iG21pTiP
- gmueckl 9y ago
  There needs to be a j in the exit condition of the inner loop of the clobber function.
drwells 9y ago
The proposed C++ tricks are useful and I think that its reasonably well-known (at least to people who do numerics) that a C++ compiler, with the right flags, can beautifully unroll loops when array sizes are known at compile-time. > brute force Cramer's rule solution This isn't really a fair comparison: trusty old DGESV will yield far more accurate results for badly conditioned problems since it does row swaps (partial pivoting). It might also rescale values; I don't remember the details. It would be much more reasonable to compare a hypothetical DGEC (double precision general matrix cramer's rule) routine.
jordigh 9y ago
D's Mir does the same sort of thing, and metaprogramming in D is much nicer, because it was designed from the ground up instead of the bug elevated to feature that metaprogramming is in C++. http://docs.mir.dlang.io/latest/index.html http://docs.mir.dlang.io/latest/index.html http://blog.mir.dlang.io/glas/benchmark/openblas/2016/09/23/glas-gemm-benchmark.html http://blog.mir.dlang.io/glas/benchmark/openblas/2016/09/23/...
- acdha 9y ago
  > instead of the bug elevated to feature that metaprogramming is in C++. That inflammatory wording doesn’t add anything to your comment. If you wanted to compare D with C++ it’d be a lot more useful to link to something which does that rather than just the documentation index.
  - jordigh 9y ago
    C++ metaprogramming was an accident. Expression templates are a huge, obscure, and very surprising hack. I don't know about you, but I find them really hard to understand. SFINAE wasn't meant to be used for compile-time branching, but when people realised that it could be used that way, it was just about the only tool for metaprogramming in C++. It certainly wasn't designed for that purpose, though. C++ has for a long time embraced metaprogramming and given new tools for it (I think there's `static if` now? Or will be soon?), but it really all started out as an accident: https://softwareengineering.stackexchange.com/questions/125314/in-c-were-sfinae-and-metaprogramming-intentional-or-just-a-byproduct-of-templ https://softwareengineering.stackexchange.com/questions/1253... > If you wanted to compare D with C++ My second link does that, compares Mir to a bunch of other numeric libraries, most in C++. Edit: (Replying to comment below by acdha because HN does not allow me to post new comments right now.) I don't have that kind of example handy. I can point you to the D Gems section of the D tour to give you a taste of what is possible in D. Click on "Gems" in the navigation menu above to see them. https://tour.dlang.org/tour/en/gems/uniform-function-call-syntax-ufcs https://tour.dlang.org/tour/en/gems/uniform-function-call-sy... If you want a detailed blog post comparing D and C++ metaprogramming, I don't have it ready yet. If you're curious to learn about it, look in more depth at the links I've provided. I'm sorry I don't have the resources to give you the detailed answer you require.
costrouc 9y ago
For those wondering the lapack guys are working on a new version of lapack to run on heterogeneous architecture including gpus. See their work at http://www.icl.utk.edu/research http://www.icl.utk.edu/research I have worked with these guys and all I can say is good luck outperforming the highly optimized routines they have written. My bet is that the the guy writing this blog used a non optimized version of lapack.
- stephencanon 9y ago
  There's no mystery here, author is only solving a 5x5 system. LAPACK will have barely finished checking function call parameters by the time a fully specialized and inlined solver is finished. Also, from a cursory inspection, author does no pivoting [edit: worse, it just uses Cramer's rule, which is insane], so their solver will be much less numerically robust (but this is admittedly less of a problem for the 5x5 case). The performance is completely expected--in fact, I would expect an optimized 5x5 solver to be faster than this.
  - stabbles 9y ago
    Cramer's rule is unstable even for 2x2 problems
    - mattnewton 9y ago
      I think gp is missing a “less” before “numerically robust”
      - stephencanon 9y ago
        Yeah, thanks. Fixed.
  - electricslpnsld 9y ago
    A better comparison might be a C++ library like Eigen that is intended for use on these small systems, that should optimize out function calls in a similar way, and that uses intrinsics.
bluescarni 9y ago
In his first factorial example, he could've avoided templates altogether with constexpr: https://godbolt.org/g/ZaoMNa https://godbolt.org/g/ZaoMNa Or just use a template specialisation: https://godbolt.org/g/TVkHMd https://godbolt.org/g/TVkHMd
- ZenoArrow 9y ago
  Thanks for sharing the godbolt site, can see myself using it for debugging.
  - andrewjw 9y ago
    the author of the site (Matt Godbolt) is in general pretty great.
1024core 9y ago
Template metaprogramming just shifts the burden of computation to compile time. Sure, I could pre-compute the first 10000 factorials, put them in an array, compile and then print the answer in constant time. But it'll work for the first 10000 numbers only.
- ska 9y ago
  It's not nearly that simple. You can use template meta-programming to make expressions both easy to read and efficient, by avoiding creation of temporary objects, fusion, lifting, etc. etc. In some domains (e.g. linear algebra) this can be very powerful because you don't have write non idiomatic expressions just to feed the compiler things it can understand.
bwasti 9y ago
in C++17 you can get generic loop unrolling templates in 8 lines of code: https://godbolt.org/g/2Di4PM https://godbolt.org/g/2Di4PM
xvilka 9y ago
Modern Fortran (especially upcoming Fortran 2018 standard) is good enough already for those kind of libraries. Moreover, big chunks of the computations are implemented in assembly in BLAS/LAPACK. What surprises me, that authors of those libraries still stuck in Fortran 90.
stabbles 9y ago
This article makes no sense at all. It compares a numerically unstable O(n!) algorithm (Cramer's rule) to a numerically stable O(n^3) algorithm (Gaussian elimination).
- sannee 9y ago
  What bothers me more is the fact that the article seems to compare LAPACKs function for solving any-sized systems of equations with a specialized function for solving 5x5 systems. I mean, it's neat trick and all, but the comparison with LAPACK does not really stand on equal ground.
blt 9y ago
Others have already mentioned Cramer's rule, but even if this code were switched to Gaussian elimination it might not beat lapack on big problems. Lapack implementations do extra tricks to optimize cache usage on big matrices. Would be interesting to include eigen in the benchmark too.
- bachmeier 9y ago
  I've read your comment several times and can't follow it. I think you're saying that even if the author continued to use Cramer's rule, he would lose to LAPACK on big problems.
  - AnimalMuppet 9y ago
    No. If he continued to use Cramer's rule, he would definitely lose to LAPACK on big problems, purely because of using Cramer's rule. But even if he didn't continue to use Cramer's rule (if he switched to Gaussian elimination), he still might lose to LAPACK, because LAPACK is pretty well optimized. At least, that's how I interpreted the comment...
    - blt 9y ago
      Correct interpretation
    - sannee 9y ago
      > he would definitely lose to LAPACK on big problems, purely because of using Cramer's rule. More importantly though, the code would probably take ages (as in, "age of the universe" ages) to compile, as Cramer's rule is O(n!) (using a naive implementation, I think one could get to O(n^4) if determinant was computed using LU factorization, but what's the point of Cramer's rule then...).
skierscott 9y ago
There’s also a closed form for the Fibonacci numbers. It runs in O(1). http://stsievert.com/blog/2015/01/31/the-mysterious-eigenvalue/ http://stsievert.com/blog/2015/01/31/the-mysterious-eigenval...

157 points by okaleniuk 9y ago

vog 9y ago

I wonder if the similar tricks work in other languages with similarily optimizing compilers, but "nicer" metaprogramming, such as Common LISP, Rust and perhaps OCaml? (with camlp4) Even more interesting would be whether this leads to the same huge performance gains in those languages as well.

fithisux 9y ago
Or even Dlang .
ufo 9y ago
These kinds of tricks certainly are more pleasant to do in languages with real metaprogramming support. For numerical applications, one interesting project that comes to mind is the Terra programming language. It uses a flexible high level language (Lua) to metaprogram a fast and predictable C-like language (Terra). The advantage compared to the examples you mentioned is that each language can focus on its strength: the low level language being similar to C allows the programmer to better reason about performance and having a high level language for metaprogramming keep things simple. https://terralang.org https://terralang.org There are also systems out there specifically designed for numeric applications that let you specify an algorithm in a high level DSL and generate highly performing C code, some of which might even tesr various versions of the code to find one that runs better on your machine (for example, the best way to iterate the data might be affected by cache sizes). I don't remember the names of these off the top of my head but the Terra PhD thesis has lots of references to them.
ska 9y ago
It was playing around with exactly these sorts of tricks in c++ (with late 90's compilers) that drove me to learn common lisp. It's a much nicer environment to do this sort of thing in, but difficult to compete head to head with a good optimizing c++ compiler on code generation in general. Fwiw, in common lisp particularly having access to to compiler macros allows for interesting approaches. I was doing primarily research code at the time though, and the overall improvement in productivity more than made up for the fractional difference in run times for me. Of course, ymmv.
ChrisRackauckas 9y ago
It does lead to performance gains. This kind of stuff is done all the time in Julia's numerical codes. Flexible macros definitely make this easier, and multiple dispatch + automatic function specialization on types + inlining + interprocedural optimizations makes a lot of what's discussed here automatic. If it's not automatic, then libraries like StaticArrays.jl make functions that use small arrays do these kinds of tricks automatically (it overloads and unrolls operations specifically to match the size of the array which is compile time information). I am sure other languages can do this too, but having it semi-automatic is pretty cool if what you do all day is write numerical codes.
- uadjet 9y ago
  Julia's flexibility in the numerical space is why it's my daily driver language for prototyping and implementation nowadays. Truly a fantastic language.
sidlls 9y ago
Just be wary of performance gains claims. Per other comments this implementation is highly optimized for a specific use and uses a numerically unstable algorithm.

barrkel 9y ago

Please check with array initialized from file or other dynamic source. With maximal inlining, there's a risk of constant propagation resulting in a lot of computation happening at compile time rather than run time.

dpwm 9y ago
As an anecdote of just how powerful this can be, I once wrote a prototype of a simulation (written in C) that used a lot of integer divisions and needed to be run many times with different parameters. As a quick and dirty hack that took about half hour to implement, I decided to just have the program as a template string and the parameters be substituted in by a python script that called GCC. To the amazement of a very sceptical friend, the result was a 5x speedup. A lot of this speedup was consistent with replacing the integer divisions with multiplicative inverses, but also merging constants in a way that would have made things quite unreadable. I was quite impressed with how far GCC had gone. Considering each instance of this simulation ran for tens of minutes, it was well worth compilation time of a few seconds. It kept the code really simple, too.
- contravariant 9y ago
  What made you choose a python script rather than C++ templates or the C preprocessor?
  - mmlb 9y ago
    The way I read GP, the python script can very well do the equivalent of gcc -DPARAM=42 -o binary-param=42 && ./binary-param-42. edit: oh nm, reread GP...
    - dpwm 9y ago
      Pretty much, except with more parameters taken from a file. The substituting in of parameters was done in Python just so that I could manually check the output for a few cases. I think if I had to do something similar again today, any script would output a shell script that used command line macro definitions. That didn't really occur to me at the time.

sdeer 9y ago

It is somewhat unclear from his writing but did he compare this metaprogramming solution against OpenBlas or Atlas or just the reference BLAS/LAPACK implementation? If I was trying to do any significant linear algebra computations, I would definitely first try an optimized BLAS/LAPACK solution first.

pletnes 9y ago
Also, it seems the author solved a triangular system. LAPACK has special routines for that. Were they used? LAPACK is typically optimized for larger systems, not 5 unknowns. Also, 5 is not a great number for vectorized operations - it might even be beneficial to zero/one pad the matrix. Optimized LAPACK is often 5-10x faster than «basic LAPACK». BLAS/LAPACK was originally written in the 70s/80s and sometimes unroll loops to the tune of 5 or 7 - it made sense then. Not so much these days. I didn’t read the article in detail, but there appear to be a lot of holes in it.
- sevensor 9y ago
  Personally, my biggest reason to prefer LAPACK in general is that its authors have already put a great deal of effort into correctness and numerical stability, so I don't have to. Even basic LAPACK is pretty fast, let alone the optimized libraries. Hand-optimizing my own special case is an absolute last resort.
  - pletnes 9y ago
    Yes, Cramer’s rule (referred to in the article) does not give numerically stable results AFAIK - which is yet another can of worms. Good catch.
    - thanatropism 9y ago
      I don't think I have heard Cramer's rule mentioned without the caveat "but don't use this for real problems".

virgilp 9y ago

[edit] I had a bug in my code; it does seem that the compiler is not picking up constants at compiletime. That is to say, it's not precomputing the final value at compiletime based on the value from those matrices. Here's how I modified the code to make sure: https://pastebin.com/iG21pTiP https://pastebin.com/iG21pTiP

gmueckl 9y ago
There needs to be a j in the exit condition of the inner loop of the clobber function.

drwells 9y ago

The proposed C++ tricks are useful and I think that its reasonably well-known (at least to people who do numerics) that a C++ compiler, with the right flags, can beautifully unroll loops when array sizes are known at compile-time. > brute force Cramer's rule solution This isn't really a fair comparison: trusty old DGESV will yield far more accurate results for badly conditioned problems since it does row swaps (partial pivoting). It might also rescale values; I don't remember the details. It would be much more reasonable to compare a hypothetical DGEC (double precision general matrix cramer's rule) routine.

jordigh 9y ago

D's Mir does the same sort of thing, and metaprogramming in D is much nicer, because it was designed from the ground up instead of the bug elevated to feature that metaprogramming is in C++. http://docs.mir.dlang.io/latest/index.html http://docs.mir.dlang.io/latest/index.html http://blog.mir.dlang.io/glas/benchmark/openblas/2016/09/23/glas-gemm-benchmark.html http://blog.mir.dlang.io/glas/benchmark/openblas/2016/09/23/...

acdha 9y ago
> instead of the bug elevated to feature that metaprogramming is in C++. That inflammatory wording doesn’t add anything to your comment. If you wanted to compare D with C++ it’d be a lot more useful to link to something which does that rather than just the documentation index.
- jordigh 9y ago
  C++ metaprogramming was an accident. Expression templates are a huge, obscure, and very surprising hack. I don't know about you, but I find them really hard to understand. SFINAE wasn't meant to be used for compile-time branching, but when people realised that it could be used that way, it was just about the only tool for metaprogramming in C++. It certainly wasn't designed for that purpose, though. C++ has for a long time embraced metaprogramming and given new tools for it (I think there's `static if` now? Or will be soon?), but it really all started out as an accident: https://softwareengineering.stackexchange.com/questions/125314/in-c-were-sfinae-and-metaprogramming-intentional-or-just-a-byproduct-of-templ https://softwareengineering.stackexchange.com/questions/1253... > If you wanted to compare D with C++ My second link does that, compares Mir to a bunch of other numeric libraries, most in C++. Edit: (Replying to comment below by acdha because HN does not allow me to post new comments right now.) I don't have that kind of example handy. I can point you to the D Gems section of the D tour to give you a taste of what is possible in D. Click on "Gems" in the navigation menu above to see them. https://tour.dlang.org/tour/en/gems/uniform-function-call-syntax-ufcs https://tour.dlang.org/tour/en/gems/uniform-function-call-sy... If you want a detailed blog post comparing D and C++ metaprogramming, I don't have it ready yet. If you're curious to learn about it, look in more depth at the links I've provided. I'm sorry I don't have the resources to give you the detailed answer you require.

costrouc 9y ago

For those wondering the lapack guys are working on a new version of lapack to run on heterogeneous architecture including gpus. See their work at http://www.icl.utk.edu/research http://www.icl.utk.edu/research I have worked with these guys and all I can say is good luck outperforming the highly optimized routines they have written. My bet is that the the guy writing this blog used a non optimized version of lapack.

stephencanon 9y ago
There's no mystery here, author is only solving a 5x5 system. LAPACK will have barely finished checking function call parameters by the time a fully specialized and inlined solver is finished. Also, from a cursory inspection, author does no pivoting [edit: worse, it just uses Cramer's rule, which is insane], so their solver will be much less numerically robust (but this is admittedly less of a problem for the 5x5 case). The performance is completely expected--in fact, I would expect an optimized 5x5 solver to be faster than this.
- stabbles 9y ago
  Cramer's rule is unstable even for 2x2 problems
  - mattnewton 9y ago
    I think gp is missing a “less” before “numerically robust”
    - stephencanon 9y ago
      Yeah, thanks. Fixed.
- electricslpnsld 9y ago
  A better comparison might be a C++ library like Eigen that is intended for use on these small systems, that should optimize out function calls in a similar way, and that uses intrinsics.

bluescarni 9y ago

In his first factorial example, he could've avoided templates altogether with constexpr: https://godbolt.org/g/ZaoMNa https://godbolt.org/g/ZaoMNa Or just use a template specialisation: https://godbolt.org/g/TVkHMd https://godbolt.org/g/TVkHMd

ZenoArrow 9y ago
Thanks for sharing the godbolt site, can see myself using it for debugging.
- andrewjw 9y ago
  the author of the site (Matt Godbolt) is in general pretty great.

1024core 9y ago

Template metaprogramming just shifts the burden of computation to compile time. Sure, I could pre-compute the first 10000 factorials, put them in an array, compile and then print the answer in constant time. But it'll work for the first 10000 numbers only.

ska 9y ago
It's not nearly that simple. You can use template meta-programming to make expressions both easy to read and efficient, by avoiding creation of temporary objects, fusion, lifting, etc. etc. In some domains (e.g. linear algebra) this can be very powerful because you don't have write non idiomatic expressions just to feed the compiler things it can understand.

bwasti 9y ago

in C++17 you can get generic loop unrolling templates in 8 lines of code: https://godbolt.org/g/2Di4PM https://godbolt.org/g/2Di4PM

xvilka 9y ago

Modern Fortran (especially upcoming Fortran 2018 standard) is good enough already for those kind of libraries. Moreover, big chunks of the computations are implemented in assembly in BLAS/LAPACK. What surprises me, that authors of those libraries still stuck in Fortran 90.

stabbles 9y ago

This article makes no sense at all. It compares a numerically unstable O(n!) algorithm (Cramer's rule) to a numerically stable O(n^3) algorithm (Gaussian elimination).

sannee 9y ago
What bothers me more is the fact that the article seems to compare LAPACKs function for solving any-sized systems of equations with a specialized function for solving 5x5 systems. I mean, it's neat trick and all, but the comparison with LAPACK does not really stand on equal ground.

blt 9y ago

Others have already mentioned Cramer's rule, but even if this code were switched to Gaussian elimination it might not beat lapack on big problems. Lapack implementations do extra tricks to optimize cache usage on big matrices. Would be interesting to include eigen in the benchmark too.

bachmeier 9y ago
I've read your comment several times and can't follow it. I think you're saying that even if the author continued to use Cramer's rule, he would lose to LAPACK on big problems.
- AnimalMuppet 9y ago
  No. If he continued to use Cramer's rule, he would definitely lose to LAPACK on big problems, purely because of using Cramer's rule. But even if he didn't continue to use Cramer's rule (if he switched to Gaussian elimination), he still might lose to LAPACK, because LAPACK is pretty well optimized. At least, that's how I interpreted the comment...
  - blt 9y ago
    Correct interpretation
  - sannee 9y ago
    > he would definitely lose to LAPACK on big problems, purely because of using Cramer's rule. More importantly though, the code would probably take ages (as in, "age of the universe" ages) to compile, as Cramer's rule is O(n!) (using a naive implementation, I think one could get to O(n^4) if determinant was computed using LU factorization, but what's the point of Cramer's rule then...).

skierscott 9y ago

There’s also a closed form for the Fibonacci numbers. It runs in O(1). http://stsievert.com/blog/2015/01/31/the-mysterious-eigenvalue/ http://stsievert.com/blog/2015/01/31/the-mysterious-eigenval...