#!/usr/bin/env python3 # -*- coding: utf-8 -*- # *************************************************************************** # * Copyright (C) 2011, Paul Lutus * # * * # * This program is free software; you can redistribute it and/or modify * # * it under the terms of the GNU General Public License as published by * # * the Free Software Foundation; either version 2 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 General Public License for more details. * # * * # * You should have received a copy of the GNU General Public License * # * along with this program; if not, write to the * # * Free Software Foundation, Inc., * # * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * # *************************************************************************** from math import * import re def f(y): return 2 * sqrt((2-y)*y) def cff(y): if(y == 0): return 0 s1 = sqrt(y) s2 = sqrt(-y+2) return pi + s2 * (y-1) * s1 - 2 * atan(s2/s1) def num_int(f,ni,a,b): if(a == b): return 0 nf = float(ni) r = b-a s = 0 di = 1/nf for i in range(ni+1): x = i*di*r+a s += f(x) return s * r / nf def gen_table(ni): nx = 10 html = True print(('Samples: %d' % ni)) for i in range(1,nx+1): x = i * 2 / float(nx) rcf = cff(x) rni = num_int(f,ni,0,x) out = '%.1f %.6f %.6f %+.6f%%' % (x,rcf,rni, 100 * ((rni-rcf)/rcf)) if(html): out = '' + re.sub('(\S+)\s*','\\1',out) + '' print(out) for n in range(2,5): gen_table(int(pow(10,n)))