Sheplyakov Alexei wrote:

+class lanczos_coeffs +{ + public: + lanczos_coeffs(); + bool sufficiently_accurate(int digits); + int get_order() const { return current_vector->size(); } + numeric calc_lanczos_A(const numeric &) const; + private: + static std::vector<numeric> coeffs_12; // Use in case Digits <= 20 + static std::vector<numeric> coeffs_30; // Use in case Digits <= 50 + static std::vector<numeric> coeffs_60; // Use in case Digits <= 100 + static std::vector<numeric> coeffs_120; // Use in case Digits <= 200

With such a code one need to break ABI (add extra class members) in order to increase the precision. May be replacing all these with static std::vector<std::vector<numeric> > would be better solution?

But if the extra class members are of a class that is not intended to be used by any user of the library and cannot be used (at least without messing around), then it is okay. Chris was careful enough to add the class to file numeric.cpp, not to any of our header files. So it is not something anybody can use when linking his/her apps. Hence, there isn't anything that can break when the library is later upgraded underneath applications using it. Regards -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>

Dear Alexei, On Wed, 13 Dec 2006, Sheplyakov Alexei wrote:

// bugme.cpp #include <ginac/ginac.h> namespace GiNaC { class lanczos_coeffs; }

void foo() { using namespace GiNaC; lanczos_coeffs bugme; // do something here }

The compiler says: "error: aggregate GiNaC::lanczos_coeffs bugme has incomplete type and cannot be defined". Hence, users would have to import the entire definition of this class from the file numeric.cpp. I think in that case they should know that they are doing something non-standard that can break at any time. If this is the only problem with it, I would not worry too much. Best wishes, Chris

Hi, Chris!

...

const numeric lgamma(const numeric &x) { - throw dunno(); + lanczos_coeffs lc; + if (lc.sufficiently_accurate(Digits)) { + numeric pi_val(cln::pi(cln::default_float_format));

What happens here if Digits > default_float_format?

I sincerely hope that it won't cause problems

... because it is impossible: /** Assign a native long to global Digits object. */ _numeric_digits& _numeric_digits::operator=(long prec) { digits = prec; cln::default_float_format = cln::float_format(prec); return *this; } Sorry for the noise.

...

All these foo.mul(bar).add(baz) are plain ugly. Any objections against s/numeric/cl_N/g (so it is possible to use natural syntax)?

If that has infix operators, I would say it would be much nicer.

Sure, it has (like [almost] all CLN types).

...

With such a code one need to break ABI (add extra class members) in order to increase the precision. May be replacing all these with static std::vector<std::vector<numeric> > would be better solution?

Richy has already answered this.

I think his answer is wrong. Leaving aside the ABI, I don't think adding more members to a class to achive better precision is a good idea anyway. Secondly, it would be nice to remove artificial limit (Digits <= 200) and make the evaluation really arbitrary-precision. To do this one need to calculate missing coefficients at run-time and store them in vector<vector <numeric> > (or vector<vector <cl_R> >).

...

+lanczos_coeffs::lanczos_coeffs() +{ if (coeffs_12[0] != 0) + return;

I think coeffs_12[0] might be not initialized at this stage (and contain random garbage).

No, that is not true. The elements of the vector have been constructed using the default constructor of numeric.

I think it is safer and simpler to add `static int count' so we don't ever run into yet another static order initailization problem. Best regards, Alexei -- All science is either physics or stamp collecting.

Hello, Jens!

@Alexei: I had to apply two of your patches manually, because the format was somehow broken. I hope I did everything correctly.

Yes, modulo minor whitespace damage :)

...

I'm little bit busy now, could you please wait till tomorrow?

okay. I will wait (as long as it takes).

I've found yet another bug (seems to be fixed in HEAD): #include <iostream> #include <stdexcept> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; int main(int argc, char** argv) { idx a(symbol("a"), 8), b(symbol("b"), 8), c(symbol("c"), 8), d(symbol("d"), 8), e(symbol("e"), 8), f(symbol("f"), 8); ex H1 = color_f(a, c, e)*color_f(b, d, e) + color_f(a, d, e)*color_f(b, c, e); ex test = (H1*H1).expand(expand_options::expand_indexed); cout << test << endl; // 2*f.e.b.d*f.e.b.c*f.e.a.c*f.e.a.d+f.symbol8.b.c*f.symbol8.a.d*f.e.b.c*f.e.a.d+f.symbol7.b.d*f.symbol7.a.c*f.e.b.d*f.e.a.c // inconsistent indices! cout << test.simplify_indexed() << endl; // prints 288, should be 216 return 0; } I wonder if it is possible to backport the fix... Best regards, Alexei -- All science is either physics or stamp collecting.

Dear Alexei, On Wed, 13 Dec 2006, Sheplyakov Alexei wrote:

ex H1 = color_f(a, c, e)*color_f(b, d, e) + color_f(a, d, e)*color_f(b, c, e);

ex test = (H1*H1).expand(expand_options::expand_indexed); cout << test << endl; // 2*f.e.b.d*f.e.b.c*f.e.a.c*f.e.a.d+f.symbol8.b.c*f.symbol8.a.d*f.e.b.c*f.e.a.d+f.symbol7.b.d*f.symbol7.a.c*f.e.b.d*f.e.a.c // inconsistent indices!

cout << test.simplify_indexed() << endl; // prints 288, should be 216

At my place this prints 144-2*f.c.a.symbol8*f.a.d.e*f.d.symbol8.b*f.c.e.b. Strange??? @dummy index renaming: yes, this has been much improved in HEAD. In principle I think it would be possible to move the class make_flat_inserter into the *.cpp files. But that would be code duplication. Not so nice. Besides that, new functions were added to indexed.cpp and these need to be called from indexed.h. Fixing all this in a binary compatible way would get way too messy IMO. Maybe we should release 1.4 and claim proudly that we now can rename dummy indices. Good thing that we don't keep an accurate NEWS file. That way we can can announce it like it is new. :-o ;-) Best wishes, Chris