/***************************************************************************  **************************************************************************                  Spherical Harmonic Transform Kit 2.7       Contact: Peter Kostelec            geelong@cs.dartmouth.edu       Copyright 1997-2003  Sean Moore, Dennis Healy,                        Dan Rockmore, Peter Kostelec       Copyright 2004  Peter Kostelec, Dan Rockmore     SpharmonicKit is free software; you can redistribute it and/or modify     it under the terms of the GNU General Public License as published by     the Free Software Foundation; either version 2 of the License, or     (at your option) any later version.       SpharmonicKit is distributed in the hope that it will be useful,     but WITHOUT ANY WARRANTY; without even the implied warranty of     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     GNU General Public License for more details.       You should have received a copy of the GNU General Public License     along with this program; if not, write to the Free Software     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.       Commercial use is absolutely prohibited.     See the accompanying LICENSE file for details.    ************************************************************************  ************************************************************************//* source code for generating cosine transforms    of Pml and Gml functions */#include <math.h>#include <string.h>   /* to declare memcpy */#include <stdlib.h>#include "newFCT.h"#include "primitive.h"#ifdef FFTPACK#include "fftpack.h"#endif#ifndef PI#define PI 3.14159265358979#endif/************************************************************************//* look up table for the power of 2 which is >= i */int Power2Ceiling(int i){  int pow_2;    pow_2 = 1;  if (i == 0)    return 0;  else    {      while ( pow_2 < i )	pow_2 *= 2;      return pow_2 ;    }}/************************************************************************//* utility functions for table management *//************************************************************************//* Computes the number of non-zero entries in a table containing   cosine series coefficients of all the P(m,l) or G(m,l) functions   necessary for computing the seminaive transform for a given   bandwidth bw and order m.  Works specifically for tables   generated by CosPmlTableGen() */int TableSize(int m,	      int bw){    return (((bw/2) * ((bw/2) + 1)) -	    ((m/2)*((m/2)+1)) - ((bw/2) * (m % 2)));}/************************************************************************//* Spharmonic_TableSize(bw) returns an integer value for   the amount of space necessary to fill out an entire spharmonic   table.  Note that in the above TableSize() formula,    you need to sum this formula over m as m ranges from 0 to   (bw-1).  The critical closed form that you need is that   \sum_{k=0}^n = \frac{(n(n+1)(2n+1)}{6}   You also need to account for integer division.   From this you should derive an upper bound on the   amount of space.    Some notes - because of integer division, you need to account for   a fudge factor - this is the additional " + bw" at the   end.  This gaurantees that you will always have slightly more   space than you need, which is clearly better than underestimating!   Also, if bw > 512, the closed form   fails because of the bw*bw*bw term (at least on Sun Sparcstations)   so the loop computation is used instead.   Also, the transpose is exactly the same size, obviously.*/int Spharmonic_TableSize(int bw){  int m, sum;    if (bw > 512)    {      sum = 0;            for (m=0; m<bw; m++)	  sum += TableSize(m,bw);            return sum;    }  else    {      return (	      (((4*(bw*bw*bw)) + (6*(bw*bw)) - (8*bw))/24)	      + bw	      );    }}/************************************************************************//* Reduced_Spharmonic_TableSize(bw,m) returns an integer value for   the amount of space necessary to fill out a spharmonic table   if interesting in using it only for orders up to (but NOT   including) order m.   This will be used in the hybrid algorithm's call of the   semi-naive algorithm (which won't need the full table ... hopefully   this'll cut down on the memory usage).   Also, the transpose is exactly the same size, obviously.   This is a "reduced" version of Spharmonic_TableSize(m).*/int Reduced_SpharmonicTableSize(int bw,				int m){    int i, sum;  sum = 0;    for (i=0; i<m; i++)      sum += TableSize(i,bw);  return sum;}/************************************************************************//* For an array containing cosine series coefficients of Pml or Gml   functions, computes the location of the first coefficient of Pml.   This supersedes the TableOffset() function.   Assumes table is generated by CosPmlTableGen()*/int NewTableOffset(int m,		   int l){    int offset;    int tm, tl;        if ((m % 2) == 1)      {	tl = l-1;	tm = m-1;      }    else      {	tl = l;	tm = m;      }    offset = ((tl/2)*((tl/2)+1)) - ((tm/2)*((tm/2)+1));    if (tl % 2 == 1)      offset += (tl/2)+1;    return offset;}/************************************************************************//* Computes the offset in the table generated by CosPmlTableGen for   the cosine series of Pml OR GML.  NOTE THAT M MUST BE EVEN -    MUST BE HANDLED BY CALLER    Superseded by NewTableOffset()*/int TableOffset(int m,		int l){    int offset;    offset = ((l/2)*((l/2)+1)) - ((m/2)*((m/2)+1));    if (l % 2 == 1)      offset += (l/2)+1;    return offset;}/************************************************************************//* generate all of the cosine series for Pmls or Gmls for a specified    value of m.  Note especially that since series are zero-striped,   all zeroes have been removed.     tablespace points to a double array of size TableSize(m,bw);   Workspace needs to be   16 * bw    Let P(m,l,j) represent the jth coefficient of the   cosine series representation of Pml.  The array   stuffed into tablespace is organized as follows:   P(m,m,0)    P(m,m,2)  ...  P(m,m,m)   P(m,m+1,1)  P(m,m+1,3)...  P(m,m+1,m+1)   P(m,m+2,0)  P(m,m+2,2) ... P(m,m+2,m+2)   etc.  Appropriate modifications are made for m odd (Gml functions).*/void CosPmlTableGen(int bw,		    int m,		    double *tablespace,		    double *workspace){    double *prev, *prevprev, *temp1, *temp2, *temp3, *temp4, *x_i, *eval_args;    double *tableptr, *cosres, *cosworkspace;    int i, j, k;    prevprev = workspace;    prev = prevprev + bw;    temp1 = prev + bw;    temp2 = temp1 + bw;    temp3 = temp2 + bw;    temp4 = temp3 + bw;    x_i = temp4 + bw;    eval_args = x_i + bw;    cosres = eval_args + bw;    cosworkspace = cosres + bw;    tableptr = tablespace;#ifdef FFTPACK    CoswSave = precomp_dct( bw );#endif        /* main loop */    /* Set the initial number of evaluation points to appropriate       amount */    /* now get the evaluation nodes */    EvalPts(bw,x_i);    ArcCosEvalPts(bw,eval_args);        /* set initial values of first two Pmls */    for (i=0; i<bw; i++)       prevprev[i] = 0.0;    if (m == 0) {	for (i=0; i<bw; i++) {	  prev[i] = 1.0 ;	}    }    else       Pmm_L2(m, eval_args, bw, prev);        if ((m % 2) == 1) { /* need to divide out sin x */	for (i=0; i<bw; i++)	  prev[i] /= sin(eval_args[i]);    }      /* set k to highest degree coefficient */    if ((m % 2) == 0)      k = m;    else      k = m-1;	    /* now compute cosine transform */#ifndef FFTPACK    kFCT(prev, cosres, cosworkspace, bw, bw, 1);#else    memcpy(cosres, prev, sizeof(double) * bw);    DCTf( cosres, bw, bw, CoswSave );#endif    for (i=0; i<=k; i+=2)      tableptr[i/2] = cosres[i];    /* update tableptr */    tableptr += k/2+1;    /* now generate remaining pmls  */    for (i=0; i<bw-m-1; i++) {	vec_mul(L2_cn(m,m+i),prevprev,temp1,bw);	vec_pt_mul(prev, x_i, temp2, bw);	vec_mul(L2_an(m,m+i), temp2, temp3, bw);	vec_add(temp3, temp1, temp4, bw); /* temp4 now contains P(m,m+i+1) */	/* compute cosine transform */#ifndef FFTPACK		kFCT(temp4, cosres, cosworkspace, bw, bw, 1);#else	memcpy(cosres, temp4, sizeof(double) * bw);	DCTf( cosres, bw, bw, CoswSave );#endif	/* update degree counter */	k++;	/* now put decimated result into table */	if ((i % 2) == 1) {	    for (j=0; j<=k; j+=2)	      tableptr[j/2] = cosres[j];	}	else {	    for (j=1; j<=k; j+=2)	      tableptr[j/2] = cosres[j];	}	/* update tableptr */	tableptr += k/2+1;	/* now update Pi and P(i+1) */	memcpy(prevprev, prev, sizeof(double) * bw);	memcpy(prev, temp4, sizeof(double) * bw);    }#ifdef FFTPACK    free( CoswSave );#endif}/****  just like above except does not compute all the  legendres for a given m, bw that semi-naive would  use, just those that semi-naive IN THE FLT_HYBRID CODE  would use  lim = degree of FIRST coefficient in the flt_hybrid        code that will NOT be computed using semi-naive  lim is defined in this strange way so that this  function will produce the SAME array as  CosPmlTableGen if lim = bw.  tablespace points to a double array of size cos_size  as computed in SplitLocales().  *****/void CosPmlTableGenLim( int bw, 			int m,			int lim,			double *tablespace,			double *workspace){    double *prev, *prevprev, *temp1, *temp2, *temp3, *temp4, *x_i, *eval_args;    double *tableptr, *cosres, *cosworkspace;    int i, j, k;    prevprev = workspace;    prev = prevprev + bw;    temp1 = prev + bw;    temp2 = temp1 + bw;    temp3 = temp2 + bw;    temp4 = temp3 + bw;    x_i = temp4 + bw;    eval_args = x_i + bw;    cosres = eval_args + bw;    cosworkspace = cosres + bw;    tableptr = tablespace;#ifdef FFTPACK    CoswSave = precomp_dct( bw );#endif    /* main loop */    /* Set the initial number of evaluation points to appropriate       amount */    /* now get the evaluation nodes */    EvalPts(bw,x_i);    ArcCosEvalPts(bw,eval_args);        /* set initial values of first two Pmls */    for (i=0; i<bw; i++)       prevprev[i] = 0.0;    if (m == 0) {	for (i=0; i<bw; i++) {	    prev[i] = 1.0;	}    }    else       Pmm_L2(m, eval_args, bw, prev);        if ((m % 2) == 1) { /* need to divide out sin x */	for (i=0; i<bw; i++)	  prev[i] /= sin(eval_args[i]);    }      /* set k to highest degree coefficient */    if ((m % 2) == 0)      k = m;    else      k = m-1;	    /* now compute cosine transform */#ifndef FFTPACK    kFCT(prev, cosres, cosworkspace, bw, bw, 1);#else    memcpy(cosres, prev, sizeof(double) * bw);    DCTf( cosres, bw, bw, CoswSave );#endif    for (i=0; i<=k; i+=2)      tableptr[i/2] = cosres[i];    /* update tableptr */    tableptr += k/2+1;    /* now generate remaining pmls  */    for (i = 0 ; i < lim - m - 1 ; i ++) {	vec_mul(L2_cn(m,m+i),prevprev,temp1,bw);	vec_pt_mul(prev, x_i, temp2, bw);	vec_mul(L2_an(m,m+i), temp2, temp3, bw);	vec_add(temp3, temp1, temp4, bw); /* temp4 now contains P(m,m+i+1) */	/* compute cosine transform */#ifndef FFTPACK		kFCT(temp4, cosres, cosworkspace, bw, bw, 1);#else	memcpy(cosres, temp4, sizeof(double) * bw);	DCTf( cosres, bw, bw, CoswSave );#endif	/* update degree counter */	k++;	/* now put decimated result into table */	if ((i % 2) == 1) {	    for (j=0; j<=k; j+=2)	      tableptr[j/2] = cosres[j];	}	else {	    for (j=1; j<=k; j+=2)	      tableptr[j/2] = cosres[j];	}	/* update tableptr */	tableptr += k/2+1;	/* now update Pi and P(i+1) */	memcpy(prevprev, prev, sizeof(double) * bw);	memcpy(prev, temp4, sizeof(double) * bw);    }#ifdef FFTPACK    free( CoswSave );#endif}/************************************************************************//* RowSize returns the number of non-zero coefficients in a row of the   cospmltable if were really in matrix form.  Helpful in transpose   computations.  It is helpful to think of the parameter l as   the row of the corresponding matrix.*/int RowSize(int m,	    int l){  if (l < m)    return 0;  else    {      if ((m % 2) == 0)	return ((l/2)+1);      else	return (((l-1)/2)+1);    }}/************************************************************************//* Transposed row size returns the number of non-zero coefficients   in the transposition of the matrix representing a cospmltable.   Used for generating arrays for inverse seminaive transform.   Unlike RowSize, need to know the bandwidth bw.  Also, in   the cospml array, the first m+1 rows are empty, but in   the transpose, all rows have non-zero entries, and the first   m+1 columns are empty.  So the input parameters are a bit different   in the you need to specify the row you want.*/int Transpose_RowSize(int row,		      int m,		      int bw){  if (row >= bw)    return 0;  else if ((m % 2) == 0)    {      if (row <= m)	return ( ((bw-m)/2) );      else	return ( ((bw-row-1)/2) + 1);    }  else    {      if (row == (bw-1))	return 0;      else if (row >= m)	return (Transpose_RowSize(row+1,m-1,bw));      else    /*  (row < m)  */	return (Transpose_RowSize(row+1,m-1,bw) - (row % 2));    }}/************************************************************************//* Inverse transform is transposition of forward transform.   Thus, need to provide transposed version of table   returned by CosPmlTableGen.  This function does that   by taking as input a cos_pml_table for a particular value   of bw and m, and loads the result as a   transposed, decimated version of it for use by an inverse    seminaive transform computation.   result needs to be of size TableSize(m,bw)*/void Transpose_CosPmlTableGen(int bw,			      int m,			      double *cos_pml_table,			      double *result){  /* recall that cospml_table has had all the zeroes     stripped out, and that if m is odd, then it is     really a Gml function, which affects indexing a bit.