Bug and updates to calc

• 3rd September 2000. Fixed an error on line 98 of lagrange.c, where I had Q0_2 instead of Q0. The error showed up in patz(67,89).
• 7th August. Changed log1(a,b,r,e) to log1(a,b,d,r,e), so as to give the user the chance to input any positive base d > 1. Updated calc.tar.gz and calc.exe.gz.
• 13th July 2000. Replaces "4" by "5" on line 72 of init.c. Updated calc.tar.gz and calc.exe.gz.
• 4th July 2000. While checking the numerical examples at the end of On the Diophantine equation u²-Dv²=±4N, Bengt Stolt, Arkiv för Matematik, 2 1-23, 1952, I uncovered three simple but important bugs in CALC: (i) in primes1.c: line 1651, *lptr=1; changed to *lptr=2;, (ii) in lagrange.c: line 56, formerly I had freed M0, then used it on line 83! Moreover I should have been using M instead of M0 in TMP2 = MULT_I(M0, j); in line 83! Updated calc.tar.gz and calc.exe.gz.
• 2nd June 2000. Fixed a bug in sqroot, which arose in the example sqroot(79,42,&[],&m,&l). (It returned 8 instead of 4). The fix was to insert an "if" statement at line 1825 of primes1.c. Could have corrupted patz(d,n). Updated calc.tar.gz.
• 29th May 2000. John Robertson raised doubts concerning about my claim that a fundamental solution would be found by examining only up to the first period in the continued fractions for the quadratic irrationalities used in patz(d,n) in the case where the period length is odd. I now use the second period in this case. The changes were minimal: an extra parameter surd_flag was added to SURD(), so that an extra period is done if and only if surd_flag = 1. Changes were also made to SURDX() and PATZ(). Updated calc.tar.gz and calc.exe.gz.
• 25th May 2000. Removed two errors in two print statements, improved printing of output and added printing to a file patz.out. Updated calc.tar.gz and calc.exe.gz.
• 24th May 2000. Removed a bug on line 93 of lagrange.c, replacing if ( tmp) by if ((Q0_2)->V[0]) % 2 == 0 && tmp). The bug had caused patz(5,19) to declare (12,5) to be a fundamental solution, whereas the correct answer is (8,3). Also removed the part of the code made redundant by my new proof of necessity in patz.pdf, which rendered unnecessary the need to look at the cases D=2 and 3, N < 0 separately. Updated calc.tar.gz and later in the evening calc.exe.gz.
• 22nd May 2000. Tightened up wrappers for patz and cornacchia. Updated calc.tar.gz and calc.exe.gz.
• 17th May 2000. Fixed a harmless memory leak of 20 bytes in CONVERGENTS, where in cases n==0 and n==1, I forgot to FREEMPI(ONE). It propagated through SURD and thence to PATZ. (This took me over 3 hours to solve!) Updated calc.tar.gz and calc.exe.gz.
• 16th May 2000. Fixed a a bug in patz(d,n) which occurred when N=1 (pell's equation). Only the fundamental solution with Norm(eta)=-1 had been printed. Updated calc.tar.gz and calc.exe.gz.
• 15th May 2000. Fixed a bug in sqroot. (The global variable sqroot2_number in SQROOT2() was not always reinitialised to a default value zero after use.
  Installed patz. This returns the positive fundamental solutions of the diophantine equation x²-d*y²=±1, d > 1, not a perfect square. The algorithm is from paper of the author and is in essence due to Lagrange 1770. Updated calc.tar.gz and calc.exe.gz.
• 6th May 2000. Fixed a bug in surd. Updated calc.tar.gz and calc.exe.gz.
• 4th May 2000. Improved surd. Also took some wrapper functiondeclarations out of fun.h and put them in wrappers.h. Changed SURD, SURDX, SURDW_W and REDUCED. Also put global MPIA variables A_SURD, U_SURD, V_SURD in nfunc.c. Updated calc.tar.gz and calc.exe.gz.
• 18th April 2000. Added cornacchia. Also fixed a bug in sqroot (on line 1532 of primes1.c e should have been e+1). Updated calc.tar.gz. Also updated calc.exe.gz.
• 15th April 2000. Updated calc.exe.gz.
• 14th April 2000. Made sqroot(a,n,&s[ ],&m) available. Updated calc.tar.gz. Also Sean Seefried fixed a bug regarding the NULL pointer (perfectpower was still misbehaving). COPYI() now copes with this and a couple of related changes to parse.y were made.
• 27th March 2000. Made log1(a,b,r,e) available. This is Algorithm 1 of (manuscript-pdf 152K). Updated calc.tar.gz and calc.exe.gz.
• 25th March 2000. Improved the algorithm for log(a,b,d,u,v,e). This is Algorithm 3 of (manuscript pdf-152K). Updated calc.tar.gz.
• 20th March 2000. Changed the algorithm for log(), as I realised the earlier algorithm, admittedly heuristic in nature, was not giving the right answers for log(2¹⁰+1,2,2,1). Updated calc.tar.gz and calc.exe.gz. (See algorithm on page 7 of On a version of Shanks' algorithm for computing the continued fraction of log_ba (pdf 120K).)
• 13th March 2000. Added log(). Updated calc.tar.gz and calc.exe.gz.
• 9th March 2000. Updated calc.exe.gz.
• 8th March 2000. Slight alteration to rootexp() so as to allow the user to enter a not necessarily squarefree polynomial. (The polynomial f(X) entered is replaced by f(X)/gcd(f(X),f'(X)).) calc.tar.gz updated.
• 18th February 2000. Sean Seefried fixed a bug in the error recovery for some functions when not all arguments were entered. This caused a subsequent crash when the function was subsequently called with the correct number of arguments. Files parse.y (and parse.c), calc.h and symbol.c were changed. calc.tar.gz and calc.exe.gz were updated.
• 8th February 2000. Fixed a bug in nprimeap which arose in the recent updating of CALC: on line 387 of wrappers.c, FREEMPI(B) should have been FREEM PI(M). I also noticed a missing "return(NULL);" after line 749 of func.c. calc.e xe.gz was also updated.
• 26th January 2000. Updated the readme in calc.exe.gz. (The description of euclid contained some inaccuracies.)
• 23rd January 2000. I updated calc.exe.gz. In so doing I realised that in the 20th January update, I had forgotten to include some files! Also Makefile_dos was incorrect. I believe that this is now fixed in the current version of calc.tar.gz. Also it is possible some bugs have slipped in during Sean's recent improvements to CALC. Please let me know if any are found.
• 20th January 2000. Sean has finished his work on CALC. calc.tar.gz has been updated. This now includes a DOS makefile based on the djgpp C compiler) and a file ytab.h=y.tab.h. calc.exe.gz will be updated soon.
• 14th January 2000. Progress report: Since mid-December 1999, vacation scholar Sean Seefried has been working on CALC to improve error recovery. This has now essentially been finished. Sean has also written a Sturm's algorithm and in a week's time will have used this to improve lagrange so as to find the simple continued fraction of all roots of a squarefree polynomial with integer coefficients with no rational roots, using the initial part of a paper of D.G. Cantor, P.H. Galyean and H.G. Zimmer.
• 26th November 1999. Some cosmetic improvements to output were made. In surd.out, the periodic part is now separated from the initial segment by a colon. Also the quadratic irrational is listed at the top of the file.
  The function convergents now outputs the arrays p[ ] and q[ ] to a file convergents.out.
  In parse.y, line 398, <= 0 was changed to < 0.
  
• 28th October 1999. Removed the dependence on external files title.dat and readme, which relied on the more command. title.dat has been incorporated into trial.c, while readme has been turned into readme.c, accessible inside CALC by typing help. A list of all CALC functions comes up and one can make a selection of a function name.
  Also now when a syntax error happens, CALC now recovers better and appears no longer to hang after a serious syntax error. All that was done was to insert a call to execerror(s, ""); in yyerror(s) in parse.c. This returns the user to the command line prompt. However one is still advised to exit CALC if a syntax error occurs.
  The MSDOS version calc.exe was also updated.
  Also files title.dat and readme are now incorporated in the executable.
• 19th October 1999. Improved euclid. It now returns the quotients, remainders and also the s_k and t_k, as well as the number of steps. Files i5I.c, parse.y, parse.c and fun.h were altered.
  MSDOS version was also updated.
• 28th September 1999. Added leastqnr(p) to CALC. This returns n_p, the least quadratic non-residue (mod p), if n_p< 65536, otherwise returns 0. MSDOS version also updated.
• 2nd September 1999. Fixed a small bug in surd() at line 742 of nfunc.c. (I had freed W a few lines too soon!) MSDOS version not yet fixed.
• 29th June 1999. Added the choice of finding all shortest solutions to axb(). Updated calc.tar.gz.
• 17th May 1999. Added power_cubicm() ( the kth power of (X,Y) on an elliptic curve over Z_p), order_cubicr() and order_cubicm(), which find the order of a point P on an elliptic curve over the rationals and Z_p, respectively. The latter is not sophisticated and simply tests each iterate kP for k=2,... calc.exe has still not been updated since 4th Feb. Updated the MSDOS version calc.exe.gz.
• 21st March 1999. Removed a bug at the end of collatz.c. I had freed the arrays m[] by carelessly first freeing the pointer to m and then the actual m[i]!. Similarly for X[]. (The bug was passed by solaris, but linux was unforgiving and dumped a core.)
• 15th March 1999. Updated calc.tar.gz (but not calc.exe). Added cycle() a cycle finding program for the generalised collatz function. (see survey (pdf)).
• 4th February 1999. Updated calc.tar.gz and calc.exe.gz. Added some hidden stuff on my gcd research for use by fellow researchers, but of no interest to the average calc user.
• 20th January 1999. Updated calc.tar.gz (but not yet calc.exe.gz). Added a shortest lattice vector program slv(). Also updated calc.exe.gz.
• 19th January 1999. Updated calc.tar.gz (but not yet calc.exe.gz) - cosmetic changes to the output of fp() and inhomfp(). Subsequently also updated calc.exe.gz.
• 18th January 1999. Updated calc.tar.gz and calc.exe.gz. I improved functions egcd(), lllgcd() and lllgcd0(), making it easier to get a complete list of all shortest multipliers. I also added an option to factor(n) whereby one can enter the number of elliptic curves used. I removed the function short() as it is no longer needed.
• 18th November 1998. Updated calc.tar.gz. Fixed a bug in i5m.c in FFM, line 763. The elements c[i] may be >= 2¹⁶, so SHIFTLB() cannot be used. In any case I don't use FFM as my implementation seems slow. Anyone who wishes to use it should uncomment the relevant code in MULTI() in i1.c. Also I am unable to free the temporary arrays B and C in FFT because of the way the recursion acts; consequently there is a buildup of nettbytes. Any suggestions on this score are welcome.
  
• 4th January 1999. Updated calc.tar.gz, so as to allow the user to enter arbitrary m1 and n1, as mentioned previously. Also updated calc.exe.gz.
• 30th December 1998. Realized I'd forgotten to warn users of LLL-based programs that m1 and n1 are less than 2¹⁶. I will render this warning unnecessary in the near future. My misuse of this in lllgcd() provided a counter-example in paper HMM, with m1/n1=100001/400000, which was in fact not a counter-example after all. (100001>655536.)
• 27th November 1998. Updated calc.exe.gz. The above bug seems to disappear when I compile with the -O2 option instead of -O3. Also updated calc.exe.gz.
• 21st November 1998. Updated calc.exe.gz. On factorizing (10⁶⁵-1)/9, calc hung in mpqsieve.c. This turned out to be due to not observing that 100000 is not less than 2¹⁶. This caused a problem on line 116 with INT0_, which is solved by introducing srange=CHANGE(SRANGE) and then using INT0 instead of INT0_. There is one bug which continues to elude me. When I factorize numbers of greater than 50 digits invoking MPQS1, while the factoring processs runs to completion, one cannot always quit gracefully and also cannot continue using calc. Thus one has to kill calc on a UNIX machine, for example. This is unsatisfactory. The problem is clearly to do with freeing and I had hoped the error lay in the 1000000 bug discovered above. But apparently not. The code is so complicated and the bug only manifests itself after 3 hours running on an Ultra Sparc machine that it seems impossible at the moment for me to detect. If any enthusiast out there has any ideas, please let me know. One example where I think it hung was in factorising the 65 digit prime formed from 64 1's with a 4 in the middle.
• 20th August 1998. Updated calc.tar.gz, removing the changes just made. Did make changes to nfunc.c, i9.c, fun.h, as regards SMITHI and SMITHI1.
• 15th August 1998. Updated calc.tar.gz, because of introduction of verbose mode in SMITH() in nfunc.c and i9.c. (Changes made to fun.h well.)
• 9th August 1998. Updated calc.tar.gz, fixing a small bug in lllhermi.c, where lines 38-39 should have been part of the body of the if(!keith97). Also updated calc.exe.gz.
• 6th August 1998. Updated calc.tar.gz. Fixed a small bug in axb(), whereby -X was returned instead of X at times.
• 5th August 1998. Updated calc.tar.gz. Added a new function axb(), which solves the linear system AX=B, where A,X,B have integer coefficients.
• 19th June 1998. Updated calc.exe.gz.
• 10th June 1998. Updated calc.tar: Deleted lines 468-469 of elliptic.c, (ie, FREEMPI(P); and FREEMPI(tmp); to avoid double freeing. Also changed j > 100 to j > 5 at line 374. Also made a slight improvement to BRENT_POLLARD as regards increasing a. calc.exe will be similarly updated soon.
• 13th May 1998. Updated calc.exe.gz.
• 9th May 1998. Updated calc.Z and calc.tar.gz. The only change being in mpqsieve.c, line 75, where I now use FBASE = 3000; SRANGE = 180000; TOLERANCE = 22; This enabled me to factorise 10¹⁰³+1 overnight, the stumbling block hitherto being a 61 digit number.
  10¹⁰³+1= 11*1237*44092859*984385009*102860539*612053256358933*
  182725114866521155647161*1471865453993855302660887614137521979
• 3rd March 1998. Updated calc.exe.gz in case my absentmindedness had extended to calc.exe.
• 26th February 1998. I had absentmindedly not updated all the relevant source files on 10th February. This has been fixed and calc.tar.gz updated.
• 11th February 1998. Updated calc.exe.gz.
• 10th February 1998. Rewrote POLLARD(). Removed segmentation error and hanging problems arising from returning (MPI *)NULL from various functions such as POLLARD, by altering parse.y and symbol.c at a few places, not using FREEMPI(N) and COPYI(N) when N = NULL.
  pollard() and perfectpower added to list of defined functions.
  Updated calc.Z and calc.tar.gz.
• 5th February 1998. Removed the array R[1280] from inside EFACTOR() and placed it in a separate file primepowers.h. This now makes elliptic.o compile much quicker.
  Also created a new function elliptic(n,m,p) to factor n using the elliptic curve method.
  Updated calc.Z, readme, krm_calc.html and calc_doc.html.
  
• 2nd February 1998. Updated calc.exe.gz.
• 1st February 1998. Updated calc.Z and calc.tar.gz. Have yet to update calc.exe.gz.
• 29th January 1998. Alex Balfour pointed out that factoring 10¹⁰⁰+1 gave rise to rubbish. The cause turned out to be five global arrays Q_[50], Q[50]; Q_PRIME[50]; K_[50], K[50] in primes1.c.
  Replacing 50 by 100 solves the problem. This bug will be fixed soon. I have also increased the Brent-Pollard step number from 2048 to 8192.
• 19th December 1997. I changed lagrange(a,n,m) to lagrange(a,&aa[],n,m). Also changed convergents(a,n,&p,&q) to convergents(a,&p[],&q[],n). Also changed the meaning of n in euclid(a,b,&c[],&n), so as to be consistent with its usage in convergents and lagrange.
• 18th December 1997. Also at the request of Edward Brisse, I added lagrange(a,n,m).
• 15th December 1997. At the request of Edward Brisse, I added euclid(a,b,&c[],&n) and convergents(a,n,&p,&q).
• 4th December 1997. Changed collatz(x) to collatz(x,e) and pell(d,&x,&y) to pell(d,e,&x,&y), so that iterates and partial quotients (respectively) are printed iff e is non zero.
• 2nd December 1997. Oops - I somehow did not comment out the call to ffm inthe MSDOS version of calc, as stated on 30th July 1997. I have just rectified this and have updated calc.exe.gz.
• 10th October 1997. Updated the MSDOS executable.
• 9th October 1997. Oops! I noticed I had absentmindedly left some print statements in lllgcd.c in LLLGCD() lines 89-123 which only pertained to the gcd of 3 integers, causing a problem in verbose mode for other than 3 integers. Fixed the source and the SUN executable.
• 25th September 1997. Upgraded the MSDOS executable.
• 24th September 1997. Upgraded the SUN executable and file primes1.c by adding PERFECT_POWER(N). This is employed in BIG_FACTOR(N) to return X if N=X^p. A corresponding declaration was added to fun.h. (Many thanks to Lew Baxter for pointing out that calc was not factoring a perfect power.)
• 30th July 1997. Upgraded the SUN and MSDOS executables on alphasun as the Fast Fourier Transform has a subtle bug. Consequently multiplication of very large numbers is now done using the brute force method.
• 14th May 1997. Upgraded the SUN and MSDOS executables on alphasun.
• 8th March 1997. Changed lllhermite1() to lllhermite(). This is the current LLL-based Hermite normal form algorithm as described off my homepage via notes.html.
• 27th February 1997. Replaced the previous LLL based Hermite normal form algorithm in SUN version on alphasun by an aesthetically more satisfying version.
• 20th February 1997. Included LLL based Hermite normal form algorithm in SUN version on alphasun.
• 13th February 1997. Sent the MSDOS update to the Archives and alphasun.
• 7th February 1997. Added a list of CALC functions that a programmer might wish to use.
• 4th February 1997. Added lllgcd() and jacobi(), which do the LLL based extended gcd algorithm and Jacobi's algorithm (1869) to the SUN version on alphasun.
• 18th September 1996. Modified sgcd() to sgcd(N) in private version
• 24th April 1996. Updated MSDOS version on alphasun. Also sent across update to Maths Archives.
• 23rd April 1996. gcda(x,n) was actually calculating the lcm. I also fixed a serious corruption of data was occurring in some functions using FACTOR() such as orderm() and lprimroot().
  (I was freeing stuff that had already been freed.) I also updated Sun version at alphasun.
• 1lth February 1996. Likewise updated MSDOS version at alphasun and Maths Archives.
• 9th February 1996. Updated Sun version at alphasun (runs faster-compiled with the O2 option.
• 11th December 1995. Updated MSDOS version at Maths Archives.
• 8th December 1995. Updated MSDOS version at alphasun.
• 7th December 1995. Updated MSDOS version at alphasun.
  Improved BRENT_POLLARD() doing something I should have done initially - ie. only computed gcd's of every accumulated 10th product.
  Also now bring in the elliptic curve algorithm after 65 digits. MSDOS version of the latter on alphasun not yet updated, nor Archives.
• 6th December 1995. Removed a bug in jacobi(n,m) which only came into effect when n < 0. Only the Sun compiled version calc stored on alphasun has been updated.
  The MSDOS version calc.exe is not yet updated on either the MSDOS Archives or alphasun.
• 5th November 1995: Removed a small error in fund(d) when d=t^2+1, t odd. Also removed a small bug in pell() resulting from using the same name pell in init.c for two functions. Sent update to Maths Archives.
• 10th October 1995: Added addcubicr() and powercubicr(). Updated the MSDOS version locally, but not at Maths Archives.
• 23rd July 1995: I fixed a bug in improvep() which only came into effect if rank(A)=number of rows of A.
• 11th June 1995: I noticed calc was not factoring 2^135+1- a trivial example. I increased the range of Brent-Pollard and the factor base in MPQS for numbers of < 20 digits. I also added an elliptic curve factorization algorithm in case MPQS yields nothing.
• 28th May 1995: Removed a small error in fund(d) when d=t^2+4.