algorithm/osp/fpsup.c
LQ15_t MUL_LQ15(LQ15_t a, LQ15_t b)
{
	/* Long multiplication principle:
	 * (From: http://www.edaboard.com/thread253439.html)
	 *
	 * x,y is the 2 multiplier (signed 64 bit)
	 * after sign valuation and (if need) sign inversion
	 * it is split into xh xl yh yl (unsigned 32 bit)
	 * the operations performed are:
	 *
	 *               xh xl  *
	 *               yh yl  =
	 *               -------
	 *               xl_yl  +
	 *            xh_yl     +
	 *            xl_yh     +
	 *         xh_yh        =
	 *         -------------
	 *  (sign) z3 z1 z2 z0
	 *
	 */

	int8_t sign = 1;

	/* Intermediates */
	uint64_t al_bl, ah_bl, al_bh, ah_bh;

	uint32_t al, ah;
	uint32_t bl, bh;

	uint64_t res0, res1;

	if (a<0) {
		a = -a;
		sign = -sign;
	}

	if (b<0) {
		b = -b;
		sign = -sign;
	}

	al = a & 0xffffffff; ah = a>>32;
	bl = b & 0xffffffff; bh = b>>32;
	if (ah != 0 || bh != 0) {
		printf("HY-DBG: OVERFLOW in %s\n", __func__);
		return 0;
	}

	al_bl = (uint64_t)al*(uint64_t)bl;
	ah_bl = (uint64_t)ah*(uint64_t)bl;
	al_bh = (uint64_t)al*(uint64_t)bh;
	ah_bh = (uint64_t)ah*(uint64_t)bh;

	res0 = al_bl & 0xffffffff;

	res1 = (al_bl>>32)+(ah_bl & 0xffffffff)+(al_bh & 0xffffffff);
	res0 |= (res1&0xffffffff) << 32;

	res1 >>= 32;              // this is 'carry out' of the previous sum
	res1 += (ah_bl>>32) + (al_bh>>32) +ah_bh;

	res0 >>= Q15_SHIFT;

	res0 &= ~(((uint64_t)(1<<15)-1) << (64-Q15_SHIFT));
	res0 |= (res1 & ((1<<15)-1))  << (64-Q15_SHIFT);

	return res0;
}
/* MUL_LQ15 */
/*
 *  LQ15 - 64 bit storage:
 *    1 signed bit
 *    63 bits of value divided into:
 *        48bits of integer
 *        15bits of fractional.
 *  Multiplication of 2 LQ15 produces 126bits of value.
 *  15bits needs to be thrown away to convert it to 15 bits of
 *  fractional. This leaves 111 bits of numbers.
 *
 * LQ15 is to support intermediates used for sqrt() Q15 guesses. So
 * the largest value needed is:
 * 	1 sign bit
 * 	2*16 integer bits
 * 	15 fractional bits
 * 	------
 * 	48bits total. Thus the upper 63 bits is always zero. (except for
 * 		sign extension).
 */
MUL_LQ15
LQ15 - 64 bit storage:
1 signed bit
63 bits of value divided into:
48bits of integer
15bits of fractional.
Multiplication of 2 LQ15 produces 126bits of value.
15bits needs to be thrown away to convert it to 15 bits of
fractional. This leaves 111 bits of numbers.

LQ15 is to support intermediates used for sqrt() Q15 guesses. So
the largest value needed is:
1 sign bit
2*16 integer bits
15 fractional bits

48bits total. Thus the upper 63 bits is always zero. (except for
sign extension).
MUL_LQ15 LQ15 - 64 bit storage : 1 signed bit 63 bits of value divided into : 48bits of integer 15bits of fractional . Multiplication of 2 LQ15 produces 126bits of value . 15bits needs to be thrown away to convert it to 15 bits of fractional . This leaves 111 bits of numbers . LQ15 is to support intermediates used for sqrt () Q15 guesses . So the largest value needed is : 1 sign bit 2 * 16 integer bits 15 fractional bits 48bits total . Thus the upper 63 bits is always zero . ( except for sign extension ) .
LQ15_t MUL_LQ15(LQ15_t a, LQ15_t b)
{
	int8_t sign = 1;
	uint64_t al_bl, ah_bl, al_bh, ah_bh;
	uint32_t al, ah;
	uint32_t bl, bh;
	uint64_t res0, res1;
	if (a<0) {
		a = -a;
		sign = -sign;
	}
	if (b<0) {
		b = -b;
		sign = -sign;
	}
	al = a & 0xffffffff; ah = a>>32;
	bl = b & 0xffffffff; bh = b>>32;
	if (ah != 0 || bh != 0) {
		printf("HY-DBG: OVERFLOW in %s\n", __func__);
		return 0;
	}
	al_bl = (uint64_t)al*(uint64_t)bl;
	ah_bl = (uint64_t)ah*(uint64_t)bl;
	al_bh = (uint64_t)al*(uint64_t)bh;
	ah_bh = (uint64_t)ah*(uint64_t)bh;
	res0 = al_bl & 0xffffffff;
	res1 = (al_bl>>32)+(ah_bl & 0xffffffff)+(al_bh & 0xffffffff);
	res0 |= (res1&0xffffffff) << 32;
	res1 >>= 32;              
	res1 += (ah_bl>>32) + (al_bh>>32) +ah_bh;
	res0 >>= Q15_SHIFT;
	res0 &= ~(((uint64_t)(1<<15)-1) << (64-Q15_SHIFT));
	res0 |= (res1 & ((1<<15)-1))  << (64-Q15_SHIFT);
	return res0;
}
LQ15_t MUL_LQ15 ( LQ15_t a , LQ15_t b ) { int8_t sign = 1 ; uint64_t al_bl , ah_bl , al_bh , ah_bh ; uint32_t al , ah ; uint32_t bl , bh ; uint64_t res0 , res1 ; if ( a < 0 ) { a = - a ; sign = - sign ; } if ( b < 0 ) { b = - b ; sign = - sign ; } al = a & 0xffffffff ; ah = a >> 32 ; bl = b & 0xffffffff ; bh = b >> 32 ; if ( ah != 0 || bh != 0 ) { printf ( " \n " , __func__ ) ; return 0 ; } al_bl = ( uint64_t ) al * ( uint64_t ) bl ; ah_bl = ( uint64_t ) ah * ( uint64_t ) bl ; al_bh = ( uint64_t ) al * ( uint64_t ) bh ; ah_bh = ( uint64_t ) ah * ( uint64_t ) bh ; res0 = al_bl & 0xffffffff ; res1 = ( al_bl >> 32 ) + ( ah_bl & 0xffffffff ) + ( al_bh & 0xffffffff ) ; res0 |= ( res1 & 0xffffffff ) << 32 ; res1 >>= 32 ; res1 += ( ah_bl >> 32 ) + ( al_bh >> 32 ) + ah_bh ; res0 >>= Q15_SHIFT ; res0 &= ~ ( ( ( uint64_t ) ( 1 << 15 ) - 1 ) << ( 64 - Q15_SHIFT ) ) ; res0 |= ( res1 & ( ( 1 << 15 ) - 1 ) ) << ( 64 - Q15_SHIFT ) ; return res0 ; }
MUL_LQ15
  LQ15 - 64 bit storage:
    1 signed bit
    63 bits of value divided into:
        48bits of integer
        15bits of fractional.
MUL_LQ15 LQ15 - 64 bit storage : 1 signed bit 63 bits of value divided into : 48bits of integer 15bits of fractional .
/* Long multiplication principle:
	 * (From: http://www.edaboard.com/thread253439.html)
	 *
	 * x,y is the 2 multiplier (signed 64 bit)
	 * after sign valuation and (if need) sign inversion
	 * it is split into xh xl yh yl (unsigned 32 bit)
	 * the operations performed are:
	 *
	 *               xh xl  *
	 *               yh yl  =
	 *               -------
	 *               xl_yl  +
	 *            xh_yl     +
	 *            xl_yh     +
	 *         xh_yh        =
	 *         -------------
	 *  (sign) z3 z1 z2 z0
	 *
	 */ /* Intermediates */ // this is 'carry out' of the previous sum