## [e6709d]: fourier / floor23.m  Maximize  Restore  History

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94``` ```function [nfft,tableout]=floor23(n) %FLOOR23 Previous number with only 2,3 factors % Usage: nceil=floor23(n); % % `floor23(n)` returns the first number less than or equal to *n*, % which can be written as a product of powers of *2* and *3*. % % The algorithm will look up the best size in a table, which is computed % the first time the function is run. If the input size is larger than the % largest value in the table, the input size will be reduced by factors of % *2*, until it is in range. % % `[nceil,table]=floor23(n)` additionally returns the table used for lookup. % % Examples: % --------- % % Return the first number smaller or equal to *26* that can be written % solely as products of powers of *2* and *3*::: % % floor23(26) % % This plot shows the behaviour of |floor23| and |ceil23| for numbers % up to 100::: % % x=1:100; % plot(x,floor23(x),x,ceil23(x)); % legend('floor23','ceil23','Location','Northwest'); % % See also: ceil23, floor235, nextfastfft % AUTHOR: Peter L. S��ndergaard persistent table; maxval=2^20; if isempty(table) % Compute the table for the first time, it is empty. l2=log(2); l3=log(3); l5=log(5); lmaxval=log(maxval); table=zeros(143,1); ii=1; prod2=1; for i2=0:floor(lmaxval/l2) prod3=prod2; for i3=0:floor((lmaxval-i2*l2)/l3) table(ii)=prod3; prod3=prod3*3; ii=ii+1; end; prod2=prod2*2; end; table=sort(table); end; % Copy input to output. This allows us to efficiently work in-place. nfft=n; % Handle input of any shape by Fortran indexing. for ii=1:numel(n) n2reduce=0; if n(ii)>maxval % Reduce by factors of 2 to get below maxval n2reduce=ceil(log2(nfft(ii)/maxval)); nfft(ii)=nfft(ii)/2^n2reduce; end; % Use a simple bisection method to find the answer in the table. from=1; to=numel(table); while from<=to mid = round((from + to)/2); diff = table(mid)-nfft(ii); if diff<0 from=mid+1; else to=mid-1; end end if nfft(ii)~=table(from) nfft(ii)=table(from-1); end; % Add back the missing factors of 2 (if any) nfft(ii)=nfft(ii)*2^n2reduce; end; tableout=table; ```