text:
# Things that are falling apart

Bergschaden (sufrace damage caused by mining)

Living in a region which saw a lot of mining over the last centuries has some advantages but also some disadvantages. One of the disadvantages is the fact that the similarities of the underground with Swiss cheese are quite astonishing. The holes close after some time but the movement caused by the closing of the underground gaps goes up and breaks houses. It is not much and repairable and the mining company pays for it but you still have all the dirt and the hassle and all the things nobody really wants.
The damage is hardly visible on the picture but you can still see at the top left of the window and the lower right a long, already plastered crack (there were more, off the picture). The walls get a rough-cast and an insulation, all done by me, so don’t expect too much work the the parser for now.
The insulation (no room for it outside, so it has to be done inside, which makes it brutally more expensive) was what it started in the first place. The kitchen has three outside walls and after insulating one (outside) we got problems with condensing water which later caused some mould to settle. Ripping off the wallpapers showed some cracks, not expected, we know where we live. One call and a day later, the cracks were inspected, okay’d and another day later repaired. Not only plastered but glued together by injecting, with high pressure a synthetic resin. That at the beginning of october and the picture is from today. What happened? The first inspection was done by a woman, well known to sit on the money but also fast to decide if the client is okay with the cheapest solution. This client was me, because the only thing I wanted was to get the cracks glued shut to be able to go on with the insulation. So it went fast and after some drying time, I just cut open the first bag of concrete when the phone rang. A colleague of the aforementioned woman wanted to look at it at the next available date. Available for him, not for me and the next date was one week later.
Ohhhhkaaay…
He finally came, looked around, and after waving some instruments at the walls and a short chat he left. I do not know exactly why he did it but I have my suspicions and I don’t think it is something nice.
One week lost.
Ohhhhkaaay…
On with the work and on with the rough cast, or so I thought.
I made the cast from some leftover concrete-mix (the only-add-water type) which was a bad idea: it was probably too old and didn’t harden.
Another week lost.
Ohhhhkaaay…
And it is not the case that this wall is the single thing I have to do, and so:
Another week lost.
And we write the 31st of October now, Halloween.
On a side-note: the final rendering for the insulation is some expensive stuff which I ordered and payed in advance and is still not available. But at least:
No week lost.
Yet.
It is unknown if we’ll have a renovated kitchen before christmas.
And how was your October?

# Adventures of a Programmer: Parser Writing Peril VII

The complete program should work on any operating system it can run on. It will be written in ISO/IEE C-99 without using any of the more exotic parts of it. The main problem is the build. There will be more than a handful of source file and dependencies to do it manually, some automation is needed. Tools for automated building that run cross-platform are rare, only two are large enough to have been shortlisted: qmake, the autotool delivered with QT (Who has bought it now, still Nokia?) which runs on several operating systems including some phones. It is very simple. Too simple, sadly. I tried it with a small subproject and put it up at Github under the name libczrational for all to see and writhe with disgust. Seems as if it has to be cmake.

The library has, for now, basic exact arithmetic implemented (+,-,*,/), basic inexact arithmetic (sqrt) and two integer functions (Bernoulli and Euler numbers). The Bernoulli and Euler number calculations are meant as examples but are already fast enough for general use (uses Brent-Harvey’s algorithm which is based on tangent and secant numbers).

The documentation looks good, but that is all that it does, for now 😉

The makefile in the packet has been generated on Debian Squeeze.

# Adventures of a Programmer: Parser Writing Peril VI

Yep, you guessed right, still no grammar 😉
But I implemented a square root for the big rationals from the last post. It is not very elegant, especially too ignorant about the exact value of the precision which it puts too high, but it works.

#ifndef MPI_LOG2
#   define MPI_LOG2 3.321928094887362347870319429489390175864831393024580612054756
#endif
#include <math.h>
static int mpq_epsilon;
void mpq_set_epsilon(int eps){
if(eps <= 0){
mpq_epsilon = 0;
}
else{
/* eps assumed to be precision in number of decimal digits */
mpq_epsilon = (int)((double)eps*MPI_LOG2) + (int)log2((double)eps);
}
}
int mpq_get_epsilon(void){
return mpq_epsilon;
}
/* sqrt() meant as a test for epsilon */
int mpq_sqrt(mp_rat *a, mp_rat *c){
if(mpq_isneg(a)){
return MP_VAL;
}
mpq_copy(a,c);
printf("epsilon = %d\n",mpq_get_epsilon());
/*
This is too brutal.
There are better ways to do it.
*/
int eps  =  mpq_get_epsilon();
mp_mul_2d(&c->numerator,  eps,&c->numerator);
mp_mul_2d(&c->denominator,eps,&c->denominator);

mp_sqrt(&c->numerator,&c->numerator);
mp_sqrt(&c->denominator,&c->denominator);
mpq_reduce(c);
return MP_OKAY;
}

The angstzuschlag of log2(eps) (the logarithm base two of eps which holds the number of decimal digits of precision) is purely arbitrary.
The problem with this approach is that any perfect square will result in unnecessary inexact results when eps is not zero.
Changing the bit-shifting to multiplying with exponents of ten will get rid of that problem, lets…seeeeee…
Segmentation Fault!
Oops!
What the…
Ah, the function mpq_reduce() uses mpq_set_int() which itself uses mpq_reduce() which resulted in a nice crash. The first problem is easy to solve:

#define mp_isone(a) (((a)->used == 1 && (a)->dp[0] == 1) ? MP_YES : MP_NO)
#define MPQ_DIVISION_BY_ZERO -997
int mpq_reduce(mp_rat *a){
mp_int gcd;
int e;
/* check if num is zero -> return 0 if it is*/
if(mp_iszero(&a->numerator)){
mp_set_int(&a->numerator,0);
mp_set_int(&a->denominator,1);
if( (e = mpq_set_sign(a,MP_ZPOS) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}
/* check if den is 1 -> return if it is*/
if(mp_isone(&a->denominator)){
mpq_normalize_sign(a);
return MP_OKAY;
}
/* check if den is zero -> return undefined if it is*/
if(mp_iszero(&a->denominator)){
mp_set_int(&a->numerator,0);
mp_set_int(&a->denominator,0);

return MPQ_DIVISION_BY_ZERO;
}
/* check if den = num -> return 1 */
if(mpq_isone(a)){
mpq_normalize_sign(a);
mp_set_int(&a->numerator,1);
mp_set_int(&a->denominator,1);

return MP_OKAY;
}
mp_init(&gcd);
/* compute gcd(num,den) */
if( (e = mp_gcd(&a->numerator,&a->denominator,&gcd) ) != MP_OKAY){
return e;
}
/* set num to num/gcd */
if( (e = mp_div(&a->numerator,&gcd,&a->numerator,NULL) ) != MP_OKAY){
return e;
}
/* set den to den/gcd */
if( (e = mp_div(&a->denominator,&gcd,&a->denominator,NULL) ) != MP_OKAY){
return e;
}
/* not to forget: */
mpq_normalize_sign(a);
return MP_OKAY;
}

The second problem is not so easy to solve because it shouldn’t have segfaulted before it hit the stack’s recursion limit. I leave that for later.
Where have I been…ah, yes, the problem with perfect squares. Using decimals as the multiplier doesn’t help at all, just give it a try.

  //mp_mul_2d(&c->numerator,  eps,&c->numerator);
//mp_mul_2d(&c->denominator,eps,&c->denominator);
mp_int tenpow;
mp_init(&tenpow);
mp_set_int(&tenpow,10);
mp_expt_d(&tenpow,eps,&tenpow);
mp_mul(&c->numerator, &tenpow ,&c->numerator);
mp_mul(&c->denominator,&tenpow,&c->denominator);
...
mp_clear(&tenpow);

It makes no difference if you expand the fraction with a power of ten or a power of two if you expand both sides of the fraction.
The actual precision gained with this method is also too low. The number of digits n of a square root of a number x with N is about N/2 depending on the number system, exactly half with base two. So we need twice the precision to get a square root of normal precision. We are using shifts, so simply double eps.
Note to self for further optimization: subtract number of digits of min(nbits(num),nbits(den)) from eps.

  /* check if eps<<1 is still smaller than the datatype used for eps!*/
mp_mul_2d(&c->numerator,  eps<<1,&c->numerator);
mp_mul_2d(&c->denominator,eps<<1,&c->denominator);

And a funny thing happens: it gives the exact result now.
Why?
Well, I do have a printf() giving me the actual value. I set epsilon to 50 which gives a result of 171 here. Which is odd, because doubling it is even curioser. Ok, ok, I’ll stop here 😉
The reason is that the number of digits of the square root of a number in base two is exactly half of the number of digits of the radicand. If the number of digits of the radicand is odd it cannot work it out exactly.
A hypothesis easily tested 😉
Two fractions given: 529/1849 ((23/43)^2), 3481/11664 ((59/108)^2), 1764/1849 ((42/43)^2), and 44521/46656 ((211/216)^2). Uh, no, that’s illegible. Let’s make a table.
With eps = 50

a/b (a/b)^2 odd eps even eps
23/43 529/1849 1258…143/235…312 23/43
59/108 3481/11664 322…847/590…109 59/108
42/43 1764/1849 229…653/235…312 42/43
211/216 44521/46656 115…927/118…219 211/216

Well, that’s what most scientist call a “failed chance to publish”. Because not many journals accept experiments with a negative outcome. That causes many already failed experiments to get duplicated because the experimentalist had no chance to know that it doesn’t work in the first playe. That heap of a male cow’s excrements wastes a lot of money and hours.
So, if you have such a failed experiment you can publish it at e.g.: arxive.org. But check if you will not get punished for it. Yes, that might happen!
But I digress.
So it is the even multiplier that counts?
Why, you ask?
Yes, that’s an interesting question.

int mpq_sqrt(mp_rat *a, mp_rat *c){
int eps;
if(mpq_isneg(a)){
return MP_VAL;
}

/* place shortcut for zero and one  here */
/* place shortcut for integers here */

mpq_copy(a,c);
eps  =  mpq_get_epsilon();
if((eps & 0x1)){
eps++;
}
eps <<= 1;

/* reduce size of eps r.t. to min(digits(num),digits(den)) here */

mp_mul_2d(&c->numerator,  eps,&c->numerator);
mp_mul_2d(&c->denominator,eps,&c->denominator);
mp_sqrt(&c->numerator,&c->numerator);
mp_sqrt(&c->denominator,&c->denominator);
mpq_reduce(c);
return MP_OKAY;
}

# Adventures of a Programmer: Parser Writing Peril V

Part five and still no grammar? Oh, my…

Before I blow the dust off of libtomfloat I decided to try it with rationals first and armed myself with my good ol’ math-textbook from primary school and D. Knuth’s second volume of TAoCP.
I also try to get a bit more structure into it this time 😉

Apropos structure: the main data structure shall hold a fraction and a fraction consists of a numerator and a denominator, as every pupil knows. What not every pupil knows (well, some do!) is that there is a sign, too. So the structure is a C struct with some aptly named inhabitants:

typedef struct  {
int sign;
mp_int numerator;
mp_int denominator;
} mp_rat;

From the data type mp_int it is a reasonable assumption that I still use Libtommath as the big-number library. The name of this structure has been chosen very carefully—the ‘r’ sounds better after a ‘p’ than after a ‘f’.

The mp_int numbers need to get initialized before and freed after use, and so do the mp_rat numbers. The implementation makes use of the respective functions of libtommath to keep it simple.
Initialisation of a single mp_rat number is the simplest, we use libtommath’s mp_init_multi() directly.

int mpq_init(mp_rat *a){
int e;
if( (e = mp_init_multi(&a->numerator,&a->denominator,NULL) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}

initializing several rationals needs a bit more work, and the inclusion of stdarg.h, too. But that was all of the unpaid overtime, we can take the rest more or less verbatim from libtommath. Including the comments.

#include <stdarg.h>
int mpq_init_multi(mp_rat *mp, ...){
mp_err res = MP_OKAY;      /* Assume ok until proven otherwise */
int n = 0;                 /* Number of ok inits */
mp_rat* cur_arg = mp;
va_list args;

va_start(args, mp);        /* init args to next argument from caller */
while (cur_arg != NULL) {
if (mpq_init(cur_arg) != MP_OKAY) {
/* Oops - error! Back-track and mp_clear what we already
succeeded in init-ing, then return error.
*/
va_list clean_args;

/* end the current list */
va_end(args);

/* now start cleaning up */
cur_arg = mp;
va_start(clean_args, mp);
while (n--) {
mpq_clear(cur_arg);
cur_arg = va_arg(clean_args, mp_rat*);
}
va_end(clean_args);
res = MP_MEM;
break;
}
n++;
cur_arg = va_arg(args, mp_rat*);
}
va_end(args);
return res;                /* Assumed ok, if error flagged above. */
}

Freeing the numbers works similar. We take the liberty to not only steal code from libtommath but the names, too.

void mpq_clear(mp_rat *a){
mp_clear_multi(&a->numerator,&a->denominator,NULL);
a->sign = MP_ZPOS;
}

void mpq_clear_multi(mp_rat *mp, ...){
mp_rat* next_mp = mp;
va_list args;
va_start(args, mp);
while (next_mp != NULL) {
mpq_clear(next_mp);
next_mp = va_arg(args, mp_rat*);
}
va_end(args);
}

One of the most useful helper-functions is the printing function, of course, so let’s implement one.

void mpq_print(mp_rat * a){
/* If ypu wonder what the zero does, try it without */
printf("%c",(a->sign == MP_ZPOS)?0:'-');
mp_fwrite(&a->numerator,10,stdout);
putchar('/');
mp_fwrite(&a->denominator,10,stdout);
}

A reading function is not implemented, yet, at least not directly but it is easy to build one from the following functions. Let’s go from left to right because that’s how fractions get written in most of the western alphabets: ±x/y
So, at first, there need to be some functions to handle the sign.

/* sets sign of mp_rat a to sign given in "sign", returns old sign */
int mpq_set_sign(mp_rat *a, int sign){
int sgn;
sgn = a->sign;
a->sign = (sign < 0 || sign == MP_NEG)?MP_NEG:MP_ZPOS;
a->numerator.sign = MP_ZPOS;
a->denominator.sign = MP_ZPOS;
return sgn;
}
/* just for symmetry: returns the sign of mp_rat a */
int mpq_get_sign(mp_rat *a){
return a->sign;
}

We keep the numerator and the denominator strictly positive but that can can change in some of the calculations, so we need something to set the sign of mp_ratfrom the signs of the numerator and denominator.

void mpq_normalize_sign(mp_rat *a){
int sign_n, sign_d;
sign_n = a->numerator.sign;
sign_d = a->denominator.sign;
a->sign = (sign_n != sign_d)?MP_NEG:MP_ZPOS;
a->numerator.sign = MP_ZPOS;
a->denominator.sign = MP_ZPOS;
}

And while we are at it, we should implement negation, too:

int mpq_neg(mp_rat *a){
int sgn;
sgn =  mpq_get_sign(a);
sgn = (sgn == MP_NEG)?MP_ZPOS:MP_NEG;
a->numerator.sign = MP_ZPOS;
a->denominator.sign = MP_ZPOS;
return mpq_set_sign(a, sgn);
}

The next one to read is the numerator and after that the denominator. Reading the slash or whatever else denotes a fraction is not part of this library. But might be in the future, who knows?

/* sets given mp_int b as the numerator in mp_rat a */
int mpq_set_num(mp_rat *a, mp_int *b){
int e;
if( (e = mp_copy(b,&a->numerator) ) != MP_OKAY){
return e;
}
if( (e = mpq_reduce(a) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}

/* sets given mp_int b as the denominator in mp_rat a */
int mpq_set_den(mp_rat *a, mp_int *b){
int e;
if( (e = mp_copy(b,&a->denominator) ) != MP_OKAY){
return e;
}
if( (e = mpq_reduce(a) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}

And now all together because everybody just loooooves syntactic sugar!

/* sets given mp_int p,q as the numerator,denominator in mp_rat a */
int mpq_set(mp_rat *a, mp_int *p, mp_int *q){
int e;

if( (e = mp_copy(p,&a->numerator) ) != MP_OKAY){
return e;
}
if( (e = mp_copy(q,&a->denominator) ) != MP_OKAY){
return e;
}
if( (e = mpq_reduce(a) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}

More sweets? Ok, Halloween is near, so here are the methods and functions to handle normal C-99 data types, unsigned long as an example:

/* sets given unsigned long b as the denominator in mp_rat a */
int mpq_set_den_int(mp_rat *a, unsigned long b){
int e;

if( (e = mp_set_int(&a->denominator,b) ) != MP_OKAY){
return e;
}
if( (e = mpq_reduce(a) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}
/* sets given unsigned longs p,q as the numerator,denominator in mp_rat a */
int mpq_set_int(mp_rat *a, unsigned long p,unsigned long q){
int e;

if( (e = mp_set_int(&a->numerator,p) ) != MP_OKAY){
return e;
}
if( (e = mp_set_int(&a->denominator,q) ) != MP_OKAY){
return e;
}
if( (e = mpq_reduce(a) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}

You might have found out that one of the functions used mpq_reduce() which has not been handled, yet. That is correct. keeping the fractions in their reduced forms is essential to working with rationals, especially if done numerically. The numbers can get very large, very fast otherwise.

But before we come it we need two more functions to handle copying.
Deep copy:

int mpq_copy(mp_rat *a, mp_rat *b){
int e;
if( (e =    mp_copy(&a->numerator,&b->numerator) ) != MP_OKAY){
return e;
}
if( (e =    mp_copy(&a->denominator,&b->denominator) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}

Pointer exchange, only:

void mpq_exch(mp_rat *a, mp_rat *b ){
mp_exch(&a->numerator,&b->numerator);
mp_exch(&a->denominator,&b->denominator);
}

Reducing a fraction is simple but computationally expensive, so we not only need some checks and balances but some shortcuts, too. With the following equations in mind, we can write some of these short cuts.
$\displaystyle{ \frac{0}{x}= 0 }$

$\displaystyle{ \frac{x}{x}= 1 }$

But $\displaystyle{ \frac{0}{0}= \text{undefined} }$ because of $\displaystyle{ \frac{x}{0}= \text{undefined} }$

$\displaystyle{ \frac{x}{1}= 1 }$

All of the above can be used more often and are worth their own functions, we take three of them.

#define mp_isone(a) (((a)->used == 1 && (a)->dp[0] == 1) ? MP_YES : MP_NO)
int mpq_iszero(mp_rat *a){
if(     mp_iszero(&a->numerator)
&&  mp_isone(&a->denominator)){
return MP_YES;
}
return MP_NO;
}
int mpq_isone(mp_rat *a){
if(mp_cmp(&a->numerator,&a->denominator) == MP_EQ){
return MP_YES;
}
return MP_NO;
}
int mpq_isinteger(mp_rat *a){
if(mp_isone(&a->denominator)){
return MP_YES;
}
return MP_NO;
}

The actual reducing is simply
$\dfrac{p/\gcd\left(p,q\right)}{q/\gcd\left(p,q\right)}$

#define MPQ_DIVISION_BY_ZERO -997
int mpq_reduce(mp_rat *a){
mp_int gcd;
int e;
/* check if num is zero -> return 0 if it is*/
if(mp_iszero(&a->numerator)){
if( (e = mpq_set_int(a,0,1) ) != MP_OKAY){
return e;
}
if( (e = mpq_set_sign(a,MP_ZPOS) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}
/* check if den is 1 -> return if it is*/
if(mp_isone(&a->denominator)){
mpq_normalize_sign(a);
return MP_OKAY;
}
/* check if den is zero -> return undefined if it is*/
if(mp_iszero(&a->denominator)){
if( (e = mpq_set_int(a,0,0) ) != MP_OKAY){
return e;
}
return MPQ_DIVISION_BY_ZERO;
}
/* check if den = num -> return 1 */
if(mpq_isone(a)){
mpq_normalize_sign(a);
if( (e = mpq_set_int(a,1,1) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}
mp_init(&gcd);
/* compute gcd(num,den) */
if( (e = mp_gcd(&a->numerator,&a->denominator,&gcd) ) != MP_OKAY){
return e;
}
/* set num to num/gcd */
if( (e = mp_div(&a->numerator,&gcd,&a->numerator,NULL) ) != MP_OKAY){
return e;
}
/* set den to den/gcd */
if( (e = mp_div(&a->denominator,&gcd,&a->denominator,NULL) ) != MP_OKAY){
return e;
}
/* not to forget: */
mpq_normalize_sign(a);
return MP_OKAY;
}

It might not be the wish of the programmers using this library to do a full stop every time they approach a division by zero, it might be more useful sometimes (I have to admit, I couldn’t come up with a single example) to keep it and take it along the calculations so it might be a good idea (I doubt it, but who am I?) to write a test for it.

int mpq_isundefined(mp_rat *a){
if(     mp_iszero(&a->numerator)
&&  mp_iszero(&a->denominator)){
return MP_YES;
}
return MP_NO;
}

Ok, that were the functions for the household chores. Now for the real stuff.
The basic computations one should be able to do are:

• Addition
• Subtraction
• Multiplication
• Division

Funnily, the simplest thing to do in integer calculation is the most complicated when handling fractions and vice versa.
Addition first.

int mpq_add(mp_rat *a, mp_rat *b, mp_rat *c){
int e;
mp_int gcd,tmp;
if(mp_iszero(&a->numerator)){
if( (e = mpq_copy(b,c) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}
if(mp_iszero(&b->numerator)){
if( (e = mpq_copy(a,c) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}
/* let libtommath take care of the signs */
(&a->numerator)->sign = a->sign;
(&b->numerator)->sign = b->sign;
if(mpq_isinteger(a) && mpq_isinteger(a) ){
if( (e = mp_add(&a->numerator,&b->numerator,&c->numerator) ) != MP_OKAY){
return e;
}
c->sign = (&c->numerator)->sign;
(&a->numerator)->sign = MP_ZPOS;
(&b->numerator)->sign = MP_ZPOS;
(&c->numerator)->sign = MP_ZPOS;
(&c->denominator)->sign = MP_ZPOS;
if( (e = mp_set_int(&c->denominator,1) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}
mp_init_multi(&gcd,&tmp,NULL);
/*gcd = gcd(b,d)*/
if( (e = mp_gcd(&a->denominator,&b->denominator,&gcd) ) != MP_OKAY){
return e;
}
/*c->numerator = ((A*D)/gcd) + ((C*B)/gcd)*/
/*c->numerator = ((A*D)/gcd) */
if( (e = mp_mul(&a->numerator,&b->denominator,&tmp) ) != MP_OKAY){
return e;
}
if( (e = mp_div(&tmp,&gcd,&c->numerator,NULL) ) != MP_OKAY){
return e;
}

/* ((C*B)/gcd)*/
if( (e = mp_mul(&b->numerator,&a->denominator,&tmp) ) != MP_OKAY){
return e;
}
if( (e = mp_div(&tmp,&gcd,&tmp,NULL) ) != MP_OKAY){
return e;
}
if( (e = mp_add(&tmp,&c->numerator,&c->numerator) ) != MP_OKAY){
return e;
}
if(mp_iszero(&c->numerator)){
if( (e = mp_set_int(&c->denominator,1) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}
/* c->denominator = (max(B,D)/gcd) * (min(B,D))*/
if( (e = mp_div(mp_max(&a->denominator, &b->denominator),&gcd,&tmp,NULL) )
!= MP_OKAY){
return e;
}
if( (e = mp_mul(&tmp, mp_min(&a->denominator, &b->denominator),
&c->denominator) ) != MP_OKAY){
return e;
}
c->sign = (&c->numerator)->sign;
if( (e = mpq_reduce(c) ) != MP_OKAY){
return e;
}
(&a->numerator)->sign = MP_ZPOS;
(&b->numerator)->sign = MP_ZPOS;
mp_clear_multi(&gcd,&tmp,NULL);
return MP_OKAY;
}

You remember it now? 😉
The two functions min() and max() respectively are additions of me to libtommath and are merged into my fork but they are short, so here they are.

/* returns a if equal, min value otherwise */
mp_int *mp_min(mp_int *a, mp_int *b){
int r;
r = mp_cmp(a,b);
if(r == MP_EQ ||  r == MP_LT){
return a;
}
return b;
}
/* returns a if equal, max value otherwise */
mp_int *mp_max(mp_int *a, mp_int *b){
int r;
r = mp_cmp(a,b);
if(r == MP_EQ || r == MP_GT ){
return a;
}
return b;
}

Subtraction, on the other side, is much simpler if we use addition instead, relying on $a- b = -a +b = a + (-b)$. We do not use deep copies of the arguments of the functions so we need not to forget to reverse every modifications we do to them as we have done so, although not commented, above.

int mpq_sub(mp_rat *a, mp_rat *b, mp_rat *c){
int e;
if(mp_iszero(&a->numerator)){
if( (e = mpq_neg(b) ) != MP_OKAY){
return e;
}
if( (e = mpq_copy(b,c) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}
if(mp_iszero(&b->numerator)){
if( (e = mpq_copy(a,c) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}

if(mpq_isinteger(a) && mpq_isinteger(a) ){
(&a->numerator)->sign = a->sign;
(&b->numerator)->sign = b->sign;
if( (e = mp_sub(&a->numerator,&b->numerator,&c->numerator) ) != MP_OKAY){
return e;
}
c->sign = (&c->numerator)->sign;
(&c->numerator)->sign = MP_ZPOS;
if( (e = mp_set_int(&c->denominator,1) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}

/* negate(b) */
if( (e = mpq_neg(b) ) != MP_OKAY){
return e;
}
/* c = a+ (-b) */
if( (e = mpq_add(a,b,c) ) != MP_OKAY){
return e;
}
/* reverse negation */
if( (e = mpq_neg(b) ) != MP_OKAY){
return e;
}
/* No reducing here, mpq_add() did it already */
return MP_OKAY;
}

Multiplication is even simpler

int mpq_mul(mp_rat *a, mp_rat *b, mp_rat *c){
int e;
if(mp_iszero(&a->numerator) || mp_iszero(&b->numerator)){
if( (e = mpq_set_int(c,0,1) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}
if( (e = mp_mul(&a->numerator,&b->numerator,&c->numerator) ) != MP_OKAY){
return e;
}
if( (e = mp_mul(&a->denominator,&b->denominator,&c->denominator) ) != MP_OKAY){
return e;
}
if( (e = mpq_reduce(c) ) != MP_OKAY){
return e;
}
c->sign = (a->sign != b->sign)?MP_NEG:MP_ZPOS;
return MP_OKAY;
}

And the simplest is, believe it or not, division.

int mpq_div(mp_rat *a, mp_rat *b, mp_rat *c){
int e;
if(mp_iszero(&a->numerator)){
if( (e = mpq_set_int(c,0,1) ) != MP_OKAY){
return e;
}
return MP_OKAY;
}
if(mp_iszero(&b->numerator)){
return MPQ_DIVISION_BY_ZERO;
}
if( (e = mpq_inverse(b) ) != MP_OKAY){
return e;
}
if( (e = mpq_mul(a,b,c) ) != MP_OKAY){
return e;
}
if( (e = mpq_inverse(b) ) != MP_OKAY){
return e;
}
if( (e = mpq_reduce(c) ) != MP_OKAY){
return e;
}
c->sign = (a->sign != b->sign)?MP_NEG:MP_ZPOS;
return MP_OKAY;
}

What I seem to have forgotten in the Tools section (don’t look for it, there is no explicitly named Tools section) is a method to compare two rationals. If comparing two big numbers is already a computationally expensive method if none of the shortcuts hit, comparing two rationals is even more so. The only way to do it for two fully reduced fractions is by way of the lowest common denominator, although the full algorithm is not needed, only the multiplication part.

Given two rationals $Q_1 = \frac{p_1}{q_1}$ and $Q_2 = \frac{p_2}{q_2}$ with $p_1 \nmid q_1$ and $p_2 \nmid q_2$ with a difference in magnitude $Q_1 - Q_2 = Q_m$ the sign of the difference ${\mathop{\mathrm{sgn}}}\; Q_m$ can be evaluated by multiplying out $P_1 = p_1q_2$ and $P_2 = p_2q_1$ and comparing $P_1$ and $P_2$.

It is should be obvious from the above that the multiplications can get quite costly if the involved parts of the fractions are very large!

Besides of the normal shortcuts by looking at the individual signs we can also use the fact that we are only interested in the magnitude of the results and not the exact values. If the magnitudes of the intermediary results are large enough in difference we don’t need to do the expensive multiplications.
I am using the number of digits here but one can go more fine-grained with the indices of the highest bits set.

int mpq_cmp(mp_rat *a, mp_rat *b){
mp_int p1,p2;
int e,wnuma,wdena,wnumb,wdenb,s1,s2;

/* a == b, actually the very same number */
if(a == b){
return MP_EQ;
}
/* a sign == MP_NEG && b->sign == MP_ZPOS){
return MP_LT;
}
/* a > b */
if(a->sign == MP_ZPOS && b->sign == MP_NEG){
return MP_GT;
}
/* some more shortcuts */
if(mpq_iszero(a)){
/* a == b (both zero) */
if(mpq_iszero(b)){
return MP_EQ;
}
/* a sign == MP_ZPOS){
return MP_LT;
}
/* a > b */
return MP_GT;
}
if(mpq_iszero(b)){
/* a > b */
if(a->sign == MP_ZPOS){
return MP_GT;
}
/* a numerator,&b->numerator);
}

/*
We can use a rough approximation of the result because we are only
interested in the magnitudes.
*/
wnuma = (&a->numerator)->used;
wdena = (&a->denominator)->used;
wnumb = (&b->numerator)->used;
wdenb = (&b->denominator)->used;
s1 = wnuma + wdenb;
s2 = wnumb + wdena;
if(s1 < s2 - 1){
return ( (a->sign == MP_NEG)?MP_ZPOS: MP_NEG);
}
if(s2 < s1 - 1){
return ( (a->sign == MP_NEG)?MP_NEG: MP_ZPOS);
}
/* This is gonna be brutal, so don't use cmp() frivolously! */
mp_init_multi(&p1,&p2,NULL);
if( (e = mp_mul(&a->numerator,&b->denominator,&p1)) != MP_OKAY){
return e;
}
if( (e = mp_mul(&b->numerator,&a->denominator,&p2)) != MP_OKAY){
mp_clear(&p2);
return e;
}

e = mp_cmp(&p1,&p2);
mp_clear_multi(&p1,&p2,NULL);
if(a->sign == MP_NEG){
return ( (e == MP_ZPOS)?MP_NEG: MP_ZPOS );
}
return e;
}

To test what we did above and if I didn’t any of my infamous C&P errors:

int main (int argc, char **argv){

mp_int num1,den1,num2,den2;
mp_rat q1,q2,q3,q4;

if (argc < 5) {
fprintf(stderr, "usage: %s integer integer integer integer \n", argv[0]);
exit(EXIT_FAILURE);
}

mp_init_multi(&num1,&den1,&num2,&den2,NULL);

mp_get_str(argv[1], &num1, 10);
mp_get_str(argv[2], &den1, 10);
mp_get_str(argv[3], &num2, 10);
mp_get_str(argv[4], &den2, 10);

printf("Numerator   1: ");mp_print(&num1);puts("");
printf("Denominator 1: ");mp_print(&den1);puts("");
printf("Numerator   2: ");mp_print(&num2);puts("");
printf("Denominator 2: ");mp_print(&den2);puts("");

mpq_init_multi(&q1,&q2,&q3,&q4,NULL);puts("000");

mpq_set(&q1,&num1,&den1);puts("111");
printf("Rational1: ");mpq_print(&q1);puts("");
mpq_set(&q2,&num2,&den2);puts("222");
printf("Rational2: ");mpq_print(&q2);puts("");

mpq_add(&q1,&q2,&q3);;
printf("R1 + R2 = ");mpq_print(&q3);puts("");
mpq_sub(&q1,&q2,&q3);
printf("R1 - R2 = ");mpq_print(&q3);puts("");
mpq_mul(&q1,&q2,&q3);
printf("R1 * R2 = ");mpq_print(&q3);puts("");
mpq_div(&q1,&q2,&q3);
printf("R1 / R2 = ");mpq_print(&q3);puts("");
printf("cmp(R1, R2) = %d\n",mpq_cmp(&q1,&q2));
printf("cmp(R2, R1) = %d\n",mpq_cmp(&q2,&q1));
printf("cmp(R1, R1) = %d\n",mpq_cmp(&q1,&q1));
exit(EXIT_SUCCESS);
}

The function mp_get_str() is a simple wrapper:

int mp_get_str(const char *str, mp_int * a, int radix)
{
return mp_read_radix(a, str, radix);
}

Seems to work here. Makes me wonder the most 😉
Some of the more interesting things like e.g.: exponentiation are left as an an exercise…no, just kidding, they are work in progress and not ready yet.

# Notes to self: GIT

Mnemonics are called mnemonics because they are hard to memorize.

The following is more or less just a reminder for myself. A good online reference is Git Reference for the basics.

To add one file (assuming /bin/sh and common GIT naming):

git add a_file.c
git commit -m "Added a_file.c with the miraculous function one_file()"

To add more than one file, commit them all and give all of them the same comment. I used this to commit the three files necessary for the most basic build.

git add tommath.h tommath_class.h makefile
git commit -m "Addition of a_file.c with function one_file()"

Push all of them to remote GIT repository

git fetch upstream
git push -v  origin master

the line git fetch upstream checks the original for intermediate changes that might ruin your daymake what you did obsoleteintersperse with what you have done. If something happened (not in libtommath, last change is a long time ago) you need to do git merge upstream/master where master is the main branch of upstream which is the remote original one. The full command merges that branch with your local branch. Uhm…you can find out the aliases (they are not fixed, the above are just the commonly used ones!) with git remote -v . Here are my results:

czurnieden@home:~/GITHUB/libtommath\$ git remote -v
origin	https://github.com/czurnieden/libtommath.git (fetch)
origin	https://github.com/czurnieden/libtommath.git (push)
upstream	https://github.com/libtom/libtommath.git (fetch)
upstream	https://github.com/libtom/libtommath.git (push)

To find out what happened in the meantime before fetching

git log upstream/master ^master

If unsure which branch you are hacking right now

git branch

Say:”Ooops!” and change to branch testing which you have made with git branch testing

git checkout testing

Check if you haven’t forgoten any file to commit:

git diff-tree --no-commit-id --name-only HEAD~20..HEAD

The number 20 is the maximum number git checks, so if you commited more, increment that number.

# On the numerical evaluation of factorials IX

## Subfactorial

The family of the factorials is larger than the mere number of permutations. On of the many offsprings is the subfactorial or derangement.
The normal implementation, using only integers can be done with a recursion or, to save some memory and leave the stack where it is, in an iterative way.
The recursion is quite complicated with libtommath (more probable is the reason that I am just too stupid to get it right 😉 )

static mp_int kill_the_stack;
static int subfactorial_error = MP_OKAY;
static mp_int * __subfactorialrecursive(unsigned long n){
mp_digit sign;
int e;
if(n==0){
if ( ( e =  mp_set_int(&kill_the_stack,1) ) != MP_OKAY){
subfactorial_error = e;
}
return &kill_the_stack;
}
if(n==1){
if ( ( e =  mp_set_int(&kill_the_stack,0) ) != MP_OKAY){
subfactorial_error = e;
}
return &kill_the_stack;
}
if(n==2){
if ( ( e =  mp_set_int(&kill_the_stack,1) ) != MP_OKAY){
subfactorial_error = e;
}
return &kill_the_stack;
}
if ( ( e =  mp_mul_d(__subfactorialrecursive(n-1),n,&kill_the_stack) ) != MP_OKAY){
subfactorial_error = e;
}
sign = (n&01)?-1:1;
if ( ( e =  mp_add_d(&kill_the_stack,sign,&kill_the_stack) ) != MP_OKAY){
subfactorial_error = e;
}
return &kill_the_stack;
}
static int subfactorial_recursive(unsigned long n, mp_int *result){
mp_init(&kill_the_stack);
subfactorial_error = MP_OKAY;
__subfactorialrecursive(n);
mp_copy(&kill_the_stack,result);
mp_clear(&kill_the_stack);
return subfactorial_error;
}

I left all o the error checks in this time to show one of the problems with recursion: making clear that an error occured. The deep copy in subfactorial_recursive() is needed to be able to clear the globali-ish variable kill_the_stack. Keeping would cause a small but significant memory-leak.

The iterative version is more…uhm…straight forward:

static int subfactorial_iterative(unsigned long n, mp_int *result){
mp_int temp1, temp2;
mp_digit k;
int e;

if(n==0){
return mp_set_int(result,1);
}
if(n==1){
return mp_set_int(result,0);
}
if(n==2){
return mp_set_int(result,1);
}
if ( ( e = mp_init_multi(&temp1, &temp2,NULL) ) != MP_OKAY){
return e;
}
if ( ( e = mp_set_int(&temp1, 0) ) != MP_OKAY){
mp_clear_multi(&temp1, &temp2,NULL);
return e;
}
if ( ( e = mp_set_int(result,1) ) != MP_OKAY){
mp_clear_multi(&temp1, &temp2,NULL);
return e;
}
for(k=3;k<=n;k++){
mp_exch ( &temp1,&temp2 );
mp_exch ( result, &temp1);
if ( ( e = mp_add(&temp1, &temp2,result))  != MP_OKAY){
mp_clear_multi(&temp1, &temp2,NULL);
return e;
}
if ( ( e = mp_mul_d(result,(k-1),result))  != MP_OKAY){
mp_clear_multi(&temp1, &temp2,NULL);
return e;
}
}
mp_clear_multi(&temp1, &temp2,NULL);
return MP_OKAY;
}

The use of libtommath’s mp_set_int() needs a restriction to the size of n. The type allowed is an unsigned long; I restrict further to the actual bit-size of mp_digit which is probably way too conservative, especially when those digits are smaller than the default of 28.
The subfactorial of 2^28 is about 1.7862e2,146,019,442, that’s a number with over two billion decimal digits. One needs a different method to calculate it exactly. This one is good up to about 100,000 which already lasts a minute with this old computer I have here, more with the iterative version.

#ifndef MP_DIGIT_SIZE
#define MP_DIGIT_SIZE (1<<DIGIT_BIT)
#endif
int mp_subfactorial(unsigned long n, mp_int *result){
if(n > MP_DIGIT_SIZE ){
return MP_RANGE;
}
#ifdef SAVE_SOME_MEMORY
return subfactorial_iterative(n, result);
#else
return  subfactorial_recursive(n, result);
#endif
}

The code above has been merged to my fork of libtommath over at github.

# On the numerical evaluation of factorials VIII

The code for the FFT multiplication is working now and had been uploaded to my fork of libtommath. It is still limited to the 28 bit digits, The Schönberg-Strassen algorithm should be able to let us get rid of that ancient limit bt it is still in the works.
BTW: for some reasons unknown to me, not all commits of the original libtommath are in my fork of it. One of those commits fixed a bug that prevented GCC from compiling libtommath, I put it in there manually wich is not the very best method to do it, but as mentioned above: I still don’t know what I made wrong when I forked libtommath.

The cut-offs and limits may need some adjusting for other processors and architectures then my poor little 1GHz AMD Duron.

Now that FFT multiplication is working we can get back to the factorial. The last state was about 12 minutes for the factorial of one million if I remember it correctly. That is definitely way too much even for such an old processor like the one in my box. To avoid linking to the last post…oops, never mind. However, the last code was:

int factorial(unsigned long n, mp_int *result){
unsigned long *p_list;
unsigned long *exp_list;
unsigned long p, i,cut;
long bit;
int shift;
mp_int temp;

unsigned long  r=0;
p_list = fill_prime_list(3, n+1, &r);
exp_list = malloc(sizeof(unsigned long)* (r+1));
if(exp_list == NULL){ return MP_MEM;}

mp_set_int(result,1);
shift = (int)prime_divisors(n,2LU);
cut = n / 2 ;
for(p = 0; p < r; p++) {
if(p_list[p] > cut) break;
exp_list[p] = prime_divisors(n,p_list[p]);
}
bit = highbit(exp_list[0]);
for( ; bit>=0; bit--) {
mp_sqr(result, result);
for(i=0; i<p; i++) {
if ((exp_list[i] & (1 << bit)) != 0) {
mp_mul_d(result,(mp_digit)p_list[i], result);
}
}
}
for(; p < r; p++){
mp_mul_d(result,(mp_digit)p_list[p], result);
}
mp_mul_2d (result, shift, result);
return MP_OKAY;
}

Anything open for optimization? What’s about the last loop over the primes with an exponent of 1 (one)? The loop adds all of the last primes to the result. What’s about summing them in a temporary variable and multiply the result with the temporary variable instead?

#define MP_FACTORIAL_BORW_PRIMORIAL_CUTOFF 200000
int factorial(unsigned long n, mp_int *result){
unsigned long *p_list;
unsigned long *exp_list;
unsigned long p, i,cut;
long bit;
int shift;
mp_int temp;

unsigned long  r=0;
p_list = fill_prime_list(3, n+1, &r);
exp_list = malloc(sizeof(unsigned long)* (r+1));
if(exp_list == NULL){ return MP_MEM;}

mp_set_int(result,1);

shift = (int)prime_divisors(n,2LU);

cut = n / 2 ;

for(p = 0; p < r; p++) {
if(p_list[p] > cut) break;
exp_list[p] = prime_divisors(n,p_list[p]);
}
mp_init(&temp);
bit = (int)highbit(exp_list[0]);
for( ; bit>=0; bit--) {
mp_sqr(result, result);
for(i=0; i<p; i++) {
if ((exp_list[i] & (1 << bit)) != 0) {
mp_mul_d(result,(mp_digit)p_list[i], result);
}
}
}
if(n<MP_FACTORIAL_BORW_PRIMORIAL_CUTOFF){
for(; p < r; p++){
mp_mul_d(result,(mp_digit)p_list[p], result);
}
}
else {
mp_init(&temp);
mp_set_int(&temp,1);
for(; p < r; p++){
mp_mul_d(&temp,(mp_digit)p_list[p], &temp);
}
mp_mul(result,&temp,result);
}
mp_mul_2d (result, shift, result);

return MP_OKAY;
}

After some tests I found a cut-off point at 200,000. The factorial of one million is now available in 6 (six) minutes. Half of the time than before!
Can we do more? Well, let’s see: the primes are consecutive, increasing and already in an array. That should enable us to use the binary splitting method. It is a bit more complicated here because we get into big-number territory in parts of it but that is good, because we can use some of the optimizations for the multiplication of large numbers. The size, or better: magnitude, number of digits, of the multiplicands is also the same, ideal for the implemented Toom-Cook 2 and 3 and FFT.

static unsigned long *primelist;
static   mp_int p_temp;

mp_int *primorial__lowlevel(unsigned long a, unsigned long b ,
unsigned long p)
{
unsigned long  c;
unsigned long long prod;

mp_init(&p_temp);

mp_int ta,tb;

if( b == a){
mp_set_int(&p_temp,p);
return  &(p_temp);
}
if( b-a > 1){
c= (b + a) >> 1;
mp_init_multi(&ta,&tb,NULL);
mp_copy(primorial__lowlevel( a   , c , primelist[a] ),&ta);
mp_copy(primorial__lowlevel( c+1 , b , primelist[b] ),&tb);
if(ta.used == 1 && tb.used == 1){
if(ta.dp[0] <= MP_DIGIT_HALF || tb.dp[0] <= MP_DIGIT_HALF ){
prod = (mp_word)ta.dp[0] * tb.dp[0];
if( prod< MP_DIGIT_SIZE){
mp_set_int(&p_temp,(mp_digit)prod);
return &p_temp;
}
}
}
mp_mul(&ta,&tb,&p_temp) ;
mp_clear_multi(&ta,&tb,NULL);
return &p_temp;
}
if(   (primelist[a]) >=  (MP_DIGIT_HALF )
|| (primelist[b]) >=  (MP_DIGIT_HALF )   ){
mp_set_int( &p_temp,primelist[a]);
mp_mul_d(&p_temp,primelist[b],&p_temp);
}
else{
c = primelist[a]*primelist[b];
mp_set_int(&p_temp,c);
}
return &p_temp;
}
void primorial(unsigned long a, unsigned long b, mp_int *result){
unsigned long  r=0;
primelist = fill_prime_list(a+1, b+1, &r);
*result = *primorial__lowlevel(0,r-1,1);
}

It even needs some quasi global (static) variables which is rarely a good idea but avoiding them wold clobber the code even more. However, let’s put it in to see if it helps.

#define MP_FACTORIAL_BORW_PRIMORIAL_CUTOFF 200000
int factorial(unsigned long n, mp_int *result){
unsigned long *p_list;
unsigned long *exp_list;
unsigned long p, i,cut;
long bit;
int shift;
mp_int temp;

unsigned long  r=0;
p_list = fill_prime_list(3, n+1, &r);
exp_list = malloc(sizeof(unsigned long)* (r+1));
if(exp_list == NULL){ return MP_MEM;}

mp_set_int(result,1);

shift = (int)prime_divisors(n,2LU);

cut = n / 2 ;

for(p = 0; p < r; p++) {
if(p_list[p] > cut) break;
exp_list[p] = prime_divisors(n,p_list[p]);
}
mp_init(&temp);
bit = (int)highbit(exp_list[0]);
for( ; bit>=0; bit--) {
mp_sqr(result, result);
for(i=0; i<p; i++) {
if ((exp_list[i] & (1 << bit)) != 0) {
mp_mul_d(result,(mp_digit)p_list[i], result);
}
}
}
if(n<MP_FACTORIAL_BORW_PRIMORIAL_CUTOFF){
for(; p < r; p++){
mp_mul_d(result,(mp_digit)p_list[p], result);
}
}
else {
mp_init(&temp);
primorial(cut,n,&temp);
mp_mul(result,&temp,result);
}
mp_mul_2d (result, shift, result);

return MP_OKAY;
}

Result: roughly the same. The theoretical complexity is sometimes different from the one of the implementation. The main culprit here is most probably the deep copy nexessary in our implementation, but it will win for much larger numbers.
But there is another one of the same in the main loop. Let’s replace that with our first algorithm:

#define MP_FACTORIAL_BORW_LOOP_CUTOFF 500000
#define MP_FACTORIAL_BORW_PRIMORIAL_CUTOFF 200000

int factorial(unsigned long n, mp_int *result){
unsigned long *p_list;
unsigned long *exp_list;
unsigned long p, i,cut;
long bit;
int shift;
mp_int temp;

unsigned long  r=0;
p_list = fill_prime_list(3, n+1, &r);
exp_list = malloc(sizeof(unsigned long)* (r+1));
if(exp_list == NULL){ return MP_MEM;}

mp_set_int(result,1);

shift = (int)prime_divisors(n,2LU);

cut = n / 2 ;

for(p = 0; p < r; p++) {
if(p_list[p] > cut) break;
exp_list[p] = prime_divisors(n,p_list[p]);
}
mp_init(&temp);
bit = (int)highbit(exp_list[0]);
if(n < MP_FACTORIAL_BORW_LOOP_CUTOFF){
for( ; bit>=0; bit--) {
mp_sqr(result, result);
for(i=0; i<p; i++) {
if ((exp_list[i] & (1 << bit)) != 0) {
mp_mul_d(result,(mp_digit)p_list[i], result);
}
}
}
}
else{
for( ; bit>=0; bit--) {
mp_sqr(result, result);
mp_set_int(&temp,1);
for(i=0; i<p; i++) {
if ((exp_list[i] & (1 << bit)) != 0) {
mp_mul_d(&temp,(mp_digit)p_list[i], &temp);
}
}
mp_mul(&temp,result,result);
}
}
if(n<MP_FACTORIAL_BORW_PRIMORIAL_CUTOFF){
for(; p < r; p++){
mp_mul_d(result,(mp_digit)p_list[p], result);
}
}
else {
primorial(cut,n,&temp);
mp_mul(result,&temp,result);
}
mp_mul_2d (result, shift, result);
return MP_OKAY;
}

After some test runs the cut-off found is something around half a million.
Run time for the factorial of one million: a mere 25 (twenty-five) seconds! Oh, and with the linear multiplication instead of the binary splitting of the final primorial it is about 30 seconds. Now we have a significant difference. It is the same difference as before but five seconds hidden in six minutes are hard to detect.

I tried ten million but with only half a gig of RAM and some memory hungry apps like the Firefox I use to write this post, it starts to swap like hell quite fast and needs about just a bit less than 13 minutes. Sounds a lot but the factorial of ten million is a large number, about 1.2024e65657052. Yes, a number with 65,657,053 decimal digits.
I tried the ten million factorial without a running Firefox and got about 10 seconds saved. Not much, even if the sound the harddrive made gave the impression of ten minutes instead 😉
The factorial of ten million is large enough to allow us a better test if the binary splitting is indeed faster. And it is: with the linear multiplication instead of the binary splitting method the run time nearly doubles and passed 24 minutes. So the theory is correct at least 😉

BTW: I just saw that libtommath shares its birthday with me 😉


text:
### For REFEREES

 A high order difference method for fractional sub-diffusion equations with the spatially variable coefficients under periodic boundary conditions Huiqin Zhang,Yuxiang Liang,Zhibo Wang Keywords:Fractional diffusion equation, Fourier method, Variable coefficients, Stability Abstract: In this paper, we propose a difference scheme with global convergence order $O(\tau^{2}+h^4)$ for a class of the Caputo fractional equation. The difficulty caused by the spatially variable coefficients is successfully handled. The unique solvability, stability and convergence of the finite difference scheme are proved by use of the Fourier method. The obtained theoretical results are supported by numerical experiments.

format: plain


text:
# Publications

## Preprints

Y. Rollin, Polyhedral approximation by Lagrangian and isotropic tori, https://arxiv.org/abs/2012.05777
• Abstract: We prove that every smoothly immersed 2-torus of $\mathbb{R}^4$ can be approximated, in the $C^0$-sense, by immersed polyhedral Lagrangian tori. In the case of a smoothly immersed (resp. embedded) Lagrangian torus of $\mathbb{R}^4$, the surface can be approximated in the $C^1$-sense by immersed (resp. embedded) polyhedral Lagrangian tori. Similar statements are proved for isotropic 2-tori of $\mathbb{R}^{2n}$.
F. Jauberteau and Y. Rollin, Polyhedral Lagrangian surfaces and moment map flow, Code attached to the paper coming up soon, https://www.math.sciences.univ-nantes.fr/~rollin/articles/dislagflow.pdf
• Abstract: We show that isotropic polyhedral surfaces of a symplectic vector space V can be understood as the vanishing locus of a moment map as- sociated to a particular Hamiltonian torus action on a finite dimensional symplectic space. This interpretation allows to introduce an associated modified moment map flow, which is an ODE whose limits are essentially the (parametrized) isotropic surfaces. This approach provides an effec- tive method for constructing Lagrangian polyhedral surfaces of $\mathbb{R}^4$ . This paper is provided with a computer program which shows the modified flow as a real-time movie, in the case of a 2-dimensional torus.
F. Jauberteau and Y. Rollin, Piecewise linear symplectomorphism and hyperKähler moment map flow, In progress.

Most of my preprints are available on www.arXiv.org

## Published

F. Jauberteau, Y. Rollin and S. Tapie, Discrete geometry and isotropic surfaces, Mém. Soc. Math. Fr. (N.S.) (161):vii+101, 2019, http://arxiv.org/abs/1802.08712
• Abstract: We consider smooth isotropic immersions from the 2-dimensional torus into R2n, for n≥2. When n=2 the image of such map is an immersed Lagrangian torus of R4. We prove that such isotropic immersions can be approximated by arbitrarily C0-close piecewise linear isotropic maps. If n≥3 the piecewise linear isotropic maps can be chosen so that they are piecewise linear isotropic immersions as well. The proofs are obtained using analogies with an infinite dimensional moment map geometry due to Donaldson. As a byproduct of these considerations, we introduce a numerical flow in finite dimension, whose limit provide, from an experimental perspective, many examples of piecewise linear Lagrangian tori in R4. The DMMF program, which is freely available, is based on the Euler method and shows the evolution equation of discrete surfaces in real time, as a movie.
E. Legendre and Y. Rollin, Hamiltonian stationary Lagrangian fibrations, J. Symplectic Geom. 17(3):753-791, 2019
• Abstract: Hamiltonian stationary Lagrangian submanifolds (HSLAG) are a natural generalization of special Lagrangian manifolds (SLAG). The latter only make sense on Calabi-Yau manifolds whereas the former are defined for any almost Kähler manifold. Special Lagrangians, and, more specificaly, fibrations by special Lagrangians play an important role in the context of the geometric mirror symmetry conjecture. However, these objects are rather scarce in nature. On the contrary, we show that HSLAG submanifolds, or fibrations, arise quite often. Many examples of HSLAG fibrations are provided by toric Kähler geometry. In this paper, we obtain a large class of examples by deforming the toric metrics into non toric almost Kähler metrics, together with HSLAG submanifolds.
V. Apostolov and Y. Rollin, ALE scalar-flat Kähler metrics on non-compact weighted projective spaces, Math. Ann. 367(3-4):1685-1726, 2017, http://arxiv.org/abs/1510.02226
• Abstract: Abstract: We construct new explicit toric scalar-flat Kähler ALE metrics on weighted projective spaces of non-compact type, which we use to obtain smooth extremal Kähler metrics on appropriate resolutions of orbifolds. In particular, we obtain new extremal metrics certain resolutions of weighted projective spaces of compact type.
Y. Rollin, K-stability and parabolic stability, Advances in Mathematics 285:1741-1766, 2015, http://arxiv.org/abs/1303.2332
• Abstract: Parabolic structures with rational weights encode certain iterated blowups of geometrically ruled surfaces. In this paper, we show that the three notions of parabolic polystability, K-polystability and existence of constant scalar curvature Kähler metrics on the iterated blowup are equivalent, for certain polarizations close to the boundary of the Kähler cone.
O. Biquard and Y. Rollin, Smoothing singular extremal Kähler surfaces and minimal Lagrangians, Advances in Mathematics 285:980 - 1024, 2015, http://arxiv.org/abs/1211.6957
• Abstract: We consider smoothings of a complex surface with singularities of class T and no nontrivial holomorphic vector field. Under an hypothesis of non degeneracy of the smoothing at each singular point, we prove that if the singular surface admits an extremal metric, then the smoothings also admit extremal metrics in nearby Kähler classes. In addition, we construct small Lagrangian stationary spheres which represent Lagrangian vanishing cycles for surfaces close to the singular one.
Y. Rollin and C. Tipler, Deformations of extremal toric manifolds, J. Geom. Anal. 24(4):1929-1958, 2014
• Abstract: Let X be a compact toric extremal Kähler manifold. Using the work of Székelyhidi, we provide a combinatorial criterion on the fan describing X to ensure the existence of complex deformations of X that carry extremal metrics. As an example, we find new CSC metrics on 4-points blow-ups of $CP^1\times CP^1$.
Y. Rollin, S.R. Simanca and C. Tipler, Deformation of extremal metrics, complex manifolds and the relative Futaki invariant, Math. Z. 273(1-2):547-568, 2013, http://arxiv.org/abs/1107.0456
• Abstract: Let $(\mathcal {X},\Omega)$ be a closed polarized complex manifold, $g$ be an extremal metric on $\mathcal X$ that represents the Kähler class $\Omega$, and $G$ be a compact connected subgroup of the isometry group $Isom(\mathcal{X},g)$. Assume that the Futaki invariant relative to $G$ is nondegenerate at $g$. Consider a smooth family $(\mathcal{M}\to B)$ of polarized complex deformations of $(\mathcal{X},\Omega)\simeq (\mathcal{M}_0,\Theta_0)$ provided with a holomorphic action of $G$. Then for every $t\in B$ sufficiently small, there exists an $h^{1,1}(\mathcal{X})$-dimensional family of extremal Kähler metrics on $\mathcal{M}_t$ whose Kähler classes are arbitrarily close to $\Theta_t$. We apply this deformation theory to analyze the Mukai-Umemura 3-fold and its complex deformations.
Y. Rollin and M. Singer, Construction of Kähler surfaces with constant scalar curvature, J. Eur. Math. Soc. (JEMS) 11(5):979-997, 2009, math.DG/0412405
• Abstract: We present new constructions of Kaehler metrics with constant scalar curvature on complex surfaces, in particular on certain del Pezzo surfaces. Some higher dimensional examples are provided as well.
O. Biquard and Y. Rollin, Wormholes in ACH Einstein manifolds, Trans. Amer. Math. Soc. 361(4):2021-2046, 2009, math.DG/0609558
• Abstract: We give a new construction of Einstein and Kaehler-Einstein manifolds which are asymptotically complex hyperbolic, inspired by the work of Mazzeo-Pacard in the real hyperbolic case. The idea is to develop a gluing theorem for 1-handle surgery at infinity, which generalizes the Klein construction for the complex hyperbolic metric.
Y. Rollin and M. Singer, Constant scalar curvature Kähler surfaces and parabolic polystability, J. Geom. Anal. 19(1):107-136, 2009, math.DG/0703212
• Abstract: A complex ruled surface admits an iterated blow-up encoded by a parabolic structure with rational weights. Under a condition of parabolic stability, one can construct a Kaehler metric of constant scalar curvature on the blow-up according to math.DG/0412405. We present a generalization of this construction to the case of parabolically polystable ruled surfaces. Thus we can produce numerous examples of Kaehler surfaces of constant scalar curvature with circle or toric symmetry.
T. Mrowka and Y. Rollin, Legendrian knots and monopoles, Algebr. Geom. Topol. 6:1-69 (electronic), 2006, math.DG/0410559, (MR)
• Abstract: We prove a generalization of Bennequin's inequality for Legendrian knots in a 3-dimensional contact manifold, under the assumption that it is the boundary of a 4-dimensional manifold M and the version of Seiberg-Witten invariants introduced by Kronheimer and Mrowka is non-vanishing. The proof requires an excision result for Seiberg-Witten moduli spaces; then, the Bennequin inequality becomes a special case of the adjunction inequality for surfaces lying inside M.
Y. Rollin and M. Singer, Non-minimal scalar-flat Kähler surfaces and parabolic stability, Invent. Math. 162(2):235-270, 2005, math.DG/0404423, (MR)
• Abstract: A new construction is presented of scalar-flat Kaehler metrics on non-minimal ruled surfaces. The method is based on the resolution of singularities of orbifold ruled surfaces which are closely related to rank-2 parabolically stable holomorphic bundles. This rather general construction is shown also to give new examples of low genus: in particular, it is shown that $\mathbb{CP}^2$ blown up at 10 suitably chosen points, admits a scalar-flat Kaehler metric; this answers a question raised by Claude LeBrun in 1986 in connection with the classification of compact self-dual 4-manifolds.
Y. Rollin, Rigidité d'Einstein du plan hyperbolique complexe, J. Reine Angew. Math. 567:175-213, 2004, math.DG/0112099, (MR)
• Abstract: We prove that every Einstein metric on the unit ball $B^4$ of $\mathbb{C}^2$, asymptotic to the Bergman metric, is equal to it up to a diffeomorphism. We need a solution of Seiberg--Witten equations in this infinite volume setting. Therefore, and more generally, if $M^4$ is a manifold with a CR-boundary at infinity, an adapted spinc-structure which has a non zero Kronheimer--Mrowka invariant and an asymptotically complex hyperbolic Einstein metric, we produce a solution of Seiberg--Witten equations with an strong exponential decay property.
Y. Rollin, Surfaces kählériennes de volume fini et équations de Seiberg-Witten, Bull. Soc. Math. France 130(3):409-456, 2002, math.DG/0106077, (MR)
• Abstract: Let M=P(E) be a ruled surface. We introduce metrics of finite volume on M whose singularities are parametrized by a parabolic structure over E. Then, we generalise results of Burns--de Bartolomeis and LeBrun, by showing that the existence of a singular Kahler metric of finite volume and constant non positive scalar curvature on M is equivalent to the parabolic polystability of E; moreover these metrics all come from finite volume quotients of $\mathbb{H}^2 \times \mathbb{CP}^1$. In order to prove the theorem, we must produce a solution of Seiberg-Witten equations for a singular metric g. We use orbifold compactifications $\overline M$ on which we approximate g by a sequence of smooth metrics; the desired solution for g is obtained as the limit of a sequence of Seiberg-Witten solutions for these smooth metrics.
Y. Rollin, Einstein rigidity of the complex hyperbolic plane, C. R. Math. Acad. Sci. Paris 334(8):671-676, 2002, http://yann.rollin.free.fr/crastein.pdf, (MR).

## Thesis

Y. Rollin, Métriques kählériennes de volume fini, uniformisation des surfaces complexes réglées et équations de Seiberg-Witten, PhD thesis, École Polytechnique, jan 2001,
http://yann.rollin.free.fr/these_rollin.pdf.

Y. Rollin, Topologie de contact et géométries à courbure spéciale via la théorie de jauge, mémoire d'habilitation , 2007, http://yann.rollin.free.fr/hab.pdf.

## Other papers

Y. Rollin, A remark on the Herzlich volume of asymptotically complex hyperbolic Einstein manifolds, 2008, http://arxiv.org/abs/0802.2474.


text:
# compute temp/deform command

## Syntax

compute ID group-ID temp/deform

• ID, group-ID are documented in compute command
• temp/deform = style name of this compute command

## Examples

compute myTemp all temp/deform

## Description

Define a computation that calculates the temperature of a group of atoms, after subtracting out a streaming velocity induced by the simulation box changing size and/or shape, for example in a non-equilibrium MD (NEMD) simulation. The size/shape change is induced by use of the fix deform command. A compute of this style is created by the fix nvt/sllod command to compute the thermal temperature of atoms for thermostatting purposes. A compute of this style can also be used by any command that computes a temperature, e.g. thermo_modify, fix temp/rescale, fix npt, etc.

The deformation fix changes the box size and/or shape over time, so each atom in the simulation box can be thought of as having a “streaming” velocity. For example, if the box is being sheared in x, relative to y, then atoms at the bottom of the box (low y) have a small x velocity, while atoms at the top of the box (hi y) have a large x velocity. This position-dependent streaming velocity is subtracted from each atom’s actual velocity to yield a thermal velocity which is used to compute the temperature.

Note

Fix deform has an option for remapping either atom coordinates or velocities to the changing simulation box. When using this compute in conjunction with a deforming box, fix deform should NOT remap atom positions, but rather should let atoms respond to the changing box by adjusting their own velocities (or let fix deform remap the atom velocities, see it’s remap option). If fix deform does remap atom positions, then they appear to move with the box but their velocity is not changed, and thus they do NOT have the streaming velocity assumed by this compute. LAMMPS will warn you if fix deform is defined and its remap setting is not consistent with this compute.

After the streaming velocity has been subtracted from each atom, the temperature is calculated by the formula KE = dim/2 N k T, where KE = total kinetic energy of the group of atoms (sum of 1/2 m v^2), dim = 2 or 3 = dimensionality of the simulation, N = number of atoms in the group, k = Boltzmann constant, and T = temperature. Note that v in the kinetic energy formula is the atom’s thermal velocity.

A kinetic energy tensor, stored as a 6-element vector, is also calculated by this compute for use in the computation of a pressure tensor. The formula for the components of the tensor is the same as the above formula, except that v^2 is replaced by vx*vy for the xy component, etc. The 6 components of the vector are ordered xx, yy, zz, xy, xz, yz.

The number of atoms contributing to the temperature is assumed to be constant for the duration of the run; use the dynamic option of the compute_modify command if this is not the case.

The removal of the box deformation velocity component by this fix is essentially computing the temperature after a “bias” has been removed from the velocity of the atoms. If this compute is used with a fix command that performs thermostatting then this bias will be subtracted from each atom, thermostatting of the remaining thermal velocity will be performed, and the bias will be added back in. Thermostatting fixes that work in this way include fix nvt, fix temp/rescale, fix temp/berendsen, and fix langevin.

Note

The temperature calculated by this compute is only accurate if the atoms are indeed moving with a stream velocity profile that matches the box deformation. If not, then the compute will subtract off an incorrect stream velocity, yielding a bogus thermal temperature. You should NOT assume that your atoms are streaming at the same rate the box is deforming. Rather, you should monitor their velocity profile, e.g. via the fix ave/chunk command. And you can compare the results of this compute to compute temp/profile, which actually calculates the stream profile before subtracting it. If the two computes do not give roughly the same temperature, then your atoms are not streaming consistent with the box deformation. See the fix deform command for more details on ways to get atoms to stream consistently with the box deformation.

This compute subtracts out degrees-of-freedom due to fixes that constrain molecular motion, such as fix shake and fix rigid. This means the temperature of groups of atoms that include these constraints will be computed correctly. If needed, the subtracted degrees-of-freedom can be altered using the extra option of the compute_modify command.

See this howto section of the manual for a discussion of different ways to compute temperature and perform thermostatting.

Output info:

This compute calculates a global scalar (the temperature) and a global vector of length 6 (KE tensor), which can be accessed by indices 1-6. These values can be used by any command that uses global scalar or vector values from a compute as input. See this section for an overview of LAMMPS output options.

The scalar value calculated by this compute is “intensive”. The vector values are “extensive”.

The scalar value will be in temperature units. The vector values will be in energy units.

none


text:
Press J to jump to the feed. Press question mark to learn the rest of the keyboard shortcuts
Community Details

687k

Subscribers

789

Online

Welcome to r/math!

This subreddit is for discussion of mathematical links and questions. Please read the FAQ and the rules below before posting.

If you're asking for help understanding something mathematical, post in the Simple Question thread or /r/learnmath. This includes reference requests - also see our lists of recommended books and free online resources. Here is a more recent thread with book recommendations.

r/math Rules
1.
No homework problems
2.
Stay on-topic
3.
Be excellent to each other
4.
No low-effort image posts
5.
No career or education related questions

Everything about X - every Wednesday

What Are You Working On? - posted Mondays

Career and Education Q&A - Every other Thursday

Simple Questions - Posted Fridays

Using LaTeX

To view LaTeX on reddit, install one of the following:

MathJax userscript (install Greasemonkey or Tampermonkey first)

[; e^{\pi i} + 1 = 0 ;]

Post the equation above like this:

[; e^{\pi i}+1=0 ;]

Useful Symbols

Basic Math Symbols

≠ ± ∓ ÷ × ∙ – √ ‰ ⊗ ⊕ ⊖ ⊘ ⊙ ≤ ≥ ≦ ≧ ≨ ≩ ≺ ≻ ≼ ≽ ⊏ ⊐ ⊑ ⊒ ² ³ °

Geometry Symbols

∠ ∟ ° ≅ ~ ‖ ⟂ ⫛

Algebra Symbols

≡ ≜ ≈ ∝ ∞ ≪ ≫ ⌊⌋ ⌈⌉ ∘∏ ∐ ∑ ⋀ ⋁ ⋂ ⋃ ⨀ ⨁ ⨂ 𝖕 𝖖 𝖗

Set Theory Symbols

∅ ∖ ∁ ↦ ↣ ∩ ∪ ⊆ ⊂ ⊄ ⊊ ⊇ ⊃ ⊅ ⊋ ⊖ ∈ ∉ ∋ ∌ ℕ ℤ ℚ ℝ ℂ ℵ ℶ ℷ ℸ 𝓟

Logic Symbols

¬ ∨ ∧ ⊕ → ← ⇒ ⇐ ↔ ⇔ ∀ ∃ ∄ ∴ ∵ ⊤ ⊥ ⊢ ⊨ ⫤ ⊣

Calculus and Analysis Symbols

∫ ∬ ∭ ∮ ∯ ∰ ∇ ∆ δ ∂ ℱ ℒ ℓ

Greek Letters

𝛢𝛼 𝛣𝛽 𝛤𝛾 𝛥𝛿 𝛦𝜀𝜖 𝛧𝜁 𝛨𝜂 𝛩𝜃𝜗 𝛪𝜄 𝛫𝜅 𝛬𝜆 𝛭𝜇 𝛮𝜈 𝛯𝜉 𝛰𝜊 𝛱𝜋 𝛲𝜌 𝛴𝜎 𝛵𝜏 𝛶𝜐 𝛷𝜙𝜑 𝛸𝜒 𝛹𝜓 𝛺𝜔

Other math subreddits
r/learnmath

96,166 subscribers

r/mathbooks

9,666 subscribers

r/cheatatmathhomework

32,460 subscribers

r/matheducation

13,949 subscribers

r/CasualMath

8,971 subscribers

r/puremathematics

10,407 subscribers

r/mathpics

17,513 subscribers

r/mathriddles

8,197 subscribers

Related subreddits
r/Mathematica

4,469 subscribers

r/matlab

17,788 subscribers

r/sagemath

685 subscribers

r/actuary

12,456 subscribers

r/algorithms

43,565 subscribers

r/compsci

437,114 subscribers

r/interdisciplinary

1,699 subscribers

r/statistics

66,984 subscribers

Moderators
u/yesmanapple
Geometry/Topology
u/inherentlyawesome
Algebraic Topology
u/amdpox
Geometric Analysis
u/AngelTC
Algebraic Geometry
u/dogdiarrhea
Dynamical Systems
u/functor7
Number Theory
16
Posted by7 days ago

## Just watched the man who knew infinity and then proceeded to research the dude and find that Ramanujan said that the number e^pi*sqrt(163) was something special. Can someone explain what’s the deal with this number?

77% Upvoted
Sort by
level 1
17 points · 7 days ago

Hmm. I don't know your level so at ELI18 level it's very closely approximated by an integer which initially looks surprising given that it looks like an arbitrary combination of transcendental and irrational constants and is itself transcendental. This is basically a consequence of the properties of a certain special function called the j-invariant and the fact that 163 is the largest value d for which Q(sqrt(-d)) has class number 1 aka unique factorisation. The other such numbers are also approximated closely (but not as closely) by integers.

level 1
Original Poster11 points · 7 days ago

https://youtu.be/DRxAVA6gYMM

In this video I still don’t get what sqrt(-163) is all about if someone could let me know

level 2
25 points · 7 days ago

-163 ist the smallest number for which the associated imaginary-quadratic number field F=Q(sqrt(-163)) has class number h(F)=1. That‘s pretty amazing and leads to a lot of other resulta regarding that number, if one likes numerical properties of some integers. :)

level 3
Original Poster5 points · 6 days ago

Without being a burden, would you mind breaking down your explanation further? ELI5? OR if you could point me to resources to understand your statement above, id really appreciate it

level 4
Number Theory4 points · 6 days ago

Basically it means that the set of all numbers in the form [;\{a+b\omega|a,b\in\mathbb{Z}\};], where [;\omega = \frac{1+\sqrt{-163}}{2};] has "unique prime factorization", similar to how [;\mathbb{Z};] does.

https://en.wikipedia.org/wiki/Unique_factorization_domain

As a simpler example, [;\mathbb{Z}[i] = \{a+b\omega|a,b\in\mathbb{Z}\};], the Gaussian integers, also have this property, but [;\mathbb{Z}[\sqrt{-5}] = \{a+b\sqrt{-5}|a,b\in\mathbb{Z}\};] does not, because (for instance) [;6;] has two distinct factorizations [;6=(2)(3) = (1+\sqrt{-5})(1-\sqrt{-5});].

The reason why this should say something about the number [;e^{\pi\sqrt{163}};] is unfortunately rather complicated. It involves elliptic curves and more specifically complex multiplication and the j-invariant.

level 4
1 point · 4 days ago

Most if the necessary math should be covered in any basic algebraic number theory book, e.g. Neukirch

But I can try to explain in a few words.

Basically, one looks at a number field F (i.e. A finite field extension of the rationals) and their ring of integers, that is, the integral closure of Z in F (i.e. the ring of all elements of F which satisfy an integral condition over Z). Lets denote this ring by O In O we can define the ideal group I(O) of all „fractional“ ideals (that‘s just the normal ideals + their „inverses“) On this group there is an equivalence relation where J1~J2 <=> there is some a from O s.t. J1*aO=J2 (in I(O)) Now we define the class group cl(O) := I(O) / ~ One can show that this group is always finite, and if the class group is the trivial group, we know that O is a PID, therefore a UFD. So we are particulary interested in fields with class number 1 (i.e. fields with |cl(O)|=1) because their rings of integers behave like the original integers.

And in your case, -163 is the smallest number for which the corresponding imaginary-quadratic field Q(x2 +163) has class number 1.

This explanation might be a little sloppy and short, for more information look at some algebraic number theory book :)

level 2
7 points · 7 days ago

It's always a delight when I find an old Numberphile video that I haven't seen yet. Even better when it's as good as this one. Reminds me of my Galois Theory class. I did alright in that class but didn't really retain much other than the flavor of it. It's such a deep topic that grad school doesn't assume that you know it, even though it's an undergrad course, because it's hard to retain it all the first time you learn it. If you're a math major, Galois Theory is must-take course. You won't regret it.

level 1
7 points · 7 days ago

As we all know, epi*sqrt(163) is exactly equal to 262537412640768744. It's used to test your computer's handling of floating-point arithmetic; most computers may give an erroneous result (e.g. 262537412640768743.99999999999925). Another such example is epi-pi, which equals 20.

Anyway, 262537412640768744 is significant because it is the smallest number that cannot nontrivially be written as the sum of three cubes.

It is also the smallest number for which the Lonely Runner Conjecture fails.

It is also one of the edges of the smallest (by volume) perfect cuboid.

All in all, a truly special number.

(yes this is a joke)

level 1
Original Poster1 point · 6 days ago


text:
# Boundary of invertibles in a normed algebra

A student and I are reading the book Introduction to Banach Spaces and Algebras, by Allan, and we're stuck. Exercise 4.5 says:

Let $A$ be a normed algebra with unit sphere $S$. Let $a\in A$. Then $a$ is a topological divisor of 0 if $$\inf\{\|ab\|+\|ba\|:b\in S \}=0.$$ Prove that every element in the frontier of $G(A)$ is a topological divisor of $0$.

Here $G(A)$ is the collection of invertible elements of $A$. I assume that the question really means to say that $A$ is a unital normed algebra. Then the book already essentially proves this result for Banach algebras (Corollary 4.13).

So if $B$ is the completion of $A$, and if $a$ is still in the frontier of $G(B)$, then we're done (the infimum obviously doesn't change if we replace $S$ by the unit sphere of $B$).

Conversely, if there is an example of $a\in\partial G(A)$ with $a\in G(B)$, then we have a counter-example to the exercise. So my question is:

If $a\in\partial G(A)$ and $B$ is the completion of $A$, then is $a\in\partial G(B)$?

Edit: Embarrassingly, I think I can now answer this!

Let $A$ be the complex polynomials, interpreted as an algebra of continuous functions on the interval $[0,1]$. A little bit of algebra shows that $G(A)$ consists of just the constant polynomials. So $G(A)$ is actually closed (not open, which would be the case if $A$ were Banach). So being careful about what "frontier" means, I guess $G(A)$ is its own frontier. But then the exercise is trivially false, as the frontier of $G(A)$ contains invertibles.

So the exercise seems wrong. But somehow my counter-example seems cheap. So a new question:

Can the frontier of $G(A)$ contain a non-invertible element which is invertible in $B$? Are there examples where $G(A)$ is open?

-
I tagged this "homework"-- it was the student's homework, but neither he nor I could do it, so now it's my homework... –  Matthew Daws Feb 29 '12 at 16:29
It still feels to me that there should be a "simpler" or "more natural" commutative counterexample - some kind of explicitly defined algebra of functions on $[0,1]$, for instance... –  user16299 Mar 1 '12 at 23:03

Let $A=\ell^1(F_2)$ sitting inside $B=C_r^*(F_2)$. There exists $a\in A$ which is self-adjoint but whose spectrum in $A$, call it $S$, is not contained in ${\mathbb R}$.
(We can take $a$ to have finite support: I forget the exact formula which works, but it can be found in Palmer Vol. II in one of the sections on "hermitian" groups, the point being that $F_2$ is not hermitian.)
Take a sequence $a_n=a-\lambda_n I$ for scalars $\lambda_n\notin S$ with $\lambda_n\to \lambda$ for some $\lambda \in S \setminus {\mathbb R}$. Then we have a sequence of invertible elements in $A$, which converge in $B$ to an element that is invertible in $B$ but not in $A$. I guess that by taking bicommutants we can get a commutative example from this one, in particular $B=C(X)$ for some $X\subseteq {\mathbb R}$.
Of course this doesn't answer your second question, where (if I understand correctly) you want $G(A)$ to be open in $A$ with respect to the norm from $B$.


text:
# 10.15: Distribution of Differences in Sample Proportions (2 of 5)

### Learning Objectives

• Draw conclusions about a difference in population proportions from a simulation.

Introduction

Recall that we are in the middle of an investigation about the difference in female and male teen depression rates. In our investigation, we are assuming that 26% of female teens and 10% of male teens are depressed. That is, we assume a 16% = 0.16 difference favoring girls.

• We saw a 0.06 gender difference in teen depression rates from the National Survey of Adolescents. Again, girls had a higher rate of depression. Does this study suggest that our assumption about a 0.16 difference in the populations is wrong?
• Or could the results have come from populations with a 0.16 difference in depression rates?

At this point, we may have a sense of the answers to these questions for samples of 64 females and 100 males. But we need to look at the long-run behavior of the differences in sample proportions. We also need to investigate the effect of sample size on our conclusion. The samples in the National Survey of Adolescents are very large.

So we continue this investigation in a Simulation WalkThrough.

A YouTube element has been excluded from this version of the text. You can view it online here: http://pb.libretexts.org/sss/?p=338

On the next page, we use the simulation shown in the WalkThrough to make inferences about a difference in population proportions. As we did in Linking Probability to Statistical Inference, we use a simulation to make observations about the sampling distribution before we develop the mathematical model that we will use in inference. The logic we use to make inferences with simulated sampling distributions is the same logic we use with mathematical models. Let’s practice that way of thinking now.

### Learn By Doing

Suppose in a study of 540 female and 475 male U.S. teens, we find that 8% of the females and 2% of the males are depressed. What does this study suggest about our assumption that the depression rate of female teens is 16% higher than that of male teens in the United States?

Here is a simulated distribution of differences for a large number of independent random samples for these sample sizes. Note that we have rescaled the axis, so the distribution may look wider than the distributions in the WalkThrough, but it actually has less variability.

https://assessments.lumenlearning.co...sessments/3632

https://assessments.lumenlearning.co...sessments/3633

https://assessments.lumenlearning.co...sessments/3634

CC licensed content, Shared previously

10.15: Distribution of Differences in Sample Proportions (2 of 5) is shared under a not declared license and was authored, remixed, and/or curated by LibreTexts.


text:
Visitor

# Hide clues in Microsoft Forms

Greetings,

I'm a teacher and using forms as a way for my students to take exam. It is a math problem so I need a number as answer. As much as possible I want it to be between a certain number for example:

The answer is 7.68374 but I want to consider an answer ranging lower or higher than that so the answer is between "6 to 8"

But the problem is when I put a restriction on form, it shows on answer box the answer must be between "6 to 8"and it gives clues. Can I hide that one? See attached photo. Thank you very much!
0 Replies


text:
# Efficient synchronization of two instances of an ordered list

What data structure or algorithm can be used to efficiently synchronize two nearly identical ordered lists? Two offline systems start with the same ordered list and each edit, insert, delete and move items in the list. Occasionally they will be online to synchronize the list. The algorithm should minimize the amount of changed items or "diffs" that need to be communicated.

For an unordered list we would simply use a modification timestamp and UUID for each item to detect edits respectively inserts/deletions. But what data structure can we use for an ordered list? Using an integer rank would require renumbering all successors after an insert, and thus result in many "diffs" when comparing items in both lists. Using a floating-point rank (having the average of the predecessor and successor ranks) would have accuracy problems after inserting many new items.

Using a doubly linked list (i.e. referencing the UUID of the predecessor and successor item) would result in two changed items and one new item for each insert, thus 3 "diffs" when comparing items in both lists.

• Do you want to recognize moves? First you say yes, then no. – reinierpost Apr 13 '12 at 9:08
• One idea is to regard the specification of the ordering as a separate item, like filesystems do, and synchronize that as well. – reinierpost Apr 13 '12 at 9:11
• Yes, I want to recognize moves. I'd prefer a solution with only one type of items, not content-items and order-items. – Jason Smith Apr 13 '12 at 15:01
• Your question is related to one I asked myself a couple of years ago: what existing theory can be applied to developing a diff that will recognize moves and file renames? But I never tried to answer that question. One idea is to look at algorithms used by software such as git. – reinierpost Apr 16 '12 at 9:55
• Mercurial (Hg) recognizes moves and renames by calculating how similar files are. If they match more than, say, 90% they are probably moved/renamed. I am not sure what algorithm they use to calculate the percentage (e.g. Levenshtein distance). That occasionally results in false positives. Since I have both lists, I can do a direct diffs and recognize changes 100% accurately. I am looking for a way to limit the number of diffs to only the items that actually changed (and not its predecessor or successor). – Jason Smith Apr 16 '12 at 15:19

For synchronizing from a master list to a slave list, it looks like a fellow named Michael Heyeck has a good, O(n) solution to this problem. Check out that blog post for an explanation and some code. Perhaps this is a good starting point to generalize to a two-way sync.

Essentially, the solution tracks both the master and slave lists in a single pass, tracking indices into each. Two data structures are then managed: a list of insertions to be replayed on the slave list, and a list of deletions.

It looks straightforward and also has the benefit of a proof of minimalism, which Heyeck followed up with in a subsequent post. The code snippet in this post is more compact, as well:

def sync_ordered_list(a, b):
x = 0; y = 0; i = []; d = []
while (x < len(a)) or (y < len(b)):
if y >= len(b): d.append(x); x += 1
elif x >= len(a): i.append((y, b[y])); y += 1
elif a[x] < b[y]: d.append(x); x += 1
elif a[x] > b[y]: i.append((y, b[y])); y += 1
else: x += 1; y += 1
return (i,d)

Again, credit to Michael Heyeck.

The standard approach is probably going to be something like this:

• Each side keeps track of the last-synchronized state of the list. In other words, at any point in time, each side knows the state $L_{t-1}$ of the list at the last synchronization point when they synchronized their lists.

• If Alice currently has list $L$ and she wants to send her updates to Bob, she computes the diff between $L_{t-1}$ and $L$ and sends that diff to Bob.

• If both Alice and Bob concurrently made changes to their list, then you need to merge their diffs to form a new list.

• The result, after sychronization is $L_t$, which is known to both parties. Now go back to the beginning.

Now you have some design decisions that you can fill in various ways:

• How should Alice compute the diff between $L_{t-1}$ and her list $L$? One way is for her to keep track of a log of the changes she has made to the list since $L_{t-1}$ and use that as the diff (or possibly to compact that diff, e.g., if she inserts an item and then deletes it, no need to include that in the log of changes). Another way is for her to compute the diff between $L_{t-1}$ and $L$ from scratch, without keeping any logs.

• How should the diff be represented compactly, to minimize network bandwidth and make it as easy as possible for Bob to apply the diff?

• How should we merge Alice's diff and Bob's diff? There are many standard algorithms here; e.g., see source code revision control systems.

I hope that this was not already obvious. Your question did not make it entirely clear how much you already know about the general approach and what specifically you were looking for in an answer.


text:
## Monday, 25 February 2008

### The Lenin of the Week

Taking example from Tommaso's column The Say of the Week, I could not resist advertising this quote:
The Higgs mechanism is just a reincarnation of the Communist Party: it controls the masses. Lenin.
I found it in a recent article of Gian Giudice, who learnt it from Luis Alvarez-Gaume who knew Lenin personally. Lenin's relationship with particle physics is in fact more intimate than you think. He was living and teaching in Geneva in the same century when the LHC construction began. For some period of time he was staying in a house at 6, rue des Plantaporrets in Jonction, the most progressive district of Geneva. By sheer coincidence (or historical necessity) me and my TH friends were often having parties in that same building. There is an inscription on the building commemorating Lenin, though the neighbours much better remember the parties.

Other famous sayings of Lenin:

Sometimes theory needs a push.
Asymptotic freedom will not satisfy the hungry masses.
It is impossible to predict the time and progress of the LHC. It is governed by its own more or less mysterious laws.
One man with a grant can control 100 without one.
A paper cited often enough becomes the truth.
Every cook should be able to run Monte Carlo.

## Saturday, 23 February 2008

### INTEGRAL

At the last Cosmo Coffee, Celine Boehm was discussing the current status of the 511 keV gamma-ray line from the Milky Way center, in light of the recent results from the INTEGRAL satellite. Photons carrying 511 keV energy arise from annihilation of positrons and electrons that are more less at rest. It is not clear which mechanism is responsible for injecting enough positrons into the interstellar medium of the galactic bulge. One hypothesis is that the positrons are scattered remnants of the ILC (Interplanetary Linear Collider) - an unfinished project of a technologically advanced civilization from that region. More contrived explanations involve black holes, radioactive nuclei from supernovae, pulsars and other fluffy toys.

Yet another, quite exciting possibility is that the positrons come from annihilation of dark matter. At first sight this seems quite natural. We are pretty confident that a lot of dark matter is present in the galactic center and its distribution should be roughly spherical, something which agreed quite well with the earlier INTEGRAL observations. Furthermore, it is likely that dark matter can annihilate into ordinary matter so that its present abundance is a thermal relic. But a dark matter particle that could explain the INTEGRAL signal must be quite peculiar. In particular, it cannot be a WIMP: its mass should be in the MeV range. If it were heavier, annihilation would yield too energetic positrons and the 511 keV line would not stand so prominently over the continuum gamma-ray radiation.

On the microscopic level, MeV dark matter can be realized as a scalar particle (so as to avoid the Lee-Weinberg bound) annihilating via an exotic heavy fermion exchange, or via an exotic heavy gauge boson. See the artist's view above. In fact, both diagrams are needed to explain the INTEGRAL signal and, at the same time, derive the correct dark matter abundance from the thermal relic density computation. The first diagram leads to the cross section that goes to a constant at small velocities and this one is relevant for the annihilation of dark matter today. The second leads to the cross section that goes like $\sigma \sim b v^2$, and it's supposed to dominate in the early hot universe.

Recently, the INTEGRAL satellite announced new results which show an asymmetry of emission (by a factor of two) with respect to the central axis of the galaxy. Morevover, the asymmetry seems to be correlated with the distribution of low mass X-ray binaries (LMXB) - systems including a neutron star or a black hole that accretes matter from its companion. LMXBs have long been one of the suspects in the positron case. The general feeling is that the new results make the dark matter explanation unlikely; Julianne on CV readily flushed MeV dark matter down her toilet. Celine, on the other hand, is more reluctant to push the button (she is, in fact, responsible for much of the stuff being thrown into the toilet). She argued that: 1) We still don't know a conventional astrophysical explanation that could account for enough positron emission, 2) While part of the gamma-ray emission maybe due to boring astrophysics, there's still a large, roughly spherical component that could be due to MeV dark matter. In fact, if only a part of the emission is assigned to dark matter, the microscopic models can more readily satisfy constraints from other experiments, for example, from measurements of the electron anomalous magnetic moment.

I guess it's fair to say that, at the moment, a conventional astrophysical explanation seems far more likely. My faith in that signal is further diminished by the fact that astrophysics provides us with too many excesses (EGRET, for example), each one having its dark matter particle that's supposed to explain it. Finally, MeV dark matter is hard to accommodate in our favourite Beyond the Standard Model schemes. But this last argument should always be taken with a grain of salt, since all existing BSM models suck. It's better to keep our eyes open...

## Monday, 18 February 2008

### Non-Gaussianity Again

I wrote not so long ago about a possible observation non-Gaussianity in the CMB spectrum. Last Wednesday Tony Riotto gave an overwiew seminar explaining what that non-Gaussianity really is. I used to think that the issue was complicated. Tony explained very clearly and pedagogically that it is bloody damn complicated.

Unfortunately, it is also important. There is an overabundance of inflationary models and several mechanisms in which these models could produce the primordial density perturbations. The number of observables to discriminate between the models is less outgrown. What we've got so far is the amplitude and the spectral index of density perturbations. Applied to the minimal slow-roll scenario, these two provide useful constraints, but at the moment they do not clearly favour a particular class of models. One way to improve the situation is to measure the spectral index with greater precision. Detecting tensor perturbations would be fantastic, but many models (and most of those sensible from the effective field theory point of view) predict tensor perturbations well below the sensitivity of near-future experiments. Non-Gaussianity of the primordial perturbation spectrum provides another complementary handle on inflationary models. A mere detection of large primordial non-Gaussianity would flush many models down the drain, including the single field slow-roll inflation.

Gaussianity of the primordial perturbation spectrum arises when the inflaton behaves approximately as a free field in a curved background. This is never quite true: the inflaton potential introduces non-linear self-interactions and, more importantly, the inflaton is coupled to gravity which is itself non-linear. In simplest inflationary scenarios, however, these non-linear effects are negligible. In such a case, the equation governing the inflaton perturbations is just the harmonic oscillator equation with a time-varying frequency. Different oscillator frequencies are decoupled and they oscillate independently, seeding the universe with Gaussian perturbations. This means that the power spectrum of the two-point correlation function (that funny thing with peaks that CMB people show) carries all the statistical information about density fluctuations. In the presence of non-Gaussianity, the three- and higher-point correlation functions are independent from the two-point one. The only way to learn about non-Gaussianity is to perform statical tests on these higher-point functions.

Non-Gaussianity is conventionally parametrized as
$\zeta(x) = \zeta_{gauss}(x) - 3/5 f_{NL} (\zeta_{gauss}^2(x) - <\zeta_{gauss}^2>)$
Here $\zeta$ is a gauge-invariant characterization of curvature perturbations. The coefficient $f_{NL}$ is not a constant in general but, for simplicity, one often presents just a number corresponding to an overall amplitude of non-Gaussian contributions over the scales of interest. The current constraint for $f_{NL}$ is that it is at most of order 100. It follows that the second non-Gaussian term is a small perturbation (recall that $\zeta \sim 10^{-5}$) and that's what is meant by saying that the primordial perturbation spectrum is approximately Gaussian.

There are several inflationary models that could lead to primordial non-Gaussianity. Tony paused on the curvaton scenario that he seems to like. In the curvaton scenario, density fluctuations are produced by a different field than the inflaton, while perturbations from the inflaton field are negligible (which may well be, for example, when the scale of inflation is low enough). There is another field - the curvaton - who has a negligible impact on the dynamics of inflation, yet it acquires isocurvature perturbations during inflation. As it decays after the end of inflation, it transmits the fluctuation to the hot plasma in the early universe. Non-gaussianity of the primordial spectrum depends on the parameter $r$, which is a ratio of the curvaton energy density to the total energy density just before the curvaton decays. Lyth et al. found $f_{NL} \sim 5/(4r)$.

One thing that makes non-Gaussianity particularly complicated is that detecting it in the CMB does not imply detecting the primordial non-Gaussianity. Many phenomena contaminate the CMB during the evolution of the universe. First of all, gravity is non-linear, so that the second order perturbation theory introduces non-Gaussianity into the CMB at the level of $f_{NL} \sim 1$. This is the minimal amount of non-Gaussianity that we expect in the CMB, if our models of density perturbations evolution are correct. For comparison, the single field slow-roll inflation predicts $f_{NL} \sim (n - 1)$, which is drowned by the gravity contribution because the spectral index n is close to 1. On top of that, there is a whole lot of dirty astrophysics that contributes to non-Gaussianity: secondary anistropies from graviational lensing or Sunyayev-Zeldovich effect, radio point sources, instrument noise etc.

Concerning the possible signal I wrote about, Tony advised not to get excited. Many effects could bias the analysis of Yadav and Wandelt, one thing pointed out was the poorly known distribution of extragalactic radio sources, another was over-substraction of the noise. On the other hand, if the signal is real, Planck will confirm it at 8 sigma confidence level.

A few details or formulas omitted here can be found in this review by Tony et al.

## Thursday, 7 February 2008

### Theta-one-three

I have been feeling dizzy all this week and I blame it on the neutrino flux. It happens to be more intense than ever - within 7 days there were 4 neutrino seminars here at CERN. I could not help but learn something about the quest for $\theta_{13}$ from the Thursday seminar by Silvia Pascoli.

Neutrino physics is the only branch of particle physics that can boast of a discovery in the last twenty years (the top quark, some may say, but that was as unexpected as election results in the former Soviet Union). The experimental progress in the last decade has been really impressive. Since the 1998 discovery of neutrino oscillations we have acquired a great deal of data concerning the neutrino masses and mixings. Yet from a theorist's perspective neutrinos are boring. It all amounts to a 3-by-3 unitary matrix that relates flavour eigenstates to mass eigenstates. The matrix is called the MNS matrix and looks like this

As usual, a 3-by-3 unitary matrix has three angles and three phases. Two of these angles were pinpointed by numerous observations of solar and atmospheric neutrino oscillations, as well as by reactor and accelerator experiments. The last angle, the notorious $\theta_{13}$, remains elusive and we only know the upper bound of 12 degrees. The other two angles are much larger than that, which gives hope to the experimenters that the last one is just behind the corner. So, while cosmologists are studying such sexy matters like dark matter or dark energy, the neutrino community is struggling to measure an angle.

One way to measure $\theta_{13}$ is by observing oscillations of muon neutrinos into electron neutrinos. The probability of this process is sensitive $\sin^2 (2 \theta_{13})$,
$P(\nu_\mu -> \nu_e) = \sin^2 \theta_{23} \sin^2 (2 \theta_{13}) \sin^2 (\Delta m_{13}^2 L/2 E)$.
At the moment, there is an ongoing search for electron neutrino appearance by the Minos experiment that shoots the beam of muon neutrinos from Fermilab to a detector in Soudan mines in Minnesota. This one is however not designed specifically for that search and no great improvement in sensitivity should be expected. In near future (2009-ish) starts another accelerator experiment called T2K (Tokai to Kamioka in Japan, where else) whose sensitivity to $\sin^2 (2 \theta_{13})$ will reach 0.01, (or 3 degrees after some complicated trigonometry). In parallel, there are reactor experiments in preparation, who will play with anti-neutrinos that escape from nuclear reactors and pollute the air. The strategy is to put two detectors close to a reactor and look for disappearance of electron anti-neutrinos on the way between the two. The next reactor experiment to start is Double Chooz, the continuation of the Chooz experiment who already put a bound on $\theta_{13}$ in the past. The Chinese also enter the race with their Daya Bay experiment. These two reactor experiments will have a similar sensitivity to $\theta_{13}$ as T2K.

The measurement of the angle $\theta_{13}$ is very important because it will open the way to measure another angle. More precisely, the CP violating phase $\delta$ in the MNS matrix. The point is that CP violating effects, that is differences in oscillation patterns between neutrinos and antineutrinos, are proportional to $\sin \delta$ but also to $\sin \theta_{13}$. If the latter is large enough, the future superbeams, beta-beams or neutrino factories may discover CP violation in the leptonic sector.

In summary, the future looks bright and gay for neutrinos, as on the picture below

No slides available.


text:
Question

# Prove that: $\frac{{\sin \theta }}{{1 - \cos \theta }} = cosec\theta + \cot \theta$

Verified
163.5k+ views
Hint: - Taking conjugate of denominator.
Given,
L.H.S $\begin{gathered} = \frac{{\sin \theta }}{{1 - \cos \theta }} \\ \\ \end{gathered}$
Multiply and Divide by $\left( {1 + \cos \theta } \right)$ , we get
$\begin{gathered} = \frac{{\sin \theta }}{{1 - \cos \theta }} \times \frac{{1 + \cos \theta }}{{1 + \cos \theta }} \\ = \frac{{\sin \theta \left( {1 + \cos \theta } \right)}}{{\left( {1 - \cos \theta } \right)\left( {1 + \cos \theta } \right)}} \\ = \frac{{\sin \theta \left( {1 + \cos \theta } \right)}}{{{{\left( 1 \right)}^2} - {{\left( {\cos \theta } \right)}^2}}} \\ = \frac{{\sin \theta \left( {1 + \cos \theta } \right)}}{{1 - {{\cos }^2}\theta }} \\ \end{gathered}$
We know that, ${\sin ^2}\theta + {\cos ^2}\theta = 1$
Or, ${\sin ^2}\theta = 1 - {\cos ^2}\theta$
Replace $\left( {1 - {{\cos }^2}\theta } \right)$ by ${\sin ^2}\theta$ , we get
L.H.S $\begin{gathered} = \frac{{\sin \theta \left( {1 + \cos \theta } \right)}}{{{{\sin }^2}\theta }} \\ = \frac{{\sin \theta \left( {1 + \cos \theta } \right)}}{{\sin \theta \times \sin \theta }} \\ \end{gathered}$
Cancel out $\sin \theta$ in numerator by $\sin \theta$ in denominator, we get
$\begin{gathered} = \frac{{1 + \cos \theta }}{{\sin \theta }} \\ = \frac{1}{{\sin \theta }} + \frac{{\cos \theta }}{{\sin \theta }} \\ \end{gathered}$
Now we can written $\frac{1}{{\sin \theta }} = \cos ec\theta$ and $\frac{{\cos \theta }}{{\sin \theta }} = \cot \theta$ , we get
L.H.S $= \cos ec\theta + \cot \theta$- (1)
Since, given R.H.S$= \cos ec\theta + \cot \theta$ - (2)
By seeing equation (1) and (2) we can tell that
L.H.S R.H.S
Hence, it proved.
Note: - These types of questions are also solve by taking R.H.S (Right Hand Side), solve it to prove L.H.S (Left Hand Side). During solving trigonometry proving we should always have basic trigonometry identities in our mind.