I wrote a program that approximates Pi by using polygons. I used the formulars in the picture beyond. In my code they are called `innerPoly` (\$c_{2n}\$) and `outerPoly` (\$C_{2n}\$). And since you can calculate the circumference of a 2n-polygon with knowing the circumference of a n-polygon you will get the circumferences \$C_8\$, \$C_{16}\$, \$C_{32}\$ etc., knowing \$C_4\$. $$ \begin{array}{l} c_{2n} =& 2 \sqrt{2n^2-n\sqrt{(2n)^2-c_n^2}} \qquad&\textrm{for the inner polygon, with}\ c_4=4\sqrt{2} \\ C_{2n} =& \frac{4 n C_n}{2n + \sqrt{(2n)^2 + C_n^2}} &\textrm{for the outer polygon, with}\ C_4=8 \end{array} $$ My thoughts are: - Would it make it anyhow better when I have a function called `void PiApproximation()` that writes my values to stdout already? I mean technically you can always put your full code into the `main()`-function, but when you have to use the same code-parts over and over again you should make a own function of it and call it, when u need it. So I guess in this case it will make no difference if I use an own function or calculate the circumferences and print it in the main-function. - What else can I improve? pi_approx.c ``` #include #include #define INNER_FOUR 4*sqrt(2); //circumference c_4 of the inner tetragon(square) #define OUTER_FOUR 8 //circumference C_4 of the outer tetragon int power(int n, int p); int main(void) { int n = 4; double innerPoly = INNER_FOUR; double outerPoly = OUTER_FOUR; printf("PI-APPROXIMATION USING POLYGONS\n"); printf("===============================\n\n"); printf(" n I c_n/2 I C_n/2 I\n"); printf("---------I-----------------I-----------------I\n"); for (int i=3; n<=8192; n=power(i,2), i++) { printf(" %4d I %1.8lf I %1.8lf I\n", n, innerPoly / 2, outerPoly / 2); innerPoly = 2 * sqrt(2 * n*n - n*sqrt(4 * n*n - innerPoly*innerPoly)); //formular c_2n outerPoly = (4 * n * outerPoly) / (2 * n + sqrt(4 * n*n + outerPoly*outerPoly)); //

My C program uses the Monte Carlo method to approximate the mathematical constant π, the ratio of a circle's circumference to its diameter (and, importantly for this code, 4 times the ratio of a circle's area to that of its bounding square). This estimation executes more slowly when parallelised using OpenMP than it does when executed in a single thread. I think that my random generator function is not thread safe, and that may be the cause of the slowdown. ``` #include #include #include #include #define DART 1000000 // Number of darts each player throws // Number of random dots in the square[-1,+1] #define MAXPLAYER 8 // Maximume number of participant players // Number of threads /** Generates a random number between two params given @param: random number scope @return: desired random number */ long double fRand(long double fMin, long double fMax) { long double f = (double)rand() / RAND_MAX; return fMin + f * (fMax - fMin); } /** Generates random dots on the board @param: playerDarts, total number of darts to throw @return: score, number of darts thrown in the circle */ int player(int playersDarts) { srand(time(NULL)); long double pi, x, y; int score = 0; for (int i = 0; i < playersDarts; i++) { x = fRand(-1.0, 1.0); y = fRand(-1.0, 1.0); if (x*x + y*y < 1.0) score++; } return score; } void main() { long double pi; long const double REAL_PI = 3.141592653589; int score = 0, playersDarts; //////////////////////////////////////////////////// // Parallel // //////////////////////////////////////////////////// // devide the total number of DARTS between players playersDarts = DART / MAXPLAYER; double beginParallel = omp_get_wtime(); #pragma omp parallel for for (int i = 1; i <= MAXPLAYER; i++) score += player(playersDarts); double endParallel = omp_get_wtime(); pi = 4.0 * ((lo

Background For a future workshop I'll have to fit arbitrary functions (independent variable is height `z`) to data from multiple sources (output of different numerical weather prediction models) in a yet unknown format (but basically gridded height/value pairs). The functions only have to interpolate the data and be differentiable. There should explicitly be no theoretical background for the type of function, but they should be smooth. The goal is to use the gridded (meaning discrete) output of the numerical weather prediction model in our pollutant dispersion model, which requires continuous functions. Workflow - choose the input model - load input data - define list of variables (not necessarily always the same) - define height ranges (for the piecewise function) - define base functions like "a0 + a1*z" for each height range and variable - optionally define weights, because some parts are more important that others - fit the piecewise functions - save the fitted functions and their derivatives as Fortran 90 free form source code (to be included in our model) I don't think 1.-6. can be automated, but the rest should be. Code ``` from __future__ import absolute_import from __future__ import division from __future__ import print_function from __future__ import unicode_literals import numpy as np import pandas as pd from scipy.optimize import curve_fit from sympy import log, ln, Piecewise, lambdify, symbols, sympify, fcode, sqrt def config(name): """Configuration of the piecewise function fitting, dependent on input name Input: name... name of experiment to fit data to, basically chooses settings Output: var_list... list of variables to fit infunc_list_dict... dictionary with var_list as keys, each having a list as value that contains strings with the sub-function to fit, from the bottom up. Only the first (lowest) may have a constant value, all others must be 0 at the height they "take over" (where their

Greg Beech says here that he asks C# candidates to produce a formula that calculates PI Given that Pi can be estimated using the function \$4 * (1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \dots)\$. I'm far from C# interview ready, but I like the challenge. Here is my attempt at the formula. Let me know if there is a better way to do it. ``` using System; class MainApp { static void Main () { double number = 0; double pi; int i = 1; do { if ((i/2) % 2 == 0) { number += (double)(1) / i; } else { number -= (double)(1) / i; } pi = 4 * number; i += 2; } while (Math.Round(pi,5) != 3.14159); Console.Clear(); Console.WriteLine("Pi can be found using the formula 4 * (1 - 1/3 + 1/5 - 1/7 + ...) with {0} iterations.\n" + "At this point 4 is multiplied by: {1} to get {2}\n" + "This Rounds up to: {3}",i,number,pi,Math.Round(pi,5)); Console.ReadKey(); } } ```

I'm trying to lighten the code review load for the maintainers of boost.math, and I was hoping you guys could help me out. I have a pull request which implements tanh-sinh quadrature, which is provably optimal for integrands in Hardy spaces. Here's my code: - Documentation - Interface - Implementation - Tests which is also reproduced below. I have a few design worries. - It is a class and not a function. This is a bit confusing; I worry that people will not recognize that the constructor is doing some one-time calculations to make integrations faster. - It takes a long time to compile. I generated the abscissas and weights to 100 digits, and then they must be lexically cast to their correct type. I could keep fewer levels of abscissas and weights in the .hpp, but then the runtime would longer for complex integrands. - For integrands in Hardy spaces, the number of correct digits always doubles on each refinement. However, we want to do just the right amount of work to deliver the requested accuracy, which is almost always overshot. Interface: ``` #ifndef BOOST_MATH_QUADRATURE_TANH_SINH_HPP #define BOOST_MATH_QUADRATURE_TANH_SINH_HPP #include #include #include #include namespace boost{ namespace math{ template class tanh_sinh { public: tanh_sinh(Real tol = sqrt(std::numeric_limits::epsilon()), size_t max_refinements = 20); template Real integrate(F f, Real a, Real b, Real* error = nullptr); private: std::shared_ptr> m_imp; }; template tanh_sinh::tanh_sinh(Real tol, size_t max_refinements) : m_imp(std::make_shared>(tol, max_refinements)) { return; } template template Real tanh_sinh::integrate(F f, Real a, Real b, Real* error) { using std::isfinite; using boost::math::constants::half; using boost::math::detail::tanh_sinh_detail; if (isfinite(a) && isfinite(b)) { if (b a.\n"); } Real avg = (a+b)*half(); Real diff = (b-a)*half(); auto u = [=](Real z) { return f(avg + d

I have written a program that implements the ADI method and Crank-Nicolson method for solving Schrodinger equations. The program is working, but it takes a very long time to run. I am looking for tips on how to improve the performance of the program. ``` #define _USE_MATH_DEFINES #include #include #include #include #include #include #include #include #include #include #include #include #define PI M_PI #define h_bar 1.0 #define M 1.0 #define L .5 #define sig .05 #define grid_point 100 #define K1 1.0E6 #define K2 1.0E6 #define omega1 std::sqrt(K1/M) #define omega2 std::sqrt(K2/M) typedef std::complex compx; #define I compx(0.,1.) #define one compx(1.,0.) #define two compx(2.,0.) class wave_function { public: wave_function(); wave_function(bool); double dt = 4.2E-06; double dx = L / grid_point; std::vector> value; void solve_triag(); double potential(int, int); double real_space(int); double density_x1[grid_point]; double density_x2[grid_point]; compx sec_space_deriv(char, int, int); void normalize(); void rho(); void alpha_beta_solver(wave_function &, compx, std::vector &mid, compx, std::vector &R, int, char); compx a = -h_bar / (two*M), b = -h_bar / (two*M); double r = dt / (dx*dx); compx A = (I*r*a / two), B = (I*r*b / two), C = I*dt / two; }; wave_function::wave_function() { value.resize(grid_point); for (int l = 0; l mid1, mid2, R1, R2; wave_function tmp; mid1.resize(grid_point); R1.resize(grid_point); for (int x2 = 0; x2 &mid, compx post_mid, std::vector &R, int coor, char x_or_y) { compx x_N, R_N;// new_mid[grid_point], new_R[grid_point]; std::vector alpha(grid_point - 2); std::vector beta1(grid_point - 1);//for x_1 std::vector beta2(grid_point - 1);//for x_2 std::vector x_1(grid_point), x_2(grid_point); alpha[0] = post_mid / mid[0]; beta1[0] = R[0] / mid[0]; beta2[0] = -b4_mid / mid[0]; //Forward run for (int l = 1; l = 0; l--) {

I'm relatively new to Python, coming from a deep C++ background. I'm mostly looking for feedback on how to make my code more idiomatic/pythonic, but I would welcome and appreciate any and all other feedback as well. This code uses the Python library PyEphem, which builds upon XEphem, an X (UNIX GUI) program that's an astrological ephemeris, which is (was historically) a table of celestial bodies and their positions in the sky. PyEphem is used to track the altitudes of astronomical bodies, which are used as an indication of when astrological aspects have or are going to occur. I'm deploying this code as an AWS Lambda microservice, thus the `lambda_handler` function. Consider also reviewing the test and deployment Makefile that corresponds to the module below. Here is the entirety of my module: ``` #!/usr/bin/env python2.7 # aspectus: aspect+prospectus. iCalendars with astrological aspect events. # Copyright (C) 2017 Frederick Eugene Aumson # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as published # by the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . """aspectus: aspect+prospectus. Generates iCalendar .ics files to populate your calendar application with events for astrological aspects (angles) between two celestial bodies. Currently supports only the trine and sextant aspects, and only between the sun and the earth. Could easily be extended to support other bodies, and perhaps with not much more work extended to support aspects between two bodies outsi

This post is the first in the MathCore series. The next post is here: Arbitrary precision π (Circular Constant) in Java - MathCore #2 Disclaimer My project is too big to be reviewed in a single question. So, I'll post one class at a time. The relevant methods needed from other classes are added as snippets. `Roots.java` The purpose of this class is to calculate the principal n-th root of a non-negative `BigDecimal` with the specified precision. The precision used internally is actually a bit more so that the answer returned is fully accurate up to the specified precision and rounding can be done with confidence. ``` package mathcore.ops; import java.math.BigDecimal; import java.math.BigInteger; import java.math.MathContext; import static java.math.BigDecimal.ONE; import static java.math.BigDecimal.ZERO; import static java.math.BigDecimal.valueOf; /** * A portion of BigMath refactored out to reduce overall complexity. * * This class handles the calculation of n-th principal roots. * * @author Subhomoy Haldar * @version 1.0 */ class Roots { // Make this class un-instantiable private Roots() { } /** * Uses the n-th root algorithm to find principal root of a verified value. * * @param a The value whose n-th root is sought. * @param n The root to find. * @param c0 The initial (unexpanded) MathContext. * @return The required principal root. */ private static BigDecimal nthRoot(final BigDecimal a, final int n, final MathContext c0) { // Obtain constants that will be used in every iteration final BigDecimal N = valueOf(n); final int n_1 = n - 1; // Increase precision by "n"; final int newPrecision = c0.getPrecision() + n; final MathContext c = BigMath.expandContext(c0, newPrecision); // The iteration limit (quadratic convergence) final int limit = n * n * (3

Recently, I was given a math assignment to calculate Gini Indexes for a table of percent distributions of aggregate income. The table takes the form of: `Year 1st Quintile 2nd Quintile 3rd Quintile 4th Quintile 5th Quintile ---- ------------ ------------ ------------ ------------ ------------ 1929 0.03 12.47 13.8 19.3 54.4 1935 4.1 9.2 14.1 20.9 51.7 ... ... ... ... ... ... ` Because of the length of the actual table I wrote a short python script to calculate the Gini Indexes. However, I'm fairly new to Python so I'd like to see what needs improvement. Method of Calculation: Calculate the log of the percentages with the quintile as the base of the logarithm (`i.e. log(0.0003)/log(0.2), log(0.1247)/log(0.4) ...`) and then average these values to find an approximate exponent for the Lorenz curve. Calculate the Gini Index by finding twice the area between `y=x` and the Lorenz curve from `0` to `1`. `import numpy as np import scipy.integrate as integrate import itertools as it def log_slope(x, y): return np.log(y)/np.log(x) # read in data from file years, data = read_values('Quintiles') # shape: [[0.03, 0.12,...], [], []..., []] accumulated_vals = [list(it.accumulate(v)) for v in data] # percentiles; remove 0.0 and 1.0 after calculating percentiles = np.linspace(0.0, 1.0, np.shape(accumulated_vals)[1] + 1)[1:-1] for j, vals in enumerate(accumulated_vals): sum = 0 for i, val in enumerate(vals[:-1]): # exclude the last accumulated value, which should be 1.0 sum += log_slope(percentiles[i], val) average = sum / (len(vals)-1) gini = 2 * integrate.quad(lambda x: x - pow(x, average), 0.0, 1.0)[0] print('{:d}: {} -> {:.5f}'.format(int(years[j]), [round(k, 4) for k in vals], gini)) ` Questions: - The code produces correct values and is alr

I'm kind of new to Java but I am writing an Android app. Right now I'm working on an async task to calculate Pi but when I run it the memory usage increases alarmingly (+5MB per second). This one piece of code may need to run thousands of times so it is important that it be very optimized. ``` protected BigDecimal doInBackground(Object... params) { digits = (int)params[0]; piView = (TextView)params[1]; boolean add = true; BigDecimal pi = new BigDecimal(0.0); for(int i=1; pi.setScale(digits, BigDecimal.ROUND_DOWN).equals(realPi) == false; i++){ if(add){ pi = pi.add(new BigDecimal(4.0 / (-1 + (i*2)))); add = false; } else { pi = pi.subtract(new BigDecimal(4.0 / (-1 + (i*2)))); add = true; } } return pi; } ``` This one piece of code may need to run thousands of times so it is important that it be very optimized.