Benchmarks for Vectorised Simpsons Rule

On a 16 Thread Ryzen 7. It came in slightly slower. My Global Warming Infrared Attention code, took, 2h 59m 23s non vectorised, and 3h 4mins and 42 seconds, vectorised.

A Vectorised Simpsons Rule for JDK 23

If you are doing numerical integration in Java, try this code on the latest JDK 23, it using the new Vectorisation methods for a faster result.

import jdk.incubator.vector.DoubleVector;

import java.util.Arrays;
import java.util.List;
import java.util.concurrent.*;
import java.util.concurrent.atomic.AtomicReference;

// Release Candidate Version, package will change in Release Version, true for JDK 23
import  jdk.incubator.vector.DoubleVector;
import jdk.incubator.vector.VectorOperators;
import jdk.incubator.vector.VectorOperators.Operator;
import jdk.incubator.vector.VectorSpecies;

import static jdk.incubator.vector.VectorOperators.ADD;


public class SimpsonsRuleVectorised {



    private static double third = 1.0/3.0;

    private static int THREADS = 16;

    public static double integrateThreaded(double a, double b, int N, DoubFunction f, double ...params) {         // precision parameter
        double h = (b - a) / (N - 1);     // step size
        ForkJoinPool pool = new ForkJoinPool(THREADS);
        double d[] = new double[N];

        for(int i=0; i {
                double mul = ii%2==0? 2*third: 4*third;
                if (ii==0) { mul = third; }
                if (ii==N-1){ mul = third; }
                double mul1 = mul;
                double x = a + h * ii;
                double fi = f.eval(x,ii, params);
                if (Double.isNaN(fi) ){
                    System.err.println(f.getClass().getName() + "IS NaN at "+x);
                }
                d[ii] = mul1*fi;
            });
        }
        try {
            pool.shutdown();
            if (!pool.awaitTermination(1, TimeUnit.HOURS)){
                pool.shutdownNow();
            }
        } catch (Exception e){}
        DoubleVector doubleVector = DoubleVector.fromArray(DoubleVector.SPECIES_64, d, 0);
        double sum = doubleVector.reduceLanes(ADD );

        return sum * h;
    }

    public static double integrate(double a, double b, int N, DoubFunction f, double ...params) {         // precision parameter
        double h = (b - a) / (N - 1);     // step size
        double fa = f.eval(a,0, params);
        double fb = f.eval(b, N-1, params);
        if (Double.isNaN(fa) ){
          System.err.println(f.getClass().getName() + "IS NaN at "+a);
        }
        if (Double.isNaN(fb)){
            System.err.println(f.getClass().getName() + "IS NaN at "+b);
        }
        // 1/3 terms
        double sum = third * (fa + fb);


        // 4/3 terms
        for (int i = 1; i < N - 1; i += 2) {
            double x = a + h * i;
            double fx = f.eval(x,i, params);
            if (Double.isNaN(fx)){
                System.err.println(f.getClass().getName() + "IS NaN at "+x);
            }
            sum += 4.0 * third * fx;
        }

        // 2/3 terms
        for (int i = 2; i < N - 1; i += 2) {
            double x = a + h * i;
            double fx = f.eval(x,i, params);
            if (Double.isNaN(fx)){
                System.err.println(f.getClass().getName() + "IS NaN at "+x);
            }
            sum += 2.0 * third * fx;
        }

        return sum * h;
    }

    public static double integrateConsecutive(double a, double b, int N, DoubFunction f, double ...params) {         // precision parameter
        double h = (b - a) / (N - 1);     // step size
        double fa = f.eval(a,0, params);
        if (Double.isNaN(fa) ){
            System.err.println(f.getClass().getName() + "IS NaN at "+a);
        }

        // 1/3 terms
        double sum = third*fa;


        // 4/3 terms
        boolean isOdd = true;
        for (int i = 1; i < N - 1; i += 1) {
            double x = a + h * i;
            double fx = f.eval(x,i, params);
            if (Double.isNaN(fx)){
                System.err.println(f.getClass().getName() + "IS NaN at "+x);
            }
            if (isOdd) {
                sum += 4.0 * third * fx;
            } else {
                sum += 2.0 * third * fx;
            }
            isOdd=!isOdd;
        }


        double fb = f.eval(b, N-1, params);
        if (Double.isNaN(fb)){
            System.err.println(f.getClass().getName() + "IS NaN at "+b);
        }
        sum = sum+ fb*third;

        return sum * h;
    }

    public static void main(String args[]){
        // Roots of polynumerial to integrate
        List in = Arrays.asList(-0.9, -0.8,-0.7, -0.6,-0.5, -0.4, -0.3, -0.2, -0.1, 0, .1,.2,.3,.4,.5, .6,.7, .8, .9 );
        DoubFunction func = new DoubFunction() {
            @Override
            double evalInner(double x, double[] params, int i) {
                return in.stream().map(y->y.doubleValue()-x).reduce(1.0,(a,b)->(a*b));
            }
        };
        double consec=0;
        long startConsec = System.currentTimeMillis();
        for(int i=1;i<1000; i++) {
            consec = integrateConsecutive(-1, 1, 100000, func);
        }
        double timeConsec = (System.currentTimeMillis() - startConsec)/1000.0;
        double standard=0;
        long startStandard = System.currentTimeMillis();
        for(int i=1; i<1000; i++) {
                standard = integrate(-1,1,100000,func);
        }
        double timeStandard = (System.currentTimeMillis() - startStandard)/1000.0;
        double threaded=0;
        long startThreaded = System.currentTimeMillis();
        for(int i=1;i<1000; i++) {
            threaded = integrateThreaded(-1, 1, 100000, func);
        }
        double timeThreaded = (System.currentTimeMillis() - startThreaded)/1000.0;
        System.out.println("Standard Integrator: "+standard+" time taken: "+timeStandard+" seconds");
        System.out.println("Consecutive Integrator: "+consec+" time taken: "+timeConsec+" seconds");
        System.out.println("threaded Integrator: "+threaded+" time taken: "+timeThreaded+" seconds");
    }

}

class AtomicDoubleLocal77 {
    private AtomicReference value = new AtomicReference(Double.valueOf(0.0));
    double getAndAdd(double delta) {
        while (true) {
            Double currentValue = value.get();
            Double newValue = Double.valueOf(currentValue.doubleValue() + delta);
            if (value.compareAndSet(currentValue, newValue))
                return currentValue.doubleValue();
        }
    }
}
  

Recent Paper has new limit on neutrino - majoron interactions from the SN1987 Supernova

The supernova of 1987, SN1987a which occurred in the large magilangic cloud, is the only supernova where neutrinos (some 24), interactions on earth have been recorded. These recorded interactions have now allow the researcher, P.I. Ballesteros and Christa Volpe, https://arxiv.org/abs/2410.11517, to limit potential interactions between neutrinos and majorons (a Majoron is a spin zero particle like a light Huggs, giving mass just to the neutrino). The limits are force contant of around 10^-7 compare. This might also limit an axial force, although that would be a spin-1 psuedovector particle, previously we look at a force contant in the range a few*10^-5 so this paper might simiilar reduce the limits on strength of how force by a factor of a hundred. The previous work on SN1987, from the year 2000, and published in Phy Rev D, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.62.023004 limited the majoron interactiom in the range, 3×10−7≲𝑔≲2×10−5 or 𝑔≳3×10−4 so left open the range 2*10^-5 to 3*10-4 allow or orignal guess of force constant (making the weak assumption that majoron and axi-photon limits would be similar.)