3
5
Fork 0
mirror of git://git.savannah.gnu.org/guix.git synced 2023-12-14 03:33:07 +01:00

gnu: Add ratpoints.

* gnu/packages/maths.scm (ratpoints): New variable.
* gnu/packages/patches/ratpoints-sturm_and_rp_private.patch: New file.
* gnu/local.mk (dist_patch_DATA): Reference patch.
This commit is contained in:
Nicolas Goaziou 2019-06-19 21:43:12 +02:00
parent 7c5f623192
commit 0c842e3a59
No known key found for this signature in database
GPG key ID: DA00B4F048E92F2D
3 changed files with 236 additions and 0 deletions

View file

@ -1233,6 +1233,7 @@ dist_patch_DATA = \
%D%/packages/patches/randomjungle-disable-static-build.patch \ %D%/packages/patches/randomjungle-disable-static-build.patch \
%D%/packages/patches/rapicorn-isnan.patch \ %D%/packages/patches/rapicorn-isnan.patch \
%D%/packages/patches/raptor2-heap-overflow.patch \ %D%/packages/patches/raptor2-heap-overflow.patch \
%D%/packages/patches/ratpoints-sturm_and_rp_private.patch \
%D%/packages/patches/ratpoison-shell.patch \ %D%/packages/patches/ratpoison-shell.patch \
%D%/packages/patches/rcs-5.9.4-noreturn.patch \ %D%/packages/patches/rcs-5.9.4-noreturn.patch \
%D%/packages/patches/rct-add-missing-headers.patch \ %D%/packages/patches/rct-add-missing-headers.patch \

View file

@ -5003,3 +5003,44 @@ command-line tools, and an Application Programming Interface (API).
This package provides the static libraries required to run programs This package provides the static libraries required to run programs
compiled against the nauty library.") compiled against the nauty library.")
(license license:asl2.0))) (license license:asl2.0)))
(define-public ratpoints
(package
(name "ratpoints")
(version "2.1.3")
(source (origin
(method url-fetch)
(uri (string-append
"http://www.mathe2.uni-bayreuth.de/stoll/programs/"
"ratpoints-" version ".tar.gz"))
(sha256
(base32
"0zhad84sfds7izyksbqjmwpfw4rvyqk63yzdjd3ysd32zss5bgf4"))
(patches
;; Taken from
;; <https://git.sagemath.org/sage.git/plain/build/pkgs/ratpoints/patches/>
(search-patches "ratpoints-sturm_and_rp_private.patch"))))
(build-system gnu-build-system)
(arguments
`(#:test-target "test"
#:make-flags
(list (string-append "INSTALL_DIR=" (assoc-ref %outputs "out")))
#:phases
(modify-phases %standard-phases
(delete 'configure) ;no configure script
(add-before 'install 'create-install-directories
(lambda* (#:key outputs #:allow-other-keys)
(let ((out (assoc-ref outputs "out")))
(mkdir-p out)
(with-directory-excursion out
(for-each (lambda (d) (mkdir-p d))
'("bin" "include" "lib"))))
#t)))))
(inputs
`(("gmp" ,gmp)))
(home-page "http://www.mathe2.uni-bayreuth.de/stoll/programs/")
(synopsis "Find rational points on hyperelliptic curves")
(description "Ratpoints tries to find all rational points within
a given height bound on a hyperelliptic curve in a very efficient way,
by using an optimized quadratic sieve algorithm.")
(license license:gpl2+)))

View file

@ -0,0 +1,194 @@
diff --git a/rp-private.h b/rp-private.h
index b4c7dad..0c7193e 100644
--- a/rp-private.h
+++ b/rp-private.h
@@ -36,7 +36,7 @@
#define LONG_SHIFT ((LONG_LENGTH == 16) ? 4 : \
(LONG_LENGTH == 32) ? 5 : \
(LONG_LENGTH == 64) ? 6 : 0)
-#define LONG_MASK (~(-1L<<LONG_SHIFT))
+#define LONG_MASK (~(-(1L<<LONG_SHIFT)))
/* Check if SSE instructions can be used.
We assume that one SSE word of 128 bit is two long's,
diff --git a/sturm.c b/sturm.c
index c78d7c6..5fd2cf5 100644
--- a/sturm.c
+++ b/sturm.c
@@ -27,7 +27,6 @@
***********************************************************************/
#include "ratpoints.h"
-
/**************************************************************************
* Arguments of _ratpoints_compute_sturm() : (from the args argument) *
* *
@@ -53,7 +52,7 @@
/* A helper function: evaluate the polynomial in cofs[] of given degree
at num/2^denexp and return the sign. */
-static long eval_sign(ratpoints_args *args, mpz_t *cofs, long degree,
+static long eval_sign(const ratpoints_args *args, const mpz_t *cofs, long degree,
long num, long denexp)
{
long n, e, s;
@@ -70,11 +69,80 @@ static long eval_sign(ratpoints_args *args, mpz_t *cofs, long degree,
return(s);
}
+static const long max = (long)(((unsigned long)(-1))>>1);
+static const long min = (long)(-(((unsigned long)(-1))>>1));
+ /* recursive helper function */
+static void iterate(long nl, long nr, long del, long der, long cleft, long cright,
+ long sl, long sr, long depth,
+ ratpoints_interval **iptr, const ratpoints_interval *ivlo,
+ const ratpoints_args *args, const long k, const long sturm_degs[],
+ const mpz_t sturm[][args->degree + 1])
+ { /* nl/2^del, nr/2^der : interval left/right endpoints,
+ cleft, cright: sign change counts at endpoints,
+ sl, sr: signs at endpoints,
+ depth: iteration depth */
+ long iter = args->sturm;
+ if(cleft == cright && sl < 0) { return; }
+ /* here we know the polynomial is negative on the interval */
+ if((cleft == cright && sl > 0) || depth >= iter)
+ /* we have to add/extend an interval if we either know that
+ the polynomial is positive on the interval (first condition)
+ or the maximal iteration depth has been reached (second condition) */
+ { double l = ((double)nl)/((double)(1<<del));
+ double u = ((double)nr)/((double)(1<<der));
+ if(*iptr == ivlo)
+ { (*iptr)->low = l; (*iptr)->up = u; (*iptr)++; }
+ else
+ { if(((*iptr)-1)->up == l) /* extend interval */
+ { ((*iptr)-1)->up = u; }
+ else /* new interval */
+ { (*iptr)->low = l; (*iptr)->up = u; (*iptr)++; }
+ }
+ return;
+ }
+ /* now we must split the interval and evaluate the sturm sequence
+ at the midpoint */
+ { long nm, dem, s0, s1, s2, s, cmid = 0, n;
+ if(nl == min)
+ { if(nr == max) { nm = 0; dem = 0; }
+ else { nm = (nr == 0) ? -1 : 2*nr; dem = 0; }
+ }
+ else
+ { if(nr == max) { nm = (nl == 0) ? 1 : 2*nl; dem = 0; }
+ else /* "normal" case */
+ { if(del == der) /* then both are zero */
+ { if(((nl+nr) & 1) == 0) { nm = (nl+nr)>>1; dem = 0; }
+ else { nm = nl+nr; dem = 1; }
+ }
+ else /* here one de* is greater */
+ { if(del > der) { nm = nl + (nr<<(del-der)); dem = del+1; }
+ else { nm = (nl<<(der-del)) + nr; dem = der+1; }
+ }
+ }
+ }
+ s0 = eval_sign(args, sturm[0], sturm_degs[0], nm, dem);
+ s1 = eval_sign(args, sturm[1], sturm_degs[1], nm, dem);
+ if(s0*s1 == -1) { cmid++; }
+ s = (s1 == 0) ? s0 : s1;
+ for(n = 2; n <= k; n++)
+ { s2 = eval_sign(args, sturm[n], sturm_degs[n], nm, dem);
+ if(s2 == -s) { cmid++; s = s2; }
+ else if(s2 != 0) { s = s2; }
+ }
+ /* now recurse */
+ iterate(nl, nm, del, dem, cleft, (s0==0) ? (cmid+1) : cmid,
+ sl, (s0==0) ? -s1 : s0, depth+1,
+ iptr, ivlo, args, k, sturm_degs, sturm);
+ iterate(nm, nr, dem, der, cmid, cright,
+ (s0==0) ? s1 : s0, sr, depth+1,
+ iptr, ivlo, args, k, sturm_degs, sturm);
+ }
+ } /* end iterate() */
+
long _ratpoints_compute_sturm(ratpoints_args *args)
{
mpz_t *cofs = args->cof;
long degree = args->degree;
- long iter = args->sturm;
ratpoints_interval *ivlist = args->domain;
long num_iv = args->num_inter;
long n, m, k, new_num;
@@ -165,75 +233,12 @@ long _ratpoints_compute_sturm(ratpoints_args *args)
/* recall: typedef struct {double low; double up;} ratpoints_interval; */
{ ratpoints_interval ivlocal[1 + (degree>>1)];
ratpoints_interval *iptr = &ivlocal[0];
- long max = (long)(((unsigned long)(-1))>>1);
- long min = -max;
long num_intervals;
long slcf = mpz_cmp_si(cofs[degree], 0);
- /* recursive helper function */
- void iterate(long nl, long nr, long del, long der, long cleft, long cright,
- long sl, long sr, long depth)
- { /* nl/2^del, nr/2^der : interval left/right endpoints,
- cleft, cright: sign change counts at endpoints,
- sl, sr: signs at endpoints,
- depth: iteration depth */
- if(cleft == cright && sl < 0) { return; }
- /* here we know the polynomial is negative on the interval */
- if((cleft == cright && sl > 0) || depth >= iter)
- /* we have to add/extend an interval if we either know that
- the polynomial is positive on the interval (first condition)
- or the maximal iteration depth has been reached (second condition) */
- { double l = ((double)nl)/((double)(1<<del));
- double u = ((double)nr)/((double)(1<<der));
- if(iptr == &ivlocal[0])
- { iptr->low = l; iptr->up = u; iptr++; }
- else
- { if((iptr-1)->up == l) /* extend interval */
- { (iptr-1)->up = u; }
- else /* new interval */
- { iptr->low = l; iptr->up = u; iptr++; }
- }
- return;
- }
- /* now we must split the interval and evaluate the sturm sequence
- at the midpoint */
- { long nm, dem, s0, s1, s2, s, cmid = 0, n;
- if(nl == min)
- { if(nr == max) { nm = 0; dem = 0; }
- else { nm = (nr == 0) ? -1 : 2*nr; dem = 0; }
- }
- else
- { if(nr == max) { nm = (nl == 0) ? 1 : 2*nl; dem = 0; }
- else /* "normal" case */
- { if(del == der) /* then both are zero */
- { if(((nl+nr) & 1) == 0) { nm = (nl+nr)>>1; dem = 0; }
- else { nm = nl+nr; dem = 1; }
- }
- else /* here one de* is greater */
- { if(del > der) { nm = nl + (nr<<(del-der)); dem = del+1; }
- else { nm = (nl<<(der-del)) + nr; dem = der+1; }
- }
- }
- }
- s0 = eval_sign(args, sturm[0], sturm_degs[0], nm, dem);
- s1 = eval_sign(args, sturm[1], sturm_degs[1], nm, dem);
- if(s0*s1 == -1) { cmid++; }
- s = (s1 == 0) ? s0 : s1;
- for(n = 2; n <= k; n++)
- { s2 = eval_sign(args, sturm[n], sturm_degs[n], nm, dem);
- if(s2 == -s) { cmid++; s = s2; }
- else if(s2 != 0) { s = s2; }
- }
- /* now recurse */
- iterate(nl, nm, del, dem, cleft, (s0==0) ? (cmid+1) : cmid,
- sl, (s0==0) ? -s1 : s0, depth+1);
- iterate(nm, nr, dem, der, cmid, cright,
- (s0==0) ? s1 : s0, sr, depth+1);
- }
- } /* end iterate() */
-
iterate(min, max, 0, 0, count2, count1,
- (degree & 1) ? -slcf : slcf, slcf, 0);
+ (degree & 1) ? -slcf : slcf, slcf, 0,
+ &iptr, &ivlocal[0], args, k, sturm_degs, sturm);
num_intervals = iptr - &ivlocal[0];
/* intersect with given intervals */
{ ratpoints_interval local_copy[num_iv];