Gpc February 2002

gpc@gnu.de

42 participants
71 discussions

Re: Complex_Arctan
by Maurice Lombardi 11 Feb '02

11 Feb '02

Frank Heckenbach <frank(a)g-n-u.de> a dit : > da Silva, Joe wrote: > > As for my comment about "tweaked" coefficients (for platforms which > > don't already provide the function), the general rule is that if you > > truncate > > an infinite series, you need to "tweak" the coefficient values to provide > > optimum accuracy for the number of terms you are using. Simply > > truncating a series may produce inaccurate results, regardless of how > > many terms you apply. That is the general rule, and I remember it was > > certainly true with trigonometric functions. However, perhaps this effect > > is not so prevalent with the "log1p" coefficients? > > As I said, I'm no expert on numerics, and I haven't heard of this > tweaking. I suppose that what Joe has in mind is the replacement of a Taylor series by a Tchebychev polynomial. The drawback of a Taylor series is that to have a given accuracy you need only one term at the origin, but a number of terms which increases when you go farther away from the origin of the development. To have the accuracy you need at the farthest point of the development you need a polynomial of too high a degree, and you get in addition problems of roundoff errors when evaluating this polynomial. A Tchebycheff polynomial have a constant accuracy over the whole range of validity, and for a given accuracy over a given range it is of lower degree than any Taylor series based polynomial in the same range. Even better is the use of a rational function (ratio of two polynomials) as is used in the CEPHES library, for the log1p function in particular. However the degree and coefficients of this best polynomial depends both on the interval sought and on the accuracy needed. It is different e.g. for single, double, long double, and you cannot simply truncate a given series to different orders. The situation is even worse for a general purpose compiler like gpc which is aimed to work with different word sizes: 64 bits for real or double for most systems, but may be 128 bits for some. Checks would be needed to select the polynomial to use. It is probably better to let this done when possible by a libc library (or an other specialized library like CEPHES) which knows the particular system for which it is built and selects the order and coefficients of the polynomial accordingly. Gpc should give only a less optimized but the most generally valid solution. Maurice -- Maurice Lombardi Laboratoire de Spectrometrie Physique, Universite Joseph Fourier de Grenoble, BP87 38402 Saint Martin d'Heres Cedex FRANCE Tel: 33 (0)4 76 51 47 51 Fax: 33 (0)4 76 63 54 95 mailto:Maurice.Lombardi@ujf-grenoble.fr

3 2

Re: Complex_Arctan
by da Silva, Joe 11 Feb '02

11 Feb '02

Hmmm ... Maybe I'm doing something wrong, but my installation of GPC (DJGPP) doesn't know what "Log1P" is, nor can I find it in the source code of the GPC unit, etc. Anyway, that's my problem, so I'll look into that further when I can. As for my comment about "tweaked" coefficients (for platforms which don't already provide the function), the general rule is that if you truncate an infinite series, you need to "tweak" the coefficient values to provide optimum accuracy for the number of terms you are using. Simply truncating a series may produce inaccurate results, regardless of how many terms you apply. That is the general rule, and I remember it was certainly true with trigonometric functions. However, perhaps this effect is not so prevalent with the "log1p" coefficients? If so, that's good - it means we don't have to go to the trouble of finding "tweaked/optimized" coefficients (or find out how to derive them). Unfortunately, I don't remember the name of the book in which I found the "tweaked" coefficients for trigonometric functions, nor whether this also had similar coefficients for other functions (it was quite a few years ago now ;-). Joe. > -----Original Message----- > From: Emil Jerabek [SMTP:ejer5183@artax.karlin.mff.cuni.cz] > Sent: Saturday, February 09, 2002 12:47 AM > To: gpc(a)gnu.de > Subject: Re: Complex_Arctan > > Joe wrote: > > > > > BTW, AFAICT, DJGPP does not have a Log1P function, > > so I wasn't able to verify the above anyway. I didn't know > > the values to use for the coefficients "Log1PCoefs[]", so > > I couldn't try the suggested Log1P code, either. Although > > ??? You don't have to guess the coefficients, they are initialized in the > code. > (If you don't want to rebuild the RTS just to test it (I didn't either), > move the > "to begin do" stuff into your main program, and import the GPC unit > (because of the > SplitReal function).) > > > Emil presented a power series which I assume represents > > these coefficients, I don't know whether these need to be > > "tweaked". If I understand this correctly, it is an infinite > > series. However, AFAIK, when you apply a finite version > > of such a series, you need to "tweak" it's values to get > > good accuracy (no, I don't know how to do this). > > > > I remember experimenting years ago with the power > > series for "sine" or "cosine" and being amazed at how > > a small finite "tweaked" power series was much more > > accurate than the "conventional" power series, no > > matter how many terms I calculated for it (after a while, > > the more terms you calculate, the greater the error > > you introduce). I presume the "log1p" power series > > also has this effect(?) > > > > Not sure what sort of "tweaking" you have in mind. Anyway, my Log1P > implementation > tries to use as few terms of the power series as possible, and it sums > them in reverse order. > I don't know if this is optimal, but it appeared to work with acceptable > precision: > the largest error I found was 2 ulp's when compiled without optimization, > and 1 ulp when > compiled with -O2, both on an i586. > > Emil >

2 1

doubt in Sets
by Anuradha 11 Feb '02

11 Feb '02

Hello, We have a small query regarding sets in GPC. Here is the problem : We have defined a set in a unit and my program is using the this defined set. But the compiler says that the identifier is undeclared. The following is my program and the unit. ========================= > program SETT; > uses > setsunit; > var > c : Myset; > begin > c:=[fs_auto]; > testSET(c); > end. > ============================ > unit setsunit; > interface > type > Myset = set of (fs_auto,fs_trans) ; > procedure testSET(x:Myset) > implementation > procedure testSET(x:myset); > begin > if x=[fs_auto] then > writeln('1'); > end; > end. > ================================= Can you kindly help us on this? Thanks and Regards Anuradha. -- S.Anuradha Generic Data Tools, Alcatel Chennai. Alcanet : 2-757-7123

2 1

Re: Complex_Arctan
by da Silva, Joe 10 Feb '02

10 Feb '02

Please see below ... Joe. > -----Original Message----- > From: Frank Heckenbach [SMTP:frank@g-n-u.de] > Sent: Wednesday, February 06, 2002 11:40 AM > To: jerabek(a)math.cas.cz; gpc(a)gnu.de > Subject: Re: Complex_Arctan > > Emil Jerabek wrote: > > > Frank Heckenbach wrote: > > > > > Emil Jerabek wrote: > > [...] > > > > Your change doesn't solve this, it just reduces the threshold > > > > to EpsReal^2 > 4 / MaxReal. (BTW, the original implementation in > math.pas also needs > > > > a similar condition.) > > > > > > > > Frank, do you _really_ think there are machines with such a small > MaxReal? > > > > > > It doesn't really matter what I think, does it? ;-) Actually, I > > > haven't looked at any non-IEEE FP implementation in detail. > > > > > > I'm applying Joe's patch (since it's certainly a little > > > improvement). If you have more ideas for improvement, let me know. > > > > > > Frank > > > > > > > So, here's another version :-) It should never overflow, AFAICS. > > > > ------------------------------- > > const Ln2 = 0.6931471805599453094172321214581765680755; > > var FunnyNumber: Real; > > > > function Complex_ArcTan (z: Complex): Complex; > > var x, y, a, b, c, d: Real; > > begin > > x := Abs (Re (z)); > > y := Abs (Im (z)); > > b := 1 - y; > > d := (1 + y) / 2; > > if b < 0 then d := -d; > > b := Abs (b); > > if x >= b then begin > > c := b / x; > > d := d * c - x / 2; > > if x <= FunnyNumber then > > b := Ln2 - Ln (4 * x * Sqrt (1 + Sqr (c)) / (1 + y)) / 2 > > else > > b := Log1P (2 * (y / x) / (b * c / 2 + x / 2)) / 4; > > c := 1; > > end else begin > > c := x / b; > > a := x * c / 2; > > d := d - a; > > if b <= FunnyNumber then > > b := Ln2 - Ln (4 * b * Sqrt (1 + Sqr (c)) / (1 + y)) / 2 > > else > > b := Log1P (2 * (y / b) / (b / 2 + a)) / 4 > > end; > > d := ArcTan2 (c, d) / 2; > > if Re (z) < 0 then d := -d; > > if Im (z) < 0 then b := -b; > > Complex_ArcTan := Cmplx (d, b) > > end; > > > > to begin do > > begin > > FunnyNumber := 2.1 / Sqrt (MaxReal); > > end; > > Well, this one is a little longer than the previous ones, so should > I use it? > [Joe da Silva] Probably. It looks more reliable than the previous version. The only concern I still have, is that the test program, emil9, does not test for accuracy with large "arctan" arguments (only up to 1e17). Since "tan" can produce huge results even with just real arguments, an "arctan" function must conversely provide accurate results with huge arguments. BTW, AFAICT, DJGPP does not have a Log1P function, so I wasn't able to verify the above anyway. I didn't know the values to use for the coefficients "Log1PCoefs[]", so I couldn't try the suggested Log1P code, either. Although Emil presented a power series which I assume represents these coefficients, I don't know whether these need to be "tweaked". If I understand this correctly, it is an infinite series. However, AFAIK, when you apply a finite version of such a series, you need to "tweak" it's values to get good accuracy (no, I don't know how to do this). I remember experimenting years ago with the power series for "sine" or "cosine" and being amazed at how a small finite "tweaked" power series was much more accurate than the "conventional" power series, no matter how many terms I calculated for it (after a while, the more terms you calculate, the greater the error you introduce). I presume the "log1p" power series also has this effect(?) ------ snip ------

3 2

Re: Complex_Arctan
by da Silva, Joe 10 Feb '02

10 Feb '02

Please see below ... Joe. > -----Original Message----- > From: Emil Jerabek [SMTP:ejer5183@artax.karlin.mff.cuni.cz] > Sent: Friday, January 25, 2002 5:23 AM > To: gpc(a)gnu.de > Subject: Complex_Arctan > > Hi all, > > I've had a look at the complex arctangent function in the RTS. It turned > out that > > - it raises a runtime error or returns an infinite imaginary part on > arguments too close > to the singularities +/- i, for example > ArcTan (Cmplx (1e6 * MinReal, +/- 1)) > > - it often returns a real number when it should return something with a > small imaginary part > instead (example: ArcTan (Cmplx (1e17, 1e17)) on IA32, YMMV), and in > general, > the imaginary part of the result is very unreliable when |Im (z)| << > |Re (z)|^2 > (because it is essentially computed as `Ln (something close to 1)'; > example: > Im (Arctan (Cmplx (1000, 0.0001))) gives 9.999989726e-11 instead of > 9.999990000e-11, > again this is machine-dependent) > > The following patch to p/rts/math.pas solves these problems (I hope it > doesn't introduce > others). > > ------------------------------ > > function Log1P (x: Real): Real; asmname 'log1p'; > > function Complex_Arctan (z: Complex): Complex; > const Ln2 = 0.6931471805599453094172321214581765680755; > var b, c, d: Real; > begin > c := Sqr (Re (z)); > b := 1 - Abs (Im (z)); > d := ((1 + Abs (Im (z))) * b - c) / 2; > [Joe da Silva] I wonder, would it be better to use the following instead, to reduce problems with overflow for large values (magnitude) of Re(z) and/or Im(z)? : c := Re(z) * 0.5 * Re(z); b := 1 - Abs(Im(z)); d := (1 + Abs(Im(z))) * 0.5 * b - c; Also, as a general question, is division more expensive than multiplication, or does the complier's optimization take care of such matters anyway? > if Abs (b) > Abs (Re (z)) then > b := Log1P (4 * (Abs (Im (z)) / b) / (b + Re (z) * (Re (z) / b))) / 4 > else if c <= MinReal then > b := (Ln2 - Ln (Abs (Re (z)))) / 2 > else > b := Log1P (4 * (Abs (Im (z)) / Re (z)) / (b * (b / Re (z)) + Re (z))) > / 4; > if Im (z) < 0 then b := -b; > Complex_Arctan := Cmplx (ArcTan2 (Re (z), d) / 2, b) > end; > ------------------------ > > > Emil Jerabek >

3 8

How to Cross Compile ?
by gerdunne＠eircom.net 08 Feb '02

08 Feb '02

Hi Can someone please tell me how do I use 'AS' utility to cross compile a Pascal prog. to the MC68000 (S format..) I've read the html file on 'AS' and tried out the -m68000 switch but i must be doing something wrong? I'm using 'GPC. Regards Ger.

2 1

Re: Constant Expressions (was: Complex_Arctan)
by da Silva, Joe 07 Feb '02

07 Feb '02

Yes, the language in both standards is horrible! :-/ AFAICS, all that stuff about the "defining point" for functions is probably just saying that you can use predefined functions only (excluding those that don't have a constant result) and not functions that you define within your program. Joe. > -----Original Message----- > From: Emil Jerabek [SMTP:ejer5183@artax.karlin.mff.cuni.cz] > Sent: Friday, February 08, 2002 4:33 AM > To: gpc(a)gnu.de > Subject: Re: Constant Expressions (was: Complex_Arctan) > > Joe wrote: > > > Please see below ... > > > > Joe. ( Feeling a bit better :-) > > > > Good news :-) > > > > > Well, actually the ISO-10206 standard (since the topic was complex > > maths, that is the applicable standard) *does* allow function calls > > in constant expressions. See sections 6.3.1, 6.3.2, 6.8.2. :-) > > > > There are some restrictions of course, mainly to ensure that the > > constant expression is indeed constant. AFAIK, only pre-defined > > functions are allowed (although I haven't looked too closely at > > that). > > > > Example from section 6.3.2 : > > > > const > > ............... > > pi = 4 * arctan(1); > > ............... > > > > You are probably right. I got somewhat confused by the language > of the standard ;-) > > From section 6.8.2: > > "...An expression shall be designated nonvarying if it does not > contain the following > ... > c) an applied occurrence of an identifier as a function-identifier > that has a defining-point contained by the program-block or that denotes > one of the required functions eof or eoln. > ..." > > Section 6.2.2.1: > > "Each identifier or label contained by the program-block shall have > a defining-point, with the exception of the identifier > of a program-heading (see 6.12)." > > From section 6.2.2.10: > > "Required identifiers that denote the required values, types, schemata, > procedures, and functions shall be used as if their defining-points > have a region enclosing the program (see 6.1.3, 6.4.2.2, 6.4.3.4, > 6.4.3.6, 6.4.3.3.3, 6.7.5, 6.7.6, and 6.10). > ..." > > > So, do those required identifiers have a defining-point or do they not? > > > Emil Jerabek

1 0

Re: Constant Expressions (was: Complex_Arctan)
by da Silva, Joe 07 Feb '02

07 Feb '02

Please see below ... Joe. ( Feeling a bit better :-) > -----Original Message----- > From: Emil Jerabek [SMTP:ejer5183@artax.karlin.mff.cuni.cz] > Sent: Wednesday, February 06, 2002 4:39 AM > To: gpc(a)gnu.de > Subject: Re: Complex_Arctan > > Frank Heckenbach wrote: > > > Emil Jerabek wrote: > [...] > > > Your change doesn't solve this, it just reduces the threshold > > > to EpsReal^2 > 4 / MaxReal. (BTW, the original implementation in > math.pas also needs > > > a similar condition.) > > > > > > Frank, do you _really_ think there are machines with such a small > MaxReal? > > > > It doesn't really matter what I think, does it? ;-) Actually, I > > haven't looked at any non-IEEE FP implementation in detail. > > > > I'm applying Joe's patch (since it's certainly a little > > improvement). If you have more ideas for improvement, let me know. > > > > Frank > > > > So, here's another version :-) It should never overflow, AFAICS. > > ------------------------------- > const Ln2 = 0.6931471805599453094172321214581765680755; > var FunnyNumber: Real; > > function Complex_ArcTan (z: Complex): Complex; > var x, y, a, b, c, d: Real; > begin > x := Abs (Re (z)); > y := Abs (Im (z)); > b := 1 - y; > d := (1 + y) / 2; > if b < 0 then d := -d; > b := Abs (b); > if x >= b then begin > c := b / x; > d := d * c - x / 2; > if x <= FunnyNumber then > b := Ln2 - Ln (4 * x * Sqrt (1 + Sqr (c)) / (1 + y)) / 2 > else > b := Log1P (2 * (y / x) / (b * c / 2 + x / 2)) / 4; > c := 1; > end else begin > c := x / b; > a := x * c / 2; > d := d - a; > if b <= FunnyNumber then > b := Ln2 - Ln (4 * b * Sqrt (1 + Sqr (c)) / (1 + y)) / 2 > else > b := Log1P (2 * (y / b) / (b / 2 + a)) / 4 > end; > d := ArcTan2 (c, d) / 2; > if Re (z) < 0 then d := -d; > if Im (z) < 0 then b := -b; > Complex_ArcTan := Cmplx (d, b) > end; > > to begin do > begin > FunnyNumber := 2.1 / Sqrt (MaxReal); > end; > ------------------------------------ > > > I encountered a compiler bug during testing :-( > Strictly speaking, the Standard doesn't allow function calls in constant > expressions, > but I expect the compiler would reject it if it were not supposed to work. > [Joe da Silva] Well, actually the ISO-10206 standard (since the topic was complex maths, that is the applicable standard) *does* allow function calls in constant expressions. See sections 6.3.1, 6.3.2, 6.8.2. :-) There are some restrictions of course, mainly to ensure that the constant expression is indeed constant. AFAIK, only pre-defined functions are allowed (although I haven't looked too closely at that). Example from section 6.3.2 : const ............... pi = 4 * arctan(1); ............... > [pas]% cat exprconst.pas > program Foo (Output); > > const > Bar = 1 / Sqrt (16); > ------ snip ------

2 1

business marketing
by jjmkt2002160715＠yahoo.com 06 Feb '02

06 Feb '02

1 0

pre-registration instructions (5598IZYG0-776@12)
by dot-biz-registrar586805＠mail.com 06 Feb '02

06 Feb '02

PUBLIC ANNOUNCEMENT: The three new Internet domain extension formats - .NAME .BIZ .INFO . have now been introduced to the Internet. However, many users are still unfamiliar with how these names work. The following are some tips on what these names are and where you can register. ############### .NAME Extension ############### - Currently in pre-registration stage. - Designed for personal use. - 2-level registration format (eg. FirstName.LastName.name) - Used for web URL (www.FirstName.Lastname.name) or email (FirstName(a)LastName.name) - Registry accepts submissions once a week, first come first serve basis. - Where to pre-register: http://www.NamesForEveryone.com. ############## .BIZ Extension ############## - Currently operational. - Designed for business use. - No registration restrictions. - Works just like a .COM name - Where to register: http://www.DotBizToday.com ############### .INFO Extension ############### - Currently operational. - Designed for information on products / organizations - No registration restrictions - Works just like a .COM name - Where to register: http://www.InfoGrab.com Over a million names of these new extensions have already been registered. Therefore, if your name is still available, please secure it immediately to avoid possible disputes. Yours truly, Domain Administrator [6760OPhl4-370lVTF8465QSdy0-002rIlB3365Cpyv7-964F@45]

1 0

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

2007

2006

2005

2004

2003

2002

2001

2000

1999

1998

1997

1996

1995

Gpc February 2002