cnum.scaffold

fn.h:

#define JMAX   100

typedef struct {
int    a;
int    b;
int    k;
double d;  // delta() for this ideal
} ideal_t;

gcd.c:

// Extended Euclidean algorithm
//
// It's assumed inputs a and b are both positive
//
// We return gcd and set ap and bp to the corresponding integer
// coefficients such that gcd(a,b) = a*ap + b*bp.  Note that one of ap or bp
// will necessarily be negative
//
int rkqc gcd(int a, int b, int *ap, int *bp)
{
int x = 0;
int y = 1;
int lx = 1;
int ly = 0;
int q, t;

// the so-called "iterative method"; this avoids recursion
//
while (b != 0)
{
q = a / b;

t = a % b;  a = b;  b = t;

t = lx - q*x;  lx = p;  x = t;
t = ly - q*y;  ly = p;  y = t;
}

*ap = x;
*bp = y;

return a;
}

product.c:

#include "fn.h"

// implements the product of two ideals (which may not be reduced)
//
// out:    the product of i1 and i2
// i1, i2: multiplier, multiplicand
// delta:  \Delta
// sdelta: sqrt of \Delta
//
void rkqc product(ideal_t *out, ideal_t *i1, ideal_t *i2, int delta, double sdelta)
{
// ignoring line 1 since I_1 and I_2 are not used subsequently
//
int kh, uh, vh, k, x, w, t;

kh = gcd(i1->a, i2->a, &uh, &vh);
k = gcd(kh, (i1->b + i2->b)/2, &x, &w);
out->a = (i1->a * i2->a / k) / k;

t = (   x * uh * i1->a * i2->b
+ x * vh * i2->a * i1->b
+ (w * i1->b * i2->b + delta)/2
);

out->b = t % (2 * out->a);
if (out->a > sdelta)
out->b -= out->a - 1;
else
out->b += sdelta  - 2*out->a + 1;

// unsure about how this is used (from g-hat and st_product)
out->k = k;
}

tau.c:

// Computes the \tau() function defined in Jozsa (and probably elsewhere)
// defined as implicitly taking \Delta (the fundamental discriminant of our
// real quadratic number field) and computing tau(b, a) = b mod 2a where
// the value is either in (\sqrt{\Delta}-2a, \sqrt{\Delta}]
// if a < \sqrt{\Delta} or in (-a, a] if a is >= \sqrt{Delta}  (see GFI pg 2)
//
// This function is used to normalize ideals to their unique
// canonical representation
//
// Note: we assume that Scaffold's % operator preserves the sign of the
//   dividend (like C99 and Java, but unlike Python and VHDL)
//
// a, b : inputs (note that b can and often will be negative)
// sdelta: sqrt{\Delta} (Delta is the overall input to the class number alg)
//
int rkqc tau(int b, int a, int sdelta)
{
int   ret;

ret = b % (2*a);

if (a > sdelta)
{
while (ret < sdelta - 2*a)
ret += 2*a;
while (ret > sdelta)
ret -= 2*a;
}
else
{
while (ret <= -a)
ret += 2*a;
while (ret > a)
ret -= 2*a;
}

return ret;
}

rho.c:

#include "fn.h"

// \rho in the Class Number algorithm steps between ideals and has various
// guarantees about how much and how little it steps, which equivalence
// class it resides in, how the delta() function reacts to it, etc
// This is all documented in the GFI
//
// \rho is defined on pg 3 of the GFI as \rho(I) = I/\gamma(I), but we
// take a simpler approach here based on the fact that 4a divides \Delta-b^2
// (GFI pg 2) and some algebra.  Here's the idea: let c = (\Delta-b^2)/4a,
// then a' = c and b' = \tau(b, a') in \rho(I).
//
// This function also updates the ideal's delta() value
// The formula is
//   \delta(\rho(I)) = ln abs(\overline(\gamma(I))/\gamma(I)) + \delta(I)
// from app A.1 in GFI;  note that this simplifies (see code)
//
// in:      the input ideal
// out:     \rho(in) is stored here, with its delta() value
// delta:   \Delta  (input to class number algorithm)
// sdelta:  \sqrt{\Delta}
//
void rkqc rho(ideal_t *out, ideal_t *in, int delta, double sdelta)
{
out->a = (abs(delta - in->b*in->b))/(4*in->a);
out->b = tau(-in->b, out->a);
out->d = log(fabs((in->b - sdelta)/(in->b + sdelta))) + in->d;
}

// rho^{-1} is sigma(rho(sigma(I))) where sigma is the algebraic conjugate
//
// The GFI does not specify how to compute rho^{-1}, so I did the algebra
// by hand.  It should probably be double-checked by someone (note: I just
// found parts of Jozsa's paper that back up part of my implementation)
//
void rkqc rho_inv(ideal_t *out, ideal_t *in, int delta, double sdelta)
{
int bs = tau(-in->b, in->a);

out->a = (delta - bs*bs)/(4*in->a);
out->b = tau(bs, out->a);
out->d = -log(fabs((out->b - sdelta)/(out->b + sdelta))) + out->d;
}

st_product.c:

#include "fn.h"

// st_product() implement the "star" product from alg A.2
// Essentially, this product finds the reduced ideal "just to the right"
// of I \cdot J
//
//
void rkqc st_product(ideal_t *res, ideal_t *i, ideal_t *j,
int delta, double sdelta)
{
// skipping line 1 of Alg A.2 beccause a and k are not defined
// (I have email pending to Brian Matt)
ideal_t tmp;
product(&tmp, i, j, delta, sdelta);

if (res->d <= tmp.d)
{
while (res->d <= tmp.d)
rho(&tmp, &tmp, delta, sdelta);
*res = tmp;
}
else
{
while (tmp->d >= tmp.d)
rho_inv(&tmp, &tmp, delta, sdelta);
rho(res, &tmp, delta, sdelta);
}
}

fn.c:

#include "fn.h"

// fn is the main entry point for the oracle function that implements f_N,
// the oracle function which enables approximation of the regulator
//
// bound: i/N + jp/L (see alg A.3 in GFI)
// delta: \Delta (the input to the class number alg)
//
//
void rkqc fn(int bound, int delta, ideal_t *j, int *dist)
{
// used a lot, so compute once
double sdelta = sqrt(delta);

// a reduced ideal is represented by a pair of integers a, b where
// I = Z + (b+sqrt(Delta))/2a Z
//
// We start with O (the ring of integers) which is an ideal of O and has
// delta(O) = 0; O is represented by a=1, b=1

// j holds all ideals as we search; j[0] is O
//
ideal_t j[JMAX];
ideal_t jstar, jtemp;

j[0].a = 1;  j[0].b = 1;  j[0].k = 1;  j[0].d = 0.0;

// run rho() on it twice, always tracking its delta
rho(&j[1], &j[0], delta, sdelta);
rho(&j[2], &j[1], delta, sdelta);

int k = 1;
while (j[k].d <= bound)
{
st_product(&jtemp, &j[k], &jstar, delta, sdelta);
k++;
}

jstar = j[k-1];
k -= 2;

while (k > 0)
{
st_product(&jtemp, &j[k], &jstar, delta, sdelta);
if (jtemp.d <= bound)
jstar = jtemp;
k--;
}

while (jstar.d < bound)
{
rho(&jstar, &jstar, delta, sdelta);
jtemp = jstar;
rho(&jstar, &jstar, delta, sdelta);
}

if (jtemp.d >= bound)
*j = jtemp;
else
*j = jstar;

*dist = floor(N * (bound - j->d));
}

ghat.c:

#include "fn.h"

// builds the ideal I, the product of the generators classically produced
// for phase II
//
// out: the output ideal I
// g[]: the input generators
// p[]: corresponding powers of g[]
// k:   number of ideals in g[]
//
void rkqc ghat(ideal_t *out, ideal_t g[], int p[], int k,
int delta, double sdelta)
{
ideal_t   a, t;
int       i, j;

a->a = 1;
a->b = 1;
a->k = 1;
a->delta = 1.0;

for (i=0; i < k; i++)
{
// compute powers first
for (j=0; j < p[i]; j++)
product(&a, &a, &g[i], delta, sdelta);
}

// now ensure a is reduced
//
while (a->k != 1)
rho(&a, &a, delta, sdelta)

*out = a;
}

fjn.c:

#include "fn.h"

// fjn is the main entry point for the oracle function that implements f_J,N,
// the oracle function which enables approximation of the regulator
//
// bound: i/N + jp/L (see alg A.3 in GFI)
// delta: \Delta (the input to the class number alg)
// sdelta: floor(sqrt(delta))
//
//
void rkqc fjn(int bound, int delta, int sdelta,
ideal_t *j, ideal_t *out, int *dist)
{
ideal_t i, k, ktemp;

i = *j;
fn(bound, delta, i, dist);
product(&k, &i, j, delta, sdelta);

while (k.k != 1)
{
rho(&k, &k, delta, sdelta);
}

if (k.d <= bound)
{
while (k.d < bound)
rho(&k, &k, delta, sdelta);
}
else
{
while (k.d >= bound)
rhoinv(&k, &k, delta, sdelta);

rho(&k, &k, delta, sdelta);
}

*out = k;
*dist = floor(N * (bound - j->d));
}

h.c:

#include "fn.h"

// h computes the h_I oracle
//
// bound: b/N + jhatp/L (see alg A.6 in GFI)
// delta: \Delta (the input to the class number alg)
//
//
void rkqc h(int a, int b, int bound, int delta, double sdelta,
ideal_t *i, ideal_t *out, int *dist)
{
ideal_t j, k;
int li, dst;

j = *i;
k.a = 1;  k.b = 1;  k.k = 1;  k.d = 0.0;
li = a;

while (k.k != 1)
{
if (i % 2 == 1)
st_product(&k, &k, &j, delta, sdelta);
st_product(&j, &j, &j, delta, sdelta);
}
fn(bound, delta, &j, &dst);
st_product(&k, &k, &j, delta, sdelta);

*out = *k;
*dist = floor(N * k.d);
}