40

Let me quote the wonderful book Numerical Recipes in C++ (but also applicable): We assume that you know enough never to evaluate a polynomial this way: p=c[0]+c[1]*x+c[2]*x*x+c[3]*x*x*x+c[4]*x*x*x*x; or (even worse!), p=c[0]+c[1]*x+c[2]*pow(x,2.0)+c[3]*pow(x,3.0)+c[4]*pow(x,4.0); Come the (computer) revolution, all persons found guilty of ...


36

In an interview it usually doesn't matter if you actually solve the problem. What is most important is the way you (try to) solve it. If they don't tell you it should be the possibly fastest solution ever you should not optimize it prematurely but instead show that you know how to write SOLID code like encapsulate propertly, make the code testable etc. show ...


30

Your code seems to be working. Here are a few things to make it better. Tuple unpacking x = point[0] y = point[1] can be written : x, y = point Remove what is not required You maintain a count variable but you already know the count : it is TIMES_TO_REPEAT. Variable name TIMES_TO_REPEAT does not seem to be such a good name. NB_POINTS would probably ...


27

You don't need to compare prev and new during each iteration. The difference between the new and the previous sum is simply the current term : $$\frac{8}{(2i+1)^2}$$ If you want this term to be smaller than error, you can solve: $$\mathrm{error} > \frac{8}{(2i+1)^2}\\ \iff (2i+1)^2 > \frac{8}{error}\\ \iff 2i+1 > \sqrt{\frac{8}{error}}\\ \iff i &...


26

I'm going to give an alternative approach to tchbot's answer There are places where SOLID prinicpals are required. I don't see this as one of them. Here, we have a simple "Are you capable of writing a loop" question - about the same level as the Fizz Buzz test. Really, the function is described in one line, and when applying mathematical formulas, I ...


23

As it turns out, a similar question was asked recently on Math.SE. Rather than reinventing-the-wheel, take advantage of built-in functionality in Python. Your norm_cdf(z) is just a numerical approximation for $$P(z) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{z} e^{-t^2/2}\ dt = \int_{-\infty}^{z} Z(t)\ dt = \frac{1}{2} \left( 1 + \mathrm{erf}\left( \frac{z}{\...


20

Flavio's answer addresses a part of something that matters greatly for this kind of problem: rounding errors will kill you. If you look at the expression: $$1-\frac13+\frac15-\frac17+...$$ This is the same as $$\frac{3-1}{1\times3} + \frac{7-5}{5\times 7} + ...$$ It is easy to see that the nth term in this sequence is $$\frac{2}{(4\times n-3)(4\times n-...


19

So, you have a polynomial sequence and want to sum its terms while they are greater than a provided tolerance (i.e.: the tolerance is lower than the current computed term). This is easily expressed using basic python constructs / builtins and the itertools module: An infinite polynomial sequence can be described using a generator expression and the ...


18

I'll only address the so-called "magic numbers" that several reviewers have mentioned. Sometimes, when you're working in pure mathematics, what seems at first glance to be a "magic number" really isn't. It may be that the numbers themselves are just part of the problem statement. I think the question boils down to this: can you come up with a name that is ...


17

If you are writing mathematical code in Python, then it is worth looking at NumPy, which is a library implementing (among other things) fast arithmetic operations on arrays of floating-point numbers. In NumPy you'd write: import numpy as np def calculate_area(f, a, b, n): """Return an approximation to the definite integral of f from a to b using ...


17

This answer uses pointer-casting for type-punning just to save space. In practice keep using your union (safe in ISO C99, and in C++ as a GNU and MSVC extension) or memcpy (safe in C and C++). This pointer-casting is only safe in MSVC, or GNU-compatible compilers with -fno-strict-aliasing Initial approximation Packed bit fields are not only unnecessary ...


16

Instead of math.pow, use the builtin ** operator. You don't need the \s at EOL, because the parentheses surrounding the expression allow it to implicitly span multiple lines. So after making both of those changes, I wind up with: def norm_cdf(z): """ Use the norm distribution functions as of Gale-Church (1993) srcfile. """ # Equation 26.2.17 from ...


15

Not sure if this helps but you could easily define a function to evaluate the value of a polynom at a given position def evaluate_polynom(pol, x): return sum(a * math.pow(x, i) for i, a in enumerate(pol)) Then (0.319381530 * t) + (-0.356563782* math.pow(t,2)) + (1.781477937 * math.pow(t,3)) + (-1.821255978* math.pow(t,4)) + (1.330274429 * math.pow(t,5)...


15

Please fix the indentation (probably just a copy-pasting issue). In Python, you don't need parenthesis around the expressions for block statements like if and while. You also don't need it around your expressions (thanks Ev. Kounis). You can use i ** 2 instead of i * i. You can use the increment operator: n += 2 instead of n = n + 2. You should probably use ...


14

Convergence testing if pi == piold: break This is not usually done, because float equality has a lot of gotchas. In this case it's possible due to the numbers being Decimal, but if you need to move away from Decimal you're going to encounter issues. Usually, convergence is measured as the absolute error decreasing below a chosen epsilon, a ...


13

First, a comment on the implementation: note that picking a random number between 1 and 10**5 is really no different than picking a random number between 0 and 1; you scale everything by the same amount, so why bother? In fact, requiring the x and y to be integers (and not including 0) reduces the accuracy of your method. Removing that factor allows the ...


13

First, this is an awesome video! Upvoted for that reason alone. :) If n=2, the program finishes within milliseconds, whereas for n=3, it takes a whopping 115 seconds. Do you know about big-O notation? Just off the top of my head after watching that video: for n=2 you're computing the number of collisions for a 1kg block and a 100**2 = 10,000kg block, ...


13

What you are approximating is $$ \pi =\sum_{m=1}^{\infty}\frac{(-1)^{s(m)}}{m},$$ where \$s(m)\$ counts the number of appearances of primes of the form \$4k+1\$ in the prime decomposition of \$m\$, compare How can we prove \$\pi =1+\frac{1}{2}+\frac{1}{3}+\frac{1}{4}-\frac{1}{5}+\frac{1}{6}+\frac{1}{7}+\cdots\,\$? on Mathematics Stack Exchange. The ...


12

As you've already realized, computational loops in pure python are often unacceptably slow. For this sort of code you're very often better off using numpy. For example: import numpy as np def pi_mc(N, seed=1234): rndm = np.random.RandomState(seed) x = rndm.uniform(size=N) y = rndm.uniform(size=N) r = x*x + y*y return 4. * np....


12

There is clearly no 'best way' to do this. From a 'Numerical Analysis' approach, you should never sum that way a series with alternating-sign terms: loss of accuracy and slow convergence speed can be expected. Use simple arithmetic to sum the terms in pairs, and then iterate by two units at a time. If you do this, you will gain a lot both in accuracy and in ...


11

The calculation would be simplified if you centered the circle at the origin. Also, by symmetry, you only need to model one quadrant. import random TIMES_TO_REPEAT = 10**5 LENGTH = 10**5 def in_circle(x, y): return x**2 + y**2 < LENGTH**2 inside_count = 0 for _ in range(TIMES_TO_REPEAT): point = random.randint(0,LENGTH), random.randint(0,...


11

I see a few things that may help you improve your program. Avoid pow for float The use of pow converts the argument to a double and returns a double. If you're casting the result to a float anyway, the use of std::powf would be more appropriate. But see the next suggestion for an even better approach. Prefer multiplication to pow with small, fixed ...


11

One could consider at least the following points: Instead of including <stdlib.h>, I'd include <cstdlib>. In getRandDart(), it might in this case be more readable to do static_cast<double>(rand()) / RAND_MAX; instead of multiplying by 1.0. In the for loop, all of x, y and d can be const, so I'd make them const. This has the potential to ...


10

The code and documentation look fine in general, so this review will focus on minor optimizations on a per-method basis. Use ThreadLocal.withInitial You can set your LOCAL_RANDOM with the following: private static final ThreadLocal<Random> LOCAL_RANDOM = ThreadLocal.withInitial(Random::new); You need to give a Supplier<Random> as argument, ...


9

There's not a lot of code here so there's not a lot to say. ^_^ Implementations like this are generally outside of the day-to-day knowledge of most programmers. You may wish to include a reference in the code to how the algorithm was derived and how it works. The chosen variable names are pretty useless. I have no idea what an or sm mean. Because I don't ...


9

FYI, in IEEE754 sqrt is a "basic" operation that's required to be correctly-rounded (rounding error <= 0.5ulp), same as + - * /. Hardware FPUs (I think) always provide sqrt if they provide the other operations, especially division which is typically implemented similarly (and with similar performance). NR iterations each involving a division are not ...


8

There doesn't seem to be a need for separate Male and Female structures as they hold the same exact data. Instead, just have one named Person and create male and female instances of it. You also don't need classes if you're not going to have any private data; just use structs. struct Person { bool isMarried; int salary; Person() : isMarried(...


8

You don't need line continuation characters inside brackets. After splitting up your reduce lines, line_profiler says most of your time is actually in (p2grid[0:n] * 3*Bgrid[0:n] / (p2 * (p2grid[0:n] + Bgrid[0:n]**2))) * dq[0:n] and (3*Bgrid[i:] / (p2grid[i:] + Bgrid[i:]**2)) * dq[i:] Extracting those outside of the loop and pulling a division out of ...


8

Comments OK I like the first comment. /******************************************************************************* * * Grant Williams * * Version 1.0.0 * August 27, 2015 * * Test Program for Brent's Method Function. * * Brent's method makes use of the bisection method, the secant method, and inverse quadratic interpolation in one algorithm. * * To ...


8

Observations and Suggestions Operate on pi directly. You can eliminate the pi = 4 * number;. You are using more operations than are necessary. I prefer decimal over double here for precision. Decimal is slower and at 6 digits of precision I don't think double would get you in trouble. More precision double could introduce errors which are not self ...


Only top voted, non community-wiki answers of a minimum length are eligible