User talk:Stanislav Sykora

WriteRealForOEIS

Hi, I just saw your utilities and thought I'd comment. I like the WriteRealForOEIS function. It's somewhat similar to my bfile function which I use daily. One thing that I learned when writing bfile is that gp is very bad at writing files: in particular, each write command opens the file, writes a single entry, then closes the file again. So writing a 10,000-term bfile with 10,000 writes takes a long time -- well, maybe 5-10 seconds, but a lot longer than it would take otherwise.

If this bothers you, the workaround I discovered was to concatenate a number of entries together so that I could write in large chunks (say, 1000 terms). This improved speed greatly. If your read/write performance is good enough, or you're sufficiently patient, you don't need this of course.

Also I thought it might be a good idea to trim off the last digit (which is likely wrong) plus all trailing 9s (together with the preceding digit of course) and trailing 0s which are also suspect. So if the value was computed to be ...3001 but the true value is ...2998 (3 ulps, not an unreasonable error for a complex calculation) then the last three terms are spurious. This way you wouldn't be tempted to leave bad digits in.

Charles R Greathouse IV 01:31, 7 November 2013 (UTC)

Hi, yes, it can be improved, of course. So far, the few seconds delay (at most) in writing the files does not bother me. As for the trailing digits, since the file I generate contains all needed data for OEIS submission and therefore needs to be manipulated manually, I thought that the concerns about precision might be left to the User to handle (e.g., I usually specify a precison of e.g. 2020 and then leave in only 2000 digits. But yes, once I see that any of these aspects should/could be automated, I will do it.

TetrationLimit

On your other file, you might find the new function lambertw to be useful here, if you're running version 2.6.1 or later (and why not?). For the convergent case e^-e < x < e^(1/e) the answer is given by

my(t=-log(x)); lambertw(t)/t

which is fast and precise.

Charles R Greathouse IV 02:01, 7 November 2013 (UTC)

Agreed; I have in fact written the same thing in the comments to my script (see the file). To tell you the truth, since lambertw(x) appeared only in PARI 2.6.1, I wanted to maintain the possibility - for myself and for the User - to use both ways and compare (maybe I sin on the side of not mistrusting new code). So TetrationLimit always does the old, sloppy work, ignoring efficiency but making sure the result is correct. It is not a replacement for Lambert W. :-)

Sounds good. I just thought I'd bring it up in case you hadn't heard about the function. (I saw that it wasn't even mentioned in the CHANGES file, so I emailed the maintainers and they added it.)

One backward-compatibility dream I have is to make a library of functions which duplicate new functions, so that a person could (say) grab a whole compatibility file or just a few functions and use programs designed for a more recent version. If the starting version is new enough you could even use lex() and version() to load only the functions needed. But then I realize it's easier just to ask people to upgrade. :)

Charles R Greathouse IV 01:14, 8 November 2013 (UTC)

Prime-related functions

I hope you're not sick of my comments yet!

You may find this useful, especially for simple functions: the final computed result is returned, so you do't need to use return() on the last line. So instead of

IsNotPrime(n) = {return (!isprime(n));}

you can use

IsNotPrime(n) = !isprime(n)

Note that NumberOfPrimeFactors and NumberOfDistinctPrimeFactors are already in GP: bigomega and omega, respectively. Similarly IsMoebiusNonZero(n) already exists: issquarefree.

Charles R Greathouse IV 00:42, 27 November 2013 (UTC)

The more you comment, the happier I am! Old, but still learning.

The duplication of intrinsic GP functions is a silly oversight. The "return" style classifies me as a C++ programmer, which I am. Since OEIS priviledges compact coding, I will do so also in the simpler functions. Will update all today.

Stanislav Sykora 09:18, 27 November 2013 (UTC)

I've duplicated at least as many GP functions by mistake, not counting the times I wrote a function which was later added. I figure anything I can do to help others avoid that, the better.

PARI is written as a C library, so the PARI/GP maintainers are (like myself) very comfortable with C. It's usually a good thing though it does make it a bit harder to adjust to the functional side of GP (e.g., apply, select, and the new | <- , operator).

Charles R Greathouse IV 17:44, 28 November 2013 (UTC)

PreviousPrime / FollowingPrime

I recommend

PreviousPrime(n)=precprime(ceil(n)-1);

FollowingPrime(n)=nextprime(floor(n)+1);

which 'do the right thing' if you happen to pass a noninteger value to them.

Charles R Greathouse IV 04:30, 15 December 2013 (UTC)

Thanks, the ceil and floor are a good idea, I will do it today. As for the names, I use on purpose the style with capitals in order not to avoid any overlap with the intrinsic PARI namespace (no need to step on other peoples' heels). The 'nextprime' function, for example, exists already, but it does >= instead of > (when the argument is a prime p, it returns p). Hence, I will use PrecPrime and SuccPrime or something like that.

May I ask you a question? As a C++ programmer, and despite eternal lack of time, my fingers are itching for a PARI/GP IDE. It is really a pity one can not do Copy and Paste, for example. I believe there was an attempt in the past, but was abandoned (?) Do you know about anybody else working on this or willing to go into it?

Stanislav Sykora 09:41, 15 December 2013 (UTC)

I use text editors with syntax highlighting, like so:

http://math.crg4.com/pari/highlight.html

You can copy and paste from gp, but if you're in Windows you need to enable quick edit mode. Right-click the gp shortcut, go to the Options tab, and click "QuickEdit Mode". (See here for a screenshot.) This lets you copy by highlighting text with the mouse (and pressing Enter?) and paste by right-clicking.

I use Linux where I can copy as described and paste with shift-control-V.

I'm not sure what other features you'd want in a PARI/GP IDE, but then again I don't use IDE's much. (Sometimes I use Geany or VisualStudio, but usually just gedit on Linux or Notepad++.) Possibly easier than writing a whole IDE would be to find some open-source IDE and add support for GP and PARI.

Charles R Greathouse IV 16:22, 15 December 2013 (UTC)

More on your functions

More comments on File:PrimesRelatedFunctions.txt and File:DigitBasedFunctions.txt since your last additions:

• LargestPrimeFactor(n) is repeated twice. Also you should use my() instead of local() essentially everywhere:

LargestPrimeFactor(n) =
{
  if (n==1,return(n));
  my(v = factor(n));
  v[#v[,1],1]
};

• SumOfPrimeFactors(n) can be simplified similarly:

SumOfPrimeFactors(n) =
{
  my(v=factor(n));
  sum(k=1,#v[,1], v[k,1]*v[k,2])
};

• The same simplification applies to ProductOfPrimeFactorsN(n).
• For Lambda(n,m) you should (as always) use my(), which also allows you to give the initial value of the variable. Further, you can use the /= operator:

Lambda(n,m) =
{
  my(result=binomial(n,m));
  forprime(p=1+m,n, result /= p);
  result
};

• DigitsToNumber(d,b) can be simplified with subst():

DigitsToNumber(d,b=10) = subst(Pol(d), 'x, b);

• A similar trick works with ReverseDigits(n,b):

ReverseDigits(n,b=10) = subst(Polrev(digits(n,b)), 'x, b);

• I don't know why you use VectorsDifferAt for IsPalindrome(n,b) when == works just as well but faster. I would do IsPalindrome(n,b=10) = n==ReverseDigits(n,b);.
• NumberOfDigits(n,b=10) = #digits(n,b);
• sumdigits() supersedes your SumOfDigits, although ProductOfDigits is still needed.
• Parity(n) = hammingweight(n)%2

Hopefully these comments are useful!

Charles R Greathouse IV 02:58, 16 May 2014 (UTC)

Charles, thanks so much for your editing help! It is time for me to go through a systematic editing cycle and merge all your comments (there are also some older ones still pending) as well as several items of mine. You basically perceive that, for lack of time and other interests, I have not fully assimilated PARI/GP and always stopped at the minimum depth that allowed me to solve a specific problem of the moment. However, since this is in public view, that approach is not quite appropriate. I will do the effort next weekend.

--Stanislav Sykora 07:14, 16 May 2014 (UTC)

With delay, I did most of the edits you recommended. There is more pruning of the scripts to be done; work in progress. In a couple of cases, I refrained from going on because of doubts about PARI (still learning). Could you please try these sequence of commands (on separate lines): precprime(2) | precprime(1) | precprime(1) | ... Why does precprime(1) return various integers; should the behaviour not be fixed, even if the argument is illogical?

--Stanislav Sykora 14:34, 6 June 2014 (UTC)

Digit-based functions

A few comments on File:DigitBasedFunctions.txt, if you don't mind.

• DigitsToNumber now duplicates the (new) function fromdigits. Nothing wrong here, just thought I'd point it out. Presumably the built-in function is faster if it matters.
• IsReversedGE could more easily be written IsReversedGE(n,b=10) = ReverseDigits(n,b) >= n; -- no need for the if.
• For SumOfDigits and ProductOfDigits there are two issues. First, you've actually created two lexically-scoped variables, the k outside the loop (never used) and the k inside the loop. To avoid confusion you shouldn't have k under my. Second, there's no need to use return here; the final line is returned, so prod(k=1,#d,d[k]) would be preferable.
• GrayDecode: Like above, there's no need for return on the last line (the developers discourage this).

Charles R Greathouse IV 03:35, 16 October 2014 (UTC)

Thanks. I will keep an eye on it. In version 2.7.1 that I dowloaded last week, fromdigits does not seem to exist. I have done the other modifications.

Since you are PARI/GP expert, may I ask you a question? So far I took it for granted that all the intrinsic functions (like, e.g., complex erfc) work up to default(realprecision), minus perhaps half a dozen digits. But is this really so? Not being quite sure, I always test extensively by, e.g., doubling precision, comparing, and then keeping just half of the digits (so far, such tests were very comforting). It seems to me that in the newest version the displayed digits are handled more dynamically, so possibly there were some second thoughts, but I can not find any explicit general guidelines.

--Stanislav Sykora 09:28, 16 October 2014 (UTC)

Yes, fromdigits() is very new.

The goal, I think, is that any single function call should be accurate to within 1 ulp. (The best you could hope for, in general, is 0.5ulp -- you don't have the precision to do better.) The complicated functions have error bounds built-in, and if they can't be confident of final digits they will truncate the result by words. So you might get an answer with just 982 digits under \p 1000, because gp didn't trust the last word. In extreme cases it might trim multiple words.

So even when you're doing something complicated, you should be good up to the last word. If gp thinks that you lost 10 digits of precision at 1000 digits, it won't truncate because you'd lose another 9 digits from the truncation.

To be more conservative you could leave on two guard words (~38 digits) for high-precision work. Doubling is pretty extreme, I wouldn't go that far (except, perhaps, in testing new features).

Of course once you get into writing your own numerical code in loops, all bets are off. If it's poorly designed you could lose a little bit of precision at each iteration, ending up with garbage. There's no easy way out here -- you'd need to know numerical analysis...!

By the way: if you find a function that returns an inaccurate result, it's very possibly a bug and should be reported.

Charles R Greathouse IV 15:29, 26 October 2014 (UTC)