# Is there a way to get BSplineFunction beyond MachinePrecision?

I have been using BSplineFunction (Mathematica 8) to generate a smooth representation of some data, which I might then do some further processing on.

As I understand it, BSplineFunction should be analytic and I should be able to pull out a value anywhere. The results of BSplineFunction however are limited to MachinePrecision, regardless of the input data's precision or accuracy. For example,

data = {{17980, -1}, {17990, -1}, {18000, 0}, {18010, 1}, {18020,
3}, {18040, 0}, {18050, -1}, {18080, 2}, {18100, 0}};

f = BSplineFunction[data, SplineDegree -> 3];


f[1/5]


will give

{18001.6, 0.158667}


Initially, I thought something like

N[f[1/5],300]


would produce an answer with 300 digits for example, but instead it only returns the same answer as above.In fact, changing the 300 to 3 also produces the same result.

Note, Precision[data] and Accuracy[data] both give $\infty$. However, Precision[f], and Accuracy[f] give MachinePrecision, and 14.9546, respectively.

PrecisionGoal and AccuracyGoal are not available options for BSplineFunction (it lets you know if you try). It will take setting WorkingPrecision to some number as an option, for example,

f = BSplineFunction[data, SplineDegree -> 3, WorkingPrecision -> 8000]


without giving any error, but the answers remain unchanged.

Am I trying to get answers out of BSplineFunction it simply isn't possible for it to give, or is there perhaps some other setting I need to explore?

Thank you for any assistance.

-
This came up recently and no one had a suggestion for precise or arbitrary precision output. I'm not saying it cannot be done however. +1 on the question. – Mr.Wizard Aug 1 '13 at 0:17

Although BSplineFunction[] is sadly limited to machine precision results, it's not too hard in this case to make a function that will give exact results for exact input. You've already given the control points, so the task is a whole lot easier than the situation in this related answer. Just as in that answer, we use the strategy of starting with BSplineBasis[], and all we have to do is make the (uniform) knots:

data = {{17980, -1}, {17990, -1}, {18000, 0}, {18010, 1}, {18020, 3},
{18040, 0}, {18050, -1}, {18080, 2}, {18100, 0}};

m = 3; n = Length[data];
knots = Join[ConstantArray[0, m], Range[0, 1, 1/(n - m)], ConstantArray[1, m]];

fnew[x_] = Table[BSplineBasis[{m, knots}, j - 1, x], {j, n}].data;


Test:

f = BSplineFunction[data, SplineDegree -> 3];
f[1/5]
{18001.6, 0.158667}

fnew[1/5]
{1350118/75, 119/750}
N[%]
{18001.6, 0.158667}

-
I have a small requirement that could you add a option WorkingPrecision in your funtion fnew like the built-in FindRoot.For instance, FindRoot[Cos[x^2] - x, {x, 1}, WorkingPrecision -> 30] gives the result {x -> 0.801070765209218366216867854087} and FindRoot[Cos[x^2] - x, {x, 1}, WorkingPrecision -> MachinePrecision] gives the result {x -> 0.801071} – Shutao TANG Jul 13 at 8:56
What's the point of that? Just do N[data, precision] and proceed as usual... – J. M. Jul 13 at 9:07
I don't agree with you. When Precison[data]!=Infinity or Precison[data]<precision , the code N[data, precision] maybe fail(just gives the result that owns Min[Precison[data],precision] precision). For example, please see here – Shutao TANG Jul 13 at 9:17
Yes, that's known; N[] can only maintain or lower the precision of a number. That has nothing to do with the algorithm I presented here; it does not care about the precision of whatever is fed to it. – J. M. Jul 13 at 9:19
...again: just feed data of the appropriate precision into the algorithm; there's nothing that needs to be modified. All the arithmetic will proceed as usual, and barring catastrophic cancellation or something, the result will have more or less the same precision. Clearly, you missed the point: for exact input, it will give exact output; for inexact input, it will (or is supposed to) give output of (nearly) the same precision. – J. M. Jul 13 at 9:26

As OP used in question

f = BSplineFunction[data, SplineDegree -> 3, WorkingPrecision -> 8000]


So WorkingPrecision is a very useful option in various numerical opertions.

In this answer, I will give a function BSplineValue that owns the WorkingPrecision option.

The definition of B-Spline curve

$$\vec{C}(u)=\sum _{i=0}^n N_{i,p}(u) \vec{P}_i \qquad (a\leq u\leq b)$$

where, $P_i$ is the control point, the $N_ {i, p} (u)$ are the pth - degree Bspline basis functions defined on the nonperiodic (and nonuniform) knot vector.

In generally, when $u_0 \in [u_i,u_{i+1})$, the B-Spline function basis $N_ {i-p, p},\cdots,N_{i, p}$ are not equal to 0.

So $$\vec{C}(u_0)=\sum _{i=0}^n N_{i,p}(u_0) \vec{P}_i=\sum _{j=i-p}^p N_{j,p}(u_0) \vec{P}_j$$

Search the index of span $[u_i,u_{i+1})$

 biSearch[knots_, {low_, high_}, u_] :=
With[{mid = Floor[(low + high)/2]},
If[u < knots[[mid]],
{low, mid}, {mid, high}]
](*Do binary search*)

searchSpan[knots_, p_, u_] :=
First@
NestWhile[
biSearch[knots, #, u] &, {p + 1, Length@knots - p},
Subtract @@ # != -1 &] - 1


### Main Implementation

Options[BSplineValue] = {WorkingPrecision -> Automatic};

BSplineValue[pts : {{_, _} ..}, deg_, u0_, opts : OptionsPattern[]] :=
Module[{workPrec, idx, knots, n = Length[pts], Pts, Knots, U0},
workPrec = OptionValue[WorkingPrecision];

knots =
Join[ConstantArray[0, deg],
Range[0, 1, 1/(n - deg)], ConstantArray[1, deg]];
idx = searchSpan[knots, deg, u0];

If[workPrec === Automatic,
workPrec =
Max[Precision@{pts, knots, u0}, \$MachinePrecision]];

(*set the precision by the built-in SetPrecision*)
{Pts, Knots, U0} = SetPrecision[{pts, knots, u0}, workPrec];

(BSplineBasis[{deg, Knots}, #, U0] & /@
Range[idx - deg, idx]).Pts[[idx - deg + 1 ;; idx + 1]]
]


### Test

data =
{{17980, -1}, {17990, -1}, {18000, 0}, {18010, 1}, {18020, 3},
{18040, 0}, {18050, -1}, {18080, 2}, {18100, 0}};

BSplineValue[data, 3, 1/5]

 {1350118/75, 119/750}

BSplineValue[data, 3, 1./5, WorkingPrecision -> 30]

 {18001.573333333333334106048558, 0.15866666666666674527045681013}

BSplineValue[data, 3, 1./5, WorkingPrecision -> 100]


-
BSplineValue[data, 3, 1./5, WorkingPrecision -> 30] seems like atypical behavior, since the third argument is a machine precision number; one would thus expect a machine precision result, even with your WorkingPrecision setting. – J. M. Jul 14 at 3:09
This now brings me to the point I was trying to hammer earlier: instead of trying to put in that option, one can just as well do BSplineValue[data, 3, N[1/5, 30]] or BSplineValue[N[data, 30], 3, 1/5] or variations thereof. – J. M. Jul 14 at 3:12
@Guesswhoitis., In this case, the data always contains Integer, so you can use BSplineValue[N[data, 30], 3, 1/5] . However, when data = {{17980, -1.2}, {17990, -1.5}, {18000, 0}, {18010, 1.5}, {18020, 3}, {18040, 0}, {18050, -1}, {18080, 2}, {18100, 0}}, then BSplineValue[N[data, 30], 3, 1/5] will fail. – Shutao TANG Jul 14 at 8:05
No, I would consider that expected behavior instead of "failure" - "garbage in, garbage out". You have the machine precision numbers in there, so the N[] does nothing, and your output is in machine precision. – J. M. Jul 14 at 8:07