What is the value of pi?
Peter C. Chapin
I'm writing a couple of short toy programs to help me learn Fortran 95.
In my latest program I need to make use of the constant pi. I could of
course define it myself (with the 'parameter' attribute?). However, I'm
wondering if it is defined in a library module to the maximum precision
recognized by the compiler. For example, the Ada standard library does
this, and I think also some other language libraries as well. I looked
around a little to see if I could find such a library definition, but I
didn't see one. Am I missing it?
Peter
George
implicit none
INTEGER,PARAMETER::dp=SELECTED_REAL_KIND(15)
real(kind=dp) :: pi
pi = 4.0_dp * atan(1.0_dp)
a1 = 4.0_dp * pi * r1**2
a1 is the surface area. Cheers,
--
George
I believe that God has planted in every heart the desire to live in
freedom.
George W. Bush
Picture of the Day http://apod.nasa.gov/apod/
Richard Maine
> On Sun, 04 Jan 2009 22:24:06 -0500, Peter C. Chapin wrote:
> > In my latest program I need to make use of the constant pi. I could of
> > course define it myself (with the 'parameter' attribute?). However, I'm
> > wondering if it is defined in a library module to the maximum precision
> > recognized by the compiler....
> INTEGER,PARAMETER::dp=SELECTED_REAL_KIND(15)
> real(kind=dp) :: pi
> pi = 4.0_dp * atan(1.0_dp)
which is not required to define pi to the maximum precision recognized
by the compiler as requested by the OP (or even the maximum precision of
kind dp). The standard does not require the intrinsics to give last-bit
accuracy. In fact, it has no accuracy requirements at all.
In practice, this has reasonable odds of being at least close to
last-bit accuracy, but that is not specified by the standard.
For my own behalf, I tend to just spell out the decimal digits of pi to
well more than the relevant number of digits. It isn't as though you
have to do this more than one time *EVER* (not even one time per
program), so it isn't exactly a burden. It took me one or two orders of
magnitude more time to write this posting than it would have taken me to
write the relavant literal. (Though it helps that I thought it amusing
to memorize 50 digits of pi in about 5th grade and the darn thing is
still there wasting neural connections that could probably be put to
better use). I have more confidence in the robust precision of compiler
decimal-to-binary conversions than in triginometric functions.
--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
Terence
Peter C. Chapin wrote:
> Hello!
>
> I'm writing a couple of short toy programs to help me learn Fortran 95.
> In my latest program I need to make use of the constant pi.
As Richard suggests, just type it:-
3.1415926535897932384626433832795029d0
(approx).
e p chandler
Richard claims to have memorized pi to 50 digits! I keep the first 20
decimal places on tap, but I have to hum "Yankee Doodle" to retrive
them. Note that there is no direct relationship between the notes and
the digits.
-- e
Peter C. Chapin
> Richard claims to have memorized pi to 50 digits! I keep the first 20
> decimal places on tap, but I have to hum "Yankee Doodle" to retrive
> them. Note that there is no direct relationship between the notes and
> the digits.
Human memory is a wonderous thing.
I find it interesting that people in this group have felt the need to
memorize pi to high precision. I guess that says something about the
sort of folks we are!
As for myself, when I was in high school a friend of mine and I
memorized pi to 100 places and used to recite it together in unison like
a chant. (Did I mention we were geeks?) At this point I can only
remember about 15 places before I start getting confused. Human memory
may be wonderous but it can also be rather unreliable.
Anyway, thanks for the replies. I'll set up my constant myself. Perhaps
a module full of such constants is in order. Hmm...
Peter
Arjen Markus
>
> Anyway, thanks for the replies. I'll set up my constant myself. Perhaps
> a module full of such constants is in order. Hmm...
>
Such a module will indeed be useful, as you can add the constants you
need
in the precision you need. A language like C may offer such macros as
M_PI,
but as soon as you divert from the set that is (sort of [*])
predefined,
you have to add the new ones yourself.
Do take care of the kind:
double precision pi = 3.1415926535897932384626433832795029
is NOT giving you pi in double precision: the number is single
precision and
is then converted to double precision losing a couple of digits.
Terence correctly added the d0.
Regards,
Arjen
[*] I am not sure if macros like M_PI are part of the C standard.
Kevin G. Rhoads
>memorized pi to 100 places
I used to have slightly over a hundred, but nowadays I trail away in the thirties or forties
3.1415926535897932384626433832795028841971693993751058+
user1
Star Trek - original series - second season
"Computer, compute to the last digit the value of pi" -- Spock (Wolf in
the Fold)
Eugene
I used to define a global module and a unit module (if necessary) to
be used in my program. For example
module class_units
real(8),parameter::pi=2*asin(1.0d0)
real(8),parameter::e=exp(1)
real(8),parameter::inch=1
real(8),parameter::ft=12.*inch
end module class_units
Paul van Delst
Highly recommended. Here's the one I use (you presumably have your own equivalent "type
kinds" module):
MODULE Fundamental_Constants
USE Type_Kinds, ONLY: fp
IMPLICIT NONE
PRIVATE
REAL(fp), PARAMETER, PRIVATE :: ONE = 1.0_fp
REAL(fp), PARAMETER, PRIVATE :: TWO = 2.0_fp
REAL(fp), PARAMETER, PUBLIC :: PI = 3.141592653589793238462643_fp
REAL(fp), PARAMETER, PUBLIC :: PI_RECIPROCAL = 0.318309886183790671537767_fp
REAL(fp), PARAMETER, PUBLIC :: PI_SQUARED = 9.869604401089358618834491_fp
REAL(fp), PARAMETER, PUBLIC :: PI_SQUARE_ROOT = 1.772453850905516027298167_fp
REAL(fp), PARAMETER, PUBLIC :: PI_LN = 1.144729885849400174143427_fp
REAL(fp), PARAMETER, PUBLIC :: PI_LOG10 = 0.497149872694133854351268_fp
REAL(fp), PARAMETER, PUBLIC :: E = 2.718281828459045235360287_fp
REAL(fp), PARAMETER, PUBLIC :: E_RECIPROCAL = 0.367879441171442321595523_fp
REAL(fp), PARAMETER, PUBLIC :: E_SQUARED = 7.389056098930650227230427_fp
REAL(fp), PARAMETER, PUBLIC :: E_LOG10 = 0.434294481903251827651129_fp
REAL(fp), PARAMETER, PUBLIC :: SPEED_OF_LIGHT = 2.99792458e+08_fp
REAL(fp), PARAMETER, PUBLIC :: PERMEABILITY = PI * 4.0e-07_fp
REAL(fp), PARAMETER, PUBLIC :: PLANCK_CONSTANT = 6.62606876e-34_fp
REAL(fp), PARAMETER, PUBLIC :: GRAVITATIONAL_CONSTANT = 6.673e-11_fp
REAL(fp), PARAMETER, PUBLIC :: ELECTRON_VOLT = 1.602176462e-19_fp
REAL(fp), PARAMETER, PUBLIC :: UNIFIED_ATOMIC_MASS_UNIT = 1.66053873e-27_fp
REAL(fp), PARAMETER, PUBLIC :: STANDARD_ATMOSPHERE = 101325.0_fp
REAL(fp), PARAMETER, PUBLIC :: STANDARD_TEMPERATURE = 273.15_fp
REAL(fp), PARAMETER, PUBLIC :: STANDARD_GRAVITY = 9.80665_fp
REAL(fp), PARAMETER, PUBLIC :: AVOGADRO_CONSTANT = 6.02214199e+23_fp
REAL(fp), PARAMETER, PUBLIC :: MOLAR_GAS_CONSTANT = 8.314472_fp
REAL(fp), PARAMETER, PUBLIC :: STEFAN_BOLTZMANN_CONSTANT = 5.670400e-08_fp
REAL(fp), PARAMETER, PUBLIC :: PERMITTIVITY = &
ONE / (PERMEABILITY * SPEED_OF_LIGHT**2)
REAL(fp), PARAMETER, PUBLIC :: BOLTZMANN_CONSTANT = &
MOLAR_GAS_CONSTANT / AVOGADRO_CONSTANT
REAL(fp), PARAMETER, PUBLIC :: C_1 = &
TWO * PLANCK_CONSTANT * SPEED_OF_LIGHT**2
REAL(fp), PARAMETER, PUBLIC :: C_2 = &
PLANCK_CONSTANT*SPEED_OF_LIGHT/BOLTZMANN_CONSTANT
REAL(fp), PARAMETER, PUBLIC :: STP_MOLAR_VOLUME = &
(MOLAR_GAS_CONSTANT * STANDARD_TEMPERATURE) / STANDARD_ATMOSPHERE
REAL(fp), PARAMETER, PUBLIC :: LOSCHMIDT_CONSTANT = &
AVOGADRO_CONSTANT / STP_MOLAR_VOLUME
END MODULE Fundamental_Constants
Every now and again I add additional constants as required.
And every now and again one has to update the fundamental physical ones (see
http://physics.nist.gov/cuu/Constants) as the metrology improves
cheers,
paulv
John Reed
"How I want a drink, alcoholic of course, after the heavy lectures involving
quantum mechanics."
John Reed
"Paul van Delst" <Paul.v...@noaa.gov> wrote in message
news:gjtgtp$q5v$1...@news.nems.noaa.gov...
Phillip Helbig---remove CLOTHES to reply
<pcc4...@gmail.com> writes:
> As for myself, when I was in high school a friend of mine and I
> memorized pi to 100 places
I made it to 150!
Gordon Sande
You seem to be mising E_LN. ;-)
James Van Buskirk
news:gjtgtp$q5v$1...@news.nems.noaa.gov...
> REAL(fp), PARAMETER, PUBLIC :: PI =
> 3.141592653589793238462643_fp
I don't like this because it doesn't work if fp is the kind type
parameter of quadruple precision real. In f03 an initialization
expression may be:
"(4) A reference to an elemental standard intrinsic function, where
each argument is an initialization expression"
(N1601.pdf, section 7.1.7) so you could use
REAL(fp), PARAMETER, PUBLIC :: PI = ATAN2(Y=0.0_fp,X=-1.0_fp)
because
"if Y = 0 and X < 0, than the result is pi if Y is positive real zero"
(N1601.pdf, section 13.7.15). Now, this isn't valid f95, but ISTR
that CVF would take a code sequence such as:
REAL(fp) PI
PI = 4*ATAN(1.0_fp)
and precompute the value of PI so that ATAN would never get
invoked at run time, provided PI could not be changed in the
subroutine in question.
Another difficulty that may be encountered is that a module such as
this only gives you one KIND of PI. This may result in inadvertently
invoking the wrong specific version of a generic procedure if one is
not careful. I.e., it's safer to use:
call wwiii(2*real(PI,pecan))
than
call wwiii(2*PI)
--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end
George
There's more "logic" in Suzanne Summers portrayal of Chrissy in Three's
Company than the stuff that came out of Nimoy's mouth.
Never mind that the only math you have to do for the away team is to
identify whichever one you've never seen before to be the one who gets it.
I find very few situations that boil down to "logic." More often, I would
say it gets into "cases." I would bet that most persons would flub the
truth value of the statement "if false, ...."
--
George
Do I think faith will be an important part of being a good president? Yes,
I do.
Gib Bogle
I obviously do not belong in this company, because I never had the
slightest desire to memorize more than about 7.
George
> [*] I am not sure if macros like M_PI are part of the C standard.
I thought I would take a look myself at 1256.pdf and found:
const struct s { int mem; } cs = { 1 };
struct s ncs; // the object ncs is modifiable
typedef int A[2][3];
const A a = {{4, 5, 6}, {7, 8, 9}}; // array of array of const int
int *pi;
const int *pci;
ncs = cs; // valid
cs = ncs; // violates modifiable lvalue constraint for =
pi = &ncs.mem; // valid
pi = &cs.mem; // violates type constraints for =
pci = &cs.mem; // valid
pi = a[0]; // invalid: a[0] has type ‘‘const int *’’
This, of course, has nothing to do with the ratio of circumference to
diameter and follows the useful C notation to prepend pointers with p.
It looks like the notation for coarrays!
--
George
America is a Nation with a mission - and that mission comes from our most
basic beliefs. We have no desire to dominate, no ambitions of empire. Our
aim is a democratic peace - a peace founded upon the dignity and rights of
every man and woman.
user1
>
> I find very few situations that boil down to "logic." More often, I would
> say it gets into "cases." I would bet that most persons would flub the
> truth value of the statement "if false, ...."
Did you ever try getting first year college students to correctly interpret
compound logic statements?
Wait a minute ... What about simple percentages and proportions ?
Then there is what we used to call story problems ... If a fish weighs ten
pounds plus half it's weight, how much does the fish weigh ?
Richard Maine
> real(8),parameter::pi=2*asin(1.0d0)
I'm not sure who "Eugene" is. I presume not the Eugene that first comes
to my mind, because I can't believe he would do that. The formula with
atan is at least numerically "sensible", even if not guaranteed to give
last-bit accuracy. But the above is just asking for numerical problems.
If you get anything even half close, that's attributable more to luck
than design. It wouldn't be completely implausible for the above to give
you an error. Consider that if the 1.0d0 argument value were even one
bit larger, then this pretty unambiguously should give you an error. And
you might want to look at a sensitivity analysis of teh behavior on the
other side of 1.0d0 as well.
Richard Maine
> Eugene <yjy...@gmail.com> wrote:
>
> > real(8),parameter::pi=2*asin(1.0d0)
>
> I'm not sure who "Eugene" is. I presume not the Eugene that first comes
> to my mind, because I can't believe he would do that. The formula with
> atan is at least numerically "sensible", even if not guaranteed to give
> last-bit accuracy. But the above is just asking for numerical problems.
> If you get anything even half close, that's attributable more to luck
> than design. It wouldn't be completely implausible for the above to give
> you an error. Consider that if the 1.0d0 argument value were even one
> bit larger, then this pretty unambiguously should give you an error. And
> you might want to look at a sensitivity analysis of teh behavior on the
> other side of 1.0d0 as well.
(Following up to myself, I know). To illustrate, I wrote the following
code and ran it (with an oldish g95 on my Mac).
program p
integer, parameter :: dp = kind(1.0d0)
real(dp) :: pi, x, eps
eps = epsilon(x)
x = 1.0d0
pi = 2*asin(x)
write (*,*) x, pi
x = 1.0d0 - eps
pi = 2*asin(x)
write (*,*) x, pi
x = 1.0d0 + eps
pi = 2*asin(x)
write (*,*) x, pi
end
The resulting output is
1. 3.141592653589793
0.9999999999999998 3.1415926114429444
1.0000000000000002 NaN
So the nominal value gives a decent result. But look at the sensitivity!
If you try an input value that is epsilon less, the input value is the
same to (if I counted correctly) 15 places, as expected. The result
differs in the 9th place. So that function looses you 6 orders of
magnitude of accuracy in that region. Us Fortran folk are supposed to be
sensitive to things like numerical instability; that's a prime example
of it.
And if you try an input value that is epsilon more, you indeed get a
NaN.
Ron Shepard
"James Van Buskirk" <not_...@comcast.net> wrote:
> > REAL(fp), PARAMETER, PUBLIC :: PI =
> > 3.141592653589793238462643_fp
>
> I don't like this because it doesn't work if fp is the kind type
> parameter of quadruple precision real.
I don't understand, why not?
$.02 -Ron Shepard
George
James has a module that gives qp the maximum precision, which would be
something that might be useful here. I wonder what Paul means with fp.
module mykinds
implicit none
integer, parameter :: dp = selected_real_kind(15,300)
integer, parameter :: qp_preferred = selected_real_kind(30,1000)
integer, parameter :: qp = (1+sign(1,qp_preferred))/2*qp_preferred+ &
(1−sign(1,qp_preferred))/2*dp
end module mykinds
Pi to quad precision is longer than 25 digits.
--
George
We know that dictators are quick to choose aggression, while free nations
strive to resolve differences in peace.
glen herrmannsfeldt
> Eugene <yjy...@gmail.com> wrote:
>> real(8),parameter::pi=2*asin(1.0d0)
> I'm not sure who "Eugene" is. I presume not the Eugene that first comes
> to my mind, because I can't believe he would do that. The formula with
> atan is at least numerically "sensible", even if not guaranteed to give
> last-bit accuracy. But the above is just asking for numerical problems.
Most do some type of argument reduction, similar to that done
by the regular trig functions. One that I know of does:
ARCSIN-ARCCOSINE FUNCTION (LONG)
1. IF X BETWEEN 0 AND 1/2, COMPUTE ARCSIN BY POLYNOMIAL
2. IF X BETWEEN 1/2 AND 1,
ARSIN(X) = PI/2-2*ARSIN(SQRT((1-X)/2))
3. IF X NEGATIVE, ARSIN(X) = -ARSIN(/X/)
4. ARCOS(X) = PI/2-ARSIN(X)
ASIN(1.D0) then depends on the internal PI/2 and the polynomial
giving zero at zero. Avoiding a discontinuity at zero likely
requires that asin(0.d0) be zero.
> If you get anything even half close, that's attributable more to luck
> than design. It wouldn't be completely implausible for the above to give
> you an error. Consider that if the 1.0d0 argument value were even one
> bit larger, then this pretty unambiguously should give you an error. And
> you might want to look at a sensitivity analysis of teh behavior on the
> other side of 1.0d0 as well.
Well, that is true but it would be a poor floating point
implementation that got 1.D0 wrong, or the comparison to
it wrong. If one tried 2.D0*asin(0.99D0+0.01D0) then
I would completely agree.
-- glen
glen herrmannsfeldt
(snip)
If you use ARSIN(X) = PI/2-2*ARSIN(SQRT((1-X)/2))
For x=1-eps, and asin(x)=x for small x,
it is PI/2-2*sqrt(eps/2), for eps around 2e-16,
the relative change should be about 1e-8
I would be surprised, though, if one didn't do either
a similar argument reduction or a polynomial expansion
around x=1, for larger arguments. An interesting question
is how they do it in hardware. It is much easier to put
in the tests and such in software than hardware.
-- glen
robert....@sun.com
> I have more confidence in the robust precision of compiler
> decimal-to-binary conversions than in triginometric functions.
I agree that spelling out constants is more likely to produce
the most accurate results than computing the values using
instrinsic functions. The case of the inverse trig functions
is especially bad because there are vendors who claim that the
Fortran standard prohibits them from producing the most accurate
representable result. They take the requirement that the
result be less than or equal to pi/2 to mean that the result
must be off by almost a full unit in the last place in some
cases. Those cases do occur in practice.
Bob Corbett
robert....@sun.com
> Pi to quad precision is longer than 25 digits.
Thirty-five digits are required if quadruple-precision means
113 significand bits and the value is rounded to nearest.
Thirty-six digits are needed to represent an arbitrary value
if there are 113 significand bits.
Bob Corbett
James Van Buskirk
news:ron-shepard-5266...@news60.forteinc.com...
The problem is that quad precision has more than 25 significant
figures. Here's a fun way to calcualte pi:
We know that cos(theta/2) = sqrt((1+cos(theta))/2). More conveniently,
let C(N) = 2*cos(theta/2**N), so C(N+1) = sqrt(2+C(N)). Then let
S(N) = sin(theta/2**N), so S(N+1) = S(N)/C(N+1). If theta = pi/2,
then C(0) = 0 and S(0) = 1, so
S(N) = S(0)/product([(C(i),i=1,N)]) = sin(pi/2**(N+1)). Now,
2**(N+1)*sin(pi/2**(N+1)) = pi - pi**3/(24*2**(2*N)) + O(1/2**(4*N)).
We may improve this a little by noting that
2**(N+1)*(4*S(N)/3-S(N-1)/6) = pi - pi**5/(480*2**(4*N)) + O(1/2**(6*N)).
This is an acceleration of the Archimedes method of calculating pi.
Finding bounds by methods available to Archimedes such that
A/2**(4*N) < 2**(N+1)*(4*sin(pi/2**(N+1))/3-sin(pi/2**N)/6) - pi < B/2**4*N
is left as an exercise to the reader. We just crunch this out to N = 30
and write out our approximation to pi and the O(1/2**(4*N)) error term.
The last two lines will be pi given by the named constant above and
pi from atan2(0.0_fp,-1.0_fp). The reader may decide for himself
whether the named constant is suitable for quad precision pi:
C:\Program Files\Microsoft Visual Studio 8\James\clf\archpi2>type
archpi2.f90
module mykinds
implicit none
integer, parameter :: fp = selected_real_kind(30,1000)
end module mykinds
program archpi2
use mykinds
implicit none
REAL(fp), PARAMETER :: PI=3.141592653589793238462643_fp
real(fp) CN
real(fp) S0
real(fp) denominator
integer n
S0 = 1
CN = 0
denominator = 6
do n = 1, 30
CN = sqrt(2+CN)
denominator = denominator*CN
write(*,*) 2.0_fp**(n+1)*S0*(8-CN)/denominator, &
real(PI**5/480/2.0_fp**(4*n))
end do
write(*,*) PI
write(*,*) atan2(0.0_fp,-1.0_fp)
end program archpi2
C:\Program Files\Microsoft Visual Studio 8\James\clf\archpi2>ifort
archpi2.f90
Intel(R) Fortran Compiler for Intel(R) EM64T-based applications, Version 9.1
Build 20061104
Copyright (C) 1985-2006 Intel Corporation. All rights reserved.
Microsoft (R) Incremental Linker Version 8.00.40310.39
Copyright (C) Microsoft Corporation. All rights reserved.
-out:archpi2.exe
-subsystem:console
archpi2.obj
C:\Program Files\Microsoft Visual Studio 8\James\clf\archpi2>archpi2
3.10456949966158679680450326455920 3.9846312E-02
3.13914757031222753256911401351779 2.4903945E-03
3.14143771670383032282085057009541 1.5564966E-04
3.14158293664190158989482476070460 9.7281036E-06
3.14159204575769079515140535096341 6.0800647E-07
3.14159261559211285331032018627372 3.8000405E-08
3.14159265121481047908401348901303 2.3750253E-09
3.14159265344135482007156179387541 1.4843908E-10
3.14159265358051580612653898017897 9.2774425E-12
3.14159265358921339845434426408714 5.7984016E-13
3.14159265358975699845451073762955 3.6240010E-14
3.14159265358979097346201612493138 2.2650006E-15
3.14159265358979309690010232075776 1.4156254E-16
3.14159265358979322961498453782697 8.8476587E-18
3.14159265358979323790966470498489 5.5297867E-19
3.14159265358979323842808221587900 3.4561167E-20
3.14159265358979323846048331031686 2.1600729E-21
3.14159265358979323846250837871934 1.3500456E-22
3.14159265358979323846263494549449 8.4377848E-24
3.14159265358979323846264285591794 5.2736155E-25
3.14159265358979323846264335031941 3.2960097E-26
3.14159265358979323846264338121950 2.0600061E-27
3.14159265358979323846264338315075 1.2875038E-28
3.14159265358979323846264338327146 8.0468987E-30
3.14159265358979323846264338327900 5.0293117E-31
3.14159265358979323846264338327947 3.1433198E-32
3.14159265358979323846264338327950 1.9645749E-33
3.14159265358979323846264338327950 1.2278593E-34
3.14159265358979323846264338327950 7.6741206E-36
3.14159265358979323846264338327950 4.7963254E-37
3.14159265358979323846264300000000
3.14159265358979323846264338327950
Peter C. Chapin
>> Anyway, thanks for the replies. I'll set up my constant myself. Perhaps
>> a module full of such constants is in order. Hmm...
>
> Highly recommended. Here's the one I use (you presumably have your own
> equivalent "type kinds" module):
>
> MODULE Fundamental_Constants...
Thanks for posting that. I'll take a closer look.
Peter
James Parsly
Go to the mall parking lot on a busy day.
Look at the license plates on the cars you find there. In my state, these
are mostly of the form:
ABC 129
For every two license plates you examine, treat the numeric parts as an
(X,Y) pair. Append a decimal point
on the front of each number. Compute x**2 + y**2. If the result is <= 1
add 1 to a tally. Also keep
track of the total number of (X,Y) pairs. (alternatively you could plot the
point on a graph on which you have
drawn a quarter circle with radius one and visually judge whether the point
lies inside the quarter circle).
Do this for a large number of license plate pairs.
When you think you have enough points, divide the smaller tally by the
bigger one. Multiply by 4. The
result is somewhere close to pi.
"James Van Buskirk" <not_...@comcast.net> wrote in message
news:gjvbdo$prs$1...@news.motzarella.org...
Craig Powers
> Here's an even more fun way to calculate pi:
>
> Go to the mall parking lot on a busy day.
>
> Look at the license plates on the cars you find there. In my state, these
> are mostly of the form:
>
> ABC 129
>
> For every two license plates you examine, treat the numeric parts as an
> (X,Y) pair. Append a decimal point
> on the front of each number. Compute x**2 + y**2. If the result is <= 1
> add 1 to a tally. Also keep
> track of the total number of (X,Y) pairs. (alternatively you could plot the
> point on a graph on which you have
> drawn a quarter circle with radius one and visually judge whether the point
> lies inside the quarter circle).
>
> Do this for a large number of license plate pairs.
>
>
> When you think you have enough points, divide the smaller tally by the
> bigger one. Multiply by 4. The
> result is somewhere close to pi.
I can foresee two problems with this approach.
The first is that Monte Carlo integration is, as I recall, notoriously
inefficient at computing the value of pi.
The second is that there is an implicit reliance on the numeric part of
the tag numbers being statistically random, and I'm not sure how
reliable that assumption is.
James Parsly
that
purpose.
I can certainly think of locations where it wouldn't work, but we live in a
fairly big city with a mall just off
the interstate. I did actually check with the state to see if all numbers
where issued. Even so, the
number 666 didn't ever show while we were sampling, and while I have seen it
since, it wouldn't
surprise me if plates with that number sometimes get exchanged and never
reissued.
Nevertheless, I find the approach amusing. Sort of a "rube-goldberg" method
for getting at pi.
"Craig Powers" <craig....@invalid.invalid> wrote in message
news:gk08bv$k6o$1...@news.motzarella.org...
Ron Shepard
"James Van Buskirk" <not_...@comcast.net> wrote:
> "Ron Shepard" <ron-s...@NOSPAM.comcast.net> wrote in message
> news:ron-shepard-5266...@news60.forteinc.com...
>
> > In article <gjtl4c$p01$1...@news.motzarella.org>,
> > "James Van Buskirk" <not_...@comcast.net> wrote:
>
> >> > REAL(fp), PARAMETER, PUBLIC :: PI =
> >> > 3.141592653589793238462643_fp
>
> >> I don't like this because it doesn't work if fp is the kind type
> >> parameter of quadruple precision real.
>
> > I don't understand, why not?
>
> The problem is that quad precision has more than 25 significant
> figures.
Ok, I thought by "doesn't work" you were making a statement about
fortran standards or syntax or something, not numerical accuracy.
> Here's a fun way to calcualte pi:
There really isn't any need for this these days. There are web sites
all over that have millions of digits of Pi. Here are the first 100:
3.14159265358979323846264338327950288419716939937510582097494459230781\
6406286208998628034825342117068
If accuracy is the issue, just use as many digits (or more) as needed.
You can fit 100 digits on a single free format line in fortran without
worrying about continuations or anything. If you need constants with
different precisions, then set the ones with the most digits first, and
then use those to set the other ones.
$.02 -Ron Shepard
glen herrmannsfeldt
(snip)
> There really isn't any need for this these days. There are web sites
> all over that have millions of digits of Pi. Here are the first 100:
> 3.14159265358979323846264338327950288419716939937510582097494459230781\
> 6406286208998628034825342117068
There is a book, "A History of PI," by Petr Beckmann, written
in 1971. In the back, it has 10000 of the 100265 digits of Pi
computed on an IBM 704 in July 1971. (Or so it says. In the
text is says that was on a 7090, but that previously 16167 digits
were computed on a 704.) Since the 704 and 7090 were commonly
programmed in Fortran, it is a good guess that is how this was
done. (Or maybe Fortran calling assembly subroutines.)
-- glen
Luka Djigas
Finally,
I was starting to wonder whether I am the only weird guy here, who
only knows 5-10 digits, depending on the state of alcohol after some
heavy lectures involving quantum mechanics :)
-- Luka
George
I don't know, Glen, I think if you're going to cook up your own pie, use
arctangent:
implicit none
INTEGER,PARAMETER::dp=SELECTED_REAL_KIND(15)
real(kind=dp) :: pi, pie, eps, x
pie = 3.141592653589793238462643_dp
pi = 4.0_dp * atan(1.0_dp)
print *, "pie = ", pie
print *, "pi = ", pi
eps = epsilon(x)
x = 1.0d0 - eps
x = 4.0_dp * atan(x)
write (*,*) "less epsilon ", x
x = 1.0d0 + eps
x = 4.0_dp * atan(x)
write (*,*) "plus epsilon ", x
print *, " 1.23456789012345"
endprogram
! g95 eps1.f03 -o p.exe
C:\MinGW\source>g95 eps1.f03 -o p.exe
C:\MinGW\source>p
pie = 3.141592653589793
pi = 3.141592653589793
less epsilon 3.1415926535897927
plus epsilon 3.1415926535897936
1.23456789012345
C:\MinGW\source>
I'm uncertain why they're different lengths. If you average the first two
values and average the second two values, the second is closer to the
second number shown to be transcendental over Q.
--
George
Natural gas is hemispheric. I like to call it hemispheric in nature because
it is a product that we can find in our neighborhoods.
James Van Buskirk
> 3.14159265358979323846264338327950288419716939937510582097494459230781\
> 6406286208998628034825342117068
How do you check a constant written out like that, except
painstakingly? It's a long winter night, come sit down by my
fireplace, I have a little story to tell you:
C:\Program Files\Microsoft Visual Studio 8\James\clf\archpi2>type gamma.f90
module mykinds
implicit none
integer, parameter :: fp = selected_real_kind(30,1000)
end module mykinds
program get_gamma
use mykinds
implicit none
integer i
real(fp) L2
real(fp) Delta(24)
integer, parameter :: L2MAX = 100
integer j
real(fp) A(11)
real(fp) f
real(fp) gamma
integer, parameter :: q = 7
real(fp) L2diff
real(fp) cleanup
real(fp), parameter :: gamma_const = 0.57721566490153386061_fp
L2 = 0
do i = 1, size(Delta)
Delta(i) = 1.0_fp/(L2MAX+i)
do j = i-1, 1, -1
Delta(j) = Delta(j+1)-Delta(j)
end do
L2 = L2-(-1)**i*Delta(1)/2.0_fp**i
end do
do i = 1, L2MAX
L2 = L2-real((-1)**i,fp)/i
end do
write(*,*) L2
A(1) = 1.0_fp/12
do i = 2, size(A), 2
A(i) = 0
do j = 1, i/2-1
A(i) = A(i)+A(j)*A(i-j)
end do
A(i) = -(2*A(i)+A(i/2)**2)/(2*i+1)
A(i+1) = 0
do j = 1, i/2
A(i+1) = A(i+1)+A(j)*A(i+1-j)
end do
A(i+1) = -2*A(i+1)/(2*i+3)
end do
f = 1
do i = 1, size(A)
A(i) = f*A(i)
f = f*(2*i)*(2*i+1)
write(*,*) A(i)*2*i
end do
gamma = 1
do i = 1, q
L2diff = 0
do j = 2**(i-1)+1, 2**i
L2diff = L2diff+1.0_fp/j
end do
L2diff = L2diff-L2
gamma = gamma+L2diff
end do
cleanup = -1/2.0_fp**(q+1)
do i = 1, size(A)
cleanup = cleanup+A(i)/2.0_fp**(2*i*q)
end do
gamma = gamma+cleanup
write(*,*) gamma
write(*,*) gamma_const
end program get_gamma
C:\Program Files\Microsoft Visual Studio 8\James\clf\archpi2>ifort gamma.f90
Intel(R) Fortran Compiler for Intel(R) EM64T-based applications, Version 9.1
Build 20061104
Copyright (C) 1985-2006 Intel Corporation. All rights reserved.
Microsoft (R) Incremental Linker Version 8.00.40310.39
Copyright (C) Microsoft Corporation. All rights reserved.
-out:gamma.exe
-subsystem:console
gamma.obj
C:\Program Files\Microsoft Visual Studio 8\James\clf\archpi2>gamma
0.693147180559945309417232121458177
0.166666666666666666666666666666667
-3.333333333333333333333333333333333E-0002
2.380952380952380952380952380952381E-0002
-3.333333333333333333333333333333333E-0002
7.575757575757575757575757575757574E-0002
-0.253113553113553113553113553113553
1.16666666666666666666666666666667
-7.09215686274509803921568627450980
54.9711779448621553884711779448621
-529.124242424242424242424242424242
6192.12318840579710144927536231884
0.577215664901532860606512090082402
0.577215664901533860610000000000000
The last line is gamma_const. Bonus points if you can tell me where
I got that from!
George
> The last line is gamma_const. Bonus points if you can tell me where
> I got that from!
There's a lot on the net about this, including a page that has 170,000
digits of this constant by one of the Borweins, who are the champs of this
stuff. Has anyone ever tried to compute pi by using the
arithmetic-geometric mean?
http://www.worldwideschool.org/library/books/sci/math/MiscellaneousMathematicalConstants/chap35.html
Unusual to have an entry in Dr. Sloane's encyclopedia too:
http://www.research.att.com/~njas/sequences/A001620
module mykinds
implicit none
integer, parameter :: dp = selected_real_kind(15,300)
integer, parameter :: qp_preferred = selected_real_kind(30,1000)
integer, parameter :: fp = (1+sign(1,qp_preferred))/2*qp_preferred+ &
(1−sign(1,qp_preferred))/2*dp
end module mykinds
program get_gamma
use mykinds
implicit none
integer i
real(fp) L2
real(fp) Delta(24)
integer, parameter :: L2MAX = 100
integer j
real(fp) A(11)
real(fp) f
real(fp) gamma
integer, parameter :: q = 7
real(fp) L2diff
real(fp) cleanup
real(fp), parameter :: gamma_const = 0.57721566490153386061_fp
L2 = 0
print *, "size delta = ", size(Delta)
do i = 1, size(Delta)
Delta(i) = 1.0_fp/(L2MAX+i)
do j = i-1, 1, -1
Delta(j) = Delta(j+1)-Delta(j)
end do
L2 = L2-(-1)**i*Delta(1)/2.0_fp**i
end do
do i = 1, L2MAX
L2 = L2-real((-1)**i,fp)/i
end do
write(*,*) "L2 is ", L2
A(1) = 1.0_fp/12
do i = 2, size(A), 2
A(i) = 0
do j = 1, i/2-1
A(i) = A(i)+A(j)*A(i-j)
end do
A(i) = -(2*A(i)+A(i/2)**2)/(2*i+1)
A(i+1) = 0
do j = 1, i/2
A(i+1) = A(i+1)+A(j)*A(i+1-j)
end do
A(i+1) = -2*A(i+1)/(2*i+3)
end do
f = 1
print *, "size of A is ", size(A)
do i = 1, size(A)
A(i) = f*A(i)
f = f*(2*i)*(2*i+1)
write(*,*) "expression is ", A(i)*2*i
end do
gamma = 1
do i = 1, q
L2diff = 0
do j = 2**(i-1)+1, 2**i
L2diff = L2diff+1.0_fp/j
end do
L2diff = L2diff-L2
gamma = gamma+L2diff
end do
cleanup = -1/2.0_fp**(q+1)
do i = 1, size(A)
cleanup = cleanup+A(i)/2.0_fp**(2*i*q)
end do
gamma = gamma+cleanup
write(*,*) "gamma is ", gamma
write(*,*) "gamma_const is ", gamma_const
end program get_gamma
! g95 gamma2.f03 -o w.exe
C:\MinGW\source>g95 gamma2.f03 -o w.exe
C:\MinGW\source>w
size delta = 24
L2 is 0.6931471805599453094172321214581767
size of A is 11
expression is 0.16666666666666666666666666666666666
expression is -0.033333333333333333333333333333333327
expression is 0.023809523809523809523809523809523808
expression is -0.033333333333333333333333333333333327
expression is 0.07575757575757575757575757575757574
expression is -0.25311355311355311355311355311355306
expression is 1.1666666666666666666666666666666665
expression is -7.092156862745098039215686274509802
expression is 54.97117794486215538847117794486214
expression is -529.1242424242424242424242424242422
expression is 6192.1231884057971014492753623188385
gamma is 0.5772156649015328606065120900824016
gamma_const is 0.57721566490153386061
C:\MinGW\source>
> The last line is gamma_const. Bonus points if you can tell me where
> I got that from!
Reverse engineering a subroutine from Numerical Recipes?
--
George
Hundreds of thousands of American servicemen and women are deployed across
the world in the war on terror. By bringing hope to the oppressed, and
delivering justice to the violent, they are making America more secure.
me...@skyway.usask.ca
Chris
Michel Olagnon
That's only 2 digits, actually. Using 3, I prefer 3 5 5 / 1 1 3
AnotherSquid
355 / 113 = 3.141592920 . The error is 0.000000266. That's seven
correct decimal digits from three! A numerical accident, but still
amazes me. According to Carl B. Boyer in "The History of Mathematics",
this fraction was first discovered by a Chinese mathematician, Tsu
Ch'ung-chih, who lived from 430 to 501 AD. He put this forth as an
"accurate but not exact" value. The West did not achieve this accuracy
until the 15th Century.
glen herrmannsfeldt
>> That's only 2 digits, actually. Using 3, I prefer 3 5 5 / 1 1 3
> 355 / 113 = 3.141592920 . The error is 0.000000266. That's seven
> correct decimal digits from three!
Well, I would include both the numerator and denominator when
counting significant digits, log10(355)+log10(113) is about 4.6,
which is still pretty good.
Richard Maine
> All you really need is 3 digits: 2 2 / 7 :-}
Oncw you have a typo on the first one, the rest are pretty much
irrelevant. :-)
Gordon Sande
> <me...@skyway.usask.ca> wrote:
>
>> All you really need is 3 digits: 2 2 / 7 :-}
>
> Oncw you have a typo on the first one, the rest are pretty much
> irrelevant. :-)
The typo is the extra blank as 3 1/7 is also known as 22 / 7.
3 1/7 is around 3.14 which is good enough for government work. ;-)
There is a well known (and often reduculed) approximation of 4
used to calculate frontages of lots with circular sections
for some municipal taxes.
Gib Bogle
Wonderful! Thanks.
George
> AnotherSquid <m...@ucar.edu> wrote:
>
>>> That's only 2 digits, actually. Using 3, I prefer 3 5 5 / 1 1 3
>
>> 355 / 113 = 3.141592920 . The error is 0.000000266. That's seven
>> correct decimal digits from three!
>
> Well, I would include both the numerator and denominator when
> counting significant digits, log10(355)+log10(113) is about 4.6,
> which is still pretty good.
I think I see your point here, Glen. When you think of representations of
pi in continued fractions, the numbers are smallish integers.
I'm certain that with the right fractions and single digits, I could give
you double precision.
How does your analysis work with continued fractions?
>
>> A numerical accident, but still
>> amazes me. According to Carl B. Boyer in "The History of Mathematics",
>> this fraction was first discovered by a Chinese mathematician, Tsu
>> Ch'ung-chih, who lived from 430 to 501 AD. He put this forth as an
>> "accurate but not exact" value. The West did not achieve this accuracy
>> until the 15th Century.
In Wolfram's mathworld, this is represented as
[ 3 7 15 1]
http://mathworld.wolfram.com/PiContinuedFraction.html
He is there identified as an astronomer.
Just looking at the sky like the rest of us.
--
George
The deliberate and deadly attacks which were carried out yesterday against
our country were more than acts of terror. They were acts of war.
George
I think the state of Indiana passed into law a means to calculate pi as a
rational. Since they voted blue this time, I think we can assume that
they've since come back to their senses.
--
George
America is a friend to the people of Iraq. Our demands are directed only at
the regime that enslaves them and threatens us. When these demands are met,
the first and greatest benefit will come to Iraqi men, women and children.
Charles Russell
> I obviously do not belong in this company, because I never had the
> slightest desire to memorize more than about 7.
How I like a drink alcoholic after solving equations ... I forget the
rest, something about magna quanta, but that much gets eight digits.
Not so long ago we had
- quote -
Re: Pi.
From:
"Michael Metcalf" <michael...@compuserve.com>
Date:
Tue, 5 Jul 2005 15:16:15 +0200
Newsgroups:
comp.lang.fortran
>> You have to define it yourself, and it will be defined in whatever
>> precision you chose....
>>
>> e.g.
>>
>>
>> INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(P=11)
>> REAL (KIND=DP), PARAMETER :: pi =
>> 3.141592653589793238462643383279502884197_DP
>>
This risks mistyping. You can get it to machine precision by:
double precision pi
pi = acos(-1.0d0)
This has to be, however, a variable rather than a named constant (in
Fortran
95 but not Fortran 2003).
Regards,
Mike Metcalf
- endquote -
Richard Maine
> This risks mistyping. You can get it to machine precision by:
>
> double precision pi
> pi = acos(-1.0d0)
You apparently missed the whole point of the prior discussion (before it
wandered off into amusements), which was that expressions like this do
*NOT* guarantee machine precision. You might happen to get it, but it is
by no means guaranteed.
James Van Buskirk
news:1it8ntw.1putxwmgffdewN%nos...@see.signature...
> Charles Russell <NOS...@bellsouth.net> wrote:
>> This risks mistyping. You can get it to machine precision by:
>> double precision pi
>> pi = acos(-1.0d0)
> You apparently missed the whole point of the prior discussion (before it
> wandered off into amusements), which was that expressions like this do
> *NOT* guarantee machine precision. You might happen to get it, but it is
> by no means guaranteed.
That wasn't the whole point. I showed that also expressions like:
>> REAL (KIND=DP), PARAMETER :: pi =
>> 3.141592653589793238462643383279502884197_DP
do not guarantee machine precision either. In my example,
http://groups.google.com/group/comp.lang.fortran/msg/b038cd6ba752c39e?hl=en
I showed that uncritically copying down a value, such as
> real(fp), parameter :: gamma_const = 0.57721566490153386061_fp
from a respected source (in this instance the 47th CRC Handbook of
Chemistry and Physics) can lead to incorrect results as well. I
see too much of a tendency for people to tell you that you should
trust authority rather than think for yourself.
George
Despite MM here, I like atan(ONE).
Richard brought up what happens to the sin(ONE) by missing a bit. Glen's
amelioration for the phase shift means that what applies to sin(ONE)
attends as well to cos(-ONE). Certainly, that's covered ground, though.
As for James' ultimate claim, my next guess is that he used an arbitrary
small number to build the gamma constant.
--
George
I understand everybody in this country doesn't agree with the decisions
I've made. And I made some tough decisions. But people know where I stand.
jfh
> On Wed, 07 Jan 2009 18:19:07 GMT, Gordon Sande wrote:
> > On 2009-01-07 12:42:05 -0400, nos...@see.signature (Richard Maine) said:
>
>
> >>> All you really need is 3 digits: 2 2 / 7 :-}
>
> >> Oncw you have a typo on the first one, the rest are pretty much
> >> irrelevant. :-)
>
> > The typo is the extra blank as 3 1/7 is also known as 22 / 7.
>
> > 3 1/7 is around 3.14 which is good enough for government work. ;-)
> > There is a well known (and often reduculed) approximation of 4
> > used to calculate frontages of lots with circular sections
> > for some municipal taxes.
>
> I think the state of Indiana passed into law a means to calculate pi as a
> rational. Since they voted blue this time, I think we can assume that
> they've since come back to their senses.
> --
> George
Indiana passed no such law. The state legislature was considering it
when an Indiana mathematics professor was visiting for another reason.
He spoke to a representative, who got the debate adjourned. So far as
I know, it has never re-started.
John Harper
George
It's funny that you, as New Zealander, know the truth of this, while I, who
grew up a half hour away in Ohio seem not to.
My family's in and from Chicago, so we'd drive through Indiana on route 6
and share the road with the Amish. It's somehow easy for me to believe
that there could exist there an irrational attachment to the rationals.
--
George
This way of life is worth defending.
Kevin G. Rhoads
>that there could exist there an irrational attachment to the rationals.
Irrational preferences are not limited to just the rationals over the irrationals,
but also to the non-transcendentals over the transcendentals (except, perhaps, in
certain areas of southern California).
Now, anyone have ideas on how one implements some reasonable facsimile of sur-real
or hyper-real numbers in finite hardware? ... Perhaps, sign-magnitude with mantissa
and infinitude/finiteness/infinitesimalness ...
jfh
Math. Intelligencer 7(2), 69-72. Anyone interested should read that.
It will be more reliable than my 24-year-old memory of what it said.
John Harper
George
The Mathematical Intelligencer is pound for pound some of the best reading
in mathematics.
--
George
You tell it that it's indicative by appending $!. That's why we made $!
such a short variable name, after all. :-)
-- Larry Wall in <1997090818...@wall.org>
Charles Coldwell
> I have more confidence in the robust precision of compiler
> decimal-to-binary conversions than in triginometric functions.
I don't.
http://sources.redhat.com/bugzilla/show_bug.cgi?id=3479
--
Charles M. "Chip" Coldwell
"Turn on, log in, tune out"
GPG Key ID: 852E052F
GPG Key Fingerprint: 77E5 2B51 4907 F08A 7E92 DE80 AFA9 9A8F 852E 052F